Zbiór fitness Przedmiotem zainteresowania będzie zbiór fitness.

(1)

Zbiór fitness

Przedmiotem zainteresowania będzie zbiór fitness.

Dane dotyczą parametrów wydolnościowych mężczyzn zmierzonych podczas biegu na 1,5 mili. Zmienną objaśnianą jest Oxygen (pobór tlenu podczas biegu).

Dla zbioru danych http://akson.sgh.waw.pl/~mz35109/RD2/LAB8/fitness.txt wykonaj następujące czynności:

Dane dotyczą parametrów wydolnościowych mężczyzn zmierzonych podczas biegu na 1,5 mili Zmienną objaśnianą jest Oxygen (pobór tlenu podczas biegu)

(2)

DRZEWO REGRESYJNE

Wywołanie komendy

Spowoduje otworzenie w oknie przeglądarki strony z dokumentacją pakietu rpart.

 co oznacza parametr cp? jaka jest jego wartość domyślna?

 co oznacza parametr minsplit? jaka jest jego wartość domyślna?

Następnie wywołujemy funkcję rysowania drzewa (tekstowo):

(3)

 Jeśli chcemy rysować graficznie:

plot(fitness.rpart,uniform=T,margin=0.1) Jego graficzna reprezentacja:

(4)

Jak widać drzewo nie ma żadnej zawartości opisowej, dzięki czemu jest w zasadzie bezużyteczne w analizie. Dopiero wywołanie dodatkowo funkcji: text(fitness.rpart)

Sprawi, ze wyrysujemy pełne drzewo:

Z kolei wywołanie funkcji: ?plot Sprawi ze w oknie przeglądarki otworzy się dokumentacja funkcji plot.

Z kolei wywołanie funkcji: ?plot.rpart Wyświetli dokumentację funkcji plot.rpart co oznaczają parametry uniform i margin?

To pozwala rozszerzyć swoją wiedzę w tym zakresie.

(5)

 dlaczego uruchamiamy teoretycznie funkcję plot() a szukamy pomocy przez ?plot.rpart?

 Czym są węzły (ang. node) ?

 Jak są numerowane węzły?

 Jak się nazywają pierwszy i ostatnie węzły?

 Czym się charakteryzują węzły nie będące liśćmi?

 Gdy warunek podziału w danym węźle jest spełniony to przechodzimy do lewego czy prawego potomka? Odp.:

do lewego

 Czym się charakteryzuje węzeł będący liściem?

Powróćmy do wydruku tekstowego drzewa: fitness.rpart

Czym są kolumny wydruku: „node”, „split” , „n”, „deviance” oraz „yval” ?

 node), - posługując się graficzną reprezentacją zwróć uwagę na sposób numeracji

 split - czy korzeń ma warunek

 n - zauważ, że suma elementów u dzieci jest równa liczbie elementów u rodzica, np. n2+n3=n1

 deviance - zauważ, że SSE u dzieci jest mniejsza niż SSE u rodzica, np. SSE2+SSE3 < SSE1

 yval - za wartość teoretyczną bierze się średnią, ale można też próbować w liściach budować regresję

 Czym są gwiazdki?

Predykcja

Zechcemy dokonać predykcji na "przeciętnego Kowalskiego", posługując się medianą cech, a w przypadku średnich taki punkt nazywałby się centroidem

oraz

Wtedy

By odpowiedzieć na następujące pytania:

 Ile jest niepowtarzalnych wartości predykcji Oxygen? Podpowiedź: Ile jest liści w drzewie?

 Czy reszty mają rozkład normalny? Oraz czy dla MNK mamy to samo ? Można wywołać następujący ciąg instrukcji R:

(6)

par(mfrow=c(2,2))

plot(predict(fitness.rpart),residuals(fitness.rpart), xlab="Predykcja Oxygen", ylab="Rezydua") qqnorm(residuals(fitness.rpart))

lines(c(-3,3),c(-3,3))

gg=lm(Oxygen ~ ., data=fitness) plot(gg$fitted.values, residuals(gg)) qqnorm(residuals(gg))

Efekt będzie następujący:

Teraz Przejdź "palcem" po wydruku graficznym drzewa i upewnij się, że predykcja jest poprawna.

