INTRODUZIONE

Lo scopo di questo progetto è lo svuluppo di un modello statistico che possa prevedere con precisione il peso dei neonati alla nascita. La base dati è stata raccolta su 2500 neonati di 3 diversi ospedali.

ANALISI PRELIMINARE

Di seguito verrà svolta una veloce analisi descrittiva delle variabili prese in considerazione.

Summary delle Variabili Numeriche

Statistiche descrittive delle variabili numeriche
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
Min. : 0.00 Min. : 0.0000 Min. :0.0000 Min. :25.00 Min. : 830 Min. :310.0 Min. :235
1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:38.00 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330
Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00 Median :3300 Median :500.0 Median :340
Mean :28.16 Mean : 0.9812 Mean :0.0416 Mean :38.98 Mean :3284 Mean :494.7 Mean :340
3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:40.00 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
Max. :46.00 Max. :12.0000 Max. :1.0000 Max. :43.00 Max. :4930 Max. :565.0 Max. :390

Dai risultati del summary ci si accorge che la variabile “Anni.madre” ha il minimo uguale a 0 (riga 1380), il che è impossibile. Andando a riordinare il dataset secondo quella variabile si nota che anche la riga 1152 presenta un valore nella variabile “Anni.madre” impossibile (cioè 1). In entrambe le righe sono riportate informazioni ragionevoli per le altre variabili; si suppone perciò ci sia stato un errore solamente nella variabile “Anni.madre”; siccome questi dati erronei potrebbero pesare sul risultato dei seguenti passaggi saranno sostituiti con il valore mediano della variabile stessa.

Statistiche descrittive delle variabili numeriche
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
Min. :13.00 Min. : 0.0000 Min. :0.0000 Min. :25.00 Min. : 830 Min. :310.0 Min. :235
1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:38.00 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330
Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00 Median :3300 Median :500.0 Median :340
Mean :28.19 Mean : 0.9812 Mean :0.0416 Mean :38.98 Mean :3284 Mean :494.7 Mean :340
3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:40.00 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
Max. :46.00 Max. :12.0000 Max. :1.0000 Max. :43.00 Max. :4930 Max. :565.0 Max. :390

Summary delle Variabili Categoriali

Statistiche descrittive delle variabili categoriali
Tipo.parto Ospedale Sesso
Ces: 728 osp1:816 F:1256
Nat:1772 osp2:849 M:1244
NA osp3:835 NA

Al fine di approfondire la distribuzione e la presenza di outlier per le variabili “Anni.madre”, “N.gravidanze”, “Gestazione”, “Peso”, “Lunghezza” e “Cranio” saranno utilizzati dei boxplots:

boxplot_var <- c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")

for (i in boxplot_var) {

  graphs <- ggplot(dati, aes_string(y = i)) +
    geom_boxplot() +
    labs(x = i, title = paste("Boxplot di", i))

  print(graphs)
}

Nei boxplot precedenti è riassunta la distribuzione delle variabili e sottolineano la presenza di molti outlier.

Prima di procedere, riassumiamo anche la distribuzione delle variabili qualitative utilizzando dei grafici a barre per dare un indicazione visiva:

barplot_var <- c("Fumatrici", "Tipo.parto", "Ospedale", "Sesso")

for (i in barplot_var) {
  cat("\nTabella di", i, "\n")
  print(table(dati[[i]]))
  cat("Percentuali:\n")
  print(round(prop.table(table(dati[[i]])) * 100, 2))
  
   graphs<- ggplot(dati, aes(x = .data[[i]])) +
    geom_bar(fill = "steelblue") +
    labs(
      x = i,
      y = "Frequenza",
      title = paste("Distribuzione di", i)
    ) +
    theme_minimal()

  print(graphs)
}
## 
## Tabella di Fumatrici 
## 
##    0    1 
## 2396  104 
## Percentuali:
## 
##     0     1 
## 95.84  4.16

## 
## Tabella di Tipo.parto 
## 
##  Ces  Nat 
##  728 1772 
## Percentuali:
## 
##   Ces   Nat 
## 29.12 70.88

## 
## Tabella di Ospedale 
## 
## osp1 osp2 osp3 
##  816  849  835 
## Percentuali:
## 
##  osp1  osp2  osp3 
## 32.64 33.96 33.40

## 
## Tabella di Sesso 
## 
##    F    M 
## 1256 1244 
## Percentuali:
## 
##     F     M 
## 50.24 49.76

Per saggiare l’ipotesi che in alcuni ospedali si fanno più parti cesarei verrà utilizzato il test chi-quadrato:

chisq.test(dati$Ospedale,dati$Tipo.parto)
## 
##  Pearson's Chi-squared test
## 
## data:  dati$Ospedale and dati$Tipo.parto
## X-squared = 1.0972, df = 2, p-value = 0.5778

Dal test chi-quadrato si ottiene un p-value maggiore di 0.05 dunque NON si rifiuta l’ipotesi nulla e perciò non c’è evidenza statistica che in alcuni ospedali ci siano più parti cesarei. Al fine di visualizzare la situazione graficamente in seguito un barplot con la % di parti cesarei rispetto all’ospedale.

perc <- prop.table(table(dati$Ospedale, dati$Tipo.parto),
                   margin = 1)[, "Ces"] * 100

df_plot <- data.frame(
  Ospedale = names(perc),
  Perc_cesarei = as.numeric(perc)
)

ggplot(df_plot, aes(x = Ospedale, y = Perc_cesarei)) +
  geom_col(fill = "steelblue") +
  geom_text(aes(label = paste0(round(Perc_cesarei, 1), "%")),
            vjust = -0.3) +
  ylim(0, 100) +
  labs(
    x = "Ospedale",
    y = "Percentuale di parti cesarei",
    title = "Percentuale di parti cesarei per ospedale"
  ) +
  theme_minimal()

Il grafico conferma il risultato del test chi-quadrato, infatti la percentale di parti cesarei per ospedale è molto simile.

Dalla letteratura, il peso medio e la lunghezza media di un neonato in Italia è di rispettivamente 3255 g e 496 mm (National, longitudinal NASCITA birth cohort study to investigate the health of Italian children and potential influencing factors; C. Pandolfini, A. Clavenna, M. Cartabia, R. Campi, M. Bonati).

Prima di tutto si controlla se i dati di “Peso”, “Lunghezza” e “Cranio” si distribuiscono secondo una distribuzione normale tramite lo Shapiro test:

variabili <- c("Peso", "Lunghezza", "Cranio")

shapiro_tab <- do.call(
  rbind,
  lapply(variabili, function(v) {
    test <- shapiro.test(dati[[v]])
    data.frame(
      Variabile = v,
      p_value = test$p.value
    )
  })
)

shapiro_tab$p_value <- format.pval(
  shapiro_tab$p_value,
  digits = 3,
  eps = 0.001
)

kable(
  shapiro_tab,
  caption = "Shapito Test per la normalità"
)
Shapito Test per la normalità
Variabile p_value
Peso <0.001
Lunghezza <0.001
Cranio <0.001

Al fine di saggiare l’ipotesi che la media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione utilizzerò il t-test poichè la distribuzione dei dati non è normale come si può vedere dai risultati dello Shapiro-Wilk test.

test_peso <- t.test(dati$Peso, mu = 3255)
test_lung <- t.test(dati$Lunghezza, mu = 496)

test_tab <- data.frame(
  Variabile = c("Peso", "Lunghezza"),
  Media = c(mean(dati$Peso, na.rm = TRUE),
            mean(dati$Lunghezza, na.rm = TRUE)),
  mu = c(3255, 496),
  p_value = c(test_peso$p.value,
              test_lung$p.value)
)

test_tab$p_value <- format.pval(
  test_tab$p_value,
  digits = 3,
  eps = 0.001
)

kable(
  test_tab,
  digits = 2,
  caption = "t-test a un campione"
)
t-test a un campione
Variabile Media mu p_value
Peso 3284.08 3255 0.00566
Lunghezza 494.69 496 0.01302

