Laboratorio 2
Erik De Luca
2022-10-21
Analisi dei dati
Importo il file TL.csv e visualizzo i dati. Le
variabili EDUCATION, NUMHH e
MARSTAT sono state trasformate in variabili fattoriali.
tl = read.csv("dataset/TL.csv") %>% tibble %>% dplyr::select(-X)
tl = tl %>% mutate_at(c( "MARSTAT"), as.factor)
tl#"NUMHH", "EDUCATION",Trasformazione variabili
Dal grafico sottostante è evidente come sia necessaria una
trasformazione delle variabili AGE e INCOME.
Attraverso un grafico del tipo scatola a baffi o un istogramma.
ggplot(tl, aes(y = FACE, x = INCOME, col = AGE)) +
geom_point()
ggplot(tl,aes(x = FACE)) +
geom_histogram(bins = 50) +
labs(title = "Istogramma della variabile FACE prima della trasformazione logaritmica")
tl = tl %>% mutate_at(c("FACE","INCOME"),log)
ggplot(tl,aes(x = FACE, y = ..density..)) +
geom_histogram(bins = 30,) +
geom_density() +
geom_density(color = "green",
fill = "green",
alpha = 0.1, # densità del colore
kernel = "gaussian",
adjust = 2) +
labs(title = "Istogramma della variabile FACE dopo della trasformazione logaritmica")
Modello
# train = sample(c(1:dim(tl)[1]), .7*dim(tl)[1], replace = F)
set.seed(1)
train = createDataPartition(tl$AGE,p = .8,list = F)
modAll = lm(FACE ~ ., data = tl, subset = train)
gtsummary::tbl_regression(modAll) | Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| AGE | -0.01 | -0.02, 0.01 | 0.5 |
| MARSTAT | |||
| not single | — | — | |
| single | -0.42 | -1.0, 0.16 | 0.2 |
| EDUCATION | 0.19 | 0.10, 0.27 | <0.001 |
| NUMHH | 0.25 | 0.09, 0.42 | 0.002 |
| INCOME | 0.42 | 0.25, 0.60 | <0.001 |
| 1 CI = Confidence Interval | |||
modFor = regsubsets(FACE ~ .,data = tl[train,],method = "forward")
modFor %>% summary## Subset selection object
## Call: regsubsets.formula(FACE ~ ., data = tl[train, ], method = "forward")
## 5 Variables (and intercept)
## Forced in Forced out
## AGE FALSE FALSE
## MARSTATsingle FALSE FALSE
## EDUCATION FALSE FALSE
## NUMHH FALSE FALSE
## INCOME FALSE FALSE
## 1 subsets of each size up to 5
## Selection Algorithm: forward
## AGE MARSTATsingle EDUCATION NUMHH INCOME
## 1 ( 1 ) " " " " " " " " "*"
## 2 ( 1 ) " " " " " " "*" "*"
## 3 ( 1 ) " " " " "*" "*" "*"
## 4 ( 1 ) " " "*" "*" "*" "*"
## 5 ( 1 ) "*" "*" "*" "*" "*"modBack = MASS::stepAIC(modAll,direction = "backward")## Start: AIC=200.54
## FACE ~ AGE + MARSTAT + EDUCATION + NUMHH + INCOME
##
## Df Sum of Sq RSS AIC
## - AGE 1 0.993 520.03 198.97
## <none> 519.04 200.54
## - MARSTAT 1 4.873 523.91 200.62
## - NUMHH 1 22.803 541.84 208.09
## - EDUCATION 1 45.512 564.55 217.20
## - INCOME 1 55.565 574.60 221.12
##
## Step: AIC=198.97
## FACE ~ MARSTAT + EDUCATION + NUMHH + INCOME
##
## Df Sum of Sq RSS AIC
## - MARSTAT 1 4.289 524.32 198.79
## <none> 520.03 198.97
## - NUMHH 1 29.582 549.61 209.25
## - EDUCATION 1 44.711 564.74 215.28
## - INCOME 1 55.019 575.05 219.29
##
## Step: AIC=198.79
## FACE ~ EDUCATION + NUMHH + INCOME
##
## Df Sum of Sq RSS AIC
## <none> 524.32 198.79
## - EDUCATION 1 44.334 568.65 214.81
## - NUMHH 1 48.977 573.30 216.62
## - INCOME 1 62.961 587.28 221.97Scelta del modello
predAll = predict.lm(modAll,newdata = tl[-train,])
predBack = predict.lm(modBack,newdata = tl[-train,])
data.frame(R2 = c(R2(predAll,tl[-train,"FACE"]),R2(predBack,tl[-train,"FACE"])),
row.names = c("Tutte le variabili", "Backward"))Validation set
testMat = model.matrix(FACE~., data = tl[-train,])
valErr = sapply(c(1:(length(tl)-1)),function(i) mean((tl$FACE[-train]-testMat[,names(coef(modFor, id = i))]%*%coef(modFor, id = i))^2))
ggplot(data.frame(MSE = valErr, numeroDiVariabili = c(1:5)),
mapping = aes(y = MSE, x = numeroDiVariabili)) +
geom_line() +
geom_point()
Utilizzando l’approccio del set di validazione risulta migliore il modello con tutte e 5 le variabili indipendenti.
Cross validation
Siccome non esiste il metodo predict per
regsubset necessitiamo allora di specificarlo prima di
poterlo chiamare in modo da semplificare i calcoli successivi. Nel punto
precedento, invece, è stata scritta l’intera procedura. Si poteva
seguire lo stesso procedimento anche in questo caso.
Il pezzo di codice che definisce la funzione predict
(copiato da web) è stato poi adattato in modo tale che potesse ricevere
in input oggetti provenienti da liste o funzioni.
predict.regsubsets = function(object, newdata, id, form = NULL, ...)
{
if(is.null(form))
form = as.formula(object$call[[2]])
mat = model.matrix(form, newdata)
coefi = coef(object, id = id)
mat[, names(coefi)] %*% coefi
}Nel seguente pezzo di codice si poteva evitare di utilizzare la
funzione map_dfc andado a calcolare il valore medio fuori
dal sapply ottenendo così un unico e valore, e
successivamente applicare un altro sapply ottenendo questa
volta un vettore anziché un dataframe.
k = 10 # numero di folds
set.seed(1)
folds = sample(1:k, nrow(tl), replace = TRUE)
cvErr = purrr::map_dfc(c(1:k),function(j) sapply(c(1:5),function(i) mean((tl$FACE[folds==j]-predict(regsubsets(FACE~., data = tl[folds!=j,], nvmax=6), tl[folds == j,], id = i,form = as.formula(FACE~.)))^2)))## New names:
## • `` -> `...1`
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`cvErrMean = apply(cvErr,2,mean)
ggplot(data.frame(MSE = cvErrMean, numeroDiVariabili = c(1:5)),
mapping = aes(y = MSE, x = numeroDiVariabili)) +
geom_line() +
geom_point()
Dal grafico soprastante si deduce che il miglior modello in questo caso è quello con 2 variabili.
Ridge e Lasso
matTl = makeX(tl[,names(tl)!="FACE"])
# per trovare lambda ottimale attraverso cross validation
cvRidge = cv.glmnet(matTl,
tl$FACE,
alpha = 0)
cvLasso = cv.glmnet(matTl,
tl$FACE,
alpha = 1)
# lambda ottimale
data.frame(Ridge = cvRidge$lambda.min,
Lasso = cvLasso$lambda.min,
row.names = "Lambda ottimale")