Un Modello Statistico per Prevedere il Peso dei Neonati

Questo documento contiene il progetto finale di Luigi Carlucci del corso “Statistica Inferenziale” per il Master in Data Science di ProfessionAI.

Il progetto ha come obiettivo finale la creazione di un modello statistico in grado di prevedere il peso dei neonati alla nascita, sulla base di un dataset contenente variabili cliniche raccolti su 2500 neonati provenienti da tre ospedali.

library(dplyr)
library(kableExtra)
library(moments) 
library(ggplot2)
library(GGally)
library(ggcorrplot)
library(patchwork)
library(car)
library(MASS)
library(Metrics)
library(lmtest)
library(plotly)
df = read.csv("neonati.csv", stringsAsFactors = T)
attach(df)
kbl(head(df)) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
26 0 0 42 3380 490 325 Nat osp3 M
21 2 0 39 3150 490 345 Nat osp1 F
34 3 0 38 3640 500 375 Nat osp2 M
28 1 0 41 3690 515 365 Nat osp2 M
20 0 0 38 3700 480 335 Nat osp3 F
32 0 0 40 3200 495 340 Nat osp2 F

1. Analisi Preliminare

Nella prima fase viene eseguita un’analisi descrittiva per esplorare le variabili, comprenderne la distribuzione e identificare eventuali anomalie.

kbl(summary(df)) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
Min. : 0.00 Min. : 0.0000 Min. :0.0000 Min. :25.00 Min. : 830 Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1256
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 Nat:1772 osp2:849 M:1244
Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00 Median :3300 Median :500.0 Median :340 NA osp3:835 NA
Mean :28.16 Mean : 0.9812 Mean :0.0416 Mean :38.98 Mean :3284 Mean :494.7 Mean :340 NA NA NA
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 NA NA NA
Max. :46.00 Max. :12.0000 Max. :1.0000 Max. :43.00 Max. :4930 Max. :565.0 Max. :390 NA NA NA

Viene verificata la normalità della variabile Peso, che verrà considerata come variabile risposta nel modello di regressione.

# Calcolo di skewness (asimmetria) e kurtosis (appiattimento)
paste("Asimmetria:", round(skewness(Peso), 3))
## [1] "Asimmetria: -0.647"
paste("Curtosi:", round(kurtosis(Peso), 3))
## [1] "Curtosi: 5.032"

La variabile Peso ha un indice di asimmetria negativo, quindi presenta una coda verso sinistra, e un indice di curtosi positivo, quindi presenta una distribuzione leptocurtica (allungata), come si può osservara anche da un grafico della densità.

ggplot(df)+
geom_density(aes(x=Peso), col="black", fill="steelblue") +
labs(title = "Densità del peso dei neonati",
x = "Peso (g)", y = "Densità")

# Test di Shapiro-Wilk per verificare la normalità
shapiro.test(Peso)
## 
##  Shapiro-Wilk normality test
## 
## data:  Peso
## W = 0.97066, p-value < 2.2e-16

Eseguendo il test di Shapiro-Wilk, viene rifiutata l’ipotesi di normalità per la variabile Peso. In questo caso, poichè la distribuzione della variabile risposta non è normale, anche i residui del modello di regressione potrebbero non avere una distribuzione normale (che rappresenta una delle assunzioni della regressione lineare), ma questo verrà verificato successivamente alla creazone del modello.

In seguito, viene creata una matrice di correlazione per visualizzare le relazioni tra le variabili qualitative continue, mostrata graficamante con la libreria ggcorrplot.

# Matrice di correlazione (solo per le variabili quantitative continue)
cor_matrix <- cor(df %>% dplyr::select(Peso, Lunghezza, Cranio, Anni.madre, Gestazione))
# Visualizzazione della matrice di correlazione con ggcorrplot
ggcorrplot(cor_matrix, hc.order = TRUE, lab = TRUE)

Tra le variabili esplicative non sembrano esserci correlazioni particolarmente alte e non dovrebbe causare problemi di multicollinearità per il modello di regressione, anche se questo verrà poi ulteriormente verificato tramite l’indicatore Variance Inflation Factor (VIF).

Le relazioni tra le variabili continue vengono visualizzate anche tramite una matrice di scatter plot (pair plot), creata con la funzione ggpairs della libreria GGally.

# Pair plot per le variabili quantitative continue

# Funzione per modificare i parametri degli scatter plot
lower_func <- function(data, mapping, ...) {
    ggplot(data, mapping) + 
      geom_point(size = 0.7, alpha = 0.5, position = "jitter") +
      geom_smooth(method = "lm", color = "red")
}

pp = ggpairs(df, 
            columns = c("Peso", "Lunghezza", "Cranio", "Anni.madre", "Gestazione"),
            lower = list(continuous=lower_func))