Ile jest niepowtarzalnych wartości predykcji Oxygen? Podpowiedź: Ile jest liści w drzewie?

Czy reszty mają rozkład normalny? Czy reszty mają rozkład normalny?

qqnorm(residuals(fitness.rpart)) lines(c(-3,3),c(-3,3))

to samo dla MNK...

gg=lm(Oxygen ~ ., data=fitness) plot(gg$fitted.values, residuals(gg)) Teraz wywołaj komendę:

> summary(fitness.rpart)

(7)

Call:

rpart(formula = Oxygen ~ ., data = fitness, cp = 0.01, minsplit = 2) n= 31

CP nsplit rel error xerror xstd 1 0.51397209 0 1.00000000 1.0827673 0.29483014 2 0.27566023 1 0.48602791 0.9650490 0.21891598 3 0.06866676 2 0.21036768 0.4065577 0.10235425 4 0.03364403 3 0.14170092 0.3755164 0.09582809 5 0.02300883 4 0.10805689 0.3932771 0.11757000 6 0.02242706 5 0.08504806 0.5139413 0.13690365 7 0.01852991 6 0.06262100 0.4590597 0.13270992 8 0.01610282 7 0.04409108 0.4680422 0.13311158 9 0.01000000 8 0.02798826 0.4942286 0.13312009

Node number 1: 31 observations, complexity param=0.5139721 mean=47.37581, MSE=27.46392

left son=2 (27 obs) right son=3 (4 obs) Primary splits:

RunTime < 8.935 to the right, improve=0.5139721, (0 missing) RestPulse < 55.5 to the right, improve=0.2291669, (0 missing) RunPulse < 171 to the right, improve=0.2002036, (0 missing) Age < 42.5 to the right, improve=0.1561956, (0 missing) MaxPulse < 190 to the right, improve=0.1210755, (0 missing)

RunTime < 12.29 to the right, improve=0.60962080, (0 missing) MaxPulse < 171 to the right, improve=0.26538490, (0 missing) RestPulse < 55.5 to the right, improve=0.14769050, (0 missing) RunPulse < 167 to the right, improve=0.12249080, (0 missing) Weight < 85.39 to the right, improve=0.08878399, (0 missing)

Age < 43 to the right, improve=0.9940683, (0 missing) RunTime < 8.64 to the right, improve=0.9940683, (0 missing) RunPulse < 161 to the left, improve=0.9940683, (0 missing) MaxPulse < 170 to the left, improve=0.9940683, (0 missing) Weight < 83.855 to the right, improve=0.3732153, (0 missing) Surrogate splits:

RunTime < 8.64 to the right, agree=1, adj=1, (0 split) RunPulse < 161 to the left, agree=1, adj=1, (0 split) MaxPulse < 170 to the left, agree=1, adj=1, (0 split)

Node number 4: 4 observations mean=38.86, MSE=0.7305915

RunTime < 10.9 to the right, improve=0.3967104, (0 missing) MaxPulse < 171 to the right, improve=0.2050189, (0 missing)

(8)

RestPulse < 56.5 to the right, improve=0.1734683, (0 missing) Age < 53 to the left, improve=0.1275460, (0 missing) Weight < 80.285 to the left, improve=0.1057539, (0 missing) Surrogate splits:

Weight < 71.33 to the left, agree=0.783, adj=0.375, (0 split) RestPulse < 56.5 to the right, agree=0.696, adj=0.125, (0 split) MaxPulse < 171 to the right, agree=0.696, adj=0.125, (0 split)

RunTime < 11.015 to the left, improve=0.7058204, (0 missing) Weight < 63.845 to the right, improve=0.3568667, (0 missing) RunPulse < 169 to the left, improve=0.2214416, (0 missing) Age < 49.5 to the right, improve=0.1609844, (0 missing) MaxPulse < 174 to the left, improve=0.1609844, (0 missing)

Age < 53 to the left, improve=0.31224620, (0 missing) Weight < 73.54 to the right, improve=0.28818430, (0 missing) RunPulse < 182.5 to the right, improve=0.18386270, (0 missing) MaxPulse < 182.5 to the right, improve=0.10608090, (0 missing) RestPulse < 46 to the right, improve=0.08776052, (0 missing)

Node number 20: 1 observations mean=40.836, MSE=0

