Wstęp

Cel prezentacji: demonstracja wykorzystania pakietu ggRandomForests do wizualizacji lasów losowych.

Uwaga: Pakiet ggRandomForests służy do wizualizacji lasów losowych stworzonych za pomocą pakietu randomForestSRC. Choć wiele wykresów tworzonych przez nowy pakiet można otrzymać przy pomocy randomForestSRC, to korzystanie z ggRandomForests daje trzy (?) udogodnienia:

Przypomnienie: lasy losowe to modyfikacja baggingu, czyli agregacji wyniku zbioru drzew dla ustalonej liczby próbek bootstrapowych, polegająca na “uniezależnieniu” od siebie drzew z tego zbioru poprzez rozpatrywanie przy każdym podziale drzewa jednynie \(m\leq p\) zmiennych wylosowanych spośród \(p\) wszystkich zmiennych.

Plan: prześledzimy przykład dla regresji zaprezentowany przez autora pakietu w jego vignette (Ehrling, 2015), bardziej szczegółowo omawiając miary prezentowane przez poszczególne wykresy i idee za nimi idące.

Dane

Będziemy korzystać ze zbioru danych “Boston” dotyczącego mieszkań w obszarach Bostonu, dostępnego m.in. w pakiecie MASS. Naszym celem jest stwierdzenie co i w jaki sposób wpływa na medianę ceny mieszkań w obszarze, przy czym zbiór zawiera następujące zmienne:

W celu stwierdzenia, które zmienne powinny wyjść “istotne” w naszej analizie rysujemy wykres każdej zmiennej ciągłej względem medv czyli tej, którą chcemy prognozować.

data(Boston, package="MASS")
Boston$chas <- as.logical(Boston$chas)
dta <- melt(Boston, id.vars=c("medv","chas"))
ggplot(dta, aes(x=medv, y=value, color=chas))+ geom_point(alpha=.4)+ geom_rug(data=dta %>% filter(is.na(value)))+ labs(y="", x="Median value of hoems ($1000s)")+ facet_wrap(~variable, scales="free_y", ncol=3)

Budujemy las:

rfsrc_Boston <- rfsrc(medv~., data = Boston, tree.err = TRUE, importance = "permute")
varsel_Boston <- var.select(rfsrc_Boston)
## minimal depth variable selection ...
## 
## 
## -----------------------------------------------------------
## family             : regr 
## var. selection     : Minimal Depth 
## conservativeness   : medium 
## x-weighting used?  : TRUE 
## dimension          : 13 
## sample size        : 506 
## ntree              : 1000 
## nsplit             : 0 
## mtry               : 5 
## nodesize           : 5 
## refitted forest    : FALSE 
## model size         : 11 
## depth threshold    : 6.749 
## PE (true OOB)      : 10.8297 
## 
## 
## Top variables:
##         depth   vimp
## rm      1.168 37.523
## lstat   1.196 57.394
## crim    2.463  8.406
## dis     2.596  5.307
## nox     2.618  8.195
## ptratio 2.945  5.830
## age     3.555  3.355
## tax     3.659  3.539
## indus   3.680  5.474
## black   3.769  1.116
## rad     6.251  1.216
## -----------------------------------------------------------
rfsrc_Boston
##                          Sample size: 506
##                      Number of trees: 1000
##           Minimum terminal node size: 5
##        Average no. of terminal nodes: 104.156
## No. of variables tried at each split: 5
##               Total no. of variables: 13
##                             Analysis: RF-R
##                               Family: regr
##                       Splitting rule: mse
##                 % variance explained: 87.2
##                           Error rate: 10.83

Błąd prognozy Out-of-Bag

Próbka Out-of-Bag (OOB): dla każdej obserwacji wyznaczamy jej prognozę z lasu biorąc pod uwagę tylko te drzewa, które uczono na próbce bez tej obserwacji. (Każda próbka bootstrapowa zawiera około 63.2% obserwacji).

Błąd prognozy OOB jest praktycznie taki sam, jak błąd obliczony za pomocą \(n\)-krotnej cross-walidacji, czyli w jednym kroku otrzymujemy zarówno las jak i oszacowanie błędu na zbiorze testowym.

