Regressione II (notebook 8)

Statistica II - 750AA

Introduzione

In questo notebook esploriamo la collinearità e multicollinearità nei modelli di regressione lineare, gli indicatori principali (correlogramma, VIF, numero di condizionamento) e le tecniche di riduzione dei modelli (stepwise) e regolarizzazione (Ridgee Lasso).

library(car)
library(ggcorrplot)

Dataset swiss

Consideriamo il dataset swiss che contiene indicatori di fertilità e socioeconomici delle province (di lingua francese) della svizzera nel 1888.

data(swiss)
?swiss
head(swiss)

Studiamo un modello di regressione lineare con la variabile Fertility (un indice di fertilità) come fattore di uscita, e tutte le altre come caratteristiche di ingresso.

# Modello iniziale
fit_swiss_full <- lm(Fertility ~ ., data = swiss)
summary(fit_swiss_full)

Call:
lm(formula = Fertility ~ ., data = swiss)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.2743  -5.2617   0.5032   4.1198  15.3213 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      66.91518   10.70604   6.250 1.91e-07 ***
Agriculture      -0.17211    0.07030  -2.448  0.01873 *  
Examination      -0.25801    0.25388  -1.016  0.31546    
Education        -0.87094    0.18303  -4.758 2.43e-05 ***
Catholic          0.10412    0.03526   2.953  0.00519 ** 
Infant.Mortality  1.07705    0.38172   2.822  0.00734 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.165 on 41 degrees of freedom
Multiple R-squared:  0.7067,    Adjusted R-squared:  0.671 
F-statistic: 19.76 on 5 and 41 DF,  p-value: 5.594e-10

Osserviamo che alcuni \(p\)-value dei coefficienti risultano molto elevati, un primo segno di collinearità, mentre il coefficiente \(R^2\) è abbastanza alto e il \(p\)-value dell’\(F\)-test è basso, segno che i fattori di ingresso nel modello spiegano una buona parte dell’uscita.

Diagnosi di collinearità

Analizziamo le correlazioni tra i fattori di ingresso mediante la matrice delle correlazioni e un correlogramma.

# Matrice di correlazione (approssimata a due cifra decimali per semplicità)
round(cor(swiss[, -1]), 2)

                 Agriculture Examination Education Catholic Infant.Mortality
Agriculture             1.00       -0.69     -0.64     0.40            -0.06
Examination            -0.69        1.00      0.70    -0.57            -0.11
Education              -0.64        0.70      1.00    -0.15            -0.10
Catholic                0.40       -0.57     -0.15     1.00             0.18
Infant.Mortality       -0.06       -0.11     -0.10     0.18             1.00

# Diagramma a nuvola di punti (non efficace: troppe variabili, eventualmente selezionarne alcune)
#pairs(mtcars, main = "Scatterplot (mtcars)")

# heatmap con ggcorrplot
library(ggcorrplot)
corr_mat_mtcars <- cor(swiss[, -1])
ggcorrplot(corr_mat_mtcars, lab = TRUE, type = "lower",
           title = "Matrice di correlazione tra predittori (mtcars)")

Abbiamo alcune correlazioni tra coppie (Examination e Education), il caso più semplice di collinearità. Studiamo gli indicatori specifici, il VIF e il condizionamento della matrice di Gram. Per il primo usiamo il comando vif() dal pacchetto car.

# Variance Inflation Factor
library(car)

vif(fit_swiss_full)

     Agriculture      Examination        Education         Catholic Infant.Mortality 
        2.284129         3.675420         2.774943         1.937160         1.107542 

Ricordiamo che valori maggiori di \(5\) indicano variabili per cui potrebbe essere presente collinearità (in particolare se maggiore di \(10\)). Notiamo che nessun VIF supera la soglia \(5\) ma Examination presenta un valore comunque grande.

Vediamo infine la matrice di Gram e il numero di condizionamento.

# estrae la matrice di design del modello (dati con una colonna di 1 per l'intercetta)
swiss_features <- model.matrix(fit_swiss_full)

# calcola la matrice di Gram
G_swiss <- t(swiss_features) %*% swiss_features

G_swiss

                 (Intercept) Agriculture Examination Education  Catholic Infant.Mortality
(Intercept)            47.00     2381.00      775.00    516.00   1933.76           937.30
Agriculture          2381.00   144347.22    33539.10  19716.10 115439.06         47298.03
Examination           775.00    33539.10    15707.00  10973.00  23120.68         15333.60
Education             516.00    19716.10    10973.00   9918.00  18392.07         10162.40
Catholic             1933.76   115439.06    23120.68  18392.07 159569.84         39544.75
Infant.Mortality      937.30    47298.03    15333.60  10162.40  39544.75         19082.41

