Modello Statistico per la Previsione del Peso Neonatale

Raccolta dei Dati e Struttura del Dataset

Per costruire il modello predittivo, abbiamo raccolto dati su 2500 neonati provenienti da tre ospedali. Le variabili raccolte includono:

-Età della madre: Misura dell’età in anni.

-Numero di gravidanze: Quante gravidanze ha avuto la madre.

-Fumo materno: Un indicatore binario (0=non fumatrice, 1=fumatrice).

-Durata della gravidanza: Numero di settimane di gestazione.

-Peso del neonato: Peso alla nascita in grammi.

-Lunghezza e diametro del cranio: Lunghezza del neonato e diametro craniale, misurabili anche durante la gravidanza tramite ecografie.

-Tipo di parto: Naturale o cesareo.

-Ospedale di nascita: Ospedale 1, 2 o 3.

-Sesso del neonato: Maschio (M) o femmina (F).

L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature.

Analisi Preliminare

Per comprendere la distribuzione delle variabili e identificare eventuali valori outlier o anomali, eseguiamo un’analisi descrittiva delle variabili elencate precedentemente.

-Età della madre, Durata della gravidanza: varibili quantitative continue

-Numero di gravidanze: variabile quantitativa discreta

-Peso del neonato, Lunghezza e diametro del cranio: variabile quantitativa continua

-Fumo materno, Sesso del neonato: variabile qualitativa nominale, rappresentata come fattore

-Tipo di parto: variabile qualitativa nominale

-Ospedale di nascita: variabile qualitativa nominale

Osserviamo i dati a disposizione

dati <- read.csv("neonati.csv", stringsAsFactors = T)
library(knitr)
attach(dati)

summary(dati)

##    Anni.madre     N.gravidanze       Fumatrici        Gestazione   
##  Min.   : 0.00   Min.   : 0.0000   Min.   :0.0000   Min.   :25.00  
##  1st Qu.:25.00   1st Qu.: 0.0000   1st Qu.:0.0000   1st Qu.:38.00  
##  Median :28.00   Median : 1.0000   Median :0.0000   Median :39.00  
##  Mean   :28.16   Mean   : 0.9812   Mean   :0.0416   Mean   :38.98  
##  3rd Qu.:32.00   3rd Qu.: 1.0000   3rd Qu.:0.0000   3rd Qu.:40.00  
##  Max.   :46.00   Max.   :12.0000   Max.   :1.0000   Max.   :43.00  
##       Peso        Lunghezza         Cranio    Tipo.parto Ospedale   Sesso   
##  Min.   : 830   Min.   :310.0   Min.   :235   Ces: 728   osp1:816   F:1256  
##  1st Qu.:2990   1st Qu.:480.0   1st Qu.:330   Nat:1772   osp2:849   M:1244  
##  Median :3300   Median :500.0   Median :340              osp3:835           
##  Mean   :3284   Mean   :494.7   Mean   :340                                 
##  3rd Qu.:3620   3rd Qu.:510.0   3rd Qu.:350                                 
##  Max.   :4930   Max.   :565.0   Max.   :390

Osservando il summary della variabile Anni.madre notiamo che il valore minimo osservato per la variabile è 0, quindi un valore imputabile. La presenza di valori non plausibili può ostacolare la qualità dell’analisi che stiamo svolgendo, pertanto osserviamo quanti e quali sono i valori imputabili.

kable(dati[Anni.madre < 10, ])

	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
1152	1	1	0	41	3250	490	350	Nat	osp2	F
1380	0	0	0	39	3060	490	330	Nat	osp3	M

A questo punto abbiamo trovato le due osservazioni i cui valori di Anni.madre sono imputabili; procediamo con la loro rimozione dal dataset

dati <- dati[!(Anni.madre <= 1), ]

summary(dati)

##    Anni.madre     N.gravidanze       Fumatrici         Gestazione   
##  Min.   :13.00   Min.   : 0.0000   Min.   :0.00000   Min.   :25.00  
##  1st Qu.:25.00   1st Qu.: 0.0000   1st Qu.:0.00000   1st Qu.:38.00  
##  Median :28.00   Median : 1.0000   Median :0.00000   Median :39.00  
##  Mean   :28.19   Mean   : 0.9816   Mean   :0.04163   Mean   :38.98  
##  3rd Qu.:32.00   3rd Qu.: 1.0000   3rd Qu.:0.00000   3rd Qu.:40.00  
##  Max.   :46.00   Max.   :12.0000   Max.   :1.00000   Max.   :43.00  
##       Peso        Lunghezza         Cranio    Tipo.parto Ospedale   Sesso   
##  Min.   : 830   Min.   :310.0   Min.   :235   Ces: 728   osp1:816   F:1255  
##  1st Qu.:2990   1st Qu.:480.0   1st Qu.:330   Nat:1770   osp2:848   M:1243  
##  Median :3300   Median :500.0   Median :340              osp3:834           
##  Mean   :3284   Mean   :494.7   Mean   :340                                 
##  3rd Qu.:3620   3rd Qu.:510.0   3rd Qu.:350                                 
##  Max.   :4930   Max.   :565.0   Max.   :390

