Dati TL
Importo i dati
Importo i dati e trasformo le variabili Face e
Income.
tl = read.csv("dataset/TL.csv") %>% tibble %>% dplyr::select(-X)
tl = tl %>% mutate_at(c( "MARSTAT"), as.factor)
colLog = c("FACE","INCOME")
tl = tl %>% mutate_at(colLog,log) %>% rename_at(colLog, function(c) paste("L",c,sep = ""))
tlCreo una partizione del dataset, l’80% dei dati saranno destinati al training del modello mentre il restante 20% verrà utilizzato per la valutazione del modello creato.
set.seed(1)
train = createDataPartition(tl$AGE,p = .7,list = F)Bagging
set.seed(1)
modBagg = randomForest(LFACE ~ .,
data = tl,
subset = train,
mtry = length(tl) - 1,
importance = T)
modBagg##
## Call:
## randomForest(formula = LFACE ~ ., data = tl, mtry = length(tl) - 1, importance = T, subset = train)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 5
##
## Mean of squared residuals: 2.311332
## % Var explained: 29.48Test
Nel grafico si possono osservare come sono distribuiti i residui
secondo la variabile risposta del modello. Essi non sono distribuiti
omogeneamente, il modello sovrastima i dati quando LFACE è
minore, mentre li sottostima quando tendono a essere sopra la media.
yPred = predict(modBagg, newdata = tl[-train,])
ggplot(mapping = aes(x = tl$LFACE[-train], y = yPred-tl$LFACE[-train])) +
geom_point() +
geom_vline(aes(xintercept = mean(tl$LFACE[-train]))) +
geom_hline(aes(yintercept = 0)) +
xlab("LFACE") +
ylab("Residui")
La media quadratica dei residui è pari a 2.76.
mean((yPred-tl$LFACE[-train])**2)## [1] 2.703698Random forest
Le variabili da usare durante la creazione di una foresta casuale sono state impostate a 3 mentre prima, nel bagging, venivano considerate tutte e 5. In questo caso la varianza spiegata è migliorata di quasi un punto percentuale.
set.seed(1)
modRFor = randomForest(LFACE ~ .,
data = tl,
subset = train,
mtry = 3,
importance = T)
modRFor ##
## Call:
## randomForest(formula = LFACE ~ ., data = tl, mtry = 3, importance = T, subset = train)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 3
##
## Mean of squared residuals: 2.247235
## % Var explained: 31.44Test
Anche in questo caso si ripresenta il problema precedente con i residui.
yPred = predict(modRFor, newdata = tl[-train,])
ggplot(mapping = aes(x = tl$LFACE[-train], y = yPred-tl$LFACE[-train])) +
geom_point() +
geom_vline(aes(xintercept = mean(tl$LFACE[-train]))) +
geom_hline(aes(yintercept = 0)) +
xlab("LFACE") +
ylab("Residui")
La media quadratica dei residui è pari a 2.61, leggermente minore a quella del bagging, ciò dimostra che il modello random forest migliore anche durante la validation con i dati di test.
mean((yPred-tl$LFACE[-train])**2)## [1] 2.618849Nel seguente grafico sono riportate le variabili con la misura in percentuale di MSE che diminuisce quando viene inclusa quella variabile e con la misura della diminuzione totale dell’impurità del nodo.
varImpPlot(modRFor)
Boosting
Per il boosting degli aberi sui dati è stato richiesto di produrre 1000 alberi con una profondità massima di 10.
modBoos = gbm(LFACE ~ .,
data = tl,
distribution = "gaussian",
n.trees = 1000,
interaction.depth = 10)
modBoos %>% summary
Nel seguente grafico è presente l’effetto marginale della variabile
LINCOME.
plot(modBoos, i = "LINCOME")
Test
Questo grafico è leggermente migliore confronto ai precedenti sulla distribuzione dei residui, inoltre se si osserva la grandezza di misura dell’asse delle ordinate si noterà quanto è diminuita confronto a prima.
yPred = predict(modBoos, newdata = tl[-train,])## Using 1000 trees...ggplot(mapping = aes(x = tl$LFACE[-train], y = yPred-tl$LFACE[-train])) +
geom_point() +
geom_vline(aes(xintercept = mean(tl$LFACE[-train]))) +
geom_hline(aes(yintercept = 0)) +
xlab("LFACE") +
ylab("Residui")
La media quadratica dei residui è pari a 0.10, all’incirca il 4% del MSE del modello random forest. Il boosting si rivela, in questo caso, il metodo migliore.
mean((yPred-tl$LFACE[-train])**2)## [1] 0.1005722Dati prostata
Importo i dati
pr = read.csv("dataset/prostate.csv") %>%
as_tibble %>% dplyr::select(-X)
pr$svi = pr$svi %>% as.logical()
prCreo una partizione del dataset, l’80% dei dati saranno destinati al training del modello mentre il restante 20% verrà utilizzato per la valutazione del modello creato.
set.seed(1)
train = createDataPartition(pr$lcavol,p = .7,list = F)Bagging
set.seed(1)
modBagg = randomForest(lpsa ~ .,
data = pr,
subset = train,
mtry = length(pr) - 1,
importance = T)
modBagg##
## Call:
## randomForest(formula = lpsa ~ ., data = pr, mtry = length(pr) - 1, importance = T, subset = train)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 9
##
## Mean of squared residuals: 0.7474968
## % Var explained: 47.98Test
Nel grafico si possono osservare come sono distribuiti i residui secondo la variabile risposta del modello.
yPred = predict(modBagg, newdata = pr[-train,])
ggplot(mapping = aes(x = pr$lpsa[-train], y = yPred-pr$lpsa[-train])) +
geom_point() +
geom_vline(aes(xintercept = mean(pr$lpsa[-train]))) +
geom_hline(aes(yintercept = 0)) +
xlab("lpsa") +
ylab("Residui")
La media quadratica dei residui è pari a 0.177.
mean((yPred-pr$lpsa[-train])**2)## [1] 0.1770565Random forest
Le variabili da usare durante la creazione di una foresta casuale sono state impostate a 4, provando dove con 4 si ottiene un modello più performante che con le altre opzioni.
set.seed(1)
modRFor = randomForest(lpsa ~ .,
data = pr,
subset = train,
mtry = 4,
importance = T)
modRFor##
## Call:
## randomForest(formula = lpsa ~ ., data = pr, mtry = 4, importance = T, subset = train)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 4
##
## Mean of squared residuals: 0.6910323
## % Var explained: 51.91Test
yPred = predict(modRFor, newdata = pr[-train,])
ggplot(mapping = aes(x = pr$lpsa[-train], y = yPred-pr$lpsa[-train])) +
geom_point() +
geom_vline(aes(xintercept = mean(pr$lpsa[-train]))) +
geom_hline(aes(yintercept = 0)) +
xlab("lpsa") +
ylab("Residui")
La media quadratica dei residui è pari a 0.198, più alta del 12% rispetto al bagging.
mean((yPred-pr$lpsa[-train])**2)## [1] 0.1981521Nel seguente grafico sono riportate le variabili con la misura in percentuale di MSE che diminuisce quando viene inclusa quella variabile e con la misura della diminuzione totale dell’impurità del nodo.
varImpPlot(modRFor)
Boosting
Per il boosting degli aberi sui dati è stato richiesto di produrre 1000 alberi con una profondità massima di 10.
modBoos = gbm(lpsa ~ .,
data = pr %>% mutate_if(is.logical, as.factor),
distribution = "gaussian",
n.trees = 1000,
interaction.depth = 10)
modBoos %>% summary
Nel seguente grafico è presente l’effetto marginale della variabile
lcavol.
plot(modBoos, i = "lcavol")
Test
yPred = predict(modBoos, newdata = pr[-train,])## Using 1000 trees...ggplot(mapping = aes(x = pr$lpsa[-train], y = yPred-pr$lpsa[-train])) +
geom_point() +
geom_vline(aes(xintercept = mean(pr$lpsa[-train]))) +
geom_hline(aes(yintercept = 0)) +
xlab("lpsa") +
ylab("Residui")
La media quadratica dei residui è pari a 0.015, all’incirca l’8% del MSE del modello bagging. Anche in questo caso il boosting si rivela il metodo migliore.
mean((yPred-pr$lpsa[-train])**2)## [1] 0.01518233Caret
Ripeto la creazione del modello random forest, questa volta tramite la libreria caret.
set.seed(1)
modCaret = train(lpsa ~ .,
data = pr,
method = "rf",
importance = T,
trControl = trainControl("oob")
)
modCaret ## Random Forest
##
## 97 samples
## 9 predictor
##
## No pre-processing
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared
## 2 0.7877793 0.5294018
## 5 0.7786405 0.5402569
## 9 0.8046302 0.5090539
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 5.yPred = predict(modCaret, newdata = pr[-train,])
ggplot(mapping = aes(x = pr$lpsa[-train], y = yPred-pr$lpsa[-train])) +
geom_point() +
geom_vline(aes(xintercept = mean(pr$lpsa[-train]))) +
geom_hline(aes(yintercept = 0)) +
xlab("lpsa") +
ylab("Residui")
### Test
La media quadratica dei residui è pari a 0.107, il 45% più bassa rispetto l’utilizzo della libreria randomForest.
mean((yPred-pr$lpsa[-train])**2)## [1] 0.1155718Controllo dei parametri
Ora ripeto la procedura ma questa volta alcuni parametri li cambio:
mtry corrisponde al numero di variabili casuali che il
modello potrà scegliere a ogni nodo tra quelle disponibili, in
trainControlo number sta per il numero di
fold nella quale dovrà suddividere i dati e
repeats è il numero di fold completi da
ripetere.
set.seed(1)
modCaretNew = train(lpsa~.,
data = pr,
method = 'rf',
metric = 'RMSE',
importance = T,
tuneGrid = expand.grid(.mtry = 7), #round(sqrt(ncol(pr)))),
trControl = trainControl(method = 'repeatedcv',
number = 10,
repeats = 3))
modCaretNew## Random Forest
##
## 97 samples
## 9 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 86, 89, 88, 86, 88, 89, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.7493033 0.6266579 0.6121985
##
## Tuning parameter 'mtry' was held constant at a value of 7La media quadratica dei residui è pari a 0.097, il 10% più bassa rispetto al precedente modello.
yPred = predict(modCaretNew, newdata = pr[-train,])
mean((yPred-pr$lpsa[-train])**2)## [1] 0.09751975