Analisi esplorativa dei dati

summary(df)
##      lcavol           lweight           age             lbph        
##  Min.   :-1.3471   Min.   :2.375   Min.   :41.00   Min.   :-1.3863  
##  1st Qu.: 0.5128   1st Qu.:3.376   1st Qu.:60.00   1st Qu.:-1.3863  
##  Median : 1.4469   Median :3.623   Median :65.00   Median : 0.3001  
##  Mean   : 1.3500   Mean   :3.629   Mean   :63.87   Mean   : 0.1004  
##  3rd Qu.: 2.1270   3rd Qu.:3.876   3rd Qu.:68.00   3rd Qu.: 1.5581  
##  Max.   : 3.8210   Max.   :4.780   Max.   :79.00   Max.   : 2.3263  
##     svi               lcp             gleason          pgg45       
##  Mode :logical   Min.   :-1.3863   Min.   :6.000   Min.   :  0.00  
##  FALSE:76        1st Qu.:-1.3863   1st Qu.:6.000   1st Qu.:  0.00  
##  TRUE :21        Median :-0.7985   Median :7.000   Median : 15.00  
##                  Mean   :-0.1794   Mean   :6.753   Mean   : 24.38  
##                  3rd Qu.: 1.1787   3rd Qu.:7.000   3rd Qu.: 40.00  
##                  Max.   : 2.9042   Max.   :9.000   Max.   :100.00  
##       lpsa           train        
##  Min.   :-0.4308   Mode :logical  
##  1st Qu.: 1.7317   FALSE:30       
##  Median : 2.5915   TRUE :67       
##  Mean   : 2.4784                  
##  3rd Qu.: 3.0564                  
##  Max.   : 5.5829

Correlogramma

Modello

Visualizzo la variabile risposta per vedere se si distribuisce come una normale (condizione necessaria per stimare un modello di regressione lineare)

ggplot(df,aes(x = lpsa)) + 
  geom_density() + 
  geom_density(color = "green",
               fill = "green",
               alpha = 0.1, # densità del colore 
               kernel = "gaussian",
               adjust = 2)

Modello lineare

modLin = lm(lpsa ~ .,data = df)
summary(modLin)
## 
## Call:
## lm(formula = lpsa ~ ., data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.76795 -0.35653 -0.00437  0.37978  1.55913 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.177306   1.338810   0.132  0.89495    
## lcavol       0.564417   0.088387   6.386 8.08e-09 ***
## lweight      0.622204   0.202179   3.077  0.00279 ** 
## age         -0.021306   0.011383  -1.872  0.06460 .  
## lbph         0.096833   0.058441   1.657  0.10113    
## sviTRUE      0.761466   0.242697   3.138  0.00233 ** 
## lcp         -0.105872   0.090661  -1.168  0.24609    
## gleason      0.049967   0.158955   0.314  0.75401    
## pgg45        0.004434   0.004485   0.989  0.32558    
## trainTRUE    0.004104   0.162772   0.025  0.97994    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7035 on 87 degrees of freedom
## Multiple R-squared:  0.6634, Adjusted R-squared:  0.6286 
## F-statistic: 19.05 on 9 and 87 DF,  p-value: < 2.2e-16

Modello polinomiale

Creazione dei modelli

Creo \(n*m\) modelli polinomiali con \(n\) numero delle variabili risposta e \(m\) grado massimo con la quale decido di valutare il modello polinomiale. Nella tabella sottostante si possono visualizzare gli \(R^2\) dei modelli.

gradoMax = 10
dfNumeric = df %>%
  select_if(is.numeric)
varIndip = names(df)[names(df)!=c("lpsa","svi","train")]
## Warning in names(df) != c("lpsa", "svi", "train"): la lunghezza più lunga
## dell'oggetto non è un multiplo della lunghezza più corta dell'oggetto
modelli = lapply(X = c(1:length(dfNumeric)), function(col) lapply(X = c(1:gradoMax),
                 function(grado) lm(lpsa ~ polym(unlist(as.vector(dfNumeric[,col])),degree = grado,raw = T),
                                    data = df)))
rDf = map_dfc(c(1:length(dfNumeric)),function(var) sapply(c(1:gradoMax),function(gr) summary(modelli[[var]][[gr]])$r.squared))
## Warning in summary.lm(modelli[[var]][[gr]]): adattamento essenzialmente
## perfetto: il riepilogo potrebbe non essere affidabile
## Warning in summary.lm(modelli[[var]][[gr]]): adattamento essenzialmente
## perfetto: il riepilogo potrebbe non essere affidabile