Adesso abbiamo 2498 osservazioni e possiamo iniziare la nostra analisi.

La variabile risposta del nostro caso di studio è la variabile Peso, che sta ad indicare il peso del neonato in grammi; verifichiamone la normalità attraverso gli indici di forma e di Shapiro-Wilk, e attraverso un’osservazione grafica.

n <- nrow(dati)

kable(moments::skewness(Peso))

x
-0.6470308

kable(moments::kurtosis(Peso)-3)

x
2.031532

shapiro.test(Peso)

## 
##  Shapiro-Wilk normality test
## 
## data:  Peso
## W = 0.97066, p-value < 2.2e-16

plot(density(Peso))

Osserviamo un p-value molto piccolo, pertanto si rifiuta l’ipotesi di normalità. Graficamente, però, osserviamo che la variabile risposta non è molto lontana da una distrbuzione normale, pertanto la costruzione di un modello di regressione classico potrebbe comunque essere una buona soluzione. Prima di iniziare la costruzione del modello, però, visualizziamo la matrice di correlazione per analizzare le relazioni tra i predittori.

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  #usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- (cor(x, y, method = "pearson"))
  txt <- format(c(r, 1), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = 1.5)
}
#correlazioni
pairs(dati,lower.panel=panel.cor, upper.panel=panel.smooth)

Osserviamo le correlazioni lineari tra la variabile risposta e le altre variabili esplicative.

Il coefficiente di correlazione lineare tra la variabile Peso e Anni.madre è negativo e pressochè nullo, questo implica che c’è poca correlazione tra le due variabili, infatti la nuvola è sparsa orizzontalmente. Situazione analoga per il numero di gravidanze che, però, presenta un coefficiente positivo ma pressocchè nullo.

Il coefficiente di correlazione lineare tra Peso e Gestazione, invece, è positivo ed è pari a 0.59. Questo implica che esiste una forte correlazione positiva tra le due variabili.

Notiamo, però, un pattern curvilineo della variabile Gestazione, che ci suggerisce un effetto quadratico che andremo a studiare durante la creazione del modello predittivo.

Situazione analoga per le variabili Lunghezza e Cranio, che presentano coefficienti positivi ancora più alti che stanno ad indicare una forte correlazione con la variabile risposta.

Per le variabili qualitative, invece, preferiamo osservare le correlazione con la variabile rispota tramite l’utilizzo di scatterplot.

Correlazione tra Peso e Fumatrici

par(mfrow=c(1,2))
boxplot(Peso)
boxplot(Peso~Fumatrici)

t.test(Peso~Fumatrici)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Fumatrici
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -45.61354 145.22674
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.153        3236.346

La mediana del peso nei neonati di madri fumatrici è leggermente più bassa rispetto a quella dei neonati di madri non fumatrici, questo ci suggerisce che il fumo possa influire negativamente sul neonato anche se non vi è una differenza così evidente. Inoltre, I neonati delle madri non fumatrici mostrano una distribuzione più ampia del peso, con alcuni valori outliers, superiori ed inferiori. I neonati delle madri fumatrici, invece, sembrano avere una distribuzione leggermente più compatta e spostata verso pesi più bassi, con meno neonati nella fascia di peso molto elevata. Questo suggerisce che il fumo materno potrebbe essere associato a un peso alla nascita più basso, con una riduzione della variabilità dei pesi.

Eseguiamo il test per saggiare l’ipotesi di ugualianza per medie di gruppi indipendenti

t.test(Peso~Fumatrici)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Fumatrici
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -45.61354 145.22674
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.153        3236.346

Il risultato del test è un p-value di 0.30 che non ci consente di rifiutare l’ipotesi nulla, quindi le medie non sono significativamente diverse.

Correlazione tra Peso e Sesso

Eseguiamo la stessa procedura, questa volta con la variabile Sesso.

par(mfrow=c(1,2))
boxplot(Peso)
boxplot(Peso~Sesso)

t.test(Peso~Sesso)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287.1051 -207.0615
## sample estimates:
## mean in group F mean in group M 
##        3161.132        3408.215

