calif <- read.table("http://www.stat.cmu.edu/~cshalizi/350/hw/06/cadata.dat", header = TRUE) ### California Housing Dataset price.deciles <- quantile(calif$MedianHouseValue, 0:10/10) cut.prices <- cut(calif$MedianHouseValue, price.deciles, include.lowest = TRUE) plot(calif$Longitude, calif$Latitude, col = grey(10:2/11)[cut.prices], pch = 20, xlab = "Longitude", ylab = "Latitude") ## gusca podrucja s visim cijenama su gradovi (LA, San Francisco...) require(tree) treefit <- tree(log(MedianHouseValue) ~ Longitude + Latitude, data = calif) plot(treefit) text(treefit, cex = 0.75) ## u listovima pisu prosjecne (log) vrijednosti u odgovarajucem elementu particije ## za novu kombinaciju (Long, Lat) jednostavno se spustamo po stablu sve dok ne dodemo do lista # te je vrijednost u tom listu upravo predikcija koju daje nas model ## velicina bridova odgovara smanjenju troska C(T) (vidi (7.4)) # buduci da imamo samo dva prediktora, mozemo prikazati rezultirajucu particiju prostora prediktora cut.prices <- cut(calif$MedianHouseValue, price.deciles, include.lowest = TRUE) plot(calif$Longitude, calif$Latitude, col = grey(10:2/11)[cut.prices], pch = 20, xlab = "Longitude", ylab = "Latitude") partition.tree(treefit, ordvars = c("Longitude", "Latitude"), add = TRUE) summary(treefit) ### Residual mean deviance (0.1662) je skoro trening MSE (C(T) u (7.3)) -- RSS se dijeli # s (n-broj listova) umjesto n ?tree ## kriterij zaustavljanja mozemo mijenjati koristeci parametre minsize, mincut i mindev ?tree.control ## smanjimo mindev -> dakle, napravit cemo neke podjele koje prvi model ne bi napravio treefit2 <- tree(log(MedianHouseValue) ~ Longitude + Latitude, data = calif, mindev = 0.001) plot(treefit2) text(treefit2, cex = 0.75) summary(treefit2) ## 68 listova cut.prices <- cut(calif$MedianHouseValue, price.deciles, include.lowest = TRUE) plot(calif$Longitude, calif$Latitude, col = grey(10:2/11)[cut.prices], pch = 20, xlab = "Longitude", ylab = "Latitude") partition.tree(treefit2, ordvars = c("Longitude", "Latitude"), add = TRUE) ## puno finija podjela u gradovima ## ukljucimo i ostale prediktore treefit3 <- tree(log(MedianHouseValue) ~ ., data = calif) plot(treefit3) text(treefit3, cex = 0.75) ## iako dopustamo samo binarne podjele (x_j < ili > od s), slozenije podjele mozemo dobiti # koristeci nekoliko binarnih -> vidi MedianIncome (cini se da je to najznacajniji prediktor) ## stabla zapravo traze (moguce) komplicirane interakcije medu varijablama -> npr. Lat i Long se koriste # samo ako je medianIncome<3.54, ali drugacije ovisno o tome je li < ili > od 2.51 itd. (npr. kod lin. regresije ili varijanti # te interekcije ce se tesko prepoznati ako jer moramo ukljuciti puno varijabli) summary(treefit3) ## 15 listova, te koristimo samo 4 od 8 prediktora treefit3 ## detalji oko izbora podjela (deviance je suma kvadrata reziduala samo u tom cvoru) ## ne mozemo crtati dobivenu particiju, ali mozemo prikazati dobivene procjene kao sto smo prikazali stvarne cijene cut.predictions <- cut(predict(treefit3), log(price.deciles), include.lowest = TRUE) plot(calif$Longitude, calif$Latitude, col = grey(10:2/11)[cut.predictions], pch = 20, xlab = "Longitude", ylab = "Latitude") cut.prices <- cut(calif$MedianHouseValue, price.deciles, include.lowest = TRUE) plot(calif$Longitude, calif$Latitude, col = grey(10:2/11)[cut.prices], pch = 20, xlab = "Longitude", ylab = "Latitude") ## obrezivanje i CV ?cv.tree ?prune.tree ## daje niz podstabala iz (7.9) pr.tree=prune.tree(treefit3) pr.tree ## size su njihove velicine, dev je RSS, a k je zapravo alfa pr.tree=prune.tree(treefit3, best=8) ## probaj 8 plot(pr.tree) text(pr.tree) summary(pr.tree) cv.calif=cv.tree(treefit3) names(cv.calif) plot(cv.calif, type="b") ## u ovom slucaju, najbolji izgleda puni model ### pogledajmo kakva je cv procjena testne greske u odnosu na stvarnu set.seed(1021) n=20640 m=10000 ## velicina skupa za trening train=sample(1:n, m) treefit4 <- tree(log(MedianHouseValue) ~ ., data = calif, subset=train) summary(treefit4) cv.calif=cv.tree(treefit4) plot(cv.calif, type="b") which.min(cv.calif$dev) y.test=log(calif$MedianHouseValue[-train]) test.er=numeric(14) for (i in 2:14) { tr=prune.tree(treefit4, best = i) yhat=predict(tr, newdata=calif[-train, ]) test.er[i]=mean((yhat-y.test)^2) } plot(cv.calif$size, cv.calif$dev/m, type="b") ## dev/m=RSS/m=MSE points(2:14, test.er[-1], type="b", col="blue") ## zbog ogromog broja podataka, tesko cemo imati overfitting ## testna greska punog modela (sa 14 cvorova) test.er[14] ## [1] 0.1404545 (gledamo logaritme cijena pa je stoga greska manjeg reda) ######## Random forrest library(randomForest) rf.cal<- randomForest(log(MedianHouseValue) ~ ., data = calif, subset = train, importance = TRUE, ntree=1000) ## treba malo duze vremena da se izvrti rf.cal # m=2 plot(1:1000, rf.cal$mse, type="l") yhat=predict(rf.cal, newdata=calif[-train, ]) mean((yhat-y.test)^2) ## [1] 0.062861 -- puno bolje nego jedno stablo varImpPlot(rf.cal, type=2) ## slicno kao za Boston podatke ## jedno stablo: summary(treefit4) plot(treefit4) text(treefit4, cex = 0.75)