Il p-value risultante da entrambi i t-test è minore di 0.05 perciò si rifiuta l’ipotesi nulla e quindi c’è evidenza statistica che la media del nostro campione è diversa da quella della popolazione (valore da letteratura). In particolare il peso medio della popolazione risulta minore all’estremo inferiore dell’intervallo di confidenza dei pesi del nostro campione mentre la lunghezza media della popolazione risulta maggiore all’estremo superiore dell’intervallo di confidenza del nostro campione.

Al fine di saggiare l’ipotesi che le misure antropometriche sono significativamente diverse tra i due sessi si utilizzerà il t-test per confronto tra gruppi indipendenti; anche se i campioni non seguono una distribuzione normale si può procedere considerando la numerosità del campione e dunque l’applicabilità del test è assicurata dal teorema dei limiti centrali.

test_peso  <- t.test(Peso ~ Sesso, data = dati)
test_lung  <- t.test(Lunghezza ~ Sesso, data = dati)
test_cranio <- t.test(Cranio ~ Sesso, data = dati)

test_tab <- data.frame(
  Variabile = c("Peso", "Lunghezza", "Cranio"),
  p_value = c(
    test_peso$p.value,
    test_lung$p.value,
    test_cranio$p.value
  )
)

test_tab$p_value <- format.pval(
  test_tab$p_value,
  digits = 3,
  eps = 0.001
)

kable(
  test_tab,
  caption = "t-test a due campioni (Sesso)"
)
t-test a due campioni (Sesso)
Variabile p_value
Peso <0.001
Lunghezza <0.001
Cranio <0.001

Tutti i t-test per campioni indipendenti riportano un p-value minore di 0.05 dunque si rifiuta l’ipotesi nulla e perciò abbiamo la conferma che le misure antropometriche sono significativamente diverse tra i due sessi.

Prima di procedere con la creazione del modello vediamo graficamente tramite alcuni scatterplot se c’è correlazione tra le variabili numeriche e la variabile “Peso”, per le stesse variabili viene fatta anche la matrice di correlazione.

plot(dati$Peso,dati$Anni.madre)

plot(dati$Peso,dati$N.gravidanze)

plot(dati$Peso,dati$Gestazione)

plot(dati$Peso,dati$Lunghezza)

plot(dati$Peso,dati$Cranio)

idx <- which(names(dati) %in% barplot_var)
dati_num <- dati[, -idx, drop = FALSE]

cor_mat <- cor(dati_num, use = "complete.obs", method = "pearson")
cor_mat <- round(cor_mat, 2)

kable(
  cor_mat,
  caption = "Matrice di correlazione di Pearson"
)
Matrice di correlazione di Pearson
Anni.madre N.gravidanze Gestazione Peso Lunghezza Cranio
Anni.madre 1.00 0.38 -0.13 -0.02 -0.06 0.02
N.gravidanze 0.38 1.00 -0.10 0.00 -0.06 0.04
Gestazione -0.13 -0.10 1.00 0.59 0.62 0.46
Peso -0.02 0.00 0.59 1.00 0.80 0.70
Lunghezza -0.06 -0.06 0.62 0.80 1.00 0.60
Cranio 0.02 0.04 0.46 0.70 0.60 1.00

Come si può notare sia dagli scatterplot che dalla matrice di correlazione le variabile “Anni.madre” e “N.gravidanze” sono indipendenti dalla variabile “Peso”; mentre per le altre variabili (“Gestazione”, “Lunghezza” e “Cranio”) si possono notare trend positivi e perciò sono correlate positivamente con coefficienti di correlazione >0.5.

Di seguito, per la valutazione delle variabili qualitative, si vedranno dei boxplot e dei test di significatività tra gruppi.

dati$Fumatrici <- factor(dati$Fumatrici, levels = c(0, 1), labels = c("No", "Sì"))

vars <- c("Fumatrici", "Tipo.parto", "Ospedale", "Sesso")

for (v in vars) {

  # Se è numerica con soli 2 valori (0/1), convertila in factor
  if (is.numeric(dati[[v]]) && length(unique(dati[[v]])) == 2) {
    dati[[v]] <- factor(dati[[v]])
  }

  p <- ggplot(dati, aes(x = .data[[v]], y = Peso)) +
    geom_boxplot() +
    labs(
      title = paste("Peso vs", v),
      x = v,
      y = "Peso"
    ) +
    theme_minimal()

  if (v == "Ospedale") {
    p <- p + theme(axis.text.x = element_text(angle = 45, hjust = 1))
  }

  print(p)
}

dati$Fumatrici <- factor(dati$Fumatrici, levels = c("No", "Sì"), labels = c(0, 1))

# dati dei test
risultati <- data.frame(
  Test = c(
    "Ospedale (ANOVA)",
    "Fumatrici (t-test)",
    "Sesso (t-test)",
    "Tipo parto (t-test)"
  ),
  p_value = c(0.183, 0.3033, 2.2e-16, 0.8968)
)

# formatta p-value leggibile
risultati$p_value <- format.pval(
  risultati$p_value,
  digits = 3,
  eps = 0.001
)

kable(
  risultati,
  caption = "Risultati dei test per la variabile Peso"
)
Risultati dei test per la variabile Peso
Test p_value
Ospedale (ANOVA) 0.183
Fumatrici (t-test) 0.303
Sesso (t-test) <0.001
Tipo parto (t-test) 0.897

Si noti che la variabile “Sesso” è l’unica in cui la differenza di peso tra i gruppi è statisticamente significativa dunque sarà l’unica, tra queste, che ci si aspetta di trovare nel modello finale.

Al fine di sviluppare un modello, si partirà dalla crezione di questo inserendo tutte le variabile anche se, come abbiamo visto, alcune feature sono ininfluenti sul nostro target o logicamente da escludere (ad esempio “Ospedale”). Fatto ciò si procederà con l’ottimizzazione del modello (si proverà ad utilizzare anche la procedura automatica del pacchetto “mass”):

mod1 <- lm(Peso ~ ., data=dati)
summary(mod1)
## 
## Call:
## lm(formula = Peso ~ ., data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1123.3  -181.2   -14.6   160.7  2612.6 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6735.1677   141.3977 -47.633  < 2e-16 ***
## Anni.madre        0.7983     1.1463   0.696   0.4862    
## N.gravidanze     11.4118     4.6665   2.445   0.0145 *  
## Fumatrici1      -30.1567    27.5396  -1.095   0.2736    
## Gestazione       32.5265     3.8179   8.520  < 2e-16 ***
## Lunghezza        10.2951     0.3007  34.237  < 2e-16 ***
## Cranio           10.4725     0.4261  24.580  < 2e-16 ***
## Tipo.partoNat    29.5027    12.0848   2.441   0.0147 *  
## Ospedaleosp2    -11.2216    13.4388  -0.835   0.4038    
## Ospedaleosp3     28.0984    13.4972   2.082   0.0375 *  
## SessoM           77.5473    11.1779   6.938 5.07e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared:  0.7289, Adjusted R-squared:  0.7278 
## F-statistic: 669.1 on 10 and 2489 DF,  p-value: < 2.2e-16
res <- residuals(mod1)

tabella_indici <- data.frame(
  Indice = c("R²", "MSE", "RMSE", "AIC", "BIC"),
  Valore = c(
    sprintf("%.3f", round(summary(mod1)$r.squared, 3)),
    sprintf("%.0f", round(mean(res^2), 0)),
    sprintf("%.0f", round(sqrt(mean(res^2)), 0)),
    sprintf("%.0f", round(AIC(mod1), 0)),
    sprintf("%.0f", round(BIC(mod1), 0))
  )
)

kable(tabella_indici,
      digits = 3,
      caption = "Metriche di valutazione del modello 1")
Metriche di valutazione del modello 1
Indice Valore
R² 0.729
MSE 74709
RMSE 273
AIC 35172
BIC 35242