## Warning in summary.lm(modelli[[var]][[gr]]): adattamento essenzialmente
## perfetto: il riepilogo potrebbe non essere affidabile

## Warning in summary.lm(modelli[[var]][[gr]]): adattamento essenzialmente
## perfetto: il riepilogo potrebbe non essere affidabile

## Warning in summary.lm(modelli[[var]][[gr]]): adattamento essenzialmente
## perfetto: il riepilogo potrebbe non essere affidabile

## Warning in summary.lm(modelli[[var]][[gr]]): adattamento essenzialmente
## perfetto: il riepilogo potrebbe non essere affidabile

## Warning in summary.lm(modelli[[var]][[gr]]): adattamento essenzialmente
## perfetto: il riepilogo potrebbe non essere affidabile

## Warning in summary.lm(modelli[[var]][[gr]]): adattamento essenzialmente
## perfetto: il riepilogo potrebbe non essere affidabile

## Warning in summary.lm(modelli[[var]][[gr]]): adattamento essenzialmente
## perfetto: il riepilogo potrebbe non essere affidabile

## Warning in summary.lm(modelli[[var]][[gr]]): adattamento essenzialmente
## perfetto: il riepilogo potrebbe non essere affidabile
## New names:
## • `` -> `...1`
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
names(rDf) = names(dfNumeric)
rownames(rDf) = c(1:gradoMax)
## Warning: Setting row names on a tibble is deprecated.
rDf
lcavol lweight age lbph lcp gleason pgg45 lpsa
0.5394320 0.1877657 0.0287617 0.0323314 0.3011959 0.1361513 0.1783507 1
0.5416070 0.1931992 0.0398406 0.0356630 0.3019325 0.2292064 0.2432055 1
0.5547384 0.2058068 0.0451996 0.0427544 0.3094325 0.2388101 0.2481146 1
0.5547416 0.2153594 0.0453597 0.0463400 0.3270213 0.2388101 0.2595939 1
0.5581873 0.2153598 0.0555537 0.0465918 0.3272028 0.2388101 0.2671779 1
0.5598660 0.2264730 0.0599538 0.0501817 0.3275041 0.2388101 0.2754673 1
0.5674491 0.2809216 0.0670901 0.0515958 0.3386659 0.2388101 0.2754677 1
0.5709401 0.2818289 0.1024416 0.0550925 0.3458153 0.2388101 0.2816128 1
0.5710294 0.2818289 0.1024416 0.0550929 0.3541366 0.2388101 0.2838632 1
0.5756956 0.2818986 0.1080621 0.0588000 0.3555225 0.2388101 0.2868540 1

Scelta dei modelli

La scelta del modello verrà fatta con l’analisi della varianza mediante la funzione ANOVA che dati dei modelli sceglie il modello meno complesso ma che allo stesso tempo spieghi sufficientemente bene i dati. Il processo verrà eseguito per ciascuna variabile risposta.

do.call(anova,modelli[[1]])
Res.Df RSS Df Sum of Sq F Pr(>F)
95 58.91478 NA NA NA NA
94 58.63656 1 0.2782223 0.4408413 0.5084928
93 56.95683 1 1.6797357 2.6615300 0.1064587
92 56.95641 1 0.0004116 0.0006522 0.9796844
91 56.51565 1 0.4407650 0.6983892 0.4056413
90 56.30091 1 0.2147369 0.3402492 0.5612120
89 55.33090 1 0.9700175 1.5369863 0.2184392
88 54.88433 1 0.4465606 0.7075723 0.4025838
87 54.87291 1 0.0114252 0.0181031 0.8932838
86 54.27602 1 0.5968844 0.9457594 0.3335285

La tabella soprastante si riferisce all’analisi della varianza per la variabile lcavol. Osservando il valore \(p\) dell statistica F si nota che i due polinomi migliori sono il modello cubico e il modello polinomiale di grado settimo.

Di seguito viene applicata la procedura a tutte le variabili risposta. La tabella rappresenta il valore \(p\) per ciascuna variabile e ciascun grado del modello polinomiale.

