Laboratorio 5
Erik De Luca
2022-11-17
knitr::opts_chunk$set(echo = TRUE)
library(readr)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## âś” ggplot2 3.4.0 âś” dplyr 1.0.9
## âś” tibble 3.1.7 âś” stringr 1.4.1
## âś” tidyr 1.2.1 âś” forcats 0.5.2
## âś” purrr 0.3.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## âś– dplyr::filter() masks stats::filter()
## âś– dplyr::lag() masks stats::lag()library(dplyr)
library(tree)
library(caret)## Caricamento del pacchetto richiesto: lattice
##
## Caricamento pacchetto: 'caret'
##
## Il seguente oggetto è mascherato da 'package:purrr':
##
## liftlibrary(rpart)
library(rpart.plot)
library(glmnet)## Caricamento del pacchetto richiesto: Matrix
##
## Caricamento pacchetto: 'Matrix'
##
## I seguenti oggetti sono mascherati da 'package:tidyr':
##
## expand, pack, unpack
##
## Loaded glmnet 4.1-4library(ggplot2)
library(coefplot)
library(earth)## Caricamento del pacchetto richiesto: Formula
## Caricamento del pacchetto richiesto: plotmo
## Caricamento del pacchetto richiesto: plotrix
## Caricamento del pacchetto richiesto: TeachingDemosImporto i dati e trasformo le variabili Face e
Income.
tl = read.csv("dataset/TL.csv") %>% tibble %>% 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)Albero
Nella raffigurazione dell’albero si nota come la variabile più
influente e importe è LINCOME, seguita da AGE.
Essa, infatti, non è solo quella che compare sul primo nodo ma è anche
la piĂą utilizzata nei nodi.
I nodi terminali sono 10 e sono stati scelti automaticamente dalla
funzione tree. In alternativa, potevano essere dati in
input il numero di nodi foglia desiderati, da un minimo di 1 (soluzione
banale), fino ad un massimo di \(n\)
nodi foglia dove \(n\) è il numero di
osservazioni utlizzate per allenare il modello. Ovviamente entrambe le
solouzioni non vanno bene, la prima perché non generebbe nessun
risultato utile in quanto tutte le osservazioni del dataset
apprterebbero alla medesima regione dello spazio in cui si è divisa la
partizione dei dati. Nel secondo caso si andrebbe incontro
all’overfitting massimo.
albero = tree(LFACE ~ ., data = tl[train,])
{
plot(albero)
text(albero, cex = 0.7)
}
Confronto l’albero con altri modelli
Regressione lineare multipla
Per scegliere il modello migliore di regressione lineare utilizzo il
criterio dell’informazione di Akaike (AIC), che va penalizzare i modelli
con troppe variabili e poche informazioni aggiuntive, come \(R_{adj}^2\). Tramite la funzione
stepAIC andiamo a selezionare il modello con il minore
AIC.
modLm = MASS::stepAIC(lm(LFACE ~ ., data = tl[train,]))## Start: AIC=140.03
## LFACE ~ AGE + MARSTAT + EDUCATION + NUMHH + LINCOME
##
## Df Sum of Sq RSS AIC
## - MARSTAT 1 2.434 378.42 139.29
## - AGE 1 3.535 379.52 139.85
## <none> 375.99 140.03
## - NUMHH 1 19.560 395.55 147.92
## - EDUCATION 1 27.044 403.03 151.57
## - LINCOME 1 77.135 453.12 174.42
##
## Step: AIC=139.29
## LFACE ~ AGE + EDUCATION + NUMHH + LINCOME
##
## Df Sum of Sq RSS AIC
## - AGE 1 2.847 381.27 138.75
## <none> 378.42 139.29
## - EDUCATION 1 25.846 404.27 150.17
## - NUMHH 1 32.756 411.18 153.47
## - LINCOME 1 87.505 465.93 177.85
##
## Step: AIC=138.75
## LFACE ~ EDUCATION + NUMHH + LINCOME
##
## Df Sum of Sq RSS AIC
## <none> 381.27 138.75
## - EDUCATION 1 25.671 406.94 149.46
## - NUMHH 1 42.432 423.70 157.33
## - LINCOME 1 84.659 465.93 175.85modLm %>% summary##
## Call:
## lm(formula = LFACE ~ EDUCATION + NUMHH + LINCOME, data = tl[train,
## ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.8398 -0.7949 0.1235 0.9196 5.2222
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.43477 1.01113 1.419 0.157534
## EDUCATION 0.15547 0.04335 3.586 0.000426 ***
## NUMHH 0.32051 0.06952 4.611 7.34e-06 ***
## LINCOME 0.65302 0.10027 6.512 6.38e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.413 on 191 degrees of freedom
## Multiple R-squared: 0.4035, Adjusted R-squared: 0.3941
## F-statistic: 43.06 on 3 and 191 DF, p-value: < 2.2e-16Lasso
modLasso = glmnet(x = makeX(tl[train,names(tl)!= "LFACE"]),
y = tl$LFACE[train],
alpha = 1)
cvModLasso = cv.glmnet(makeX(tl[train,names(tl)!= "LFACE"]),
tl$LFACE[train],
alpha = 1)
{
plot(cvModLasso)
abline(v = log(cvModLasso$lambda.min), col = "blue", lty = 2)
}
modLasso = glmnet(x = makeX(tl[train,names(tl)!= "LFACE"]),
y = tl$LFACE[train],
alpha = 1,
lambda = cvModLasso$lambda.min,
standardize = TRUE,
family = "gaussian")
coefplot(modLasso)
coef.glmnet(modLasso, s = cvModLasso$lambda.min)## 7 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 1.94394349
## AGE -0.01146647
## MARSTATnot single 0.30992159
## MARSTATsingle .
## EDUCATION 0.15876524
## NUMHH 0.25386762
## LINCOME 0.64806305Ridge
modRidge = glmnet(x = makeX(tl[train,names(tl)!= "LFACE"]),
y = tl$LFACE[train],
alpha = 0)
cvModRidge = cv.glmnet(makeX(tl[train,names(tl)!= "LFACE"]),
tl$LFACE[train],
alpha = 0)
{
plot(cvModRidge)
abline(v = log(cvModRidge$lambda.min), col = "blue", lty = 2)
}
modRidge = glmnet(x = makeX(tl[train,names(tl)!= "LFACE"]),
y = tl$LFACE[train],
alpha = 0,
lambda = cvModRidge$lambda.min,
standardize = TRUE,
family = "gaussian")
coefplot(modRidge)
Confronto modelli
Come deducibile la variabile LINCOME è la più influente
su tutti gli modelli. Interessante notare come i modelli ridge e lasso
contengano la variabile MARSTAT mentre il modello di
regressione lineare no, oltretutto il coefficiente è elvato ma
giustificatodal fatto che la variabile può assumere solo valori 0 e
1.
Potatura
Produco un albero di regressione e stampo il grafico delle variabili piĂą influenti nel modello.
alberoRpart = rpart::rpart(LFACE ~ .,
data = tl,
subset = train,
method = "anova",
control = rpart.control(minsplit = 10, cp = 0.003))
vip::vip(alberoRpart)
Di seguito viene raffigurato il grafico ad albero con un elevato numero
di nodi. Successivamente è presente un grafico rappresentate l’errore
relativo dell’albero al variare del cp (cost complexity
pruning) e si noti come inizialmente la curva decresce per poi
tornare a crescere.
rpart.plot(alberoRpart)
{
plotcp(alberoRpart)
points(x = which.min(alberoRpart$cptable[,"xerror"]),
y = min(alberoRpart$cptable[,"xerror"]),
col = "red",
pch = 18)
}
Il cp ottimale risulta essere 0.025. Stampo, allora, il nuovo albero potato. In questo caso l’albero avrà solamente 4 nodi foglia e verranno utilizzate 2 variabili.
rpart.plot(rpart(LFACE ~ .,
data = tl,
subset = train,
method = "anova",
control = rpart.control(cp = alberoRpart$cptable[which.min(alberoRpart$cptable[,"xerror"]),"CP"])))
Potatura tramite cross-validation
Ora eseguo la stessa operazione ma utilizzando la cross-validation come metodo per potare l’albero.
set.seed(1)
alberoTree = tree(LFACE ~ .,
data = tl,
subset = train,
control = tree.control(nobs = length(train),
minsize = 5))
cvAlbero = cv.tree(alberoTree)
cv = data.frame(size = cvAlbero$size, dev = cvAlbero$dev)
ggplot(cv, mapping = aes(x = size,y = dev), size = 10) +
geom_line() +
geom_vline(mapping = aes(xintercept = which.min(dev)), color = "red", linetype = 2)
L’albero migliore risulta essere quello con 10 nodi terminali.
alberoPotato = prune.tree(alberoTree, best = which.min(cv$dev))
{
plot(alberoPotato)
text(alberoPotato)
}
MSE
set.seed(1)
calcolaMSE = function(nodiFinali, dati, varDipendente, subsTrain, subsTest)
{
t = prune.tree(tree(get(varDipendente) ~ .,
data = dati,
subset = subsTrain,
control = tree.control(nobs = length(subsTrain),minsize = nodiFinali)),
best = nodiFinali)
return(sqrt(mean((unlist(dati[subsTest,colnames(dati) == varDipendente] - predict(t,newdata = dati[subsTest,])))^2)))
}
mseTrain = sapply(c(2:10), function(X) calcolaMSE(X,dati = tl,varDipendente = "LFACE",subsTrain = train, subsTest = train))## Warning in prune.tree(tree(get(varDipendente) ~ ., data = dati, subset =
## subsTrain, : best is bigger than tree size
## Warning in prune.tree(tree(get(varDipendente) ~ ., data = dati, subset =
## subsTrain, : best is bigger than tree sizemseTest = sapply(c(2:10), function(X) calcolaMSE(X,dati = tl,varDipendente = "LFACE",subsTrain = c(1:length(tl$AGE))[-train], subsTest = c(1:length(tl$AGE))[-train]))## Warning in prune.tree(tree(get(varDipendente) ~ ., data = dati, subset =
## subsTrain, : best is bigger than tree size
## Warning in prune.tree(tree(get(varDipendente) ~ ., data = dati, subset =
## subsTrain, : best is bigger than tree sizemseCV = sapply(c(2:10), function(X) calcolaMSE(X,dati = tl,varDipendente = "LFACE",subsTrain = c(1:length(tl$AGE))[train], subsTest = c(1:length(tl$AGE))[-train]))## Warning in prune.tree(tree(get(varDipendente) ~ ., data = dati, subset =
## subsTrain, : best is bigger than tree size
## Warning in prune.tree(tree(get(varDipendente) ~ ., data = dati, subset =
## subsTrain, : best is bigger than tree sizeNel grafico seguente sono rappresentati gli errori quadratici medi sugli alberi con diversi nodi e diversi dati. La curva train corrisponde ai modelli addestrati sui dati train e validati con i medesimi dati, infatti è la curva con l’MSE minore delle altre. La curva test, invece, corrisponde ai modelli addesstrati sui dati test e validati con i medesimi dati, questa curva è molto simile alla precedente ma tende a essere maggiore in alcuni punti in quanto la dimensione del dataset test è minore di quella del dataset train. La curva CV, corrisponde ai modelli addestrati con i dati train ma validati sui dati test, essa risulta essere la più alta ma anche la più veritiera (meno distorta).
df = tibble(nodi = c(2:10),mseTrain,mseTest,mseCV) %>% na.omit()
ggplot(df,mapping = aes(x = nodi)) +
geom_line(aes(y = mseTrain,
color = "Train")) +
geom_line(aes(y = mseTest,
color = "Test")) +
geom_line(aes(y = mseCV,
color = "CV")) +
labs(y = "MSE", x = "Nodi",color = "Legenda") 
MARS
Stimo un primo modello MARS.
modMars = earth(LFACE ~ ., data = tl, nprune = 7, degree = 1)
modMars %>% summary## Call: earth(formula=LFACE~., data=tl, degree=1, nprune=7)
##
## coefficients
## (Intercept) 11.5522454
## h(EDUCATION-14) -1.6308494
## h(EDUCATION-15) 4.1403612
## h(EDUCATION-16) -2.4230637
## h(4-NUMHH) -0.4417816
## h(LINCOME-10.6919) 0.7900094
##
## Selected 6 of 15 terms, and 3 of 5 predictors (nprune=7)
## Termination condition: Reached nk 21
## Importance: LINCOME, NUMHH, EDUCATION, AGE-unused, MARSTATsingle-unused
## Number of terms at each degree of interaction: 1 5 (additive model)
## GCV 2.194151 RSS 556.0856 GRSq 0.3753119 RSq 0.4200775plot(modMars)
Provo a vedere se includendo delle interazioni migliori il modello. Il miglior modello, secondo \(R^2\), è il modello con 3 gradi d’interazione.
gradi = which.max(sapply(c(1:5), function(d) earth(LFACE ~ ., data = tl, nprune = 7, degree = d) %>% summary %>% .$rsq))
modMars = earth(LFACE ~ ., data = tl, nprune = 7, degree = gradi)
modMars %>% summary## Call: earth(formula=LFACE~., data=tl, degree=gradi, nprune=7)
##
## coefficients
## (Intercept) 10.7741142
## h(EDUCATION-15) 2.8382491
## h(EDUCATION-16) -1.8479160
## h(4-NUMHH) -0.4034998
## h(LINCOME-10.1266) 1.3528800
## h(EDUCATION-14) * h(LINCOME-9.39266) -0.7683679
## h(EDUCATION-15) * h(LINCOME-10.7144) 0.8783589
##
## Selected 7 of 19 terms, and 3 of 5 predictors (nprune=7)
## Termination condition: Reached nk 21
## Importance: LINCOME, NUMHH, EDUCATION, AGE-unused, MARSTATsingle-unused
## Number of terms at each degree of interaction: 1 4 2
## GCV 2.130717 RSS 519.7476 GRSq 0.393372 RSq 0.4579731La variabile AGE non viene utilzzita, mentre vengono
considerate alcune interazioni tra le variabili EDUCATION e
LINCOME.
vip::vip(modMars)