pp[5,1] = pp[5,1] + scale_x_continuous(breaks = c(1500,3000,4500))
pp[5,2] = pp[5,2] + scale_x_continuous(breaks = c(350,450,550))
pp

Viene poi esplorato l’impatto delle variabili qualitative Sesso, Fumatrici (fumo materno) e Tipo.parto sulla variabile Peso, tramite dei boxlot.

# Boxplot per esplorare la varabile Peso rispetto a variabili qualitative
bp1 = ggplot(df, aes(x = "", y = Peso)) + 
  geom_boxplot(fill="steelblue") +
  ggtitle("Distribuzione del Peso (tutti i neonati)")

bp2 = ggplot(df, aes(x = Sesso, y = Peso, fill = Sesso)) + 
  geom_boxplot() +
  ggtitle("Distribuzione del Peso per Sesso")

bp3 = ggplot(df, aes(x = as.factor(Fumatrici), y = Peso, fill = as.factor(Fumatrici))) + 
  geom_boxplot() +
  ggtitle("Distribuzione del Peso per Fumatrici") + 
  scale_x_discrete(name = "Fumatrici", labels = c("0" = "No", "1" = "Si")) +
  scale_fill_discrete(name = "Fumatrici", labels = c("No", "Si"))

bp4 = ggplot(df, aes(x = Tipo.parto, y = Peso, fill = Tipo.parto)) + 
  geom_boxplot() +
  ggtitle("Distribuzione del Peso per Tipo di parto")

bp1 + bp2 + bp3 + bp4

Vengono ora saggiate le seguenti ipotesi con i relativi test statistici adatti:

  1. In alcuni ospedali si fanno più parti cesarei

  2. Le medie di Peso e Lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione

  3. Le misure antropometriche (Peso, Lunghezza, Cranio) sono significativamente diverse tra i due sessi

Nel primo punto, viene verificato se la distribuzione dei parti cesarei varia significativamente tra gli ospedali. Viene quindi creata una tabella di contingenza tra le variabili Ospedale e Tipo.parto e viene eseguito un test Chi-quadrato, adatto per saggiare l’ipotesi di indipendenza tra variabili qualitative.

# Creazione della tabella di contingenza tra Ospedale e Tipo.parto
tab_parti <- table(df$Ospedale, df$Tipo.parto)
kbl(tab_parti) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
Ces Nat
osp1 242 574
osp2 254 595
osp3 232 603 
# Esecuzione del test Chi-quadrato
chi2_test <- chisq.test(tab_parti)
print(chi2_test)
## 
##  Pearson's Chi-squared test
## 
## data:  tab_parti
## X-squared = 1.0972, df = 2, p-value = 0.5778

Dal test Chi-quadrato è risultato un p-value di \(0.577\). Pertanto, non si può rifiutare l’ipotesi di indipendenza tra le variabili Ospedale e Tipo.parto, cioè non vi è una differenza significativa nel tipo di parto tra i tre ospedali. Quindi in nessun ospedale si fanno significativamente più parti cesarei.

Nel secondo punto, vengono confrontate le medie campionarie delle variabili Peso e Lunghezza con quelle della popolazione (Peso = 3300 g, Lunghezza = 500 mm, fonte: Ospedale Bambino Gesu) per verificare che siano significativamente uguali. In questo caso è adeguato utilizzare un T-test a campione singolo se i dati sono normali o un Test di Wilcoxon se non sono normali.

Era già stato verificato, tramite il test di Shapiro-Wilk per la normalità, che i dati della variabile Peso non hanno una distribuzione normale, per cui viene eseguito il test di Wilcoxon.

# Esecuzione del test di Wilcoxon per la variabile Peso
wilcox.test(Peso, mu = 3300)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  Peso
## V = 1495594, p-value = 0.9612
## alternative hypothesis: true location is not equal to 3300

Dal test di Wilcoxon è risultato un p-value di \(0.961\). Pertanto, non si può rifiutare l’ipotesi che la media del peso di questo campione di neonati sia significativamente uguale a quella della popolazione.

Viene ora eseguito il test di Shapiro-Wilk per verificare la normalità della variabile Lunghezza.

shapiro.test(Lunghezza)
## 
##  Shapiro-Wilk normality test
## 
## data:  Lunghezza
## W = 0.90941, p-value < 2.2e-16

Anche in questo caso è risultato che la variabile Lunghezza non ha una distribuzione normale, per cui viene eseguito il test di Wilcoxon.

# Esecuzione del test di Wilcoxon per la variabile Lunghezza
wilcox.test(Lunghezza, mu = 500)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  Lunghezza
## V = 877236, p-value < 2.2e-16
## alternative hypothesis: true location is not equal to 500

Con un p-value \(< 2.2\times10^{-16}\) viene qui rifiutata l’ipotesi che la media della lunghezza di questo campione di neonati sia significativamente uguale a quella della popolazione. Quindi si può affermare che la media della lunghezza calcolata per questo campione, pari a 494.7, è significativamente inferiore alla media della popolazione.