Notiamo che la mediana del peso dei neonati di sesso maschile è notevolemente più alta rispetto a quella dei neonati di sesso femminile. In entrambi i casi abbiamo dei valori dispersi con la presenza di outliers, quindi deduciamo che i neonati di sesso maschile pesino mediamente di più rispetto ai neonati di sesso femminile.

t.test(Peso~Sesso)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287.1051 -207.0615
## sample estimates:
## mean in group F mean in group M 
##        3161.132        3408.215

In questo caso, invece, il p-value risulta essere molto piccolo quindi possiamo rifiutare l’ipotesi nulla di ugualianza e affermare che le medie sono significativamente diverse.

Correlazione tra Peso e Tipo.parto

par(mfrow=c(1,2))
boxplot(Peso)
boxplot(Peso~Tipo.parto)

t.test(Peso~Tipo.parto)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Tipo.parto
## t = -0.12968, df = 1493, p-value = 0.8968
## alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
## 95 percent confidence interval:
##  -46.27992  40.54037
## sample estimates:
## mean in group Ces mean in group Nat 
##          3282.047          3284.916

La tipologia di parto non sembra influenzare significativamente il peso medio, ma il parto naturale mostra una maggiore dispersione dei pesi, con più casi di neonati molto piccoli o molto grandi rispetto a neonati nati con parto cesareo.

t.test(Peso~Tipo.parto)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Tipo.parto
## t = -0.12968, df = 1493, p-value = 0.8968
## alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
## 95 percent confidence interval:
##  -46.27992  40.54037
## sample estimates:
## mean in group Ces mean in group Nat 
##          3282.047          3284.916

In questo caso il p-value ci porta a non rifiutare l’ipotesi nulla di ugualianza, infatti possiamo notare che le medie sono pressochè uguali, quindi non significativamente diverse.

Correlazione tra Peso e Ospeadale

par(mfrow=c(1,2))
boxplot(Peso)
boxplot(Peso~Ospedale)

Notiamo che l’ospedale in cui nasce il neonato non è significativo per indicarne il peso. Infatti notiamo che le mediane sono all’incirca allo stesso livello, la dispersione dei valori sembra abbastanza simile, con presenza per la maggior parte di valori outliers inferiori, solo alcuni di valori sono superiori alla media.

Da quanto emerso da questa analisi delle variabili e delle loro correlazioni con la variabile risposta, possiamo affermare che le variabili Ospedale e Tipo.parto non sono significativamente utili per l’analisi che andremo a svolgere e neanche logicamente significative in ottica previsionale per perseguire l’obiettivo della nostra analisi.

Saggiare le ipotesi

1. In alcuni ospedali si fanno più parti cesarei

Per testare l’indipendenza o associazione delle due variabili categoriali usiamo il test del Chi-quadro. Il test del Chi-quadro testa l’indipendenza delle due variabili e prevede la creazione di una tabella di contingenza per incrociare le modalità delle due variabili e mostrare le frequenze assolute per ogni incrocio tra le modalità.

Iniziamo con la creazione della tabella di contingenza

tabella_contingenza <- table(Ospedale, Tipo.parto)
tabella_contingenza

##         Tipo.parto
## Ospedale Ces Nat
##     osp1 242 574
##     osp2 254 595
##     osp3 232 603

Cerchiamo un riscontro grafico rispetto alla tabella creata

#costruiamo il grafico baloon plot
ggpubr::ggballoonplot(data=as.data.frame(tabella_contingenza),
                      fill = "blue")

Dal ballonplot le proporzioni sembrano rispettarsi. Per una conferma, effettuiamo il test del Chi-quadro per saggiare l’ipotesi di indipendenza e fissiamo un valore soglia alfa dell’1%

test.indipendenza <- chisq.test(tabella_contingenza)
test.indipendenza

## 
##  Pearson's Chi-squared test
## 
## data:  tabella_contingenza
## X-squared = 1.0972, df = 2, p-value = 0.5778

La statistica calcolata è di 1.1 circa e fa riferimento ad una distribuzione di chi-quadro con 2 gradi di libertà e porta ad un p-valore di 0.54 pertanto non si rifiuta l’ipotesi di indipendenza.

Disegniamo la Chi-quadro con 2 gradi di libertà.

plot(density(rchisq(2500,2)), xlim = c(0,130))
#individuiamo il valore soglia sulla coda destra della distribuzione, corrispondente a un livello di significatività alfa del 1% 
abline(v=qchisq(0.99,2),col = 2)
#individuiamo la statistica test
points(x = 1.097, 0, cex = 3, col = 4, pch = 20)

Notiamo appunto che la chi-quadro ricade proprio nella zona di non rifiuto, quindi possiamo dire che, come avevamo immaginato dal ballonplot, le due variabili non sembrano essere associate in qualche modo.