Per il numero di condizionamento usiamo la funzione norm() per calcolare (una) norma della matrice \(G\) e solve() per l’inversa.

kappa_G <- norm(G_swiss, type="2") * norm(solve(G_swiss, type="2"))

kappa_G

[1] 681041.3

Un valore maggiore di \(1000\) indica che \(G\) è mal condizionata (segno di collinearità). Vediamo per confronto un semplice esempio generato.

# genera una matrice quasi non-invertibile
eps = 10^{-3}

G = matrix( c(1, 1-eps, 1-eps, 1 ), byrow=TRUE, nrow=2, ncol=2)

G

      [,1]  [,2]
[1,] 1.000 0.999
[2,] 0.999 1.000

kappa_G = norm(G)*norm(solve(G))
print(paste("Numero di condizionamento:", round(kappa_G, 0)))

[1] "Numero di condizionamento: 1999"

Un indicatore alternativo è il determinante della matrice di correlazione dei fattori di ingresso.

round(cor(swiss), 2)

                 Fertility Agriculture Examination Education Catholic Infant.Mortality
Fertility             1.00        0.35       -0.65     -0.66     0.46             0.42
Agriculture           0.35        1.00       -0.69     -0.64     0.40            -0.06
Examination          -0.65       -0.69        1.00      0.70    -0.57            -0.11
Education            -0.66       -0.64        0.70      1.00    -0.15            -0.10
Catholic              0.46        0.40       -0.57     -0.15     1.00             0.18
Infant.Mortality      0.42       -0.06       -0.11     -0.10     0.18             1.00

det(cor(swiss))

[1] 0.03496976

Riduzione del modello

Per eliminare eventuali variabili collineari, consideriamo gli approcci basati sui test di significatività dei coefficienti e sugli indici AIC/BIC.

Riprendiamo i \(p\)-value del modello completo.

summary(fit_swiss_full)

Call:
lm(formula = Fertility ~ ., data = swiss)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.2743  -5.2617   0.5032   4.1198  15.3213 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      66.91518   10.70604   6.250 1.91e-07 ***
Agriculture      -0.17211    0.07030  -2.448  0.01873 *  
Examination      -0.25801    0.25388  -1.016  0.31546    
Education        -0.87094    0.18303  -4.758 2.43e-05 ***
Catholic          0.10412    0.03526   2.953  0.00519 ** 
Infant.Mortality  1.07705    0.38172   2.822  0.00734 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.165 on 41 degrees of freedom
Multiple R-squared:  0.7067,    Adjusted R-squared:  0.671 
F-statistic: 19.76 on 5 and 41 DF,  p-value: 5.594e-10

Una soluzione rapida (ma potenzialmente rischiosa!) è di eliminare direttamente le variabili con \(p\)-value maggiore di una determinata soglia (ad esempio 0.1).

fit_swiss_reduced <- lm(Fertility ~ .- Agriculture - Examination, data = swiss)
summary(fit_swiss_reduced)

Call:
lm(formula = Fertility ~ . - Agriculture - Examination, data = swiss)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.4781  -5.4403  -0.5143   4.1568  15.1187 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      48.67707    7.91908   6.147 2.24e-07 ***
Education        -0.75925    0.11680  -6.501 6.83e-08 ***
Catholic          0.09607    0.02722   3.530  0.00101 ** 
Infant.Mortality  1.29615    0.38699   3.349  0.00169 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.505 on 43 degrees of freedom
Multiple R-squared:  0.6625,    Adjusted R-squared:  0.639 
F-statistic: 28.14 on 3 and 43 DF,  p-value: 3.15e-10

Il coefficiente \(R^2\) (o meglio quello Adjusted che tiene conto del numero di variabili) è sceso leggermente ma non drasticamente. Confrontiamo con il modello in cui eliminiamo solo Examination.

fit_swiss_reduced_one_step <- lm(Fertility ~ .- Examination, data = swiss)
summary(fit_swiss_reduced_one_step)

Call:
lm(formula = Fertility ~ . - Examination, data = swiss)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.6765  -6.0522   0.7514   3.1664  16.1422 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      62.10131    9.60489   6.466 8.49e-08 ***
Agriculture      -0.15462    0.06819  -2.267  0.02857 *  
Education        -0.98026    0.14814  -6.617 5.14e-08 ***
Catholic          0.12467    0.02889   4.315 9.50e-05 ***
Infant.Mortality  1.07844    0.38187   2.824  0.00722 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.168 on 42 degrees of freedom
Multiple R-squared:  0.6993,    Adjusted R-squared:  0.6707 
F-statistic: 24.42 on 4 and 42 DF,  p-value: 1.717e-10