Nel terzo punto vengono confrontate le distribuzioni di Peso, Lunghezza e Cranio tra neonati maschi e femmine per verificare che siano significativamente diverse tra i due sessi. In questo caso viene utilizzato un T-test per campioni indipendenti (definiti in base alla variabile categorica Sesso) se i dati sono normali o un Test di Wilcoxon se non sono normali. Era già stato verificato, tramite il test di Shapiro-Wilk per la normalità, che i dati delle variabili Peso e Lunghezza non hanno una distribuzione normale, per cui viene eseguito il test di Wilcoxon.

# Esecuzione del test di Wilcoxon per la variabile Peso
wilcox.test(Peso ~ Sesso)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Peso by Sesso
## W = 538641, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

Con un p-value \(< 2.2\times10^{-16}\) viene rifiutata l’ipotesi che la media del peso sia significativamente uguale tra i due sessi.

# Esecuzione del test di Wilcoxon per la variabile Lunghezza
wilcox.test(Lunghezza ~ Sesso)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Lunghezza by Sesso
## W = 594455, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

Anche in questo caso, con un p-value \(< 2.2\times10^{-16}\) viene rifiutata l’ipotesi che la media della lunghezza sia significativamente uguale tra i due sessi.

Viene ora eseguito il test di Shapiro-Wilk per verificare la normalità della variabile Cranio.

shapiro.test(Cranio)
## 
##  Shapiro-Wilk normality test
## 
## data:  Cranio
## W = 0.96357, p-value < 2.2e-16

Poichè anche la variabile Cranio non ha una distribuzione normale, viene eseguito il test di Wilcoxon.

# Esecuzione del test di Wilcoxon per la variabile Lunghezza
wilcox.test(Cranio ~ Sesso)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Cranio by Sesso
## W = 641638, p-value = 9.633e-15
## alternative hypothesis: true location shift is not equal to 0

Con un p-value di \(9.633 \times 10^{-15}\) viene rifiutata l’ipotesi che la media del diametro del cranio sia significativamente uguale tra i due sessi.

2. Creazione del Modello di Regressione

In questa fase viene sviluppato un modello di regressione lineare multipla in cui Peso è la variabile risposta e che includa tutte le variabili del dataset come variabili esplicative. In questo modo, è possibile quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato.

# Costruzione del modello di regressione lineare multipla
mod1 <- lm(Peso ~ ., data = df)
# Riassunto del modello
summary(mod1)
## 
## Call:
## lm(formula = Peso ~ ., data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1124.40  -181.66   -14.42   160.91  2611.89 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6738.4762   141.3087 -47.686  < 2e-16 ***
## Anni.madre        0.8921     1.1323   0.788   0.4308    
## N.gravidanze     11.2665     4.6608   2.417   0.0157 *  
## Fumatrici       -30.1631    27.5386  -1.095   0.2735    
## Gestazione       32.5696     3.8187   8.529  < 2e-16 ***
## Lunghezza        10.2945     0.3007  34.236  < 2e-16 ***
## Cranio           10.4707     0.4260  24.578  < 2e-16 ***
## Tipo.partoNat    29.5254    12.0844   2.443   0.0146 *  
## Ospedaleosp2    -11.2095    13.4379  -0.834   0.4043    
## Ospedaleosp3     28.0958    13.4957   2.082   0.0375 *  
## SessoM           77.5409    11.1776   6.937 5.08e-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.2 on 10 and 2489 DF,  p-value: < 2.2e-16

Il Summary del modello mostra le stime dei coefficienti di regressione (“Estimate”) per tutte le variabili, che in questo caso rappresentano gli effetti marginali delle singole variabili esplicative sulla variabile risposta Peso. Tuttavia, nonostante tutte le variabili (con l’esclusione della variabile Anni.madre) risultano avere un effetto sulla variabile Peso, non tutte le variabili risultano essere significative, come indicato dai valori alti del P-value.

Questo suggerisce la possibilità di sviluppare un modello migliore, escludendo alcune delle variabili esplicative inserite in questo modello iniziale. In particolare, è probabile che eliminando le variabili Anni.madre, Fumatrici e Ospedale si otterrebbe un medello migliore e più parsimonioso.

Viene anche calcolato il Variance Inflation Factor (VIF) per controllare la multicollinearità tra le variabili indipendenti inserite nel modello.

vif(mod1)
##                  GVIF Df GVIF^(1/(2*Df))
## Anni.madre   1.187454  1        1.089704
## N.gravidanze 1.186428  1        1.089233
## Fumatrici    1.007392  1        1.003689
## Gestazione   1.695810  1        1.302233
## Lunghezza    2.085755  1        1.444214
## Cranio       1.630796  1        1.277026
## Tipo.parto   1.004242  1        1.002119
## Ospedale     1.004071  2        1.001016
## Sesso        1.040643  1        1.020119