2. La media del peso e della lunghezza di questo campione di neonati non sono significativamente diversi a quelle della popolazione

Per saggiare queste ipotesi utilizziamo il one-simple t-test, ovvero il t-test che confronta la media del campione (ovvero la nostra variabile Peso) con la media di un valore noto, ovvero la media della popolazione. Nonostante il peso sia un valore non lineare possiamo utilizzare il t-test in quanto abbiamo una numerosità campionaria elevata di 2500 osservazioni.

Dall’ Organizzazione Mondiale della Sanità (https://www.who.int) sappiamo che il peso medio dei neonati della popolazione è di 3300 grammi e che la lunghezza media è di 500mm, mentre la media del peso dei neonati del campione è di 3284 grammi.

Saggiamo l’ipotesi che la media del peso della popolazione (3300) non sia significativamente diversa dalla media del campione (3284) utilizzando il one-simple t-test, fissando un valore soglia alfa pari al 5%:

t.test(Peso,
       mu = 3300,
       conf.level = 0.95,
       alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  Peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.490 3304.672
## sample estimates:
## mean of x 
##  3284.081

Osserviamo che il p-value è significativamente superiore al valore soglia di 0.05 fissato, pertanto non rifiutiamo l’ipotesi nulla di ugualianza delle medie del peso.

Effettuiamo la stessa procedura con la variabile lunghezza e fissando un valore soglia alfa pari al 5%:

t.test(Lunghezza,
       mu = 500,
       conf.level = 0.95,
       alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  Lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6598 495.7242
## sample estimates:
## mean of x 
##   494.692

Osserviamo che il p-value è significativamente inferiore al valore soglia di 0.05 fissato, pertanto rifiutiamo l’ipotesi nulla di ugualianza delle medie della lunghezza.

Concludiamo quindi che la media del peso di questo campione di neonati non è significativamente diversa alla media del peso della popolazione; stessa cosa non può dirsi della media della lunghezza del campione che risulta significativamente diversa dalla media della lunghezza della popolazione.

3. Le misure antropometriche sono significativamente diverse tra i due sessi

Per saggiare questa ipotesi utilizziamo il T-test per gruppi indipendenti. Questo test confronta le medie di due gruppi indipendenti e verifica se la differenza tra di esse è significativamente diversa.

t.test(Peso ~ Sesso, data = dati, var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.115, df = 2488.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287.4841 -207.3844
## sample estimates:
## mean in group F mean in group M 
##        3161.061        3408.496

t.test(Lunghezza ~ Sesso, data = dati, var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  Lunghezza by Sesso
## t = -9.5823, df = 2457.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -11.939001  -7.882672
## sample estimates:
## mean in group F mean in group M 
##        489.7641        499.6750

t.test(Cranio ~ Sesso, data = dati, var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  Cranio by Sesso
## t = -7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.110504 -3.560417
## sample estimates:
## mean in group F mean in group M 
##        337.6231        342.4586

I test effettuati mostrano dei valori di p-value molto bassi, pertanto ci sentiamo confidenti nel rifiutare l’ipotesi nulla di indipendenza: le misure antropometriche nei neonati sono significativamente differenti tra maschi e femmine.

Creazione del modello

Come anticipato, la variabile risposta non si discosta molto da una distribuzione normale, pertanto creiamo il primo modello di regressione lineare multipla. Dal punto di vista logico per l’obiettivo dell’analisi, togliamo a priori le variabili Ospedale e Tipo.parto dal nostro modello.

mod1 <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + Sesso, data=dati)
summary(mod1)

## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Sesso, data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1160.6  -181.3   -15.7   163.6  2630.7 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6712.2405   141.3339 -47.492  < 2e-16 ***
## Anni.madre       0.8803     1.1491   0.766    0.444    
## N.gravidanze    11.3789     4.6767   2.433    0.015 *  
## Fumatrici      -30.3958    27.6080  -1.101    0.271    
## Gestazione      32.9472     3.8288   8.605  < 2e-16 ***
## Lunghezza       10.2316     0.3011  33.979  < 2e-16 ***
## Cranio          10.5198     0.4271  24.633  < 2e-16 ***
## SessoM          78.0787    11.2132   6.963 4.24e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2490 degrees of freedom
## Multiple R-squared:  0.7272, Adjusted R-squared:  0.7264 
## F-statistic: 948.3 on 7 and 2490 DF,  p-value: < 2.2e-16

Notiamo che:

• la variabile Gestazione è una variabile molto significativa, ed indica che per ogni settimana in più di gestazione il peso del neonato aumenta di 32.95g.

• la variabile Lunghezza è significativa per determinare il peso del neonato, infatti ci dice che per ogni cm in più di lunghezza il peso del neonato aumenta di 10.23g

• la variabile Cranio è anch’essa significativa e indica che per ogni mm di diametro del cranio, il peso del neonato aumenta di 10.52g

• la variabile Sesso è molto significativa per determinare il peso del neonato in quanto indica che, tenendo fissate le altre variabili, i maschi pesano 78.08g in più delle femmine.

Inoltre l’R quadro aggiustato è di 0.73 circa, quindi una variabilità spiegata del 73%.

Tra le variabili meno significative, invece, troviamo la variabile Anni.madre che abbiamo studiato essere poco correlata alla variabile Peso. Proviamo a rimuoverla dal nostro modello

mod2 <- update(mod1, ~.-Anni.madre)
summary(mod2)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1150.24  -181.32   -15.73   162.95  2635.69 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6682.2637   135.7983 -49.207  < 2e-16 ***
## N.gravidanze    12.6996     4.3470   2.921  0.00352 ** 
## Fumatrici      -30.5728    27.6048  -1.108  0.26818    
## Gestazione      32.6437     3.8079   8.573  < 2e-16 ***
## Lunghezza       10.2309     0.3011  33.979  < 2e-16 ***
## Cranio          10.5366     0.4265  24.707  < 2e-16 ***
## SessoM          78.1596    11.2117   6.971 4.01e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2491 degrees of freedom
## Multiple R-squared:  0.7271, Adjusted R-squared:  0.7265 
## F-statistic:  1106 on 6 and 2491 DF,  p-value: < 2.2e-16

Vediamo che l’R quadro aggiustato è rimasto pressochè invariato, così come la media delle stime è rimasta invariata. Proviamo a confrontare i due modelli creati

anova(mod2, mod1)

## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Sesso
## Model 2: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1   2491 187949505                           
## 2   2490 187905214  1     44292 0.5869 0.4437

Dal test Anova, osserviamo che il p-value è superiore al livello di significatività 0.05, pertanto l’aggiunta della variabile “Anni.madre” al modello non comporta un miglioramento significativo nella spiegazione della variabilità del peso del neonato, pertanto siamo confidenti nel rimuovere la variabile. Utilizziamo il criterio di informazione di Bayes per confrontare i due modelli e scegliere il migliore.

BIC(mod2, mod1)

##      df      BIC
## mod2  8 35200.24
## mod1  9 35207.48

Il BIC del modello 2 è inferiore al BIC del modello 1 pertanto continuiamo a preferire il modello semplificato. Infine eseguiamo un test per verificare la collinearità delle variabili esplicative.

car::vif(mod2)

## N.gravidanze    Fumatrici   Gestazione    Lunghezza       Cranio        Sesso 
##     1.026101     1.006621     1.676199     2.079728     1.624705     1.040400

Tutti i VIF sono inferiori a 5 pertanto non vi è rischio di collinearità. Per il rasoio di Occam scegliamo il modello più semplice per continuare la nostra analisi. Dall’ultimo modello creato e dall’analisi svolta in precedenza notiamo che la varibile fumatrici non sembra essere significativamente correlata alla variabile risposta. Proviamo a rimuoverla dal nostro modello per osservarne i risultati

mod3 <- update(mod2, ~.-Fumatrici)
summary(mod3)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.37  -180.98   -15.57   163.69  2639.09 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.7251   135.8036 -49.201  < 2e-16 ***
## N.gravidanze    12.4554     4.3416   2.869  0.00415 ** 
## Gestazione      32.3827     3.8008   8.520  < 2e-16 ***
## Lunghezza       10.2455     0.3008  34.059  < 2e-16 ***
## Cranio          10.5410     0.4265  24.717  < 2e-16 ***
## SessoM          77.9807    11.2111   6.956 4.47e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1327 on 5 and 2492 DF,  p-value: < 2.2e-16

Rimuovendo la variabile Fumatrici, otteniamo un R quadro aggiustato identico a quello del modello precedente; in media le statistiche significative non hanno perso significatività. Eseguiamo il test Anova

anova(mod3, mod2)

## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Sesso
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1   2492 188042054                           
## 2   2491 187949505  1     92548 1.2266 0.2682

Il risultato è 0.27 quindi superiore al livello di significatività di 0.05, quindi la presenza della variabile “Fumatrici” all’interno del modello non comporta un miglioramento significativo nella spiegazione della variabilità del peso del neonato, pertanto siamo confidenti nel rimuovere la variabile. Eseguiamo il criterio BIC

BIC(mod3, mod2, mod1)

##      df      BIC
## mod3  7 35193.65
## mod2  8 35200.24
## mod1  9 35207.48

Per il criterio BIC il modello migliore è il modello 3. Verifichiamo, infine, la presenza di collinearità o meno

car::vif(mod3)

## N.gravidanze   Gestazione    Lunghezza       Cranio        Sesso 
##     1.023462     1.669779     2.075747     1.624568     1.040184

Tutti i valori sono inferiori a 5 pertanto non vi è presenza di collinearità.

Osserviamo adesso alcuni effetti secondari notati durante l’osservazione degli scatterplot. Durante l’analisi abbiamo rilevato un pattern curvilineo della variabile Gestazione, suggerendo un effetto quadratico.

plot(Gestazione,Peso,pch=20)

Pertanto, potrebbe essere utile inserire all’interno del modello una nuova variabile, ovvero Gestazione^2.

mod4 <- update(mod3, ~.+I(Gestazione^2))
summary(mod4)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2), data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1144.0  -181.5   -12.9   165.8  2661.9 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -4646.7158   898.6322  -5.171 2.52e-07 ***
## N.gravidanze       12.5489     4.3381   2.893  0.00385 ** 
## Gestazione        -81.2309    49.7402  -1.633  0.10257    
## Lunghezza          10.3502     0.3040  34.045  < 2e-16 ***
## Cranio             10.6376     0.4282  24.843  < 2e-16 ***
## SessoM             75.7563    11.2435   6.738 1.99e-11 ***
## I(Gestazione^2)     1.5168     0.6621   2.291  0.02206 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.5 on 2491 degrees of freedom
## Multiple R-squared:  0.7276, Adjusted R-squared:  0.7269 
## F-statistic:  1109 on 6 and 2491 DF,  p-value: < 2.2e-16

La variabile Gestazione^2 è statisticamente significativa con p-value = 0.022, quindi porta informazione utile al modello. La variabile Gestazione però perde mediamente di significatività nel modello quadratico, forse a causa della collinearità tra le due variabili. L’R quadro aggiustato aumenta ma di pochissimo. Confrontiamo i due modelli tra di loro

BIC(mod4, mod3)

##      df      BIC
## mod4  8 35196.21
## mod3  7 35193.65

Il modello 3 risulta essere il modello migliore, quindi per il rasoio di Occam scegliamo il modello più semplice.

Analisi della qualità del modello

Procediamo, adesso, con l’analisi della qualità dei residui attraverso delle visualizzazioni grafiche

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

Il primo grafico mostra una nuvola relativamente orizzontale, con un leggero pattern curvo quindi maggiore dispersione per valori bassi e alti dei fitted con la presenza di valori outlier che spiccano: osservazioni 1551, 155 e 1306. Quindi globalmente è accettabile ma ci sono dei sospetti di eteroschedasticità.

Nel Q-Q plot i residui si posizionano sulla bisettrice del grafico ma la coda superiore è un po al di sopra della linea teorica, questo indica sicuramente la presenza di outliers; ad ogni modo la loro distribuzione non è troppo lontana da una distribuzione normale.

Nel grafico Scale-Location si vede una leggera tendenza crescente, segno che la varianza cresce con i valori previsti e questo indica la presenza di eteroschedasticità.

Nell’ultimo grafico notiamo qualche valore di leva ma non sembra superare molto la distanza di Cook, quindi ci sono pochi punti potenzialmente influenti, ma niente di troppo preoccupante.

Dopo aver effettuato un’analisi grafica, procediamo con i test statistici:

shapiro.test(residuals(mod3))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod3)
## W = 0.97414, p-value < 2.2e-16

lmtest::bptest(mod3) 

## 
##  studentized Breusch-Pagan test
## 
## data:  mod3
## BP = 90.297, df = 5, p-value < 2.2e-16

lmtest::dwtest(mod3)

## 
##  Durbin-Watson test
## 
## data:  mod3
## DW = 1.9532, p-value = 0.1209
## alternative hypothesis: true autocorrelation is greater than 0

plot(density(residuals(mod3)))

Come era stato preannunciato dall’analisi preliminare, anche i residui non si distribuiscono normalmente, inoltre il test di Breusch-Pagan mette in evidenza la presenza di eteroschedasticità, quindi la varianza dei residui non è costante. Invece il test di Durbin-Watson ci dice che i residui non sono correlati. Il modello nel complesso è ragionevole.

Indaghiamo adesso la presenza di valori leverage o outliers, che superano la distanza di cock, guardando i grafici

lev <- hatvalues(mod3)
plot(lev)
p = sum (lev)
soglia = 2 * p/n
abline(h=soglia, col=2)

kable(sum(lev>soglia))

x
152

Notiamo la presenza di molti valori di leva (152) che si trovano lontani dallo spazio dei regressori.

Visualizziamo, adesso, la presenza degli outliers

plot(rstudent(mod3))
abline(h=c(-2,2))

Notiamo qualche valore che supera i valori soglia. Effettuiamo il test statistico applicando la correzione di Bonferroni del p-value.

car::outlierTest(mod3)

##       rstudent unadjusted p-value Bonferroni p
## 1551 10.046230         2.6345e-23   6.5810e-20
## 155   5.025345         5.3818e-07   1.3444e-03
## 1306  4.824963         1.4848e-06   3.7092e-03

Come anticipato, sono tre le osservazioni più problematiche. Tra queste, solo l’osservazione 1551 supera di poco la distanza di Cook

cook <- cooks.distance(mod3)
plot(cook)

kable(max(cook))

x
0.8297645

Quindi, il punto 1551 è sicuramente critico, è sia un outlier per i residui sia molto influente secondo Cook, anche i punti 155 e 1306 meritano attenzione in quanto sono outlier di residuo, ma meno influenti rispetto all’osservazione 1551.

Possiamo concludere l’analisi della qualità del modello e dei residui dicendo che si tratta sicuramente di un modello accettabile, ma è importante porre l’attenzione sul fatto che:

• Ci sono tre punti potenzialmente influenti (1551, 155, 1306), con particolare nota all’osservazione 1551

• Il modello mostra qualche traccia di eteroschedasticità, probabilmente causata dalla presenza di valori outliers e dal fatto che i residui non siano perfettamente normali

• I residui non sono distribuiti normalmente, ma dato l’alto numero di osservazioni è accettabile.

Previsioni e Risultati

Stimiamo il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.

nuovo <- data.frame(
  Anni.madre = mean(Anni.madre, na.rm = TRUE),
  N.gravidanze = 3,
  Fumatrici = 0,  # metti quello che preferisci o media di categoria
  Gestazione = 39,
  Lunghezza = mean(Lunghezza, na.rm = TRUE),
  Cranio = mean(Cranio, na.rm = TRUE),
  Tipo.parto = "Nat",
  Ospedale = "osp1",
  Sesso = "F"
)

kable(predict(mod4, newdata = nuovo, interval = "prediction"))

fit	lwr	upr
3267.238	2728.524	3805.952

La stima del peso alla nascita per una neonata, figlia di una madre alla terza gravidanza e con una gestazione di 39 settimane, è di circa 3.267 grammi. L’intervallo di predizione al 95% è compreso tra 2.728 grammi e 3.805 grammi.

Visualizzazioni

Proseguiamo con la visualizzazione grafica per osservare i risultati del modello e mostrare le relazioni più significative tra le variabili.

Studiamo la relazione tra la variabile peso e la variabile gestazione: dalle analisi precedenti era emerso come all’aumentare delle settimane di gestazione della madre aumentasse anche il peso del neonato.

Creiamo un’evidenza grafica

library(ggplot2)
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")+
    labs(title = "Peso del neonato in funzione della gestazione e distinzione per sesso")+
  scale_color_discrete(name = "Sesso", labels = c("Femmine", "Maschi"))

## `geom_smooth()` using formula = 'y ~ x'

La variabile Gestazione è indicata con valori discreti in settimane, infatti sul grafico visualizziamo i valori allineati sui valori della variabile. Per dare un’impressione di continuità, modifichiamo il grafico per dare più credibilità alla variabile temporale Gestazione

ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col = Sesso), position= "jitter")+
  geom_smooth(aes(x=Gestazione,
                 y=Peso,
                 col = Sesso), se=F, method = "lm")+
  labs(title = "Peso del neonato in funzione della gestazione e distinzione per sesso")+
  scale_color_discrete(name = "Sesso", labels = c("Femmine", "Maschi"))

## `geom_smooth()` using formula = 'y ~ x'

Dal grafico ottenuto viene mostrata la relazione tra Peso e Gestazione divisa dal Sesso del neonato, infatti notiamo due rette di regressione.

Le due rette relative ai sessi mostrano una pendenza molto simile, indicando che l’effetto della gestazione sul peso è analogo sia per i maschi che per le femmine. Tuttavia, la retta dei maschi si colloca sistematicamente al di sopra di quella delle femmine, suggerendo che, a parità di settimane di gestazione, i neonati maschi tendono ad avere un peso medio più elevato rispetto alle femmine.

Aggiungiamo al grafico una retta di regressione generale non distinta per il sesso, in maniera tale da osservare il comportamento complessivo della relazione tra peso e gestazione.

ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col = Sesso), position= "jitter")+
  geom_smooth(aes(x=Gestazione,
                 y=Peso,
                 col = Sesso), se=F, method = "lm")+
  geom_smooth(aes(x=Gestazione,
                 y=Peso,), col = "black" , se=F, method = "lm")+
  labs(title = "Peso del neonato in funzione della gestazione e distinzione per sesso")+
  scale_color_discrete(name = "Sesso", labels = c("Femmine", "Maschi"))

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Notiamo che nelle prime settimane di gestazione, la retta generale è molto vicina alla retta delle femmine, il che potrebbe indicare che nella fase iniziale del periodo gestazionale il peso medio non differisce sostanzialmente tra maschi e femmine, o che la distribuzione dei dati è più densa per le femmine in queste settimane.