Błąd ten jest liczony dla każdego drzewa przy budowie lasu (opcja tree.err = TRUE funkcji rfsrc()) i można go łatwo narysować za pomocą funkcji plot.gg_error(), która właściwie tylko wrzuca to do ggplota.

gg_e <- gg_error(rfsrc_Boston)
plot(gg_e)

Wniosek: Jak widać błąd OOB w tym przypadku stabilizuje się już przy 250 drzewach, ale dla celów analizy ważnych zmiennych lepiej jest mieć ich więcej, żeby każda zmienna była uwzględniona w dostatecznej ich ilości.

Prognoza lasu losowego

Funkcja gg_rfsrc() zwraca prognozy lasu (domyślnie OOB, mogą być też zwykłe), które można przedstawić na wykresie, choć to nie jest zbyt ciekawe…

plot(gg_rfsrc(rfsrc_Boston), alpha=.5)

Ważność zmiennych

Problem: ponieważ las losowy nie zakłada żadnej postaci funkcyjnej nie daje też możliwości testowania istotności zmiennych. Warto pamiętać, że to jest zarazem największa zaleta lasu (drzew), gdyż pozwala na dowolną nieliniowość zależności między zmiennymi.

Co można zrobić? Analizować strukturę lasu, czyli drzew wchodzących w jego skład. Ponieważ jest ich dużo wszelkie miary muszą być łatwe w agregacji. Dwa podstawowe podejścia, to:

Variable importance

VIMP liczymy porównując błąd predykcji OOB przed i po zakłóceniu danej zmiennej.

Jak zakłócać zmienne?

  • permutacja zmiennej (permute)

  • predykcja OOB jest wyznaczana losowo na węzłach dzielących populację ze względu na wartości danej zmiennej (random)

  • \(x\) is assigned to the opposite node whenever a split on \(x\) is encountered in the in-bag tree” (anti)

Funkcja vimp() z pakietu randomForestSRC liczy VIMP według wybranej metody zakłócania (do wyboru: "permute", "random", "anti", "permute.ensemble", "random.ensemble", "anti.ensemble", opcje z ensemble są szybsze). W nowym pakiecie dostajemy to za pomocą plot.gg_vimp(), która tworzy wykres ggplot:

plot(gg_vimp(rfsrc_Boston, importance = "permute"))

Uwaga: VIMP może być ujemne - oznacza to, że szum jest bardziej informatywny od danej zmiennej. Nie ma żadnej referencyjnej wartości VIMP, do której można by porównywać otrzymane żeby stwierdzić czy są “wystarczająco duże”.

Minimal depth

Węzły w drzewie numerujemy od 0 (korzeń), czyli według ich odległości od korzenia. Minimal depth to średnia odległość od korzenia pierwszych węzłów dzielących populację według wybranej zmiennej, więc mniejsze wartości tej miary wskazują na zmienne częściej biorące udział w dzieleniu znacznej części populacji, czyli mające większy wpływ na sortowanie obserwacji.

Pakiet randomForestSRC dostarcza rankingu zmiennych pod względem minimal depth (funkcja var.select()), które oblicza znajdując maksymalne poddrzewa względem każdej ze zmiennych, czyli poddrzewa z korzeniem dzielącym próbę ze względu na tę zmienną.

Nowy pakiet zamienia wynik var.select() na ramkę danych (funkcja gg_minimal_depth()) i tworzy wykres ggplot za pomocą plot().

gg_md <- gg_minimal_depth(varsel_Boston)
plot(gg_md)

W przeciwieństwie do VIMP, dla minimal depth można wyznaczyć wartość progową (pionowa linia przerywana), do której porównuje się pozostałe wartości, żeby stwierdzić, czy są duże, czy nie. W tym przypadku wzięto średnią z rozkładu minimal depth.

Za pomocą plot.gg_minimal_vimp() można porównać pozycję zmiennych w rankingu VIMP oraz minimal depth:

plot(gg_minimal_vimp(gg_md))