In questo caso l’Adjusted \(R^2\) è rimasto invariato.

Problema: discutere cosa accade se invece si rimuove la variabile Education.

Ricordiamo che AIC e BIC cercano di indicare la performance del modello in termini del compromesso tra bias (numero di parametri) e varianza (valore della verosimiglianza). Per calcolarli è sufficiente usare la funzione AIC()

AIC(fit_swiss_full, fit_swiss_reduced, fit_swiss_reduced_one_step)

Vediamo quindi che il modello in cui abbiamo eliminato solo la variabile Examination è preferibile in termini di AIC. Similmente per il BIC:

BIC(fit_swiss_full, fit_swiss_reduced, fit_swiss_reduced_one_step)

Ovviamente non si può procedere a mano con la riduzione del modello, specialmente se ci sono molte caratteristiche. Per questo le procedure stepwise sono automatizzate. Usiamo il comando step (già caricato).

Selezione all’indietro.

# Backward
fit_back <- step(fit_swiss_full, direction = "backward", trace = 0) # mettere trace = 1 o TRUE per avere informazioni sui passi del'algoritmo

summary(fit_back)

Call:
lm(formula = Fertility ~ Agriculture + Education + Catholic + 
    Infant.Mortality, data = swiss)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.6765  -6.0522   0.7514   3.1664  16.1422 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      62.10131    9.60489   6.466 8.49e-08 ***
Agriculture      -0.15462    0.06819  -2.267  0.02857 *  
Education        -0.98026    0.14814  -6.617 5.14e-08 ***
Catholic          0.12467    0.02889   4.315 9.50e-05 ***
Infant.Mortality  1.07844    0.38187   2.824  0.00722 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.168 on 42 degrees of freedom
Multiple R-squared:  0.6993,    Adjusted R-squared:  0.6707 
F-statistic: 24.42 on 4 and 42 DF,  p-value: 1.717e-10

Selezione in avanti.

# Forward
# iniziamo con il modello costante (banale)
fit_swiss_trivial <- lm(Fertility ~ 1, data = swiss)


fit_forw <- step(fit_swiss_trivial,
                 direction = "forward",
                 scope = formula(fit_swiss_full), #specifica un modello `massimo`
                 trace = FALSE)
summary(fit_forw)

Call:
lm(formula = Fertility ~ Education + Catholic + Infant.Mortality + 
    Agriculture, data = swiss)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.6765  -6.0522   0.7514   3.1664  16.1422 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      62.10131    9.60489   6.466 8.49e-08 ***
Education        -0.98026    0.14814  -6.617 5.14e-08 ***
Catholic          0.12467    0.02889   4.315 9.50e-05 ***
Infant.Mortality  1.07844    0.38187   2.824  0.00722 ** 
Agriculture      -0.15462    0.06819  -2.267  0.02857 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.168 on 42 degrees of freedom
Multiple R-squared:  0.6993,    Adjusted R-squared:  0.6707 
F-statistic: 24.42 on 4 and 42 DF,  p-value: 1.717e-10

E infine una selezione in entrambe le direzioni.

# Both directions
fit_both <- step(fit_swiss_full, direction = "both", trace = FALSE, k=log(nrow(swiss))) #in questo caso impostiamo BIC

summary(fit_both)

Call:
lm(formula = Fertility ~ Agriculture + Education + Catholic + 
    Infant.Mortality, data = swiss)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.6765  -6.0522   0.7514   3.1664  16.1422 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      62.10131    9.60489   6.466 8.49e-08 ***
Agriculture      -0.15462    0.06819  -2.267  0.02857 *  
Education        -0.98026    0.14814  -6.617 5.14e-08 ***
Catholic          0.12467    0.02889   4.315 9.50e-05 ***
Infant.Mortality  1.07844    0.38187   2.824  0.00722 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.168 on 42 degrees of freedom
Multiple R-squared:  0.6993,    Adjusted R-squared:  0.6707 
F-statistic: 24.42 on 4 and 42 DF,  p-value: 1.717e-10

Attenzione: Le procedure stepwise ottimizzano un criterio (AIC, BIC), ma non sempre portano a un modello interpretabile o stabile. È buona pratica consideriare gli indicatori (\(p\)-value, VIF) anche dopo la riduzione automatica.

vif(fit_both)

     Agriculture        Education         Catholic Infant.Mortality 
        2.147153         1.816361         1.299916         1.107528 

There were 50 or more warnings (use warnings() to see the first 50)

PCR

Un approccio differente consiste nell’applicare la PCA sulle caratteristiche in input. Dividiamo al solito in training e test set con il comando createDataPartition.