Con l’aumentare delle settimane di gestazione, però, la retta generale si colloca progressivamente in una posizione intermedia ma più vicina alla retta di regressione dei maschi, indicando che nelle settimane gestazionali più avanzate, i maschi tendono a pesare di più e a influenzare maggiormente la media generale.

Possiamo dire quindi che il sesso del neonato ha un impatto significativo sul peso alla nascita, con i maschi mediamente più pesanti rispetto alle femmine. Tuttavia, l’effetto della gestazione sul peso è consistente tra i sessi: entrambi mostrano un incremento simile del peso con l’aumentare delle settimane di gestazione.

Proviamo adesso a costruire un grafico differente, analizzando l’impatto del numero di settimane di gestazione e del fumo sul peso del neonato.

ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col = factor(Fumatrici)), position= "jitter")+
  geom_smooth(aes(x=Gestazione,
                 y=Peso,
                 col = factor(Fumatrici)), se=F, method = "lm")+
  
  labs(title = "Peso del neonato in funzione della gestazione e distinzione per fumatrici")+
  scale_color_discrete(name = "Fumatrici", labels = c("Non fumatrici", "Fumatrici"))

## `geom_smooth()` using formula = 'y ~ x'

summary(dati$Gestazione[dati$Fumatrici == 1])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   33.00   38.00   40.00   39.27   40.00   42.00

