Analisi iniziale dei dati

Prendiamo in esame il dataset TL che raccoglie dati di tipo assicurativo ed iniziamo ad esplorarne le caratteristiche.

data<-read.csv("TL.csv",header=TRUE,row.names=1) 
plot(data) #-1 così tolgo la colonna X che mi serve da indice

summary(data)
##       AGE          MARSTAT            EDUCATION         NUMHH     
##  Min.   :20.00   Length:275         Min.   : 2.00   Min.   :1.00  
##  1st Qu.:39.00   Class :character   1st Qu.:13.00   1st Qu.:2.00  
##  Median :47.00   Mode  :character   Median :16.00   Median :3.00  
##  Mean   :47.83                      Mean   :14.52   Mean   :2.96  
##  3rd Qu.:57.00                      3rd Qu.:16.50   3rd Qu.:4.00  
##  Max.   :82.00                      Max.   :17.00   Max.   :9.00  
##      INCOME              FACE         
##  Min.   :     260   Min.   :     800  
##  1st Qu.:   38000   1st Qu.:   50000  
##  Median :   65000   Median :  150000  
##  Mean   :  208975   Mean   :  747581  
##  3rd Qu.:  120000   3rd Qu.:  590000  
##  Max.   :10000000   Max.   :14000000

Notiamo come non siano tutte di tipo int ma abbiamo anche una variabile (MARSTAT nello specifico) di tipo character. inoltre possiamo notare come le 2 variabili FACE e INCOME siano fortemente concentrate sui primi valori, presentando una skewness positiva.

hist(data$FACE, breaks=20) 

hist(data$INCOME, breaks=20)

Per rimediare a ciò possiamo applicare la traformazione logaritmica a queste due variabili e vedere se migliora la loro distribuzione.

data$FACE<-log(data$FACE)
data$INCOME<-log(data$INCOME)

hist(data$FACE, breaks=20)

hist(data$INCOME,breaks=20)

Oltre a ricordare meglio una distribuzione normale, i 2 grafici presentano anche più simmetria tra i valori.

Trasformiamo la variabile MARSTAT in un fattoriale in modo da poter creare un primo grafico di correlazione tra le variabili indipendenti, per avere un’anteprima di cosa aspettarsi con i modelli che eseguiremo successivamente.

library(corrplot)
## corrplot 0.92 loaded
data$MARSTAT<-as.factor(data$MARSTAT)
data$MARSTAT<-as.numeric(data$MARSTAT) #1 married,2 single
corrplot(cor(data),method="number")

Essendo FACE la nostra variabile risposta ci interessa soprasttutto vedere le correlazioni che ha con le altre variabili, in particolare si può notare una correlazione positiva più o meno forte con tutte le variabili tranne che per MARSTAT con la quale è presenta una correlazione negativa (la correlazione con AGE è abbastanza piccola da poter essere omessa).

Modelli best subset e stepwise

Attraverso la funzione regsubset() del pacchetto leaps costruiamo i 3 modelli:

library(leaps)
mfor<-regsubsets(FACE~.,data,method = "forward")
mback<-regsubsets(FACE~.,data,method = "backward")
mbest<-regsubsets(FACE~.,data,method = "exhaustive")

back.summary<-summary(mback)
forw.summary<-summary(mfor)
best.summary<-summary(mbest)

summary(mfor)
## Subset selection object
## Call: regsubsets.formula(FACE ~ ., data, method = "forward")
## 5 Variables  (and intercept)
##           Forced in Forced out
## AGE           FALSE      FALSE
## MARSTAT       FALSE      FALSE
## EDUCATION     FALSE      FALSE
## NUMHH         FALSE      FALSE
## INCOME        FALSE      FALSE
## 1 subsets of each size up to 5
## Selection Algorithm: forward
##          AGE MARSTAT EDUCATION NUMHH INCOME
## 1  ( 1 ) " " " "     " "       " "   "*"   
## 2  ( 1 ) " " " "     "*"       " "   "*"   
## 3  ( 1 ) " " " "     "*"       "*"   "*"   
## 4  ( 1 ) " " "*"     "*"       "*"   "*"   
## 5  ( 1 ) "*" "*"     "*"       "*"   "*"

All’interno del summary è possibile vedere l’ordine di scelta delle variabili con il metodo stepwise, nel caso del metodo forward possiamo vedere che la prima variabile ad essere inserita è EDUCATION, mentre l’ultima è AGE.

Ora facciamo creiamo i plot per le variabili predittive per metterle a confronto.

par(mfrow=c(1,4))
plot(best.summary$rss,type="l")
plot(best.summary$cp,type="l")
plot(best.summary$adjr2,type="l")
plot(best.summary$bic,type="l")