Weight < 73.54 to the right, improve=0.3829283, (0 missing) Age < 50 to the right, improve=0.3319523, (0 missing) RestPulse < 46 to the right, improve=0.2220905, (0 missing) MaxPulse < 167 to the left, improve=0.1815867, (0 missing) RunTime < 9.515 to the right, improve=0.1733780, (0 missing)

RunTime < 9.515 to the right, improve=0.5395235, (0 missing)

(9)

Weight < 80.17 to the left, improve=0.4061913, (0 missing) Age < 50 to the right, improve=0.2825593, (0 missing) MaxPulse < 169 to the left, improve=0.1877790, (0 missing) RestPulse < 45.5 to the right, improve=0.1317747, (0 missing) Surrogate splits:

RunPulse < 174 to the left, agree=0.818, adj=0.333, (0 split) MaxPulse < 176 to the left, agree=0.818, adj=0.333, (0 split)

plot(fitness.rpart,uniform=TRUE,margin=0.1) text(fitness.rpart)

Zagadnienie braków danych (ang. missing values)

załóżmy, że brakuje nam informacji o zmiennej Age.Czy predykcji MNK jest możliwa?

co wtedy można zrobić? czy predykcja przy pomocy drzewa jest możliwa? Sprawdźmy...

summary(fitness.rpart)

co oznaczają zmienne 'surrogate' w węzłach?

predict(fitness.rpart,newdata=data.frame(Age=NaN, Weight=69, RunTime=NaN, RestPulse=57,RunPulse=168, MaxPulse=172))

plot(fitness.rpart,uniform=TRUE,margin=0.1) text(fitness.rpart)

(10)

OPTYMALNA WIELKOŚĆ DRZEWA

Uzasadnienie:

Dysponujemy niewielką liczbą obserwacji (n=31), a zbudowaliśmy relatywnie duże drzewo, w których mamy 16 kryteriów podziału (a więc 16 parametrów), stąd istnieje duże ryzyko przeuczenia modelu, a więc nadmiernego dopasowania do tych konkretnych danych i zgubieniu z oczu prawdziwego procesu generowania danych (ang. DGP), które za tym stoi . Ponieważ będziemy się posługiwać metodą CV, ustawmy wspólną wartość ziarna

> set.seed(289)

> fitness2.rpart <- rpart(Oxygen ~ ., data=fitness,cp=0.001, minsplit=2)

> printcp(fitness2.rpart) Regression tree:

rpart(formula = Oxygen ~ ., data = fitness, cp = 0.001, minsplit = 2)

Variables actually used in tree construction:

[1] Age RunPulse RunTime Weight Root node error: 851.38/31 = 27.464 n= 31

CP nsplit rel error xerror xstd

1 0.5139721 0 1.0000000 1.12998 0.301569 2 0.2756602 1 0.4860279 0.95172 0.204636 3 0.0686668 2 0.2103677 0.42355 0.108958 4 0.0336440 3 0.1417009 0.36754 0.096277 5 0.0230088 4 0.1080569 0.33281 0.083356 6 0.0224271 5 0.0850481 0.37549 0.085346 7 0.0185299 6 0.0626210 0.35721 0.115581 8 0.0161028 7 0.0440911 0.38512 0.117331 9 0.0083371 8 0.0279883 0.39605 0.117071 10 0.0075833 9 0.0196512 0.41225 0.119227 11 0.0033934 10 0.0120679 0.41317 0.121580 12 0.0023192 11 0.0086746 0.40981 0.121878 13 0.0020579 12 0.0063553 0.42362 0.123989 14 0.0011531 13 0.0042974 0.42088 0.122729 15 0.0010078 14 0.0031443 0.43253 0.122540 16 0.0010000 15 0.0021364 0.43253 0.122540

>

Wywołaj teraz komendę printcp printcp(fitness.rpart)

Regression tree:

rpart(formula = Oxygen ~ ., data = fitness, cp = 0.01, minsplit = 2)

Variables actually used in tree construction:

[1] Age RunTime Weight

Root node error: 851.38/31 = 27.464 n= 31

CP nsplit rel error xerror xstd

(11)