Maksymalne poddrzewa

Ishwaran et al. (2010) wprowadzają pojęcie maksymalnych poddrzew dla zmiennej \(v\):

Definicja: Dla każdej zmiennej \(v\), \(T_v\) nazywamy \(v\)-poddrzewem \(T\) jeśli korzeń \(T_v\) jest dzielony względem zmiennej \(v\). \(T_v\) nazywamy maksymalnym \(v\)-poddrzewem jeśli \(T_v\) nie jest poddrzewem żadnego większego \(v\)-drzewa.

Uwaga: maksymalne \(v\)-poddrzewo może nie istnieć lub może ich być kilka.

Maksymalne poddrzewa są wartościową koncepcją, ponieważ:

  • pozwalają mierzyć ważność zmiennej - rzadko zdarza się, że mało informatywne zmienne stanowią kryterium podziału węzłów blisko korzenia, więc to znika przy uśrednianiu,

  • można za ich pomocą identyfikować ważne interakcje zmiennych, np. poprzez analizę \((v,w)\)-maksymalnych poddrzew, czyli maksymalnych \(w\)-poddrzew wewnątrz maksymalnych \(v\)-poddrzew,

  • dotyczą wyłącznie struktury drzewa (np. zasad podziału), nie są związane z błędem predykcji, który nie zawsze jest najważniejszy w analizie (np. w analizie przeżycia nie jest jasne co jest dobrą miarą błędu predykcji),

  • w przeciwieństwie do VIMP, które opierają się na losowaniu, maksymalne poddrzewa można dokładnie analizować.

Definicja: Niech \(D_v\) będzie odległością od korzenia całego drzewa do korzenia najbliższego maksymalnego \(v\)-poddrzewa dla ustalonego \(v\). Wtedy \(D_v\) jest nieujemną zmienną losową przyjmującą wartości ze zbioru \(\{0,\ldots,D(T)\}\), gdzie \(D(T)\) jest wysokością drzewa. \(D_v\) to minimal depth \(v\).

Twierdzenie: Niech \(D(T)\geq 1\) będzie ustalone, a liczba węzłów o głębokości d, \(l_d\), wynosi \(2^d\). Załóżmy, że \(\pi_{v,j}(t)\) i \(\theta_{v,j}(t)\) zależą jedynie od głębokości węzła \(t\), wtedy: \[\mathbb{P}(D_v=d)= \left[\prod_{j=0}^{d-1}(1-\pi_{v,j}\theta_{v,j})^{l_j}\right][1- (1-\pi_{v,d}\theta_{v,d})^{l_d}],\quad 0\leq d\leq D(T)-1,\]

gdzie:

  • \(\pi_{v,j}(t)\) oznacza prawdopodobieństwo, że \(v\) jest kandydatką (jedną z \(m\leq p\)) do podziału na węźle \(t\) o głębokości \(j\) przy założeniu, że nie istnieje maksymalne \(v\)-poddrzewo na głębokości mniejszej niż \(j\),

  • \(\theta_{v,j}(t)\) oznacza prawdopodobieństwo, że \(v\) dzieli węzeł \(t\) o głębokości \(j\) przy założeniu, że \(v\) jest kandydatką do podziału na węźle \(t\) oraz nie istnieje maksymalne \(v\)-poddrzewo na głębokości mniejszej niż \(j\).

Dowód: Z definicji \(\pi_{v,j}\) i \(\theta_{v,j}\) jeśli \(t\) jest węzłem o głębokości \(j\leq d\), to \[\mathbb{P}(v \text{ nie dzieli } t|D_v\geq j)=1-\pi_{v,j}\theta_{v,j}.\]