Dal summary del modello generato possiamo vedere che le variabili quantitative significative sono “Lunghezza”, “Cranio” e “Gestazione” (come si era già visto dagli scatterplot e dalla matrice di correlazione) e la variabile qualitativa “Sesso” (come visto dal boxplot e dal t-test). Si ottiene un R2 complessivo del 73% ed un RMSE di 273; inoltre i valori di AIC e BIC sono rispettivamente 35172 e 35242. Come punto di partenza non è male, ma nella fase di miglioramento del modello si procederà col rimuovere le variabili superflue al fine di ridurre la complessità del modello mantenendone alta la performance poichè al momento le variabili considerate sono 9.

Si procede perciò con la rimozione delle variabili “Anni.madre”, “Fumatrici” e “Ospedale”. Si lasciano momentaneamente le variabili “N.Gravidanze” e “Tipo.parto”.

mod2 <- update(mod1,.~.-Fumatrici - Anni.madre -Ospedale)
summary(mod2)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1129.31  -181.70   -16.31   161.07  2638.85 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6707.2971   135.9911 -49.322  < 2e-16 ***
## N.gravidanze     12.7558     4.3366   2.941   0.0033 ** 
## Gestazione       32.2713     3.7941   8.506  < 2e-16 ***
## Lunghezza        10.2864     0.3007  34.207  < 2e-16 ***
## Cranio           10.5057     0.4260  24.659  < 2e-16 ***
## Tipo.partoNat    30.0342    12.0969   2.483   0.0131 *  
## SessoM           77.9285    11.1905   6.964 4.22e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.3 on 2493 degrees of freedom
## Multiple R-squared:  0.7277, Adjusted R-squared:  0.727 
## F-statistic:  1110 on 6 and 2493 DF,  p-value: < 2.2e-16
res <- residuals(mod2)

tabella_indici <- data.frame(
  Indice = c("R²", "MSE", "RMSE", "AIC", "BIC"),
  Valore = c(
    sprintf("%.3f", round(summary(mod2)$r.squared, 3)),
    sprintf("%.0f", round(mean(res^2), 0)),
    sprintf("%.0f", round(sqrt(mean(res^2)), 0)),
    sprintf("%.0f", round(AIC(mod2), 0)),
    sprintf("%.0f", round(BIC(mod2), 0))
  )
)

kable(tabella_indici,
      digits = 3,
      caption = "Metriche di valutazione del modello 2")
Metriche di valutazione del modello 2
Indice Valore
R² 0.728
MSE 75041
RMSE 274
AIC 35175
BIC 35222

Il nuovo modello mantiene un R2 del 73% e un RMSE di 274 però con meno feature e perciò risultando meno complesso; è stato deciso di mantenere la variabile “Tipo.parto” per vedere il suo comportamente con la rimozione delle altre variabili. Inoltre, non si nota una grossa differenza tra il modello 1 e il modello 2: nel caso dell’AIC è leggermente più performante il modello 1 mentre nel caso del BIC il contrario, in questo caso si tende a preferire il modello meno complesso perciò il modello 2.

Di seguito una valutazione della multicollinearità si nota che nessun valore supera il valore soglia di 5 perciò non c’è multicollinearità.

car::vif(mod2)
## N.gravidanze   Gestazione    Lunghezza       Cranio   Tipo.parto        Sesso 
##     1.024171     1.669258     2.080039     1.626199     1.003438     1.040060

Si continua con l’ottimizzazione: prima togliendo la variabile “Tipo.parto” e valutando il modello e poi togliendo la variabile “N.gravidanze” e valutando il modello.

mod3 <- update(mod2,.~.-Tipo.parto)
summary(mod3)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.44  -180.81   -15.58   163.64  2639.72 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.1445   135.7229 -49.226  < 2e-16 ***
## N.gravidanze    12.4750     4.3396   2.875  0.00408 ** 
## Gestazione      32.3321     3.7980   8.513  < 2e-16 ***
## Lunghezza       10.2486     0.3006  34.090  < 2e-16 ***
## Cranio          10.5402     0.4262  24.728  < 2e-16 ***
## SessoM          77.9927    11.2021   6.962 4.26e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1328 on 5 and 2494 DF,  p-value: < 2.2e-16
res <- residuals(mod3)

tabella_indici <- data.frame(
  Indice = c("R²", "MSE", "RMSE", "AIC", "BIC"),
  Valore = c(
    sprintf("%.3f", round(summary(mod3)$r.squared, 3)),
    sprintf("%.0f", round(mean(res^2), 0)),
    sprintf("%.0f", round(sqrt(mean(res^2)), 0)),
    sprintf("%.0f", round(AIC(mod3), 0)),
    sprintf("%.0f", round(BIC(mod3), 0))
  )
)

kable(tabella_indici,
      digits = 3,
      caption = "Metriche di valutazione del modello 3")
Metriche di valutazione del modello 3
Indice Valore
R² 0.727
MSE 75226
RMSE 274
AIC 35179
BIC 35220

R2 continua a rimanere vicino al 73% e l’RMSE 274, inoltre anche i valori di AIC e BIC continuano a rimanere stabili; dunque si preferisce il modello 3 al 2 poichè meno complesso.

mod4 <- update(mod3,.~.-N.gravidanze)
summary(mod4)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso, 
##     data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.2  -184.3   -17.6   163.3  2627.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6651.1188   135.5172 -49.080  < 2e-16 ***
## Gestazione     31.2737     3.7856   8.261 2.31e-16 ***
## Lunghezza      10.2054     0.3007  33.939  < 2e-16 ***
## Cranio         10.6704     0.4245  25.139  < 2e-16 ***
## SessoM         79.1049    11.2117   7.056 2.22e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared:  0.7261, Adjusted R-squared:  0.7257 
## F-statistic:  1654 on 4 and 2495 DF,  p-value: < 2.2e-16
res <- residuals(mod4)

tabella_indici <- data.frame(
  Indice = c("R²", "MSE", "RMSE", "AIC", "BIC"),
  Valore = c(
    sprintf("%.3f", round(summary(mod4)$r.squared, 3)),
    sprintf("%.0f", round(mean(res^2), 0)),
    sprintf("%.0f", round(sqrt(mean(res^2)), 0)),
    sprintf("%.0f", round(AIC(mod4), 0)),
    sprintf("%.0f", round(BIC(mod4), 0))
  )
)

kable(tabella_indici,
      digits = 3,
      caption = "Metriche di valutazione del modello 4")
Metriche di valutazione del modello 4
Indice Valore
R² 0.726
MSE 75475
RMSE 275
AIC 35186
BIC 35221

Il modello 4 continua a mantenere pressocchè invariati i valori di R2 intorno al 73%, RMSE di 275, AIC e BIC. Si preferisce dunque il modello 4 al 3 essendo ulteriormente semplificato.

Si mettono a confronto ora le perfomance dei 4 modelli.

modelli <- list(mod1, mod2, mod3, mod4)  # inserisci i tuoi modelli lm qui
nomi_modelli <- c("Modello 1", "Modello 2", "Modello 3", "Modello 4")  # puoi cambiare i nomi

num_var_modelli <- c(9, 6, 5, 4)  # modifica qui con il numero di predittori che vuoi per ciascun modello

estrai_indici <- function(mod, n_var) {
  R2   <- summary(mod)$r.squared
  RSE  <- summary(mod)$sigma
  AICv <- AIC(mod)
  BICv <- BIC(mod)
  
  data.frame(
    R2           = round(R2, 3),
    RMSE         = round(RSE, 0),
    AIC          = round(AICv, 0),
    BIC          = round(BICv, 0),
    Num_Variabili = n_var
  )
}

tabella_modelli <- do.call(
  rbind,
  Map(estrai_indici, modelli, num_var_modelli)
)

tabella_modelli$Modello <- nomi_modelli
tabella_modelli <- tabella_modelli[, c("Modello", "Num_Variabili", "R2", "RMSE", "AIC", "BIC")]

kable(tabella_modelli,
      caption = "Confronto tra modelli lineari",
      align = "c")
Confronto tra modelli lineari
Modello Num_Variabili R2 RMSE AIC BIC
Modello 1 9 0.729 274 35172 35242
Modello 2 6 0.728 274 35175 35222
Modello 3 5 0.727 275 35179 35220
Modello 4 4 0.726 275 35186 35221