Poiche gli indicatori VIF sono inferiori a 5 per tutte le variabili, non vi sono problemi di multicollinerarità, come suggerito anche in precedenza osservando la matrice di correlazione. Pertanto, non è necessario rimuovere variabili per questo motivo.

3. Selezione del Modello Ottimale

In questa fase viene selezionato il modello più parsimonioso, quindi eliminando le variabili non significative. Questa procedura viene eseguita utilizzare la funzione stepAIC della libreria “MASS”, che applica una selezione stepwise per minimizzare il criterio di informazione del modello. In questo caso, viene utilizzato come criterio il BIC (Bayesian Information Criterion), che tende ad essere più parsimonioso nella selezione, mantenendo solo variabili strettamente rilevanti.

# Selezione del modello ottimale con procedura stepwise basata su BIC
n = nrow(df)
stepwise.mod = stepAIC(mod1, direction = "both", k = log(n))
## Start:  AIC=28139.32
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     46578 186809099 28132
## - Fumatrici     1     90019 186852540 28133
## - Ospedale      2    685979 187448501 28133
## - N.gravidanze  1    438452 187200974 28137
## - Tipo.parto    1    447929 187210450 28138
## <none>                      186762521 28139
## - Sesso         1   3611021 190373542 28179
## - Gestazione    1   5458403 192220925 28204
## - Cranio        1  45326172 232088693 28675
## - Lunghezza     1  87951062 274713583 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     46578 186762521 28139
## - 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     47456 186852540 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     54392 187547285 28126
## - 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     54816 188010731 28125
## - Sesso         1   3655292 191720838 28158
## - Gestazione    1   5464853 193530399 28181
## - Cranio        1  46108583 234174130 28658
## - Lunghezza     1  87632762 275698308 29066
summary(stepwise.mod)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = df)
## 
## 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

In seguito alla procedura di selezione, il modello migliore risultante include le variabili: N.gravidanze, Gestazione, Lunghezza, Cranio, Sesso.

# Confronto del BIC
BIC(mod1, stepwise.mod)
##              df      BIC
## mod1         12 35241.84
## stepwise.mod  7 35220.10

Rispetto al modello iniziale, il modello selezionato ha un valore del BIC minore, risultando quindi un medello migliore da un punto di vista di semplicità e interpretazione.

Vengono testati ora anche modelli con interazioni tra alcune variabili. Viene considerata prima l’interazione tra le variabili N.gravidanze e Gestazione:

# Creazione di un modello con interazione tra le variabili N.gravidanze e Gestazione
mod_int_1 = lm(formula = Peso ~ N.gravidanze * Gestazione + Lunghezza + Cranio + 
    Sesso, data = df)
summary(mod_int_1)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze * Gestazione + Lunghezza + Cranio + 
##     Sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1148.72  -180.77   -14.14   163.93  2638.53 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             -6588.5062   158.7705 -41.497  < 2e-16 ***
## N.gravidanze              -66.1971    70.1095  -0.944    0.345    
## Gestazione                 29.9109     4.3658   6.851 9.20e-12 ***
## Lunghezza                  10.2482     0.3006  34.091  < 2e-16 ***
## Cranio                     10.5463     0.4263  24.741  < 2e-16 ***
## SessoM                     77.9137    11.2017   6.956 4.47e-12 ***
## N.gravidanze:Gestazione     2.0278     1.8036   1.124    0.261    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2493 degrees of freedom
## Multiple R-squared:  0.7271, Adjusted R-squared:  0.7265 
## F-statistic:  1107 on 6 and 2493 DF,  p-value: < 2.2e-16

In questo caso, il termine di interazione (N.gravidanze:Gestazione) non risulta significativo (p-value = 0.261).

BIC(stepwise.mod, mod_int_1)
##              df      BIC
## stepwise.mod  7 35220.10
## mod_int_1     8 35226.66

Anche il BIC del nuovo modello risulta maggiore, quindi si può concludere che non è necessario aggiungere questa interazione al modello.

Viene ora considerata l’interazione tra le variabili Lunghezza e Cranio:

# Creazione di un modello con interazione tra le variabili Lunghezza e Cranio
mod_int_2 = lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza * Cranio + 
    Sesso, data = df)