anovaDf = map_dfc(c(1:length(dfNumeric)),function(var) do.call(anova,modelli[[var]])$`Pr(>F)`)
## New names:
## • `` -> `...1`
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
names(anovaDf) = names(dfNumeric)
rownames(anovaDf) = c(1:gradoMax)
## Warning: Setting row names on a tibble is deprecated.
anovaDf
lcavol lweight age lbph lcp gleason pgg45 lpsa
NA NA NA NA NA NA NA NA
0.5084928 0.4193802 0.3014353 0.5825591 0.7546450 0.0010893 0.0063702 0.0059259
0.1064587 0.2198235 0.4716249 0.4230622 0.3199202 0.2815133 0.4437571 0.0438298
0.9796844 0.2849954 0.9008267 0.5685518 0.1291893 NA 0.2426141 0.6583964
0.4056413 0.9939129 0.3214538 0.8798074 0.8766941 NA 0.3415817 0.1340191
0.5612120 0.2490835 0.5141153 0.5683229 0.8415357 NA 0.3202040 0.0141919
0.2184392 0.0119231 0.4063878 0.7201286 0.2256398 NA 0.9946898 0.0028084
0.4025838 0.7410318 0.0667046 0.5733782 0.3314362 NA 0.3917165 0.0208734
0.8932838 NA NA 0.9954645 0.2949435 NA 0.6037439 0.1475784
0.3335285 0.9270268 0.4610377 0.5620847 0.6682415 NA 0.5497148 0.3377160

Stima modello polinomiale

Sono stati stimati i due modelli polinomiali presi in considerazione precedentemente: quello di grado 3 e quello di grado 7. Rappresentandoli su un grafico si nota come il modello polinomiale di grado 7 si adatti troppo ai dati, specialmente sugli outlier. Sarà preferibile, quindi, il modello polinomiale di grado 3 in modo tale che dia delle stime più affidabili.

ggplot(df,aes(y = lpsa, x = lcavol)) +
  geom_point() +
  geom_smooth(method = "lm",
              formula = y ~ poly(x,3),
              color = "paleturquoise",
              fill = "paleturquoise",
              alpha = 0.3) +
  geom_smooth(method = "lm",
              formula = y ~ poly(x,7),
              color = "orchid",
              fill = "orchid",
              alpha = 0.1)

Grafico dei coefficienti

mod = lm(lpsa ~ poly(lcavol,3), data = df)
arm::coefplot(mod)

Funzione a gradino

La funzione a gradini è composta da 4 intervalli, otteremo così 3 variabili dummy più l’intercetta e le relative stime dei loro coefficienti. L’intercetta corrisponderà al valore stimato della variabile dipendente quando la variabile indipendente appartiene all’intevallo inferiore. Mentre le stime degli altri coefficienti vanno sommate al valore dell’intercetta per ottenere la stima della variabile risposta.

coef(lm(lpsa ~ cut(lcavol,breaks = c(-Inf,0.2,1.5,2.8,Inf)), data = df))
##                                                 (Intercept) 
##                                                   1.2088323 
##  cut(lcavol, breaks = c(-Inf, 0.2, 1.5, 2.8, Inf))(0.2,1.5] 
##                                                   0.9714319 
##  cut(lcavol, breaks = c(-Inf, 0.2, 1.5, 2.8, Inf))(1.5,2.8] 
##                                                   1.6443988 
## cut(lcavol, breaks = c(-Inf, 0.2, 1.5, 2.8, Inf))(2.8, Inf] 
##                                                   2.9944228
ggplot(df,aes(y = lpsa, x = lcavol)) +
  geom_point() +
  geom_smooth(method = "lm",
              formula = y ~ cut(x,breaks = c(-1.1,0.2,1.5,2.8,5)),
              color = "paleturquoise",
              fill = "paleturquoise",
              alpha = 0.3)
## Warning: Removed 4 rows containing missing values (geom_smooth).

##Regressione Lasso

Dal grafico delle stime dei coefficienti risulta evidente quali sono state le variabili portate a zero, ovvero le variabili poco influenti sulla variabile risposta.

modLasso = glmnet(makeX(df[, names(df) != "lpsa"]),
                  df$lpsa,
                  alpha = 1
                    # cv.glmnet(makeX(df[, names(df) != "lpsa"]),
                    #                 df$lpsa,
                    #                 alpha = 1
                    #                 )$lambda.1se
                  )
coefplot::coefplot(modLasso)

Regressione locale

Un valore troppo basso di span risulta addatarsi troppo ai dati andando incontro a overfitting.

modLocale = loess(lpsa ~ lcavol,
                  data = df,
                  span = 0.2
                  )
ggplot(df,aes(y = lpsa, x = lcavol)) +
  geom_point() +
  geom_smooth(method = "loess",
              formula = y ~ x,
              span = 0.2,
              aes(color = "0.2"),
              fill = "paleturquoise",
              alpha = 0.3) +
  geom_smooth(method = "loess",
              formula = y ~ x,
              span = 0.6,
              aes(color = "0.6"),
              fill = "orchid",
              alpha = 0.25) + 
  labs(color = "Span") +
  scale_color_manual(values = c("0.6" = "orchid", "0.2" = "paleturquoise"))