Dalla tabella riassuntiva si nota come i valori di R2, RMSE, AIC e BIC varino molto poco in confronto alla grande riduzione di complessività del modello ottenuta rimuovendo le variabili superflue. Dunque si preferisce il modello più semplice cioè il modello 4.

Al fine di confermare viene provata la procedura stepwise automatica del pacchetto “MASS” sia con metodo AIC (k=2), sia con metodo BIC (k=log(2500)).

mod_aic <- MASS::stepAIC(mod1, direction = "both", k = 2)
## Start:  AIC=28075.39
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     36392 186809099 28074
## - Fumatrici     1     89979 186862686 28075
## <none>                      186772707 28075
## - Tipo.parto    1    447233 187219939 28079
## - N.gravidanze  1    448762 187221469 28079
## - Ospedale      2    686397 187459103 28081
## - Sesso         1   3611588 190384294 28121
## - Gestazione    1   5446558 192219264 28145
## - Cranio        1  45338051 232110758 28617
## - Lunghezza     1  87959790 274732497 29038
## 
## Step:  AIC=28073.88
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     90897 186899996 28073
## <none>                      186809099 28074
## + Anni.madre    1     36392 186772707 28075
## - Tipo.parto    1    448222 187257321 28078
## - Ospedale      2    692738 187501837 28079
## - N.gravidanze  1    633756 187442855 28080
## - Sesso         1   3618736 190427835 28120
## - Gestazione    1   5412879 192221978 28143
## - Cranio        1  45588236 232397335 28618
## - Lunghezza     1  87950050 274759149 29036
## 
## Step:  AIC=28073.1
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      186899996 28073
## + Fumatrici     1     90897 186809099 28074
## + Anni.madre    1     37311 186862686 28075
## - Tipo.parto    1    440684 187340680 28077
## - Ospedale      2    701680 187601677 28079
## - N.gravidanze  1    610840 187510837 28079
## - Sesso         1   3602797 190502794 28119
## - Gestazione    1   5346781 192246777 28142
## - Cranio        1  45632149 232532146 28617
## - Lunghezza     1  88355030 275255027 29039
summary(mod_aic)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1113.18  -181.16   -16.58   161.01  2620.19 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6707.4293   135.9438 -49.340  < 2e-16 ***
## N.gravidanze     12.3619     4.3325   2.853  0.00436 ** 
## Gestazione       31.9909     3.7896   8.442  < 2e-16 ***
## Lunghezza        10.3086     0.3004  34.316  < 2e-16 ***
## Cranio           10.4922     0.4254  24.661  < 2e-16 ***
## Tipo.partoNat    29.2803    12.0817   2.424  0.01544 *  
## Ospedaleosp2    -11.0227    13.4363  -0.820  0.41209    
## Ospedaleosp3     28.6408    13.4886   2.123  0.03382 *  
## SessoM           77.4412    11.1756   6.930 5.36e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared:  0.7287, Adjusted R-squared:  0.7278 
## F-statistic: 836.3 on 8 and 2491 DF,  p-value: < 2.2e-16
mod_bic <- MASS::stepAIC(mod1, direction = "both", k = log(2500))
## Start:  AIC=28139.46
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     36392 186809099 28132
## - Fumatrici     1     89979 186862686 28133
## - Ospedale      2    686397 187459103 28133
## - Tipo.parto    1    447233 187219939 28138
## - N.gravidanze  1    448762 187221469 28138
## <none>                      186772707 28140
## - Sesso         1   3611588 190384294 28180
## - Gestazione    1   5446558 192219264 28204
## - Cranio        1  45338051 232110758 28675
## - Lunghezza     1  87959790 274732497 29096
## 
## Step:  AIC=28132.12
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     90897 186899996 28126
## - Ospedale      2    692738 187501837 28126
## - Tipo.parto    1    448222 187257321 28130
## <none>                      186809099 28132
## - N.gravidanze  1    633756 187442855 28133
## + Anni.madre    1     36392 186772707 28140
## - Sesso         1   3618736 190427835 28172
## - Gestazione    1   5412879 192221978 28196
## - Cranio        1  45588236 232397335 28670
## - Lunghezza     1  87950050 274759149 29089
## 
## Step:  AIC=28125.51
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Ospedale      2    701680 187601677 28119
## - Tipo.parto    1    440684 187340680 28124
## <none>                      186899996 28126
## - N.gravidanze  1    610840 187510837 28126
## + Fumatrici     1     90897 186809099 28132
## + Anni.madre    1     37311 186862686 28133
## - Sesso         1   3602797 190502794 28165
## - Gestazione    1   5346781 192246777 28188
## - Cranio        1  45632149 232532146 28664
## - Lunghezza     1  88355030 275255027 29086
## 
## Step:  AIC=28119.23
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Tipo.parto    1    463870 188065546 28118
## <none>                      187601677 28119
## - N.gravidanze  1    651066 188252743 28120
## + Ospedale      2    701680 186899996 28126
## + Fumatrici     1     99840 187501837 28126
## + Anni.madre    1     43845 187557831 28127
## - Sesso         1   3649259 191250936 28160
## - Gestazione    1   5444109 193045786 28183
## - Cranio        1  45758101 233359778 28657
## - Lunghezza     1  88054432 275656108 29074
## 
## Step:  AIC=28117.58
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      188065546 28118
## - N.gravidanze  1    623141 188688687 28118
## + Tipo.parto    1    463870 187601677 28119
## + Ospedale      2    724866 187340680 28124
## + Fumatrici     1     91892 187973654 28124
## + Anni.madre    1     45044 188020502 28125
## - Sesso         1   3655292 191720838 28158
## - Gestazione    1   5464853 193530399 28181
## - Cranio        1  46108583 234174130 28658
## - Lunghezza     1  87632762 275698308 29066
summary(mod_bic)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.44  -180.81   -15.58   163.64  2639.72 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.1445   135.7229 -49.226  < 2e-16 ***
## N.gravidanze    12.4750     4.3396   2.875  0.00408 ** 
## Gestazione      32.3321     3.7980   8.513  < 2e-16 ***
## Lunghezza       10.2486     0.3006  34.090  < 2e-16 ***
## Cranio          10.5402     0.4262  24.728  < 2e-16 ***
## SessoM          77.9927    11.2021   6.962 4.26e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1328 on 5 and 2494 DF,  p-value: < 2.2e-16
modelli <- list(mod_aic, mod_bic, mod4)
nomi_modelli <- c("Step AIC", "Step BIC", "Modello 4")  # nomi da visualizzare in tabella

num_var_modelli <- c(7, 5, 4)  # inserisci i numeri corretti per ciascun modello

estrai_indici <- function(mod, n_var) {
  R2   <- summary(mod)$r.squared
  RSE  <- summary(mod)$sigma
  AICv <- AIC(mod)
  BICv <- BIC(mod)
  
  data.frame(
    R2           = round(R2, 3),
    RMSE         = round(RSE, 0),
    AIC          = round(AICv, 0),
    BIC          = round(BICv, 0),
    Num_Variabili = n_var
  )
}

tabella_modelli <- do.call(
  rbind,
  Map(estrai_indici, modelli, num_var_modelli)
)

tabella_modelli$Modello <- nomi_modelli
tabella_modelli <- tabella_modelli[, c("Modello", "Num_Variabili", "R2", "RMSE", "AIC", "BIC")]

kable(tabella_modelli,
      caption = "Confronto tra modelli: Step AIC, Step BIC e Modello 4",
      align = "c")
Confronto tra modelli: Step AIC, Step BIC e Modello 4
Modello Num_Variabili R2 RMSE AIC BIC
Step AIC 7 0.729 274 35170 35228
Step BIC 5 0.727 275 35179 35220
Modello 4 4 0.726 275 35186 35221

In entrambi i casi non si arriva ad una conclusione uguale al modello 4, ma si può notare che il benificio in termini di R2, RMSE, AIC e BIC è trascurabile in confronto alla semplicità del modello 4 rispetto agli altri.