Il grafico mostra la relazione tra la durata della gestazione e il peso alla nascita, distinguendo tra madri fumatrici e non fumatrici. Anche in questo caso osserviamo due rette di regressione distinte. La retta relativa alle non fumatrici si colloca al di sopra di quella delle fumatrici, indicando che i neonati di madri non fumatrici tendono ad avere un peso superiore rispetto a quelli di madri fumatrici, a parità di settimana gestazionale. Ci sono, però, alcune implicazioni. Infatti il fumo materno durante la gravidanza appare associato a un ridotto peso alla nascita. Questo effetto è visibile lungo tutto l’intervallo di settimane di gestazione rappresentato, dalla 33esima settimana in poi, confermando il possibile impatto negativo del fumo sullo sviluppo fetale.

Conclusioni

L’analisi approfondita dei dati clinici ci ha permesso di identificare diversi fattori determinanti del peso neonatale, confermando l’importanza di questo progetto nell’ottica di Neonatal Health Solutions per il miglioramento della cura prenatale e l’ottimizzazione delle risorse cliniche.

Gestazione

La durata della gestazione si conferma come il principale predittore del peso alla nascita, infatti all’aumentare delle settimane di gestazione, si osserva un chiaro incremento del peso neonatale. Il risultato conferma ciò che ci si aspetta dal punto di vista fisiologico: i neonati nati a termine tendono ad avere un peso maggiore e, di conseguenza, un rischio inferiore di problemi alla nascita. Per questo motivo, dal punto di vista clinico, sarebbe utile concentrarsi su interventi che prevengano i parti prematuri, in modo da migliorare direttamente la salute dei neonati.

Sesso

I neonati di sesso maschile tendono a pesare di più rispetto alle femmine a parità di condizioni, questa differenza è statisticamente significativa e riflette tendenze biologiche osservate. Dal punto di vista clinico, conoscere il sesso del neonato può essere di fondamentale importanza per l’accuratezza nella stima del peso fetale durante il monitoraggio prenatale.

Fumo materno

Le madri fumatrici presentano neonati con un peso medio inferiore rispetto alle non fumatrici. Nel grafico elaborato, le rette di regressione separate per fumatrici e non fumatrici mostrano un trend chiaro: pur mantenendo una crescita con la gestazione, la curva delle fumatrici si posiziona sistematicamente più in basso. Inoltre, le gravidanze di madri fumatrici sembrano più concentrate nelle settimane più avanzate di gestazione. Dal punto di vista clinico sarebbe di grande impatto la promozione di programmi di cessazione del fumo in gravidanza come una priorità strategica. La sensibilizzazione delle pazienti sull’impatto negativo del fumo sul peso neonatale può contribuire significativamente alla prevenzione di complicanze neonatali.