par(mfrow=c(1,4))
plot(forw.summary$rss,type="l")
plot(forw.summary$cp,type="l")
plot(forw.summary$adjr2,type="l")
plot(forw.summary$bic,type="l")

par(mfrow=c(1,4))
plot(back.summary$rss,type="l")
plot(back.summary$cp,type="l")
plot(back.summary$adjr2,type="l")
plot(back.summary$bic,type="l")

I grafici sono sostanzialmente uguali per i 3 metodi quidni possiamo fare un’analisi unica: notiamo come l’RSS sia l’unico che non penalizza per il numero di parametri, infatti dà come miglior modello quello con il maggior numero di parametri,cioè 5; mentre gli altri 3 applicano una penalizzazzione ma mentre Rquadro aggiustato e il CP sono concordi nell’indicare 4 come numero ottimale di parametri, il BIC è più parsimonioso e indica come migliore il modello con soli 3 parametri.

Scelta del modello tramite errori sul set di validazione e cross-validation

Proviamo ora a scegliere il miglior modello tramite un calcolo diretto dell’errore di valutazione.

Il primo metodo consisterà nel dividere i nostri dati in training set e test set: sul primo fitteremo il modello e sul secondo calcoleremo l’errore di valutazione. Calcoleremo l’errore per ogni numero di predittori possibili, così da poter indicare non solo il migliore tra il metodo di subset ma anche tra i numeri di predittori.

: regsubset() non ha e non eredita la funzione predict(), per ovviare a questo problema si calcoleranno a mano i valori previsti tramite il prodotto tra il vettore dei coefficienti per ogni modello stimato e la matrice dei valori per il set test, per fare ciò verra creata una matrice modello in modo tale da poter calcolare la matrice prodotto tra quest’ultima e il vettore dei coefficienti.

Iniziamo con il bestsubset:

set.seed(3)
data<-read.csv("TL.csv",header=TRUE) 
v<-sample(data[,1],200)
data<-data[,-1]
data$FACE<-log(data$FACE)
data$INCOME<-log(data$INCOME)
data$MARSTAT<-as.factor(data$MARSTAT)
data$MARSTAT<-as.numeric(data$MARSTAT) #1 married,2 single
train<-data[v,]
test<-data[-v,]

mbest<-regsubsets(FACE~.,train,method = "exhaustive")
ValError<-c()
test.mat <- model.matrix (FACE~., data = test)

for(i in 1:5){
  coefi<-coef(mbest, i)
  pred<-test.mat[, names(coefi)] %*% coefi
  ValError[i] <- mean((test[,1] - pred)^2)
}
which.min(ValError)
## [1] 3
a<-ValError

Ora lo stepwise foward:

mfor<-regsubsets(FACE~.,train,method = "forward")
for(i in 1:5){
  coefi<-coef(mfor, i)
  pred<-test.mat[, names(coefi)] %*% coefi
  ValError[i] <- mean((test[,1] - pred)^2)
}
which.min(ValError)
## [1] 3
b<-ValError

Ed infine lo stepwise backward:

mback<-regsubsets(FACE~.,train,method = "backward")
for(i in 1:5){
  coefi<-coef(mback, i)
  pred<-test.mat[, names(coefi)] %*% coefi
  ValError[i] <- mean((test[,1] - pred)^2)
}
which.min(ValError)
## [1] 3
c<-ValError

Tutte e 3 indicano come modello ottimale quello che ha al suo interno 3 predittori, possiamo vedere graficamente che sono sostanzialmente uguali.

par(mfrow=c(1,3))
plot(a,type="l",xlab="K",ylab="MSE-Best subset")
plot(b,type="l",xlab="K",ylab="MSE-Stepwise forward")
plot(c,type="l",xlab="K",ylab="MSE-Stepwise backward")

L’analisi appena fatta però è troppo suscettibile alla particolare partizione fatta, prendendo un’altra divisione dei nostri dati potremmo trovare risultati diversi da quelli trovati in questo momento. La cross-validation evita questo problema, sia che si applichi in modalità Leave One Out Cross Validation sia che la si applicchi con i fold.

Proveremo inizialmente con un f=10:

prima creo un vettore in cui ripeto i numeri da 1 a 10 per n volte, poi con sample lo ordino in maniera casuale

f<-10
n<-nrow(data)
folds <- sample ( rep (1:f, length = n))

Creiamo la matrice nella quale salvare i risultati delle varie validazioni, nello specifico avremo una matrice 10x5, una riga per ogni fold e una colonna per ogni numero di variabili possibili del modello

cv.errors <- matrix (NA, nrow=f, ncol=5,dimnames = list (NULL , paste (1:5))) #in di,names sepcifico solo il nome delle colonne perchè alle righe viene assegnato in automatico il numero , mentre per le colonne mette V1,V2...

Implementiamo il ciclo di validazione per ogni j-esimo fold