Ponadto, \(\mathbb{P}(D_v=0)=\pi_{v,0}\theta_{v,0}\). Stąd, prawdopodobieństwo, że nie istnieje maksymalne \(v\)-poddrzewo na głębokości mniejszej niż \(d\geq 1\) wynosi \[\mathbb{P}(D_v\geq d)=\mathbb{P}(D_v\geq 1) \prod_{j=1}^{d-1} \mathbb{P}(v\text{ nie dzieli na głębokości } j|D_v\geq j)=\] \[=[1-\mathbb{P}(D_v=0)]\prod_{j=1}^{d-1} \mathbb{P}(v\text{ nie dzieli na głębokości } j|D_v\geq j) = \prod_{j=0}^{d-1} \mathbb{P}(1-\pi_{v,j}\theta_{v,j})^{l_j},\] gdzie \(l_j=2^j\) jest równe całkowitej liczbie węzłów na głębokości \(j\). Ponieważ \(D_v\geq d\), prawdopodobieństwo, że \(v\) dzieli węzeł na głębokości \(d\) wynosi 1 minus prawdopodobieństwo, że każdy węzeł na głębokości \(d\) dzielony jest inną zmienną niż \(v\) czyli \(1-(1-\pi_{v,d}\theta_{v,d})^{l_d}\).

Ostateczny wynik dostajemy ze wzoru \[\mathbb{P}(D_v=d)=\mathbb{P}(v\text{ dzieli na głębokości } d|D_v\geq d)\cdot \mathbb{P}(D_v\geq d).\]

Zależność między prognozą i zmiennymi

VIMP i minimal depth mogą być podstawą do wyboru mniejszego podzbioru z rozpatrywanych zmiennych. W tym celu chcemy zbadać jak poszczególne zmienne wpływają na predykcję lasu.

Uwaga: choć las losowy jest często nazywany “czarną skrzynką”, to jego prognoza jest przecież (skomplikowaną) funkcją rozważanych zmiennych: \[\hat{f}_{rf}=f(x).\] Do analizy powyższej zależności najlepiej wykorzystać metody graficzne.

Variable dependence

Wykresy variable dependence pokazują zależność prognozy od określonej zmiennej. Każdy punkt na wykresie to jedna obserwacja, a ponieważ wartość prognozy jest zależna od wszystkich zmiennych jednocześnie, to interpretację wykresu należy przeprowadzać z uwagą.

Funkcja gg_variable() wyciąga zmienne wykorzystane do uczenia lasu oraz prognozę OOB z lasu oraz obiektu typu predict z pakietu randomForestSRC. Natępnie, funkcja plot.gg_variable() zwraca listę wykresów ggplot dla zmiennych zadeklarowanych jako xvar lub jeden wykres jeśli wybrano opcję panel=TRUE.

W naszym przykładzie dla 10 zmiennych o najniższym minimal depth dostajemy:

gg_v <- gg_variable(rfsrc_Boston)
xvar <- gg_md$topvars[1:10]
plot(gg_v, xvar=xvar, panel=TRUE, alpha=.4)+geom_smooth(se=TRUE, level=.95, span=1.2)+ labs(y="Median value of homes ($1000s)", x="")

Jak widać, wykresy te są podobne do tych, które otrzymaliśmy na początku przy eksploracji danych, czego można się spodziewać jeśli dany las losowy dobrze prognozuje.

Uwaga: dla zmiennych jakościowych odpowiednim wykresem jest boxplot, stąd dla nich rysujemy osobny wykres (w tym przypadku mamy tylko jedną taką zmienną, w dodatku mało ważną według wszelkich miar).

plot(gg_v, xvar="chas", points=FALSE, alpha=.4)+ geom_boxplot(notch=TRUE, na.rm=TRUE)+ labs(y="Median value of homes ($1000s)")

Partial dependence

Wykresy partial dependence pokazują zależność predykcji od określonej zmiennej po wyłączniu zależności od pozostałych zmiennych, co liczy się następująco: dla wybranej zmiennej \(x\) wybierane są równoodległe punkty wzdłuż jej rozkładu, aby następnie dla każdego z nich obliczyć średnią predykcję lasu dla wszystkich pozostałych zmiennych, czyli: \[\tilde{f}(x)=\frac{1}{n}\sum_{i=1}^n \hat{f}(x,x_{i,o}),\] gdzie \(\hat{f}\) jest prognozą lasu, a \(x_{i,o}\) to wartości wszystkich zmiennych poza \(x\) dla obserwacji \(i\).