summary(mod_int_2)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza * Cranio + 
##     Sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1150.74  -180.48   -13.62   165.82  2866.17 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.801e+03  1.018e+03  -1.770  0.07685 .  
## N.gravidanze      1.295e+01  4.321e+00   2.996  0.00276 ** 
## Gestazione        3.810e+01  3.964e+00   9.610  < 2e-16 ***
## Lunghezza        -3.047e-01  2.202e+00  -0.138  0.88996    
## Cranio           -4.759e+00  3.191e+00  -1.491  0.13601    
## SessoM            7.327e+01  1.119e+01   6.545 7.20e-11 ***
## Lunghezza:Cranio  3.158e-02  6.528e-03   4.838 1.39e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.4 on 2493 degrees of freedom
## Multiple R-squared:  0.7295, Adjusted R-squared:  0.7289 
## F-statistic:  1121 on 6 and 2493 DF,  p-value: < 2.2e-16

Il termine di interazione (Lunghezza:Cranio) risulta qui significativo (p-value = \(1.39 \times 10^{-6}\)).

BIC(stepwise.mod, mod_int_2)
##              df      BIC
## stepwise.mod  7 35220.10
## mod_int_2     8 35204.57

In questo caso, inoltre, il BIC risulta minore del modello senza interazione. Viene eseguito anche un test ANOVA, che prende in considerazione le varianze spiegate dai due modelli:

anova(mod_int_2, stepwise.mod)
## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza * Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1   2493 186316621                                  
## 2   2494 188065546 -1  -1748926 23.401 1.395e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Il risultato del test ANOVA conferma che vi è un aumento significativo della varianza spegata dal modello che include l’interazione, che quindi può essere considerato un modello migliore.

Vengono ora testati possibili effetti non lineari. In particolare, osservando gli scatter plot creati in precedenza, viene prima aggiunto al modello l’effetto quadratico della variabile Lunghezza:

# Creazione del modello con Lunghezza^2
mod_quad_1 = update(mod_int_2, ~.+I(Lunghezza^2))
summary(mod_quad_1)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Lunghezza^2) + Lunghezza:Cranio, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1176.69  -179.20   -11.78   165.68  1306.79 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.262e+03  1.003e+03  -2.256 0.024135 *  
## N.gravidanze      1.425e+01  4.254e+00   3.350 0.000821 ***
## Gestazione        4.042e+01  3.909e+00  10.339  < 2e-16 ***
## Lunghezza        -2.137e+01  3.169e+00  -6.744 1.91e-11 ***
## Cranio            2.741e+01  4.725e+00   5.800 7.46e-09 ***
## SessoM            7.186e+01  1.102e+01   6.523 8.29e-11 ***
## I(Lunghezza^2)    4.476e-02  4.914e-03   9.109  < 2e-16 ***
## Lunghezza:Cranio -3.449e-02  9.689e-03  -3.560 0.000378 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269 on 2492 degrees of freedom
## Multiple R-squared:  0.7383, Adjusted R-squared:  0.7375 
## F-statistic:  1004 on 7 and 2492 DF,  p-value: < 2.2e-16
BIC(mod_int_2, mod_quad_1)
##            df      BIC
## mod_int_2   8 35204.57
## mod_quad_1  9 35130.51

Aggiungendo al modello il termine quadratico della variabile Lunghezza, questo risulta essere molto significativo (p-value \(< 2\times10^{-16}\) ). Inoltre, vi è un’ulteriore diminuzione del BIC e anche un aumento di un punto del coefficiente adjusted R-squared.

Viene aggiunto ora anche l’effetto quadratico della variabile Gestazione:

# Creazione del modello con Gestazione^2
mod_quad_2 = update(mod_quad_1, ~.+I(Gestazione^2))
summary(mod_quad_2)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Lunghezza^2) + I(Gestazione^2) + Lunghezza:Cranio, 
##     data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1191.64  -181.49   -12.11   165.26  1325.43 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -3.400e+03  1.048e+03  -3.244 0.001193 ** 
## N.gravidanze      1.451e+01  4.245e+00   3.418 0.000640 ***
## Gestazione        2.862e+02  6.769e+01   4.228 2.44e-05 ***
## Lunghezza        -3.082e+01  4.091e+00  -7.532 6.93e-14 ***
## Cranio            2.042e+01  5.090e+00   4.012 6.19e-05 ***
## SessoM            7.328e+01  1.100e+01   6.664 3.28e-11 ***
## I(Lunghezza^2)    4.947e-02  5.070e-03   9.757  < 2e-16 ***
## I(Gestazione^2)  -3.227e+00  8.872e-01  -3.637 0.000281 ***
## Lunghezza:Cranio -2.046e-02  1.041e-02  -1.966 0.049358 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 268.3 on 2491 degrees of freedom
## Multiple R-squared:  0.7396, Adjusted R-squared:  0.7388 
## F-statistic: 884.6 on 8 and 2491 DF,  p-value: < 2.2e-16
BIC(mod_quad_1, mod_quad_2)
##            df      BIC
## mod_quad_1  9 35130.51
## mod_quad_2 10 35125.09