1 0.513972 0 1.000000 1.08277 0.294830 2 0.275660 1 0.486028 0.96505 0.218916 3 0.068667 2 0.210368 0.40656 0.102354 4 0.033644 3 0.141701 0.37552 0.095828 5 0.023009 4 0.108057 0.39328 0.117570 6 0.022427 5 0.085048 0.51394 0.136904 7 0.018530 6 0.062621 0.45906 0.132710 8 0.016103 7 0.044091 0.46804 0.133112 9 0.010000 8 0.027988 0.49423 0.133120

>

jaka jest interpretacja parameteru cp w odniesieniu do SSE?

Węzły (i towarzyszące im splity) są uporządkowane zgodnie z malejącym "potencjałem" redukcji relatywnego SSE, a więc z malejącą wartością parametru regularyzacji cp. co oznacza kolumna rel error? ile wynosi redukcja rel terror między sąsiednimi wierszami? czy zależy nam na minimalizacji rel error? co oznacza kolumna xerror? jaka jest optymalna wartość parametru regularyzacji CP ? co oznacza kolumna xstd?

> plotcp(fitness2.rpart)

co jest na osiach? co oznacza pojedyńczy punkt na wykresie? co oznaczają pionowe "wąsy" przy punktach? co oznacza pionowa przerywana kreska? jaka jest jego wartość? jaka jest optymalna wartość CP przy zastosowaniu reguły jednego odchylenia standardowego 1SE? optymalne CP zgodnie z 1SE: jest to taka maksymalna wartość CP, której błąd nie jest większy niż minimalny błąd plus jego odchylenie standardowe.

Budujemy optymalne drzewo:

fitness3.rpart <- prune(fitness2.rpart, cp=0.03)

> print(fitness3.rpart) n= 31

node), split, n, deviance, yval

(12)

* denotes terminal node

1) root 31 851.381500 47.37581

2) RunTime>=8.935 27 384.980400 45.92970 4) RunTime>=12.29 4 2.922366 38.86000 * 5) RunTime< 12.29 23 147.366000 47.15922 10) RunTime>=10.9 8 27.753940 44.97612 * 11) RunTime< 10.9 15 61.150420 48.32353 * 3) RunTime< 8.935 4 28.814820 57.13700 6) Age>=43 2 0.053792 54.46100 * 7) Age< 43 2 0.117128 59.81300 *

>

posługujemy się cp niewiele mniejszą niż optymalne cp, czy to dobrze? co by było gdybyśmy użyli cp=0.0337 ? Porównanie różnych modeli na zbiorze treningowym:

Porównanie różnych modeli na zbiorze treningowym:

(13)

Graficzna reprezentacja predykcji drzewa regresyjnego

Aby móc to zilustrować, posłużymy się tylko dwie "najważniejszymi" zmiennymi

> fitness5.rpart <- rpart(Oxygen ~ RunTime + Age, data=fitness,cp=0.00001, minsplit=2)

> print(fitness5.rpart) n= 31

node), split, n, deviance, yval * denotes terminal node

1) root 31 851.38150000 47.37581

2) RunTime>=8.935 27 384.98040000 45.92970 4) RunTime>=12.29 4 2.92236600 38.86000 8) RunTime>=13.555 1 0.00000000 37.38800 * 9) RunTime< 13.555 3 0.03332067 39.35067 18) Age>=49 2 0.02080800 39.30500

36) RunTime>=12.755 1 0.00000000 39.20300 * 37) RunTime< 12.755 1 0.00000000 39.40700 * 19) Age< 49 1 0.00000000 39.44200 *

5) RunTime< 12.29 23 147.36600000 47.15922 10) RunTime>=10.9 8 27.75394000 44.97612 20) RunTime< 11.015 1 0.00000000 40.83600 * 21) RunTime>=11.015 7 8.16464200 45.56757 42) Age< 47.5 4 0.70762270 44.96375 84) RunTime< 11.79 3 0.02169267 44.72467 168) Age< 44.5 1 0.00000000 44.60900 * 169) Age>=44.5 2 0.00162450 44.78250 * 85) RunTime>=11.79 1 0.00000000 45.68100 * 43) Age>=47.5 3 4.05408300 46.37267