Si procede ora col valutare le interazione tra variabili ed eventuali effetti quadratici; in seguito all’osservazione degli scatterplot sembra esserci qualche interazione quadratica, più precisamente lo scatterplot delle variabili “Lunghezza” e “Cranio” assomiglia all’andamento di una radice quadratica; di seguito proverò comunque varie combinazioni:

mod5 <- update(mod4,.~.+ I(Lunghezza^2))
summary(mod4)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso, 
##     data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.2  -184.3   -17.6   163.3  2627.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6651.1188   135.5172 -49.080  < 2e-16 ***
## Gestazione     31.2737     3.7856   8.261 2.31e-16 ***
## Lunghezza      10.2054     0.3007  33.939  < 2e-16 ***
## Cranio         10.6704     0.4245  25.139  < 2e-16 ***
## SessoM         79.1049    11.2117   7.056 2.22e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared:  0.7261, Adjusted R-squared:  0.7257 
## F-statistic:  1654 on 4 and 2495 DF,  p-value: < 2.2e-16
summary(mod5)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Lunghezza^2), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1156.75  -182.00   -12.52   166.67  1783.83 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    156.673587 724.783055   0.216    0.829    
## Gestazione      41.174410   3.860512  10.666  < 2e-16 ***
## Lunghezza      -19.913124   3.165788  -6.290 3.73e-10 ***
## Cranio          10.795854   0.417182  25.878  < 2e-16 ***
## SessoM          71.369249  11.043854   6.462 1.24e-10 ***
## I(Lunghezza^2)   0.031241   0.003269   9.555  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 270.2 on 2494 degrees of freedom
## Multiple R-squared:  0.7358, Adjusted R-squared:  0.7352 
## F-statistic:  1389 on 5 and 2494 DF,  p-value: < 2.2e-16
mod6 <- update(mod4,.~.+ I(Cranio^2))
summary(mod6)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Cranio^2), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1127.21  -183.53   -14.88   163.98  2610.62 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   76.70852 1153.12221   0.067    0.947    
## Gestazione    37.73144    3.91778   9.631  < 2e-16 ***
## Lunghezza     10.44511    0.30147  34.647  < 2e-16 ***
## Cranio       -31.41831    7.17690  -4.378 1.25e-05 ***
## SessoM        74.30205   11.16712   6.654 3.50e-11 ***
## I(Cranio^2)    0.06226    0.01060   5.875 4.80e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.2 on 2494 degrees of freedom
## Multiple R-squared:  0.7298, Adjusted R-squared:  0.7293 
## F-statistic:  1347 on 5 and 2494 DF,  p-value: < 2.2e-16
mod7 <- update(mod6,.~.+ I(Lunghezza^2))
summary(mod7)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Cranio^2) + I(Lunghezza^2), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1157.40  -181.85   -11.74   166.59  1772.46 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.150e+01  1.141e+03   0.010    0.992    
## Gestazione      4.108e+01  3.901e+00  10.531  < 2e-16 ***
## Lunghezza      -2.035e+01  4.125e+00  -4.933 8.62e-07 ***
## Cranio          1.231e+01  9.194e+00   1.339    0.181    
## SessoM          7.143e+01  1.105e+01   6.463 1.23e-10 ***
## I(Cranio^2)    -2.237e-03  1.357e-02  -0.165    0.869    
## I(Lunghezza^2)  3.168e-02  4.233e-03   7.485 9.86e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 270.2 on 2493 degrees of freedom
## Multiple R-squared:  0.7358, Adjusted R-squared:  0.7351 
## F-statistic:  1157 on 6 and 2493 DF,  p-value: < 2.2e-16
mod8 <- update(mod4,.~.+ sqrt(Lunghezza))
summary(mod8)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     sqrt(Lunghezza), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1156.88  -182.35   -11.26   166.34  1697.87 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     19654.0735  2662.6953   7.381 2.12e-13 ***
## Gestazione         41.0694     3.8439  10.684  < 2e-16 ***
## Lunghezza          66.0294     5.6513  11.684  < 2e-16 ***
## Cranio             10.7914     0.4166  25.902  < 2e-16 ***
## SessoM             71.0596    11.0303   6.442 1.41e-10 ***
## sqrt(Lunghezza) -2444.0637   247.0873  -9.891  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.8 on 2494 degrees of freedom
## Multiple R-squared:  0.7364, Adjusted R-squared:  0.7359 
## F-statistic:  1394 on 5 and 2494 DF,  p-value: < 2.2e-16
mod9 <- update(mod4,.~.+ sqrt(Cranio))
summary(mod9)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     sqrt(Cranio), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1128.10  -182.32   -14.34   162.75  2615.52 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  20141.8679  4414.6940   4.562 5.30e-06 ***
## Gestazione      37.9841     3.9178   9.695  < 2e-16 ***
## Lunghezza       10.4536     0.3013  34.690  < 2e-16 ***
## Cranio          91.3317    13.2911   6.872 7.99e-12 ***
## SessoM          74.0594    11.1629   6.634 3.98e-11 ***
## sqrt(Cranio) -2961.9675   487.8181  -6.072 1.46e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273 on 2494 degrees of freedom
## Multiple R-squared:  0.7301, Adjusted R-squared:  0.7295 
## F-statistic:  1349 on 5 and 2494 DF,  p-value: < 2.2e-16
mod10 <- update(mod8,.~.+ sqrt(Cranio))
summary(mod10)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     sqrt(Lunghezza) + sqrt(Cranio), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1158.26  -182.58   -11.34   165.84  1669.34 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     18382.784   4369.099   4.207 2.67e-05 ***
## Gestazione         40.853      3.890  10.503  < 2e-16 ***
## Lunghezza          67.779      7.394   9.166  < 2e-16 ***
## Cranio              4.458     17.260   0.258    0.796    
## SessoM             71.201     11.039   6.450 1.34e-10 ***
## sqrt(Lunghezza) -2521.529    324.987  -7.759 1.24e-14 ***
## sqrt(Cranio)      232.712    634.022   0.367    0.714    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.9 on 2493 degrees of freedom
## Multiple R-squared:  0.7365, Adjusted R-squared:  0.7358 
## F-statistic:  1161 on 6 and 2493 DF,  p-value: < 2.2e-16
mod11 <- update(mod4,.~.+ Lunghezza*Cranio)
summary(mod11)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     Lunghezza:Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1139.05  -181.44   -15.46   163.32  2850.10 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.839e+03  1.019e+03  -1.804   0.0714 .  
## Gestazione        3.692e+01  3.951e+00   9.344  < 2e-16 ***
## Lunghezza        -2.013e-01  2.205e+00  -0.091   0.9273    
## Cranio           -4.409e+00  3.194e+00  -1.380   0.1676    
## SessoM            7.449e+01  1.121e+01   6.648 3.64e-11 ***
## Lunghezza:Cranio  3.113e-02  6.537e-03   4.763 2.01e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.8 on 2494 degrees of freedom
## Multiple R-squared:  0.7286, Adjusted R-squared:  0.728 
## F-statistic:  1339 on 5 and 2494 DF,  p-value: < 2.2e-16
# Lista modelli
modelli <- list(
  "Modello 4" = mod4,
  "Modello 5" = mod5,
  "Modello 6" = mod6,
  "Modello 7" = mod7,
  "Modello 8" = mod8,
  "Modello 9" = mod9,
  "Modello 10" = mod10,
  "Modello 11" = mod11
)

# Funzione per estrarre metriche (p-value globale modello + R2)
estrai_metriche <- function(mod) {
  rmse <- sqrt(mean(residuals(mod)^2))
  f_stat <- summary(mod)$fstatistic
  p_val <- pf(f_stat[1], f_stat[2], f_stat[3], lower.tail = FALSE)
  r2 <- summary(mod)$r.squared
  
  data.frame(
    RMSE = rmse,
    AIC = AIC(mod),
    BIC = BIC(mod),
    p_value = p_val,
    R2 = r2
  )
}