Anche il termine quadratico della variabile Gestazione risulta essere molto significativo (p-value = \(0.0002\)). Inoltre, vi è un’ulteriore diminuzione del BIC e un leggero aumento del coefficiente adjusted R-squared.

anova(mod_quad_2, mod_int_2)
## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Lunghezza^2) + I(Gestazione^2) + Lunghezza:Cranio
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza * Cranio + Sesso
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1   2491 179360498                                  
## 2   2493 186316621 -2  -6956123 48.304 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Il risultato del test ANOVA conferma che vi è un aumento significativo della varianza spegata dal modello con i termini quadratici di Lunghezza e Gestazione. Questo modello viene quindi considerato il modello migliore tra quelli creati.

4. Analisi della Qualità del Modello

In questa fase, viene prima confrontata la capacità predittiva del modello finale rispetto ai i modelli precedenti, utilizzando tre metriche:

  • BIC

  • adjusted R-squared

  • Root Mean Squared Error (RMSE)

Il RMSE è una metrica che rappresenta la media delle distanze dei valori predetti dal modelli rispetto ai valori osservati. In questo caso, un valore più piccolo del RMSE indica una migliore capacità del modello di prevedere i dati.

# Creazione della lista dei modelli
modelli <- list(
  "Base" = mod1,
  "Stepwise" = stepwise.mod,
  "Interazione" = mod_int_2,
  "Quadratico" = mod_quad_2
)
# Calcolo delle metriche per ciascun modello
metriche <- data.frame(
  BIC = sapply(modelli, BIC),
  Adjusted_R2 = sapply(modelli, function(m) summary(m)$adj.r.squared),
  RMSE = sapply(modelli, function(m) rmse(df$Peso, predict(m)))
)
kbl(metriche) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
BIC Adjusted_R2 RMSE
Base 35241.84 0.7278038 273.3222
Stepwise 35220.10 0.7264542 274.2740
Interazione 35204.57 0.7288893 272.9957
Quadratico 35125.09 0.7388017 267.8511

Come si può osservare dalla tabella riassuntiva, tutte le metriche sono migliori per il modello finale, che include i termini quadratici.

Viene ora effettuata un’analisi dei residui del modello.

par(mfrow=c(2,2))
plot(mod_quad_2, pch=20)

Dall’analisi grafica si osserva che i residui mostrano in buona parte una distribuzione normale, ma vi sono alcuni valori, in particolare nella coda superiore, che si allontanano visibilmente. Inoltre, i residui non sembrano distribuiti uniformemente intorno alla media di 0.

Infine, osservando l’ultimo grafico, che mostra possibili valori influenti sulle stime di regressione, si nota come una delle osservazioni superi entrambe le soglie problematiche di 0.5 e 1.

Vengono ora eseguiti dei test sui residui.

# Test di normalità (Shapiro-Wilk)
shapiro.test(residuals(mod_quad_2))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_quad_2)
## W = 0.99039, p-value = 6.954e-12
ggplot()+
geom_density(aes(x=residuals(mod_quad_2)), col="black", fill="steelblue") +
labs(title = "Densità dei residui",
x = "Residui", y = "Densità")

Il test di Shapiro-Wilk conferma che i residui non appaiano perfettamente normali. Come si osserva dal grafico della densità, questo è principalmente dovuto ad alcuni valori estremi nelle due code.

# Test di omoschedasticità (Breusch-Pagan)
bptest(mod_quad_2)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod_quad_2
## BP = 134.53, df = 8, p-value < 2.2e-16

Anche il test di Breusch-Pagan per l’omoschedasticità non viene superato, per cui la varianza dei residui non risulta costante, come ipotizzato precedentemente.

# Test di incorrelazione (Durbin-Watson)
dwtest(mod_quad_2)
## 
##  Durbin-Watson test
## 
## data:  mod_quad_2
## DW = 1.9471, p-value = 0.09295
## alternative hypothesis: true autocorrelation is greater than 0

Dal test di Durbin-Watson per l’incorrelazione, l’ipotesi nulla non viene rifiutata, per cui i residui non sono autocorrelati.

Per identificare precisamente possibili valori influenti, vengono prima calcolati i valori di leva (leverage), cioè osservazioni che si trovano lontano nello spazio delle variabili esplicative.

# Visualizzazione dei valori di leva
leverage = hatvalues(mod_quad_2)
p = sum(leverage)
soglia = 2 * p/n
plot(leverage, pch=20, main="Valori di leva")
abline(h=soglia, col=2)

Dal grafico risultano presenti diversi valori di leva che superano la soglia (2*p/n), di cui uno abbondantemente.

Vengono individuati anche i valori estremi della variabile risposta (outliers), tramite l’analisi dei residui studentizzati.

# Visualizzazione dei residui studentizzati
plot(rstudent(mod_quad_2), pch=20, main="Residui studentizzati", ylab="Studentized residuals")