86) RunTime< 11.335 2 0.46272200 45.59900 172) RunTime< 11.125 1 0.00000000 45.11800 * 173) RunTime>=11.125 1 0.00000000 46.08000 * 87) RunTime>=11.335 1 0.00000000 47.92000 * 11) RunTime< 10.9 15 61.15042000 48.32353 22) Age< 53 13 41.19838000 47.88100 44) Age>=50 4 2.49102900 46.34250

88) RunTime< 10.485 3 0.80502870 45.96767 176) RunTime< 9.815 1 0.00000000 45.44100 * 177) RunTime>=9.815 2 0.38896200 46.23100 354) RunTime>=10.235 1 0.00000000 45.79000 * 355) RunTime< 10.235 1 0.00000000 46.67200 * 89) RunTime>=10.485 1 0.00000000 47.46700 * 45) Age< 50 9 25.03146000 48.56478

90) RunTime>=10.19 3 2.97420500 47.71267 180) RunTime< 10.725 2 0.12450050 47.02350 360) RunTime< 10.425 1 0.00000000 46.77400 * 361) RunTime>=10.425 1 0.00000000 47.27300 * 181) RunTime>=10.725 1 0.00000000 49.09100 * 91) RunTime< 10.19 6 18.78983000 48.99083 182) RunTime< 10.075 4 12.26305000 48.25400 364) RunTime>=9.735 1 0.00000000 45.31300 * 365) RunTime< 9.735 3 0.73040470 49.23433 730) Age>=43.5 2 0.11664450 48.91450

1460) RunTime>=9.175 1 0.00000000 48.67300 * 1461) RunTime< 9.175 1 0.00000000 49.15600 * 731) Age< 43.5 1 0.00000000 49.87400 *

183) RunTime>=10.075 2 0.01170450 50.46450

(14)

366) RunTime< 10.105 1 0.00000000 50.38800 * 367) RunTime>=10.105 1 0.00000000 50.54100 * 23) Age>=53 2 0.85805000 51.20000

46) RunTime< 10.13 1 0.00000000 50.54500 * 47) RunTime>=10.13 1 0.00000000 51.85500 * 3) RunTime< 8.935 4 28.81482000 57.13700 6) RunTime>=8.64 2 0.05379200 54.46100 12) RunTime< 8.785 1 0.00000000 54.29700 * 13) RunTime>=8.785 1 0.00000000 54.62500 * 7) RunTime< 8.64 2 0.11712800 59.81300 14) RunTime< 8.4 1 0.00000000 59.57100 * 15) RunTime>=8.4 1 0.00000000 60.05500 *

plot(fitness5.rpart, uniform=TRUE, margin=0.1) text(fitness5.rpart)

(15)

Obszary różnicące się predykcją

plot(fitness$RunTime, fitness$Age, xlim=c(8.0,14.0), ylim=c(37,58), xlab="RunTime", ylab="Age") lines(c(8.935,8.935), c(37,58))

lines(c(12.29,12.29), c(37,58)) lines(c(8.64,8.64), c(37,58)) lines(c(12.29,12.29), c(37,58)) lines(c(10.9,10.9), c(37,58)) lines(c(11.015,11.015), c(37,58)) lines(c(8.935,10.9), c(53,53))

Aby wykres był do końca czytelny, powinniśmy jeszcze w pola wprowadzić wartość teoretyczną poszczególnego liścia len=60

xp <- seq(8.0, 15.0, length=len) yp <- seq(35, 60, length=len)

(16)

siatka <- expand.grid(RunTime=xp, Age=yp) zp <- predict(fitness5.rpart, siatka)

persp(xp, yp, matrix(zp, len), theta = 150, phi = 30, col = "lightblue", zlim=c(37,60), xlab = "RunTime", ylab = "Age", zlab =

"Oxygen")

# Wykres 3D

Drzewa klasyfikacyjne: dane Cars93

library(MASS) data(Cars93) dim(Cars93) names(Cars93)

Stosując funkcje dim oraz names poznaj licznościowo zbiór Cars93 I nazwy kolumn opisujących dane.

Będziemy modelować zmienną kategoryczną "typ samochodu"

table(Cars93$Type)

(17)

Zechcemy połączyć pewne kategorie samochodów do tego samego "koszyka", w ten sposób aby mieć cztery kategorie:

# D - duży: Large, Van

# SR - średni: Midsize, Compact