library(caret)

# Dividiamo in training e test set il dataset ciqual_reduced

train_index <- createDataPartition(swiss$Fertility,  p = 0.8, list = FALSE)

swiss_train <- swiss[train_index, ]
swiss_test <- swiss[-train_index, ]

Effetuiamo la PCA sulle variabili di ingresso (escludiamo la prima colonna, la variabile di uscita) con il comando preProcess() dal pacchetto caret. In questo modo sarà più semplice fare previsioni con il modello (altrimenti bisognerebbe trasformare ogni volta le caratteristiche di input di un punto di test)

# possiamo specificare quante componenti principali, in questo caso pcaComp=3. È anche possibile impostare una soglia sulla percentuale di varianza cumulata con l'opzione thresh=0.95
pca_model <- preProcess(swiss_train[, -1], method = c("center", "scale", "pca"), pcaComp=3)

Possiamo ora calcolare i fattori di ingresso proiettati sulle componenti principali con il comando predict (sia sul training che sul test set).

pca_train <- predict(pca_model, swiss_train[, -1])
head(pca_train)

pca_test <- predict(pca_model, swiss_test[, -1])
head(pca_test)

Effettuiamo quindi la regressione usando come fattori di ingresso le PC selezionate.

fit_pcr <- lm(swiss_train$Fertility~., data=pca_train)
summary(fit_pcr)

Call:
lm(formula = swiss_train$Fertility ~ ., data = pca_train)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.8587  -5.9633   0.1315   6.2244  16.1768 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  69.7692     1.3886  50.243  < 2e-16 ***
PC1           5.2321     0.8567   6.107 5.58e-07 ***
PC2           4.1389     1.3819   2.995  0.00501 ** 
PC3          -2.3331     1.5193  -1.536  0.13361    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.672 on 35 degrees of freedom
Multiple R-squared:  0.5815,    Adjusted R-squared:  0.5456 
F-statistic: 16.21 on 3 and 35 DF,  p-value: 8.996e-07

vif(fit_pcr)

PC1 PC2 PC3 
  1   1   1 

Le prime 2 componenti spiegano già parte della variabilità e riducono (necessariamente) la collinearità.

Regolarizzazione: Ridge e LASSO

Usiamo il pacchetto caret per usare Ridge, Lasso (anche altri modelli regolarizzati sarebbero disponibili) e determinare i parametri usando cross-validation. Un pacchetto alternativo più specifico è glmnet.

Definiamo prima il metodo di cross-validation.

train_control <- trainControl(method = "cv", number = 10) #10-fold cv

Per la regressione Ridge, specifichiamo il metodo “glmnet” con il parametro \(\alpha=0\) (glmnet è un modello Elastic Net che permette di interpolare tra Ridge e LASSO).

ridge_model <- train(Fertility ~ ., data = swiss_train,  
                        method = "glmnet", 
                        trControl = train_control, 
                        tuneGrid = expand.grid(alpha = 0, lambda = seq(0, 1, length = 100)))

Plottiamo l’errore di CV al variare del parametro di regolarizzazione \(\lambda\).

plot(ridge_model, type='l')

Vediamo quindi che il parametro ottimale (scelto tramite cross-validation) è 10.7070707. Vediamo i parametri determinati.

coef(ridge_model$finalModel, ridge_model$bestTune$lambda)

6 x 1 sparse Matrix of class "dgCMatrix"
                 s=0.8585859
(Intercept)      63.07000721
Agriculture      -0.12330908
Examination      -0.25350028
Education        -0.74633617
Catholic          0.08684621
Infant.Mortality  1.10689452

Per la regressione LASSO, specifichiamo il metodo “glmnet” con il parametro \(\alpha=1\).

lasso_model <- train(Fertility ~ ., data = swiss_train,  
                        method = "glmnet", 
                        trControl = train_control, 
                        tuneGrid = expand.grid(alpha = 1, lambda = seq(0, 10, length = 100)))

Plottiamo l’errore di CV al variare del parametro di regolarizzazione \(\lambda\).

plot(lasso_model, type='l')

lasso_model$bestTune$lambda

[1] 0.1414141

coef(lasso_model$finalModel, lasso_model$bestTune$lambda)

6 x 1 sparse Matrix of class "dgCMatrix"
                 s=0.1414141
(Intercept)       63.2699451
Agriculture       -0.1469023
Examination       -0.1271760
Education         -0.8786019
Catholic           0.1050504
Infant.Mortality   1.0893602

In alternativa possiamo usare il comando cv.glm() dal pacchetto glmnet.

library(glmnet)
X <- as.matrix(swiss_train[, -1])
y <- swiss_train$Fertility
cv_ridge <- cv.glmnet(X, y, alpha = 0)
plot(cv_ridge)