# Costruiamo la tabella
tabella_modelli <- do.call(rbind, lapply(modelli, estrai_metriche))

# Descrizioni manuali
descrizioni <- c(
  "Gestazione + Lunghezza + Cranio + Sesso",
  "Modello 4 + Lunghezza^2",
  "Modello 4 + Cranio^2",
  "Modello 4 + Lunghezza^2 + Cranio^2",
  "Modello 4 + sqrt(Lunghezza)",
  "Modello 4 + sqrt(Cranio)",
  "Modello 4 + sqrt(Lunghezza) + sqrt(Cranio)",
  "Modello 4 + l'interazione tra Lunghezza e Cranio"
)

# Aggiungiamo la colonna Descrizione e rimuoviamo la colonna Modello ridondante (nomi riga)
tabella_modelli <- cbind(
  Descrizione = descrizioni,
  tabella_modelli
)

# Gestione p_value molto piccolo
tabella_modelli$p_value <- ifelse(
  tabella_modelli$p_value < 1e-10,
  "<1e-10",
  signif(tabella_modelli$p_value, 3)
)

# Visualizziamo con kable o print
if(requireNamespace("knitr", quietly = TRUE)) {
  knitr::kable(tabella_modelli, digits = 3, caption = "Confronto tra modelli")
} else {
  print(tabella_modelli)
}
Confronto tra modelli
Descrizione RMSE AIC BIC p_value R2
Modello 4 Gestazione + Lunghezza + Cranio + Sesso 274.728 35185.60 35220.54 <1e-10 0.726
Modello 5 Modello 4 + Lunghezza^2 269.833 35097.71 35138.48 <1e-10 0.736
Modello 6 Modello 4 + Cranio^2 272.847 35153.24 35194.01 <1e-10 0.730
Modello 7 Modello 4 + Lunghezza^2 + Cranio^2 269.832 35099.68 35146.28 <1e-10 0.736
Modello 8 Modello 4 + sqrt(Lunghezza) 269.493 35091.40 35132.17 <1e-10 0.736
Modello 9 Modello 4 + sqrt(Cranio) 272.720 35150.91 35191.68 <1e-10 0.730
Modello 10 Modello 4 + sqrt(Lunghezza) + sqrt(Cranio) 269.485 35093.26 35139.86 <1e-10 0.736
Modello 11 Modello 4 + l’interazione tra Lunghezza e Cranio 273.487 35164.96 35205.73 <1e-10 0.729

Si può notare che il beneficio su R2, AIC e BIC non è rilevante perciò si opterà per tenere il modello 4 con un R2 del 72.6% (RMSE=274.73 g) avendo una performance simile agli altri modelli ottenuti ma una semplicità maggiore. Si sottolinea comunque il fatto che un RMSE di 275 grammi circa significa un errore potenziale di più o meno questa cifra che comparate con il peso medio (3284 g) significa un errore di più o meno 8.4% che non è poco. Tuttavia l’RMSE degli altri modelli testati è molto simile perciò si continua a preferire la semplicità del modello 4.

mod10 <- update(mod4,.~.-Gestazione -Sesso)
summary(mod4)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso, 
##     data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.2  -184.3   -17.6   163.3  2627.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6651.1188   135.5172 -49.080  < 2e-16 ***
## Gestazione     31.2737     3.7856   8.261 2.31e-16 ***
## Lunghezza      10.2054     0.3007  33.939  < 2e-16 ***
## Cranio         10.6704     0.4245  25.139  < 2e-16 ***
## SessoM         79.1049    11.2117   7.056 2.22e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared:  0.7261, Adjusted R-squared:  0.7257 
## F-statistic:  1654 on 4 and 2495 DF,  p-value: < 2.2e-16
summary(mod10)
## 
## Call:
## lm(formula = Peso ~ Lunghezza + Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1295.09  -184.36   -11.95   162.85  2792.61 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6306.9134   124.8791  -50.50   <2e-16 ***
## Lunghezza      11.6312     0.2682   43.37   <2e-16 ***
## Cranio         11.2847     0.4298   26.26   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 281.4 on 2497 degrees of freedom
## Multiple R-squared:  0.7129, Adjusted R-squared:  0.7127 
## F-statistic:  3101 on 2 and 2497 DF,  p-value: < 2.2e-16

Provando a rimuovere invece le variabili rimaste si ottengono i seguenti R2:

  • Rimuovendo Cranio: 65.5%
  • Rimuovendo Lunghezza: 60%
  • Rimuovendo Gestazione: 71.8%
  • Rimuovendo Sesso: 72%

Di conseguenza “Cranio” e “Lunghezza” sono le features più impattanti ed importanti mentre “Gestazione” e “Sesso” si potrebbero togliere ma si possono considerare variabili di controllo importanti.

Si procede con un analisi dei residui del modello finale (modello 4) e valutazione della presenza di leverage o outliers (come visto inizialmente probabilmente ce ne saranno):

par(mfrow=c(2,2))
plot(mod4)

# Test diagnostici

# Normalità residui
shapiro_res <- shapiro.test(residuals(mod4))

# Eteroschedasticità
bp_res <- lmtest::bptest(mod4)

# Autocorrelazione residui
dw_res <- lmtest::dwtest(mod4)

# Tabella
tabella_test <- data.frame(
  Test = c("Shapiro-Wilk", "Breusch-Pagan", "Durbin-Watson"),
  p_value = c(
    round(shapiro_res$p.value, 2),
    round(bp_res$p.value, 2),
    round(dw_res$p.value, 2)
  )
)

kable(tabella_test,
      caption = "Risultati dei test diagnostici del modello mod4",
      align = "c")
Risultati dei test diagnostici del modello mod4
Test p_value
Shapiro-Wilk 0.00
Breusch-Pagan 0.00
Durbin-Watson 0.13 
n=length(dati$Peso)

# Leverage
lev <- hatvalues(mod4)

# Plot leverage
plot(lev, main = "Leverage per osservazione", ylab = "Leverage", xlab = "Riga")
p <- sum(lev)
soglia <- 2 * p / n
abline(h = soglia, col = 2, lty = 2)

# Outliers
plot(rstudent(mod4))
abline(h=c(-2,2))

# Distanza di Cook
cook <- cooks.distance(mod4)
plot(cook)

In seguito allo Shapiro test si rifiuta l’ipotesi di normalità perciò i residui non sono distribuiti normalmente ed il risultato del Breusch-Pagan test ci fa rifiutare l’ipotesi di omoschedasticità. Il Durbin-Watson test ci conferma l’incorrelazione e questo è positivo. Si nota nel grafico “Residual vs Leverage” un punto che supera la distanza di Cook (il punto 1551). Si è dunque approfondita la presenza di punti di leverage e/o outliers individuandone rispettivamente 135 e 3, i quali potrebbero essere la causa, o parte di essa, dei risultati negativi ottenuti durante l’analisi dei residui. In ogni caso, con una quantità grande di dati, lo Shapiro test non è abbastanza per rifiutare il modello e graficamente (grafico Residuals vs Fitted) si può notare che i residui sono distribuiti intorno allo 0 (positivo) e non si notano pattern inoltre (grafico Q-Q Residuals) seguono la diagonale (positivo). Il grafico Scale-Location (anche se il Breusch-Pagan test è fallito) non mostra trend e i residui sono distribuiti in una nuvola (casuali) il che è positivo; si notano però (anche nel grafico Residuals vs Fitted) dei punti iniziali fuori dalla nuvola il che è dovuto ai vari valori di leverage e/o outliers. Si procede con la visualizzazione di questi valori di leverage e outliers in modo da capire se sono dati credibili o se contengono errori.

#Tabella leverage

idx_lev <- which(lev > soglia)

tabella_leverage <- cbind(
  Leverage = round(lev[idx_lev], 3),  # arrotondato a 3 cifre
  dati[idx_lev, ]
)

kable(tabella_leverage,
      digits = 3,
      caption = "Osservazioni con leverage elevato")