# M - mały: small

# SP - sportowy: sporty

typ <- ifelse(Cars93$Type=="Large" | Cars93$Type=="Van", "D", "SR") typ2 <- ifelse(Cars93$Type=="Small", "M", typ)

typ3 <- ifelse(Cars93$Type=="Sporty", "SP", typ2)

typ3 <- ifelse(Cars93$Type=="Sporty", "SP", ifelse(Cars93$Type=="Small","M", ifelse(Cars93$Type=="Large" | Cars93$Type=="Van", "D", "SR")))

table(typ3)

Cars93 = data.frame(Cars93, typ=typ3) dim(Cars93)

names(Cars93)

ponieważ posłużymy się CV użyjmy wspólnego ziarna set.seed(222)

Z wielu potencjalnych zmiennych objaśniających wybierzmy tylko 5

cars.tree <- rpart(typ ~ Length + Weight + EngineSize + Horsepower + RPM, data=Cars93, cp=0.0001, minsplit=5) niewielka wartość cp i minsplit, powoduje, że mamy "duże" drzewo

(18)

par(mfrow=c(1,1), mar=c(0,0,0,0), pty="m")

plot(cars.tree,compress=T,uniform=TRUE, branch=0.4, margin=0.2) text(cars.tree, use.n=T, cex=0.7)

co oznaczają liczby w liściach?

Wydruk tekstowy drzewa?

print(summary(cars.tree), digits=4)

(19)

co oznacza wartość loss?

co oznacza wartość loss przy jednostkowej macierzy kosztu?

co oznacza yprob? co oznacza yval? i jaką przyjęto regułę decyzyjną przy jego wyborze?

(20)

Optymalne drzewo klasyfikacyjne

printcp(cars.tree)

jakie jest optymalne drzewo, kierując się minimalizacją błędu kroswalidacyjnego? jakie jest optymalne drzewo, kierując się regułą 1SE?

plotcp(cars.tree)

Z czego wynika brak monotoniczności tej funkcji błędu? Jak jest uzasadniania reguła 1SE ?

Wizualizacja drzewa klasyfikacyjnego

w tym celu wybierzemy tylko dwie zmienne

par(mfrow=c(1,1), mar=c(0,0,0,0), pty="m")

plot(cars2zmienne.tree,compress=T,uniform=TRUE,branch=0.4,margin=0.2) text(cars2zmienne.tree, use.n=TRUE, cex=0.7)

print(summary(cars2zmienne.tree), digits=4)

(21)

Rysujemy obszary klasyfikacji:

plot(c(1,5.70), c(1695,4105), type="p", xlab="Weight", ylab="EngineSize") text(Cars93$EngineSize, Cars93$Weight, Cars93$typ, cex=0.6)

lines(c(1.95,1.95), c(1695,4105)) lines(c(1.95,5.70), c(3460,3460)) lines(c(3.2,3.2), c(1695,3460)) lines(c(1.95,5.70), c(3700,3700)) lines(c(3.25,3.25), c(3460,3700)) lines(c(3.25,5.70), c(3630,3630))

(22)

Metody łączenia klasyfikatorów

Uzasadnienie: Załóżmy, że dysponujemy n (np. n=1000) niezależnymi klasyfikatorami. Każdy z nich charakteryzuje się prawdopodobieństwem p (np. p = 60%) poprawnej klasyfikacji.

Załóżmy bardzo słaby klasyfikator p=0.51

dysponujemy 1000 niezależnych klasyfikatorów n=1000

Ile klasyfiaktorów na 1000 wskaże poprawną odpowiedź rbinom(1,n,p)

rbinom(1,n,p) rbinom(1,n,p) rbinom(1,n,p) rbinom(1,n,p)

Jeśli będziemy się posługiwali głosowaniem, tj. klasyfikowali do tej kategorie,do której większość klasyfikatorów z 1000 klasyfikowało, to jak często będziemy poprawnie klasyfikować?

rbinom(1,n,p)>500 rbinom(1,n,p)>500 rbinom(1,n,p)>500 rbinom(1,n,p)>500 rbinom(1,n,p)>500

prop.table(table(rbinom(1000000,n,p)>500))

?rbinom

Wniosek: z poprawności 51% awansowaliśmy na 72%