# confronta i lambda trovati tramite i due metodi
print(paste("lambda_min glmnet:", cv_ridge$lambda.min))

[1] "lambda_min glmnet: 0.8545232653727"

print(paste("lambda_min_caret:", ridge_model$bestTune$lambda))

[1] "lambda_min_caret: 0.858585858585859"

# confronto dei coefficienti del modello trovati per i due metodi
coef(cv_ridge)

6 x 1 sparse Matrix of class "dgCMatrix"
                  lambda.1se
(Intercept)      60.63706601
Agriculture       0.01045061
Examination      -0.28175926
Education        -0.32409486
Catholic          0.04008528
Infant.Mortality  0.76982050

coef(ridge_model$finalModel, ridge_model$bestTune$lambda)

6 x 1 sparse Matrix of class "dgCMatrix"
                 s=0.8585859
(Intercept)      63.07000721
Agriculture      -0.12330908
Examination      -0.25350028
Education        -0.74633617
Catholic          0.08684621
Infant.Mortality  1.10689452

Operiamo similmente con LASSO.

cv_lasso <- cv.glmnet(X, y, alpha = 1)
plot(cv_lasso)

# confronta i lambda trovati tramite i due metodi
print(paste("lambda_min glmnet:", cv_lasso$lambda.min))

[1] "lambda_min glmnet: 0.226968880540014"

print(paste("lambda_min_caret:", lasso_model$bestTune$lambda))

[1] "lambda_min_caret: 0.141414141414141"

# confronto dei coefficienti del modello trovati per i due metodi
coef(cv_lasso)

6 x 1 sparse Matrix of class "dgCMatrix"
                  lambda.1se
(Intercept)      57.12305060
Agriculture       .         
Examination      -0.08853346
Education        -0.58404464
Catholic          0.05533704
Infant.Mortality  0.92894431

coef(lasso_model$finalModel, lasso_model$bestTune$lambda)

6 x 1 sparse Matrix of class "dgCMatrix"
                 s=0.1414141
(Intercept)       63.2699451
Agriculture       -0.1469023
Examination       -0.1271760
Education         -0.8786019
Catholic           0.1050504
Infant.Mortality   1.0893602

Dataset CIQUAL

Riprendiamo il dataset CIQUAL sulle caratteristiche nutritive degli alimenti. Carichiamo la versione .csv già ripulita nel corso della lezione sulla PCA.

ciqual_tidy <- read.csv('datasets/ciqual_tidy.csv', header = TRUE)

head(ciqual_tidy)

Eliminiamo le prime due colonne (con le classi di alimenti) e usiamo la terza come colonna di nomi (in modo da non includerla nella regressione).

ciqual_reduced <- ciqual_tidy[, -(1:3)]

row.names(ciqual_reduced) <- ciqual_tidy[,3]

head(ciqual_reduced)

Poniamoci come obiettivo di calcolare l’energia (in kcal) a partire dalle altre informazioni nutritive. Rimuoviamo quindi le colonne (praticamente duplicate) dell’energia (teniamo solo la prima).

ciqual_reduced <- ciqual_reduced[,-(2:4)]

head(ciqual_reduced)

Consideriamo il modello di regressione lineare completo.

ciqual_fit <- lm( ciqual_reduced$Energy..Regulation.EU.No.1169.2011..kJ.100g. ~., data=ciqual_reduced)

summary(ciqual_fit)

Visualizziamo le significatività e i VIF

# plot delle significatività (-log(p-value))
library(ggplot2)
ciqual_coefs <- summary(ciqual_fit)$coefficients