# Test degli outliers
outlierTest(mod_quad_2)
##       rstudent unadjusted p-value Bonferroni p
## 1551  6.249419         4.8311e-10   1.2078e-06
## 1306  4.968338         7.2110e-07   1.8028e-03
## 1399 -4.463879         8.4073e-06   2.1018e-02
## 155   4.387990         1.1917e-05   2.9792e-02
## 1694  4.290034         1.8546e-05   4.6366e-02

Effettuando un test dei residui studentizzati, risulta che in questo caso 5 valori vengono classificati come outliers, quindi valori potenzialmente influenti.

Viene ora considerato l’effetto contemporaneo di leverage e outliers, calcolando e visualizzando le distanze di Cook:

cook = cooks.distance(mod_quad_2)
plot(cook, pch=20, main="Distanza di Cook")
abline(h=0.5, col=2)

# Valori superiori alla soglia della distanza di Cook di 0.5.
cook[cook>0.5]
##     1551 
## 5.311502

Come notato in precedenza, soltanto un valore supera la soglia problematica della distanza di Cook (0.5), e può quindi essere considerato come influente. Tuttavia, assumendo che l’osservazione non sia un errore di misurazione, il modello viene considerato sufficientemente adeguato.

5. Previsioni e Risultati

In questa fase il modello viene utilizzato per fare previsioni pratiche, stimando il peso di un neonato considerando diverse combinazioni di valori delle variabili esplicative. Per semplificare la procedura viene creata una funzione che prende come argomenti i valori delle variabili esplicative del modello (i valori mediani del dataset sono impostati come valori di default) e restituisce la previsione della variabile Peso.

# Creazione della funzione per la previsione del peso del neonato
previsione_peso = function(N.gravidanze=1, 
                           Gestazione=39, 
                           Lunghezza=500, 
                           Cranio=340, 
                           Sesso="M") {
  
  nuovi_dati <- data.frame(N.gravidanze=c(N.gravidanze), 
                           Gestazione=c(Gestazione), 
                           Lunghezza=c(Lunghezza),
                           Cranio=c(Cranio),
                           Sesso=as.factor(c(Sesso)))
  
  peso_predetto = predict(mod_quad_2, newdata=nuovi_dati)
  
  return(cat("Per i seguenti valori:\n\nNumero di gravidanze:", N.gravidanze,
             "\nSettimane di gestazione:", Gestazione,
             "\nLunghezza del neonato (mm):", Lunghezza,
             "\nDiametro craniale del neonato (mm):", Cranio, 
             "\nSesso del neonato:", Sesso, 
             "\n\nIl peso previsto del neonato (g) è:", peso_predetto))
}

Non inserendo alcun argomento, la funzione restituisce la previsione con i valori di default:

previsione_peso()
## Per i seguenti valori:
## 
## Numero di gravidanze: 1 
## Settimane di gestazione: 39 
## Lunghezza del neonato (mm): 500 
## Diametro craniale del neonato (mm): 340 
## Sesso del neonato: M 
## 
## Il peso previsto del neonato (g) è: 3364.558

Impostando la variabile Sesso come "F", la funzione restituisce la previsione del peso di una neonata:

previsione_peso(Sesso="F")
## Per i seguenti valori:
## 
## Numero di gravidanze: 1 
## Settimane di gestazione: 39 
## Lunghezza del neonato (mm): 500 
## Diametro craniale del neonato (mm): 340 
## Sesso del neonato: F 
## 
## Il peso previsto del neonato (g) è: 3291.283

Si possono inserire diverse combinazioni di argomenti (es. N.gravidanze, Gestazione, Sesso) per modificare i valori delle variabili e restituire la previsione corrispondente:

previsione_peso(N.gravidanze=3)
## Per i seguenti valori:
## 
## Numero di gravidanze: 3 
## Settimane di gestazione: 39 
## Lunghezza del neonato (mm): 500 
## Diametro craniale del neonato (mm): 340 
## Sesso del neonato: M 
## 
## Il peso previsto del neonato (g) è: 3393.577
previsione_peso(N.gravidanze=3, Sesso="F")
## Per i seguenti valori:
## 
## Numero di gravidanze: 3 
## Settimane di gestazione: 39 
## Lunghezza del neonato (mm): 500 
## Diametro craniale del neonato (mm): 340 
## Sesso del neonato: F 
## 
## Il peso previsto del neonato (g) è: 3320.302
previsione_peso(Gestazione=28)
## Per i seguenti valori:
## 
## Numero di gravidanze: 1 
## Settimane di gestazione: 28 
## Lunghezza del neonato (mm): 500 
## Diametro craniale del neonato (mm): 340 
## Sesso del neonato: M 
## 
## Il peso previsto del neonato (g) è: 2594.537
previsione_peso(Gestazione=28, Sesso="F")
## Per i seguenti valori:
## 
## Numero di gravidanze: 1 
## Settimane di gestazione: 28 
## Lunghezza del neonato (mm): 500 
## Diametro craniale del neonato (mm): 340 
## Sesso del neonato: F 
## 
## Il peso previsto del neonato (g) è: 2521.262