Idea: uśredniamy predykcje każdej obserwacji ze zbioru uczącego dla każdej z wybranego zbioru wartości \(X=x\).

Funkcja gg_partial() tworzy listę obiektów (po jednym dla każdej zmiennej) korzystając z wyniku funkcji plot.variable() z pakietu randomForestSRC, co następnie można narysować korzystając z plot.gg_partial().

Uwaga: wykresy cząstkowe są znacznie bardziej wymagające obliczeniowo, co należy mieć na uwadze wybierając liczbę punktów z rozkładu każdej ze zmiennych, które mają być rozpatrzone (opcja npts funkcji plot.variable()).

partial_Boston <- plot.variable(rfsrc_Boston, xvar=xvar, partial=TRUE, sorted=FALSE, show.plots = FALSE)
gg_p <- gg_partial(partial_Boston)
plot(gg_p, panel=TRUE)+ labs(y="Median value of homes ($1000s)", x="") + geom_smooth(se=FALSE)

Wniosek: teraz wyraźnie widać, że zmienne lstat i rm odgrywają dużo większą rolę w prognozowaniu niż pozostałe zmienne (patrz: płaskie wykresy); jest to zależność nieliniowa!

Analiza interakcji

Interakcje minimal depth

Najbardziej podstawowa koncepcja analizy struktury lasów, minimal depth, pozwala również analizować interakcje między zmiennymi. Zamiast \(v\)-maksymalnych poddrzew wystarczy rozpatrywać \((v,w)\)-maksymalne poddrzewa, czyli maksymalne \(w\)-poddrzewa wewnątrz maksymalnych \(v\)-poddrzew.

Funkcja randomForestSRC::find.interaction dostarcza macierzy \(p\times p\) z uogólnioną miarą minimal depth (opcja method="maxsubtree") dla każdej pary zmiennych, przy czym:

  • wyrazy na diagonali \((i,i)\) to unormowane względem korzenia (wielkości drzewa) minimal depth zmiennej \(i\),

  • wyrazy poza diagonalą \((i,j)\) to minimal depth zmiennej \(j\) względem maksymalnego \(i\)-poddrzewa unormowane względem rozmiaru tego maksymalnego poddrzewa.

W nowym pakiecie funkcja plot.gg_interaction() daje odpowiedni wykres:

interaction_Boston <- find.interaction(rfsrc_Boston, method = "maxsubtree", sorted = TRUE)
## 
##                               Method: maxsubtree
##                     No. of variables: 13
##   Variables sorted by minimal depth?: TRUE
## 
##           rm lstat crim  dis  nox ptratio  age  tax indus black  rad   zn
## rm      0.08  0.15 0.20 0.19 0.23    0.25 0.23 0.25  0.31  0.25 0.43 0.57
## lstat   0.15  0.08 0.18 0.16 0.20    0.26 0.22 0.25  0.30  0.24 0.44 0.59
## crim    0.25  0.21 0.16 0.26 0.32    0.47 0.29 0.48  0.49  0.32 0.63 0.78
## dis     0.21  0.21 0.28 0.17 0.33    0.38 0.27 0.38  0.41  0.31 0.57 0.72
## nox     0.29  0.25 0.27 0.29 0.17    0.51 0.34 0.50  0.54  0.36 0.66 0.79
## ptratio 0.26  0.27 0.34 0.33 0.43    0.19 0.33 0.43  0.47  0.38 0.62 0.73
## age     0.27  0.29 0.31 0.32 0.41    0.45 0.23 0.45  0.44  0.36 0.65 0.78
## tax     0.32  0.36 0.42 0.40 0.50    0.51 0.41 0.24  0.54  0.46 0.66 0.75
## indus   0.35  0.36 0.43 0.41 0.51    0.52 0.41 0.53  0.24  0.46 0.67 0.77
## black   0.37  0.34 0.39 0.39 0.52    0.59 0.39 0.61  0.58  0.24 0.73 0.87
## rad     0.62  0.65 0.69 0.69 0.79    0.77 0.68 0.77  0.79  0.72 0.40 0.91
## zn      0.80  0.84 0.84 0.84 0.90    0.88 0.83 0.88  0.89  0.85 0.93 0.54
## chas    0.82  0.80 0.83 0.83 0.88    0.91 0.83 0.91  0.92  0.85 0.94 0.98
##         chas
## rm      0.71
## lstat   0.68
## crim    0.79
## dis     0.81
## nox     0.82
## ptratio 0.85
## age     0.84
## tax     0.88
## indus   0.85
## black   0.88
## rad     0.95
## zn      0.99
## chas    0.63
plot(gg_interaction(interaction_Boston), xvar=xvar, panel = TRUE)