ciqual_coefs_df <- data.frame(
  Term = rownames(ciqual_coefs)[-1],  # escludi l'intercetta
  NegLogPValue = -log10(ciqual_coefs[-1, "Pr(>|t|)"]),  # escludi l'intercetta
  VIF = vif(ciqual_fit)
)
ggplot(ciqual_coefs_df, aes(x = reorder(Term, NegLogPValue), y = NegLogPValue)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  coord_flip() +
  labs(title = "-log10(p-value) dei coefficienti del modello CIQUAL", x = "Termine", y = "-log10(p-value)") +
  theme_minimal()

ggplot(ciqual_coefs_df, aes(x = reorder(Term, VIF), y = VIF)) +
  geom_bar(stat = "identity", fill = "darkorange") +
  coord_flip() +
  labs(title = "Variance Inflation Factor (VIF) dei coefficienti del modello CIQUAL", x = "Termine", y = "VIF") +
  theme_minimal()

Molte caratteristiche hanno un \(p\)-value alto, ma potrebbe esserci collinearità (ad esempio tra gli zuccheri).

Riduzione stepwise

Procediamo con una riduzione all’indietro.

ciqual_reduced_stepwise <- step(ciqual_fit, direction = "backward", trace=0)

summary(ciqual_reduced_stepwise)

Plottiamo al solito significatività e VIF.

ciqual_coefs <- summary(ciqual_reduced_stepwise)$coefficients

ciqual_coefs_df <- data.frame(
  Term = rownames(ciqual_coefs)[-1],  # escludi l'intercetta
  NegLogPValue = -log10(ciqual_coefs[-1, "Pr(>|t|)"]),  # escludi l'intercetta
  VIF = vif(ciqual_reduced_stepwise)
)
ggplot(ciqual_coefs_df, aes(x = reorder(Term, NegLogPValue), y = NegLogPValue)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  coord_flip() +
  labs(title = "-log10(p-value) dei coefficienti del modello CIQUAL", x = "Termine", y = "-log10(p-value)") +
  theme_minimal()

ggplot(ciqual_coefs_df, aes(x = reorder(Term, VIF), y = VIF)) +
  geom_bar(stat = "identity", fill = "darkorange") +
  coord_flip() +
  labs(title = "Variance Inflation Factor (VIF) dei coefficienti del modello CIQUAL", x = "Termine", y = "VIF") +
  theme_minimal()

Procediamo in alternativa con una selezione in avanti

ciqual_reduced_stepwise_forward <- step(lm(ciqual_reduced$Energy..Regulation.EU.No.1169.2011..kJ.100g.~ 1, data=ciqual_reduced), 
                                        scope=formula(ciqual_fit), direction="forward", trace=0)

summary(ciqual_reduced_stepwise_forward)

ciqual_coefs <- summary(ciqual_reduced_stepwise)$coefficients

ciqual_coefs_df <- data.frame(
  Term = rownames(ciqual_coefs)[-1],  # escludi l'intercetta
  NegLogPValue = -log10(ciqual_coefs[-1, "Pr(>|t|)"]),  # escludi l'intercetta
  VIF = vif(ciqual_reduced_stepwise)
)
ggplot(ciqual_coefs_df, aes(x = reorder(Term, NegLogPValue), y = NegLogPValue)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  coord_flip() +
  labs(title = "-log10(p-value) CIQUAL ridotto", x = "Termine", y = "-log10(p-value)") +
  theme_minimal()

ggplot(ciqual_coefs_df, aes(x = reorder(Term, VIF), y = VIF)) +
  geom_bar(stat = "identity", fill = "darkorange") +
  coord_flip() +
  labs(title = "Variance Inflation Factor (VIF)  CIQUAL ridotto", x = "Termine", y = "VIF") +
  theme_minimal()

Vediamo che in entrambi i casi il modello ridotto ha comunque un gran numero di fattori di ingresso. Possiamo però forzare il numero massimo di passi in modo da ridurre le variabili di input.

ciqual_reduced_stepwise_forward <- step(lm(ciqual_reduced$Energy..Regulation.EU.No.1169.2011..kJ.100g.~ 1, data=ciqual_reduced), 
                                        scope=formula(ciqual_fit), direction="forward", steps=5, trace=0)

summary(ciqual_reduced_stepwise_forward)

Verifichiamone la collinearità.

vif(ciqual_reduced_stepwise_forward)

Vediamo che c’è una grande collinearità tra le due variabili relative alle proteine (abbastanza ovvio, sono due colonne duplicate). Quindi possiamo ridurre ulteriormente eliminandone una.

# rimuovi la variabile dalla formula usando update()
ciqual_minimal <- lm( update(ciqual_reduced_stepwise_forward, . ~ . - Protein..crude..N.x.6.25..g.100g.), data=ciqual_reduced)

summary(ciqual_minimal)
vif(ciqual_minimal)

Ovviamente gli indicatori come l’\(R^2\) sono diminuiti, ma abbiamo un modello decisamente semplificato, intepretabile e senza collinearità.

PCR

Dividiamo al solito in training e test set con il comando createDataPartition.

library(caret)

# Dividiamo in training e test set il dataset ciqual_reduced

train_index <- createDataPartition(ciqual_reduced$Energy..Regulation.EU.No.1169.2011..kJ.100g., p = 0.8, list = FALSE)

ciqual_reduced_train <- ciqual_reduced[train_index, ]
ciqual_reduced_test <- ciqual_reduced[-train_index, ]

Effetuiamo la PCA sulle variabili di ingresso (escludiamo la prima colonna, la variabile di uscita) con il comando preProcess() dal pacchetto caret. In questo modo sarà più semplice fare previsioni con il modello (altrimenti bisognerebbe trasformare ogni volta le caratteristiche di input di un punto di test)

# specifichiamo una soglia sulla percentuale di varianza cumulata con l'opzione thresh
pca_model <- preProcess(ciqual_reduced_train[, -1], method = c("center", "scale", "pca"), thresh =0.95)

Possiamo ora calcolare i fattori di ingresso proiettati sulle componenti principali con il comando predict (sia sul training che sul test set).

pca_ciqual_train <- predict(pca_model, ciqual_reduced_train[, -1])
head(pca_ciqual_train)

pca_ciqual_test <- predict(pca_model, ciqual_reduced_test[, -1])
head(pca_ciqual_test)

Effettuiamo quindi la regressione usando come fattori di ingresso le PC selezionate.

fit_pcr <- lm(ciqual_reduced_train$Energy..Regulation.EU.No.1169.2011..kJ.100g. ~., data=pca_ciqual_train)
summary(fit_pcr)
vif(fit_pcr)

Problema: come cambia la performance del modello se specifichiamo un numero basso di PC?

Regolarizzazione LASSO

Usiamo il pacchetto glmnet per la regolarizzazione con LASSO.

X = as.matrix(ciqual_reduced_train[,-1])
y = as.matrix(ciqual_reduced_train[, 1])
cv_lasso <- cv.glmnet(X, y, alpha = 1)
plot(cv_lasso)

print(paste("lambda_min glmnet:", cv_lasso$lambda.min))

[1] "lambda_min glmnet: 9.36755294235116"

# confronto dei coefficienti del modello trovati per i due metodi
coef(cv_lasso)

50 x 1 sparse Matrix of class "dgCMatrix"
                                          lambda.1se
(Intercept)                              149.9657571
Water..g.100g.                            -0.6032222
Protein..g.100g.                          13.6781669
Protein..crude..N.x.6.25..g.100g.          .        
Carbohydrate..g.100g.                      4.9960658
Fat..g.100g.                              21.1724823
Sugars..g.100g.                           -2.6449047
fructose..g.100g.                          .        
galactose..g.100g.                         .        
glucose..g.100g.                           1.6691924
lactose..g.100g.                           .        
maltose..g.100g.                           .        
sucrose..g.100g.                          -5.6103001
Starch..g.100g.                            7.0525990
Fibres..g.100g.                            .        
Polyols..g.100g.                           1.1199149
Ash..g.100g.                               .        
Alcohol..g.100g.                          15.8331160
Organic.acids..g.100g.                     .        
FA.saturated..g.100g.                      9.6775364
FA.mono..g.100g.                           2.3824688
FA.poly..g.100g.                           .        
Cholesterol..mg.100g.                      .        
Salt..g.100g.                              .        
Calcium..mg.100g.                          .        
Chloride..mg.100g.                         .        
Copper..mg.100g.                           .        
Iron..mg.100g.                             3.3243452
Iodine..µg.100g.                           .        
Magnesium..mg.100g.                        .        
Manganese..mg.100g.                        .        
Phosphorus..mg.100g.                       .        
Potassium..mg.100g.                        .        
Selenium..µg.100g.                         .        
Sodium..mg.100g.                           .        
Zinc..mg.100g.                             .        
Retinol..µg.100g.                          .        
Beta.carotene..µg.100g.                    .        
Vitamin.D..µg.100g.                        0.1470879
Vitamin.E..mg.100g.                        9.1144274
Vitamin.K1..µg.100g.                       .        
Vitamin.K2..µg.100g.                       .        
Vitamin.C..mg.100g.                        .        
Vitamin.B1.or.Thiamin..mg.100g.            .        
Vitamin.B2.or.Riboflavin..mg.100g.         .        
Vitamin.B3.or.Niacin..mg.100g.             .        
Vitamin.B5.or.Pantothenic.acid..mg.100g.   .        
Vitamin.B6..mg.100g.                       .        
Vitamin.B9.or.Folate..µg.100g.             .        
Vitamin.B12..µg.100g.                      .        

vif(cv_lasso)

Error in UseMethod("vcov") : 
  no applicable method for 'vcov' applied to an object of class "cv.glmnet"

Confrontiamo con il risultato di caret in cross-validation.

train_control <- trainControl(method = "cv", number = 20)


lasso_model <- train(Energy..Regulation.EU.No.1169.2011..kJ.100g. ~ ., data = ciqual_reduced_train,  
                        method = "glmnet", 
                        trControl = train_control, 
                        tuneGrid = expand.grid(alpha = 1, lambda = seq(0, 20, length = 100)))

plot(lasso_model)

Confrontiamo i coefficienti selezionati.

lasso_model$bestTune$lambda

[1] 10.70707

cbind( coef(lasso_model$finalModel, lasso_model$bestTune$lambda), coef(cv_lasso))

50 x 2 sparse Matrix of class "dgCMatrix"
                                            s=10.70707  lambda.1se
(Intercept)                               7.472305e+01 149.9657571
Water..g.100g.                           -4.185668e-01  -0.6032222
Protein..g.100g.                          1.536153e+01  13.6781669
Protein..crude..N.x.6.25..g.100g.         .              .        
Carbohydrate..g.100g.                     9.716083e+00   4.9960658
Fat..g.100g.                              2.087898e+01  21.1724823
Sugars..g.100g.                          -1.080517e+01  -2.6449047
fructose..g.100g.                         4.098578e+00   .        
galactose..g.100g.                        .              .        
glucose..g.100g.                          1.341604e+01   1.6691924
lactose..g.100g.                          1.700698e+01   .        
maltose..g.100g.                          8.277471e+00   .        
sucrose..g.100g.                         -6.037565e+00  -5.6103001
Starch..g.100g.                           4.672916e+00   7.0525990
Fibres..g.100g.                           .              .        
Polyols..g.100g.                          1.941517e+00   1.1199149
Ash..g.100g.                              4.099386e-02   .        
Alcohol..g.100g.                          2.319661e+01  15.8331160
Organic.acids..g.100g.                    4.130024e+01   .        
FA.saturated..g.100g.                     1.232813e+01   9.6775364
FA.mono..g.100g.                          4.030925e+00   2.3824688
FA.poly..g.100g.                          .              .        
Cholesterol..mg.100g.                     4.800914e-02   .        
Salt..g.100g.                             .              .        
Calcium..mg.100g.                         .              .        
Chloride..mg.100g.                        .              .        
Copper..mg.100g.                         -2.794927e+01   .        
Iron..mg.100g.                            6.204059e+00   3.3243452
Iodine..µg.100g.                          2.964907e-04   .        
Magnesium..mg.100g.                      -1.063540e-01   .        
Manganese..mg.100g.                       6.147647e-01   .        
Phosphorus..mg.100g.                     -1.356200e-02   .        
Potassium..mg.100g.                       .              .        
Selenium..µg.100g.                        .              .        
Sodium..mg.100g.                          .              .        
Zinc..mg.100g.                            .              .        
Retinol..µg.100g.                         .              .        
Beta.carotene..µg.100g.                   .              .        
Vitamin.D..µg.100g.                       2.044689e+00   0.1470879
Vitamin.E..mg.100g.                       1.186171e+01   9.1144274
Vitamin.K1..µg.100g.                      .              .        
Vitamin.K2..µg.100g.                      .              .        
Vitamin.C..mg.100g.                       .              .        
Vitamin.B1.or.Thiamin..mg.100g.          -1.176463e+02   .        
Vitamin.B2.or.Riboflavin..mg.100g.        .              .        
Vitamin.B3.or.Niacin..mg.100g.            8.112183e+00   .        
Vitamin.B5.or.Pantothenic.acid..mg.100g.  .              .        
Vitamin.B6..mg.100g.                      .              .        
Vitamin.B9.or.Folate..µg.100g.            .              .        
Vitamin.B12..µg.100g.                     .              .        

Esercizi suggeriti

1. Studia la multicollinearità nel dataset mtcars usando VIF e numero di condizionamento della matrice di Gram. Applica una riduzione del modello usando stepwise selection (sia forward che backward) e confronta i modelli risultanti in termini di AIC, BIC, VIF e interpretabilità.

2. Utilizza il dataset mtcars e applica Ridge e Lasso regression per prevedere il consumo di carburante (mpg) in base al resto delle variabili. Confronta le prestazioni dei due modelli utilizzando RMSE e visualizza i coefficienti stimati. Discuti quale modello offre una migliore interpretabilità.

3. Discuti la multicollinearità in un modello di regressione polinomiale (anche su dati generati).

4. Genera un dataset in cui la PCR (su un certo numero \(k\) di variabili ruotate) peggiora la performance rispetto a OLS.

5. Utilizza il dataset diamonds per eseguire un’analisi delle componenti principali. Seleziona un numero ridotto di componenti che spiegano una certa percentuale della varianza totale e utilizza queste componenti in un modello di regressione lineare per prevedere il prezzo. Valuta le prestazioni rispetto a un modello di regressione lineare standard.

6. Utilizza bootstrap per stimare gli intervalli di confidenza dei coefficienti di un modello di regressione lineare regolarizzato su un dataset a scelta (ad esempio mtcars o swiss).