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:
wykresy powstają przy użyciu ggplot2
, przy czym można je łatwo modyfikować, gdyż funkcje zwracają obiekty klasy ggplot
,
otrzymane dane/wykresy stanowią samodzielne obiekty, co pozwala na łatwą modyfikację,
rozdzielenie danych od wykresów: pakiet dostarcza funkcji operujących na obiektach powstałych przy pomocy randomForestSRC
i zwracających dane, które później można przekształcić na wykresy za pomocą innych funkcji z ggRandomForests
.
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.
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
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.
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)
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) opiera się na pomyśle, że jeśli zmienna jest ważna, to jej zakłócenie powinno zwiększyć błąd predykcji.
Minimal depth opiera się na pomyśle, że ważne zmienne występują w drzewie blisko korzenia, ponieważ wtedy na podstawie ich wartości dzielona jest duża część populacji.
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”.
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))
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).\]
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.
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)")
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!
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.
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 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.
Nie taki box czarny: dwie główne miary ważności zmiennych to minimal depth i variable importance, przy czym dla pierwszej z nich możemy obliczyć średnią co daje pewien punkt odniesienia.
W dość przystępny sposób możemy badać dziesiątki interakcji między zmiennymi (wystarczy jeden wykres).
Pakiet ggRandomForests
dostarcza prostych funkcji generujących ładne wykresy w ggplocie.
Ehrling, J., 2015. “ggRandomForests: Random Forests for Regression”
Ishwaran, H., U. B. Kogular, E. Z. Gorodeski, A. J. Minn, M. S. Lauer, 2010. “High-dimensional variable selection for survival data”, Journal of American Statistical Association, 105, 205-217.
dokumentacja pakietu ggRandomForests
dokumentacja pakietu randomForestSRC