Interpretacja: podobnie jak dla zwyczajnego minimal depth szukamy małych wartości.

Interakcje VIMP

W pakiecie randomForestSRC można również wyznaczyć miarę ważności interakcji opartą na VIMP, która nie jest wykorzystana w pakiecie ggRandomForests.

Idea: dla każdej pary zmiennych liczymy ich łączny VIMP (paired VIMP) oraz sumę ich indywidualnych VIMP (additive VIMP). Duże różnice między tymi dwoma miarami sugerują ważne interakcje (zaburzenie obu zmiennych naraz ma inny skutek niż zaburzenie tych zmiennych osobno).

interaction_Boston <- find.interaction(rfsrc_Boston, importance = "permute", method = "vimp", sorted = TRUE)
## Pairing lstat with rm 
## Pairing lstat with crim 
## Pairing lstat with nox 
## Pairing lstat with ptratio 
## Pairing lstat with indus 
## Pairing lstat with dis 
## Pairing lstat with tax 
## Pairing lstat with age 
## Pairing lstat with rad 
## Pairing lstat with black 
## Pairing lstat with chas 
## Pairing lstat with zn 
## Pairing rm with crim 
## Pairing rm with nox 
## Pairing rm with ptratio 
## Pairing rm with indus 
## Pairing rm with dis 
## Pairing rm with tax 
## Pairing rm with age 
## Pairing rm with rad 
## Pairing rm with black 
## Pairing rm with chas 
## Pairing rm with zn 
## Pairing crim with nox 
## Pairing crim with ptratio 
## Pairing crim with indus 
## Pairing crim with dis 
## Pairing crim with tax 
## Pairing crim with age 
## Pairing crim with rad 
## Pairing crim with black 
## Pairing crim with chas 
## Pairing crim with zn 
## Pairing nox with ptratio 
## Pairing nox with indus 
## Pairing nox with dis 
## Pairing nox with tax 
## Pairing nox with age 
## Pairing nox with rad 
## Pairing nox with black 
## Pairing nox with chas 
## Pairing nox with zn 
## Pairing ptratio with indus 
## Pairing ptratio with dis 
## Pairing ptratio with tax 
## Pairing ptratio with age 
## Pairing ptratio with rad 
## Pairing ptratio with black 
## Pairing ptratio with chas 
## Pairing ptratio with zn 
## Pairing indus with dis 
## Pairing indus with tax 
## Pairing indus with age 
## Pairing indus with rad 
## Pairing indus with black 
## Pairing indus with chas 
## Pairing indus with zn 
## Pairing dis with tax 
## Pairing dis with age 
## Pairing dis with rad 
## Pairing dis with black 
## Pairing dis with chas 
## Pairing dis with zn 
## Pairing tax with age 
## Pairing tax with rad 
## Pairing tax with black 
## Pairing tax with chas 
## Pairing tax with zn 
## Pairing age with rad 
## Pairing age with black 
## Pairing age with chas 
## Pairing age with zn 
## Pairing rad with black 
## Pairing rad with chas 
## Pairing rad with zn 
## Pairing black with chas 
## Pairing black with zn 
## Pairing chas with zn 
## 
##                               Method: vimp
##                     No. of variables: 13
##            Variables sorted by VIMP?: TRUE
##    No. of variables used for pairing: 13
##     Total no. of paired interactions: 78
##             Monte Carlo replications: 1
##     Type of noising up used for VIMP: permute
## 
##                 Var 1   Var 2  Paired Additive Difference
## lstat:rm      56.8039 36.9035 96.8927  93.7074     3.1852
## lstat:crim    56.8039  8.1213 59.6208  64.9252    -5.3044
## lstat:nox     56.8039  8.2821 64.6194  65.0860    -0.4666
## lstat:ptratio 56.8039  5.9608 63.2144  62.7647     0.4497
## lstat:indus   56.8039  5.5189 61.4704  62.3228    -0.8524
## lstat:dis     56.8039  5.3315 59.5169  62.1355    -2.6185
## lstat:tax     56.8039  3.5309 59.2245  60.3348    -1.1104
## lstat:age     56.8039  3.3444 57.4446  60.1483    -2.7037
## lstat:rad     56.8039  1.1802 56.9792  57.9842    -1.0050
## lstat:black   56.8039  1.0915 58.1270  57.8954     0.2316
## lstat:chas    56.8039  0.4194 57.3192  57.2233     0.0959
## lstat:zn      56.8039  0.3105 57.0955  57.1144    -0.0190
## rm:crim       36.9871  8.1213 44.4437  45.1084    -0.6647
## rm:nox        36.9871  8.2821 45.7231  45.2692     0.4538
## rm:ptratio    36.9871  5.9608 43.3398  42.9479     0.3918
## rm:indus      36.9871  5.5189 42.5835  42.5060     0.0774
## rm:dis        36.9871  5.3315 41.6360  42.3187    -0.6827
## rm:tax        36.9871  3.5309 41.1085  40.5180     0.5905
## rm:age        36.9871  3.3444 39.8959  40.3315    -0.4357
## rm:rad        36.9871  1.1802 37.8334  38.1674    -0.3340
## rm:black      36.9871  1.0915 38.1124  38.0786     0.0338
## rm:chas       36.9871  0.4194 37.4370  37.4065     0.0304
## rm:zn         36.9871  0.3105 37.3305  37.2976     0.0329
## crim:nox       8.1916  8.2821 15.6469  16.4738    -0.8269
## crim:ptratio   8.1916  5.9608 13.9146  14.1525    -0.2378
## crim:indus     8.1916  5.5189 13.2738  13.7105    -0.4368
## crim:dis       8.1916  5.3315 13.3594  13.5232    -0.1638
## crim:tax       8.1916  3.5309 11.6122  11.7226    -0.1103
## crim:age       8.1916  3.3444 11.9025  11.5361     0.3664
## crim:rad       8.1916  1.1802  9.2566   9.3719    -0.1152
## crim:black     8.1916  1.0915  9.4164   9.2832     0.1333
## crim:chas      8.1916  0.4194  8.5155   8.6111    -0.0956
## crim:zn        8.1916  0.3105  8.5202   8.5022     0.0181
## nox:ptratio    8.2864  5.9608 14.1916  14.2473    -0.0557
## nox:indus      8.2864  5.5189 13.5124  13.8053    -0.2929
## nox:dis        8.2864  5.3315 12.7904  13.6180    -0.8275
## nox:tax        8.2864  3.5309 11.4727  11.8174    -0.3447
## nox:age        8.2864  3.3444 11.3393  11.6309    -0.2915
## nox:rad        8.2864  1.1802  9.3334   9.4667    -0.1333
## nox:black      8.2864  1.0915  9.3753   9.3780    -0.0027
## nox:chas       8.2864  0.4194  8.6757   8.7059    -0.0302
## nox:zn         8.2864  0.3105  8.4846   8.5969    -0.1123
## ptratio:indus  5.7937  5.5189 11.3393  11.3126     0.0267
## ptratio:dis    5.7937  5.3315 10.8397  11.1252    -0.2855
## ptratio:tax    5.7937  3.5309  9.0704   9.3246    -0.2542
## ptratio:age    5.7937  3.3444  9.1567   9.1381     0.0186
## ptratio:rad    5.7937  1.1802  6.8342   6.9739    -0.1397
## ptratio:black  5.7937  1.0915  6.8853   6.8852     0.0001
## ptratio:chas   5.7937  0.4194  6.1559   6.2131    -0.0572
## ptratio:zn     5.7937  0.3105  6.1367   6.1042     0.0325
## indus:dis      5.3929  5.3315 10.0540  10.7245    -0.6705
## indus:tax      5.3929  3.5309  8.7810   8.9238    -0.1428
## indus:age      5.3929  3.3444  8.7228   8.7374    -0.0145
## indus:rad      5.3929  1.1802  6.5678   6.5732    -0.0053
## indus:black    5.3929  1.0915  6.4417   6.4845    -0.0428
## indus:chas     5.3929  0.4194  5.7930   5.8124    -0.0194
## indus:zn       5.3929  0.3105  5.7123   5.7034     0.0089
## dis:tax        5.2460  3.5309  8.5468   8.7769    -0.2301
## dis:age        5.2460  3.3444  8.4659   8.5904    -0.1245
## dis:rad        5.2460  1.1802  6.4085   6.4263    -0.0177
## dis:black      5.2460  1.0915  6.2908   6.3375    -0.0467
## dis:chas       5.2460  0.4194  5.6870   5.6654     0.0215
## dis:zn         5.2460  0.3105  5.5293   5.5565    -0.0272
## tax:age        3.5739  3.3444  6.8829   6.9183    -0.0355
## tax:rad        3.5739  1.1802  4.6664   4.7542    -0.0878
## tax:black      3.5739  1.0915  4.6105   4.6654    -0.0549
## tax:chas       3.5739  0.4194  3.9642   3.9933    -0.0292
## tax:zn         3.5739  0.3105  3.8621   3.8844    -0.0224
## age:rad        3.2516  1.1802  4.4635   4.4319     0.0316
## age:black      3.2516  1.0915  4.3945   4.3431     0.0513
## age:chas       3.2516  0.4194  3.6624   3.6710    -0.0086
## age:zn         3.2516  0.3105  3.5767   3.5621     0.0145
## rad:black      1.2450  1.0915  2.3329   2.3365    -0.0036
## rad:chas       1.2450  0.4194  1.6695   1.6644     0.0051
## rad:zn         1.2450  0.3105  1.5534   1.5555    -0.0021
## black:chas     1.0791  0.4194  1.5228   1.4985     0.0243
## black:zn       1.0791  0.3105  1.3888   1.3896    -0.0008
## chas:zn        0.4685  0.3105  0.7791   0.7790     0.0002