Inserendo i valori di tutte le varibili presenti nel modello, viene stimato il peso previsto con la massima precisione:

previsione_peso(N.gravidanze=2, Gestazione=33, Lunghezza=512, Cranio=345, Sesso="F")
## Per i seguenti valori:
## 
## Numero di gravidanze: 2 
## Settimane di gestazione: 33 
## Lunghezza del neonato (mm): 512 
## Diametro craniale del neonato (mm): 345 
## Sesso del neonato: F 
## 
## Il peso previsto del neonato (g) è: 3179.727

6. Visualizzazioni

Infine, vengono creati alcuni grafici per visualizzare alcune relazione più significative tra le variabili incluse nel modello.

Ad esempio, viene visualizzata tramite uno scatter plot la relazione tra la variabile risposta Peso e la variabile esplicativa Gestazione, utilizzando colori diversi per la variabile Sesso (mostrando anche le relative rette di regressione) e diverse dimensioni dei punti per la variabile Lunghezza:

ggplot(data=df) + 
  geom_point(aes(x=Gestazione, y=Peso, colour=Sesso, size=Lunghezza),
             alpha=0.2,
             position = "jitter") +
  geom_smooth(aes(x=Gestazione, y=Peso, colour=Sesso),
              se=F, method="lm") +
  labs(title="Impatto di Gestazione, Sesso e Lunghezza sul Peso")

In seguito viene invece visualizzata la relazione tra Peso e Lunghezza, utilizzando la dimensione dei punti per la variabile Gestazione:

ggplot(data=df) + 
  geom_point(aes(x=Lunghezza, y=Peso, colour=Sesso, size=Gestazione),
             alpha=0.2,
             position = "jitter") +
  geom_smooth(aes(x=Lunghezza, y=Peso, colour=Sesso),
              se=F, method="lm") +
  labs(title="Impatto di Lunghezza, Sesso e Gestazione sul peso")

Creando uno scatter plot in tre dimensioni utilizzando la libreria plotly, viene visualizzata la relazione della variabile risposta Peso con entrambe le variabili esplicative Gestazione e Lunghezza, mantenendo due colori diversi per la variabile Sesso:

fig <- plot_ly(data = df, 
               x = ~Gestazione, y = ~Lunghezza, z = ~Peso, 
               color = ~Sesso,
               opacity=0.5,
               marker=list(size=5)) %>%   
  layout(title = "Impatto di Gestazione, Lunghezza e Sesso sul Peso", legend = list(title = list(text = "Sesso")))

fig

Inoltre, in questo modo è possibile utilizzare la dimensione dei punti per inserire anche una quinta variabile, ad esempio N.gravidanze:

fig <- plot_ly(data = df, 
               x = ~Gestazione, y = ~Lunghezza, z = ~Peso, 
               color = ~Sesso, size = ~N.gravidanze, sizes = c(50, 500),
               text = ~paste('N gravidanze:', N.gravidanze)) %>%   
  layout(title = "Impatto di Gestazione, Lunghezza, Sesso e Num. gravidenze sul Peso", legend = list(title = list(text = "Sesso")))

fig

7. Conclusioni

Questa analisi inferenziale su dati clinici di 2500 neonati ha permesso di identificare quali variabili sono più predittive del peso alla nascita. Inoltre, attraverso la creazione di un modello predittivo, è possibile stimare con precisione il peso previsto sulla base di valori misurabili delle variabili, consentendo di identificare rapidamente eventuali anomalie e pianificare in anticipo l’utilizzo di cure intensive.

L’analisi ha permesso anche di rispondere ad alcuni quesiti preliminari. In particolare si è concluso che non vi è una maggiore incidenza di parti cesarei in determinati ospedali considerati nell’analisi. Inoltre, è stato provato che non vi è una differenza significativa tra la media del peso dei neonati nei tre ospedali coinvolti rispetto a quella della popolazione, mentre la lunghezza è risultata significativamente inferiore. Infine, è stata confermata l’ipotesi che le misure antropometriche dei neonati (peso, lunghezza, diametro del cranio) sono significativamente diverse tra i due sessi.

Alcune conclusioni interessanti sono emerse dall’analisi di regressione multipla. Ad esempio, è risultato che fattori come il fumo materno e l’età della madre non influiscono significativamente sul peso alla nascita. Tuttavia, occorre sottolineare che ciò potrebbe essere dovuto alle caratteristiche di questo campione specifico, e l’inclusione di un maggior numero di ospedali nell’analisi potrebbe portare a risultati diversi e più generalizzabili.