for (j in 1:f){
  test.mat <- model.matrix (FACE~., data = data[folds==j,])
  mbest<-regsubsets(FACE~.,data=data[folds != j,])
  for (i in 1:5){
    coefi<-coef(mbest, i)
    pred<-test.mat[, names(coefi)] %*% coefi
    cv.errors[j, i] <- mean((data$FACE[folds == j]- pred)^2)
  }
}

Facciamo la media rispetto alle colonne,quindi rispetto ai folds, in modo tale da trovare quale modello ha il minor MSE test

mean.cv.errors<-c()
for (i in 1:5){
  mean.cv.errors[i]<- mean(cv.errors[,i])
}
mean.cv.errors
## [1] 2.747175 2.581845 2.413218 2.408639 2.412599
which.min(mean.cv.errors)
## [1] 4
plot(mean.cv.errors,type="b",xlab="K",ylab="Mean MSE Test")

Il grafico ci indica che lMSE test minore lo sitrova per un K=4, questo risultato ci viene indicato per una cross-validation che utilizza il metodo Best Subset, per completezza rifacciamo il tutto utilizzando il metodo Stepwise forward in modo tale da poter confrontare i dati

cv.errors <- matrix (NA, nrow=f, ncol=5,dimnames = list (NULL , paste (1:5)))
for (j in 1:f){
  test.mat <- model.matrix (FACE~., data = data[folds==j,])
  mforw<-regsubsets(FACE~.,data=data[folds != j,],method = "forward")
  for (i in 1:5){
    coefi<-coef(mforw, i)
    pred<-test.mat[, names(coefi)] %*% coefi
    cv.errors[j, i] <- mean((data$FACE[folds == j]- pred)^2)
  }
}  

mean.cv.errors<-c()
for (i in 1:5){
  mean.cv.errors[i]<- mean(cv.errors[,i])
}
mean.cv.errors
## [1] 2.747175 2.581845 2.413218 2.408639 2.412599
which.min(mean.cv.errors)
## [1] 4
plot(mean.cv.errors,type="b",xlab="K",ylab="Mean MSE Test")

Il risultato è identico.

Ridge e Lasso

Useremo il packages glmnet per eseguire entrambi i metodi, nello specifico useremo la funzione glmnet() per creare il modello e cv.glmnet per trovare il valore ottimale di lambda attraverso una cross-validation in maniera automatica senza doverla scrivere a mano.

A differenza di lm() e regsubsets(), la sintassi è differente e in più viene richiesto che la var.risposta sia un vettore e le variabili indipendenti una matrice.

library(glmnet)
## Caricamento del pacchetto richiesto: Matrix
## Loaded glmnet 4.1-6
x<-model.matrix(FACE~.,data = train)
y<-train$FACE
xtest<-model.matrix(FACE~.,data = test)
ytest<-test$FACE
lambda<-cv.glmnet(x,y,alpha=0)
minlambda<-lambda$lambda.min
rmod<-glmnet (x,y, alpha = 0,lambda = minlambda)
ridg.pred<-predict(rmod, s=minlambda, newx=xtest)
mean((ridg.pred - ytest)^2)
## [1] 1.733046
predict(rmod , type = "coefficients", s =minlambda)
## 7 x 1 sparse Matrix of class "dgCMatrix"
##                       s1
## (Intercept)  4.654078829
## (Intercept)  .          
## AGE         -0.001400292
## MARSTAT     -0.438235669
## EDUCATION    0.167316762
## NUMHH        0.254459567
## INCOME       0.424364064

Ovviamente il metodo Ridge non manda nessuna variabile a 0, per questo risultano tutte le variabili indipendenti significative.

Per implementare il Lasso basata ripetere il codice precedente cambiando il parametro alpha da 0 a 1

x<-model.matrix(FACE~.,data = train)
y<-train$FACE
xtest<-model.matrix(FACE~.,data = test)
ytest<-test$FACE
lambda<-cv.glmnet(x,y,alpha=1)
minlambda<-lambda$lambda.min
rmod<-glmnet (x,y, alpha = 1,lambda = minlambda)
ridg.pred<-predict(rmod, s=minlambda, newx=xtest)
mean((ridg.pred - ytest)^2)
## [1] 1.735878
predict(rmod , type = "coefficients", s =minlambda)
## 7 x 1 sparse Matrix of class "dgCMatrix"
##                     s1
## (Intercept)  3.8340325
## (Intercept)  .        
## AGE          .        
## MARSTAT     -0.3397008
## EDUCATION    0.1765092
## NUMHH        0.2864429
## INCOME       0.4607279

Il metodo Lasso manda a 0 solo la variabile AGE per il lambda dato (che abbiamo scelto tramite una cross-validation, quindi il migliore), le altre variabili risultano tutte significative per il modello.