Coplots

Coplots czyli conditioning plots służą do wizualizacji zależności prognozy od jednej zmiennej pod warunkiem drugiej (lub kilku innych). Obserwacje są grupowane według wartości zmiennych z warunku po czym dla każdej grupy rysujemy wykres variable dependence.

Uwaga: grupowanie najłatwiej przeprowadzić dla zmiennych jakościowych. Dla ilościowych, jak w tym przypadku, często trzeba grupować dość arbitralnie, np. według kwantyli.

W naszym przykładzie najbardziej interesująca jest interakcja między zmiennymi lstat i rm, wykres wygląda następująco:

rm_pts <- quantile_pts(rfsrc_Boston$xvar$rm, groups=6, intervals=TRUE) # cathegorize continuous variables
rm_grp <- cut(rfsrc_Boston$xvar$rm, breaks=rm_pts) # create 6 factor intervals
gg_v$rm_grp <- rm_grp # add the group factor to the gg_variable object
levels(gg_v$rm_grp) <- paste("rm in ", levels(gg_v$rm_grp), sep="") # modify labels
plot(gg_v, xvar = "lstat", smooth = TRUE, alpha = .5)+ geom_smooth(se=FALSE, span=1.5)+ labs(y = "Median value of homes ($1000s).", x="Lower status of the population (percent)")+ theme(legend.position = "none")+ facet_wrap(~rm_grp)

Wniosek: niezależnie od rm mediana ceny mieszkań maleje wraz ze wzrostem udziału w populacji osób o niskim statusie, przy czym zależność ta jest znacznie silniejsza dla większej średniej liczby pokoi w mieszkaniu.

Podsumowanie

Bibliografia