Osservazioni con leverage elevato
Leverage Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
15 0.006 33 3 0 34 2400 470 298 Ces osp3 M
34 0.006 27 0 0 39 3150 480 382 Nat osp1 F
42 0.004 28 1 0 41 2720 500 310 Ces osp2 M
61 0.005 34 0 0 39 3620 545 334 Nat osp1 F
67 0.005 29 0 0 33 2400 450 320 Nat osp2 M
96 0.005 39 3 0 37 3640 510 376 Ces osp2 M
101 0.007 31 0 0 34 1370 390 287 Nat osp2 F
106 0.013 29 4 0 30 1340 400 273 Ces osp1 M
117 0.005 34 1 0 34 2160 435 303 Ces osp3 M
131 0.007 30 0 0 34 2290 450 285 Nat osp2 M
151 0.011 20 0 0 41 2280 450 280 Ces osp3 M
155 0.007 30 0 0 36 3610 410 330 Nat osp1 M
190 0.005 26 2 0 39 4050 525 390 Nat osp1 M
205 0.005 45 2 0 38 3850 505 384 Nat osp3 M
206 0.009 39 1 0 31 1500 405 295 Nat osp3 M
220 0.007 23 1 0 40 3520 445 363 Nat osp1 F
249 0.005 36 1 0 34 2500 440 338 Nat osp3 F
295 0.004 18 0 0 40 1850 460 305 Nat osp3 F
304 0.004 36 0 0 37 2060 420 326 Nat osp2 F
305 0.005 41 2 0 33 2300 450 323 Nat osp1 M
310 0.029 40 3 0 28 1560 420 379 Nat osp3 F
312 0.013 26 1 0 32 1280 360 276 Nat osp2 M
315 0.005 33 2 0 35 2910 450 355 Nat osp3 F
348 0.004 32 0 0 38 3560 460 360 Nat osp3 M
378 0.015 27 0 0 28 1285 400 274 Nat osp1 F
383 0.004 22 1 0 39 3600 550 346 Nat osp2 F
445 0.007 27 0 0 32 1550 410 289 Nat osp1 F
471 0.004 29 0 0 40 3680 560 346 Nat osp2 M
486 0.005 22 0 0 33 2250 445 310 Nat osp3 F
492 0.008 34 2 0 33 1410 380 295 Nat osp2 F
565 0.005 24 1 0 40 4250 520 386 Nat osp2 F
587 0.008 16 1 0 31 1900 440 300 Nat osp2 F
592 0.006 30 1 0 32 2260 440 322 Ces osp3 F
615 0.005 27 1 0 35 2100 440 345 Ces osp2 F
629 0.004 23 0 0 34 2450 475 329 Nat osp1 F
638 0.006 25 0 0 33 1720 420 300 Nat osp1 M
656 0.005 38 3 0 41 2320 450 330 Nat osp2 F
666 0.004 32 1 0 40 4240 555 345 Nat osp2 F
684 0.009 30 1 0 39 3000 475 390 Ces osp2 F
697 0.006 30 0 0 39 2820 510 300 Ces osp3 F
702 0.005 25 0 0 33 2220 450 320 Nat osp1 F
726 0.004 30 0 0 35 2350 446 344 Nat osp1 F
748 0.008 35 0 0 33 1390 390 277 Nat osp1 F
750 0.007 24 0 0 35 1450 405 280 Nat osp1 F
765 0.006 26 1 0 33 1970 445 300 Nat osp3 M
805 0.014 30 2 0 29 1190 360 272 Nat osp2 F
821 0.004 22 0 0 41 3050 495 310 Nat osp1 F
895 0.005 30 1 0 34 2580 470 347 Nat osp2 M
928 0.022 25 0 0 28 830 310 254 Nat osp1 F
947 0.008 34 3 0 32 1615 390 297 Nat osp3 F
956 0.008 25 0 0 41 2210 430 310 Nat osp3 F
964 0.005 26 0 0 40 4010 550 335 Nat osp2 F
968 0.004 23 1 0 39 2900 470 366 Nat osp1 F
991 0.004 35 2 0 40 4050 550 385 Nat osp1 M
1014 0.008 17 0 0 37 2050 390 295 Nat osp2 F
1067 0.008 26 3 0 31 1960 420 300 Nat osp2 F
1091 0.009 30 1 0 33 1770 410 275 Nat osp3 M
1096 0.004 37 0 0 34 1750 420 306 Ces osp3 F
1130 0.007 33 11 0 43 3400 475 360 Nat osp1 M
1157 0.004 26 0 0 40 4060 505 380 Ces osp1 F
1166 0.004 36 3 0 39 2950 505 307 Nat osp1 F
1181 0.006 30 1 0 36 4070 500 373 Nat osp2 M
1188 0.006 21 0 0 40 4140 550 320 Ces osp1 M
1200 0.005 21 0 0 40 3200 525 310 Nat osp1 F
1238 0.005 19 0 0 33 2270 440 315 Nat osp1 M
1248 0.015 26 1 0 30 1170 370 266 Nat osp2 M
1273 0.007 32 1 0 33 2040 480 307 Ces osp1 F
1283 0.004 29 0 0 40 4250 550 376 Ces osp3 F
1293 0.006 30 3 0 38 4600 485 380 Nat osp1 M
1294 0.005 24 1 0 40 3250 460 360 Nat osp2 F
1356 0.005 32 1 0 40 4090 525 390 Ces osp3 F
1357 0.007 22 0 0 32 2340 445 304 Nat osp1 F
1361 0.004 25 0 0 39 3570 520 315 Nat osp1 F
1385 0.012 33 0 0 29 1620 410 292 Nat osp3 F
1395 0.005 30 2 0 39 3790 505 304 Ces osp3 M
1400 0.006 22 0 0 34 2590 485 314 Nat osp2 M
1402 0.005 35 1 0 39 2660 460 364 Ces osp1 F
1420 0.005 32 1 0 38 3770 530 380 Nat osp3 F
1428 0.008 30 1 0 36 1280 385 292 Nat osp2 F
1429 0.021 24 4 0 29 1280 390 355 Nat osp1 F
1551 0.049 35 1 0 38 4370 315 374 Nat osp3 F
1553 0.007 30 4 0 35 4520 520 360 Nat osp2 F
1556 0.006 37 0 0 41 2420 490 300 Ces osp1 M
1560 0.005 17 0 0 41 2800 455 328 Nat osp3 M
1593 0.005 41 3 0 35 1500 420 304 Nat osp1 M
1606 0.005 28 0 0 41 3050 460 352 Nat osp3 F
1610 0.008 37 3 0 33 2000 470 293 Ces osp1 F
1619 0.015 31 0 0 31 990 340 278 Ces osp2 F
1628 0.005 31 0 0 35 2120 410 312 Nat osp1 F
1634 0.005 32 1 0 38 2500 430 328 Nat osp1 M
1686 0.009 27 0 0 31 1800 430 308 Ces osp3 M
1692 0.004 15 0 0 35 2700 470 300 Nat osp1 F
1693 0.005 25 0 0 41 3100 500 305 Nat osp2 F
1701 0.010 22 0 0 32 1430 380 301 Nat osp1 M
1712 0.007 28 0 0 39 3800 520 300 Nat osp3 F
1735 0.004 15 0 0 37 2750 440 345 Nat osp2 M
1780 0.026 25 2 0 25 900 325 253 Nat osp3 F
1802 0.004 27 2 0 39 3600 540 330 Nat osp1 M
1806 0.004 42 2 0 37 3290 455 355 Nat osp1 M
1809 0.008 35 0 0 32 1780 420 277 Ces osp1 F
1827 0.006 32 1 0 33 2100 420 310 Nat osp1 M
1858 0.004 32 1 0 37 3860 520 371 Nat osp1 M
1868 0.005 29 0 0 40 3470 525 390 Nat osp1 M
1893 0.004 22 1 0 34 3030 470 312 Nat osp2 F
1911 0.004 25 1 0 39 3480 520 312 Nat osp1 M
1920 0.004 26 0 1 39 4930 550 350 Ces osp2 F
1977 0.005 39 4 0 34 2970 480 350 Ces osp2 F
2037 0.004 42 3 0 37 4020 525 368 Ces osp1 M
2040 0.011 27 1 0 38 3240 410 359 Ces osp1 F
2066 0.004 30 2 0 41 3540 470 355 Nat osp1 M
2089 0.006 32 1 1 33 1780 400 305 Ces osp1 F
2114 0.013 36 0 0 31 1180 355 270 Nat osp3 F
2115 0.012 35 1 0 32 1890 500 309 Nat osp2 F
2120 0.018 32 0 0 27 1140 370 267 Nat osp3 F
2140 0.006 30 2 0 33 1600 410 290 Ces osp1 F
2149 0.013 39 3 0 30 1300 380 276 Nat osp1 M
2175 0.021 37 8 0 28 930 355 235 Nat osp1 F
2200 0.011 33 0 0 30 1750 410 294 Nat osp2 M
2215 0.005 29 1 0 40 2440 465 298 Nat osp2 F
2216 0.008 22 0 0 32 2580 470 330 Nat osp1 F
2220 0.004 23 3 1 41 2620 475 313 Nat osp3 M
2224 0.006 41 1 0 33 2000 425 312 Ces osp3 M
2225 0.005 27 0 0 35 3140 465 290 Nat osp2 F
2257 0.006 40 0 0 34 1580 400 300 Nat osp2 F
2307 0.014 26 1 0 30 1170 370 273 Nat osp3 M
2318 0.005 17 0 0 41 3100 500 307 Ces osp1 M
2337 0.005 31 3 1 37 3440 460 362 Ces osp1 M
2359 0.005 25 6 0 33 2230 430 313 Nat osp3 F
2391 0.004 36 1 0 41 4400 565 366 Nat osp1 F
2408 0.010 37 2 0 31 1690 405 290 Nat osp2 M
2437 0.024 28 1 0 27 980 320 265 Nat osp1 M
2452 0.023 28 0 0 26 930 345 245 Ces osp3 F
2458 0.008 31 0 0 31 1730 430 300 Nat osp3 F
2476 0.004 19 0 0 40 2910 520 315 Nat osp1 F
2478 0.006 32 1 0 33 2740 475 324 Ces osp2 F 
# Tabella outliers

out_test <- car::outlierTest(mod4)

if (!is.null(out_test)) {

  idx_out <- as.numeric(names(out_test$rstudent))
  rstud_val <- out_test$rstudent

  tabella_outlier <- cbind(
    rstudent = rstud_val,
    dati[idx_out, ]
  )

  kable(tabella_outlier,
        digits = 3,
        caption = "Osservazioni outlier (outlierTest di Bonferroni)")

} else {
  cat("Nessun outlier significativo individuato.\n")
}
Osservazioni outlier (outlierTest di Bonferroni)
rstudent Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
1551 9.986 35 1 0 38 4370 315 374 Nat osp3 F
155 4.951 30 0 0 36 3610 410 330 Nat osp1 M
1306 4.781 23 0 0 41 4900 510 352 Nat osp2 F

In seguito alla visualizzazione delle righe dei valori di leverage e outliers non sono stati trovati dati dei quali si potesse affermare con certezza l’erroneità dei valori presenti perciò il modello finale sarà il modello 4 con i limiti risultati dell’analisi dei residui e dalla presenza di parecchi valori tra leverage e outliers di cui uno sulla soglia d’allarme della distanza di Cook. Sicuramente questi valori contribuiscono anche all’RMSE che porta ad un errore di 8.4% rispetto al peso medio e all’R2 di 72.6% che non è male ma neanche ottimo. In linea di massima il modello si potrà utilizzare se l’errore sul peso predetto è accettabile, nel caso serva più precisione si può utilizzare il modello per avere una stima ma va approfondito a parte.

Di seguito la previsione del peso di un neonato, rispettivamente di sesso femminile e maschile, utilizzando come lunghezza e diametro del cranio la media dei valori fissati i valori di terza gravidanza per la variabile “N.gravidanze” e 39° settimana per la variabile “Gestazione”:

dati_filtrati <- subset(dati,
                         N.gravidanze == 3 & Gestazione == 39)

print("Peso previsto di un neonato di sesso femminile nato da madre alla terza gravidanza e alla 39° settimana di gestazione:")
## [1] "Peso previsto di un neonato di sesso femminile nato da madre alla terza gravidanza e alla 39° settimana di gestazione:"
as.numeric(predict(mod4, newdata = data.frame(Lunghezza=mean(dati_filtrati$Lunghezza, na.rm = TRUE), Gestazione=39, Cranio=mean(dati_filtrati$Cranio, na.rm = TRUE), Sesso="F")))
## [1] 3305.478

Di seguito altre previsioni:

## [1] "Peso previsto di un neonato di sesso femminile, di lunghezza 496 mm, con un diametro del cranio di 350 mm e alla 39° settimana di gestazione:"
## [1] 3365.072
## [1] "Peso previsto di un neonato di sesso maschile, di lunghezza 496 mm, con un diametro del cranio di 350 mm e alla 39° settimana di gestazione:"
## [1] 3444.177
## [1] "Peso previsto di un neonato di sesso femminile, di lunghezza 496 mm, con un diametro del cranio di 380 mm e alla 39° settimana di gestazione:"
## [1] 3685.183

Si conclude con alcune visualizzazioni grafiche di come varia la nostra variabile target (Peso) in funzione delle variabili ritenute significative per il modello:

ggplot(data=dati)+
  geom_point(aes(x=Lunghezza,
                 y=Peso,
                 col=Sesso))+
  geom_smooth(aes(x=Lunghezza,
                 y=Peso,
                 col=Sesso), se=F, method="lm")
## `geom_smooth()` using formula = 'y ~ x'

ggplot(data=dati)+
  geom_point(aes(x=Cranio,
                 y=Peso,
                 col=Sesso))+
  geom_smooth(aes(x=Cranio,
                 y=Peso,
                 col=Sesso), se=F, method="lm")
## `geom_smooth()` using formula = 'y ~ x'

Sia nella rappresentazione grafica del peso rispetto alla lughezza che rispetto al diametro del cranio si nota come le rette (sesso maschile e sesso femminile) ottenute siano molto vicine, ulteriore conferma del fatto che la variabile “Sesso” potrebbe essere rimossa al fine di semplificare il modello (è stato scelto infatti di tenerla come variabile di controllo).

ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col=Sesso))+
  geom_smooth(aes(x=Gestazione,
                 y=Peso,
                 col=Sesso), se=F, method="lm")
## `geom_smooth()` using formula = 'y ~ x'

X_num <- dati[, c("Lunghezza", "Cranio", "Gestazione")]
cor(X_num, use = "complete.obs")
##            Lunghezza    Cranio Gestazione
## Lunghezza  1.0000000 0.6033410  0.6189204
## Cranio     0.6033410 1.0000000  0.4608289
## Gestazione 0.6189204 0.4608289  1.0000000

Nel grafico del peso rispetto alle settimane di gestazione si nota anche qui la poca influenza del sesso del bambino viste le rette molto vicine; c’è da sottolineare però che, nel caso del diametro del cranio e delle settimane di gestazione, le rette sono si vicine ma parallele con la retta del sesso maschile sopra in entrambi i casi a quella del sesso femminile, dunque è una maggiore conferma della scelta di tenere la variabile “Sesso” come variabile di controllo.

Un ulteriore ricontrollo è stato fatto, tramite la matrice di correlazione, sulla possibile dipendenza tra le variabili (dal momento che logicamente lunghezza, diametro del cranio e settimane di gestazione sicuramente hanno una dipendenza anche se piccola): il risultato, come ci si aspettava, è un fattore di correlazione basso tra diametro del cranio e settimane di gestazione ma un fattore di correlazione del 60% tra lunghezza-diametro del cranio e lunghezza-settimane di gestazione. Si mantiene comunque il modello 4 come modello finale considerando “Lunghezza” e “Cranio” vere e proprio variabili del modello mentre “Sesso” e “Gestazione” variabili di controllo del modello.