Elaborazione di un modello statistico predittivo per il peso neonatale.

Otteniamo dalla Neonatal Health Solutions un dataset, popolato di dati a disposizione di tre ospedali gestiti dalla compagnia, per l’elaborazione di un modello predittivo che riesca a stimare il peso neonatale in base alle variabili prese in considerazione nel dataset stesso.

Analisi preliminare del dataset e studio descrittivo delle variabili.

Carichiamo il dataset in un data frame e stampiamo le prime righe per avere un’idea degli oggetti e delle variabili presenti.

neonatalcare_data = read.csv("neonati.csv")

str(neonatalcare_data)
## 'data.frame':    2500 obs. of  10 variables:
##  $ Anni.madre  : int  26 21 34 28 20 32 26 25 22 23 ...
##  $ N.gravidanze: int  0 2 3 1 0 0 1 0 1 0 ...
##  $ Fumatrici   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Gestazione  : int  42 39 38 41 38 40 39 40 40 41 ...
##  $ Peso        : int  3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
##  $ Lunghezza   : int  490 490 500 515 480 495 480 510 500 510 ...
##  $ Cranio      : int  325 345 375 365 335 340 345 349 335 362 ...
##  $ Tipo.parto  : chr  "Nat" "Nat" "Nat" "Nat" ...
##  $ Ospedale    : chr  "osp3" "osp1" "osp2" "osp2" ...
##  $ Sesso       : chr  "M" "F" "M" "M" ...
head(neonatalcare_data)
##   Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1         26            0         0         42 3380       490    325        Nat
## 2         21            2         0         39 3150       490    345        Nat
## 3         34            3         0         38 3640       500    375        Nat
## 4         28            1         0         41 3690       515    365        Nat
## 5         20            0         0         38 3700       480    335        Nat
## 6         32            0         0         40 3200       495    340        Nat
##   Ospedale Sesso
## 1     osp3     M
## 2     osp1     F
## 3     osp2     M
## 4     osp2     M
## 5     osp3     F
## 6     osp2     F

Il dataframe è stato popolato da 2500 osservazioni per 10 variabili differenti. Esse riguardano le rilevazioni inerenti le nascite avvenute nei tre ospedali della compagnia, sia per il nascituro che per le madri.

In particolare abbiamo a disposizione informazioni grezze per:

  • Anni.Madre - Rilevazione dell’età della madre al momento del parto - Variabile quantitativa continua su scala di rapporti;

  • N.gravidanze - Numero di gravidanze avuto dalla madre precedentemente al parto in oggetto - Variabile quantitativa discreta su scala di rapporti;

  • Fumatrici - Osservazione riguardante lo stato di fumatrice della madre. Rilevazione codificata secondo lo schema No/Si -> 0/1 - Variabile qualitativa codificata su scala nominale;

  • Gestazione - Tempo di gestazione del feto espresso in settimane - Variabile quantitativa continua su scala di rapporti;

  • Peso - Peso del neonato, rilevato alla nascita, espresso in grammi - Variabile quantitativa continua su scala di rapporti;

  • Lunghezza - Lunghezza del neonato alla nascita, rilevata in millimetri - Variabile quantitativa continua su scala di rapporti;

  • Cranio - diametro del cranio del neonato rilevato alla nascita, espresso in millimetri - Variabile quantitativa continua su scala di rapporti;

  • Tipo.parto - Rilevazione della tipologia di parto del neonato, Naturale o Cesareo, indicato con “Nat” o “Ces” - Variabile qualitativa su scala nominale;

  • Ospedale - Rilevazione dell’ospedale in cui è avvenuto il parto. Codificato con relativa numerazione “osp1”, “osp2”, “osp3” - Variabile qualitativa su scala nominale;

  • Sesso - Rilevazione del sesso del neonato - codificata con “M” per i maschi e “F” per le femmine - Variabile qualitativa su scala nominale;

Sulla base di questa analisi preliminare sulla tipologia delle variabili, identifichiamo in Peso la nostra variabile risposta, mentre le altre saranno variabili esplicative nel nostro modello.

Analisi descrittiva delle variabili.

Anni.madre:

attach(neonatalcare_data)

agen_distribution = table(Anni.madre)
agen_frtable = agen_distribution/sum(agen_distribution)
age_cumsum = cumsum(agen_distribution)
age_cumrel = age_cumsum/sum(agen_distribution)

final_agetable = cbind(ni=agen_distribution, fi=agen_frtable, Ni=age_cumsum, Fi=age_cumrel)
colnames(final_agetable) = c("Fr.Assolute", "Fr.Relative", "Fr. Assolute Cumulate", "Fr. relative cumulate")
knitr::kable(final_agetable)
Fr.Assolute Fr.Relative Fr. Assolute Cumulate Fr. relative cumulate
0 1 0.0004 1 0.0004
1 1 0.0004 2 0.0008
13 1 0.0004 3 0.0012
14 2 0.0008 5 0.0020
15 6 0.0024 11 0.0044
16 13 0.0052 24 0.0096
17 18 0.0072 42 0.0168
18 24 0.0096 66 0.0264
19 45 0.0180 111 0.0444
20 66 0.0264 177 0.0708
21 74 0.0296 251 0.1004
22 100 0.0400 351 0.1404
23 115 0.0460 466 0.1864
24 131 0.0524 597 0.2388
25 180 0.0720 777 0.3108
26 184 0.0736 961 0.3844
27 197 0.0788 1158 0.4632
28 172 0.0688 1330 0.5320
29 174 0.0696 1504 0.6016
30 200 0.0800 1704 0.6816
31 147 0.0588 1851 0.7404
32 159 0.0636 2010 0.8040
33 110 0.0440 2120 0.8480
34 96 0.0384 2216 0.8864
35 66 0.0264 2282 0.9128
36 64 0.0256 2346 0.9384
37 41 0.0164 2387 0.9548
38 38 0.0152 2425 0.9700
39 27 0.0108 2452 0.9808
40 19 0.0076 2471 0.9884
41 13 0.0052 2484 0.9936
42 8 0.0032 2492 0.9968
43 2 0.0008 2494 0.9976
44 4 0.0016 2498 0.9992
45 1 0.0004 2499 0.9996
46 1 0.0004 2500 1.0000 
summary(Anni.madre)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   25.00   28.00   28.16   32.00   46.00
library("moments")
mother_age_data = data.frame(
  Asimmetria = round(skewness(Anni.madre), 2), 
  Curtosi = round(kurtosis(Anni.madre)-3, 2),
  Deviazione_Standard = round(sd(Anni.madre), 2)
)

knitr::kable((mother_age_data), "pipe")
Asimmetria Curtosi Deviazione_Standard
0.04 0.38 5.27

Osservando un riepilogo degli indici di posizione, variabilità e forma di Anni.madre, nonchè una tabella di frequenza delle osservazioni, notiamo subito come fra le osservazioni della variabile siano state segnate due nascite anomale, probabilmente frutto di un errore di trascrizione, relativamente a madri partorienti a 0 e 1 anno di età.

Non considerando questi due dati, quindi, avremmo avuto un min per la variabile Anni.Madre di 13, invece che 0. Lo scostamento della media dalla mediana (28,16 contro 28) non è comunque tale da farci pensare che i due valori outlier possano aver influito sugli indici di posizione in maniera significativa. Come mostra la tabella di frequenza, 0 e 1 hanno delle frequenze relative minori dello 0.000 sulle n osservazioni.

Procediamo col visualizzare la curva di densità della distribuzione per Anni.Madre:

library("ggplot2")
ggplot()+
  geom_density(aes(x=Anni.madre), col="black", fill="turquoise")+
  geom_vline(xintercept = mean(Anni.madre), linewidth=1.5, color = "lightcoral")+
  labs(x="Età della partoriente", y="Densità", title = "Variabile Anni.madre", subtitle = "Curva di densità della distribuzione")+
  geom_text(aes(x = mean(Anni.madre) + 5, y = 0.075, label = sprintf("Media\n %.3f", mean(Anni.madre))), inherit.aes = FALSE) +
  ylim(0, 0.08)+
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

Con un’asimmetria di 0.04 la curva di densità della distribuzione della variabile Anni.Madre presenta una tendenza quasi nulla all’asimmetria positiva, con tuttavia una curtosi a 0.38 che indica una curva leptocurtica.

N.gravidanze:

birthn_distribution = table(N.gravidanze)
birthn_frtable = birthn_distribution/sum(birthn_distribution)
birth_cumsum = cumsum(birthn_distribution)
birth_cumrel = (birth_cumsum/sum(birthn_distribution))

final_birthtable = cbind(ni=birthn_distribution, fi=birthn_frtable, Ni=birth_cumsum, Fi=birth_cumrel)
colnames(final_birthtable) = c("Fr.Assolute", "Fr.Relative", "Fr. Assolute Cumulate", "Fr. relative cumulate")
knitr::kable(final_birthtable)
Fr.Assolute Fr.Relative Fr. Assolute Cumulate Fr. relative cumulate
0 1096 0.4384 1096 0.4384
1 818 0.3272 1914 0.7656
2 340 0.1360 2254 0.9016
3 150 0.0600 2404 0.9616
4 48 0.0192 2452 0.9808
5 21 0.0084 2473 0.9892
6 11 0.0044 2484 0.9936
7 1 0.0004 2485 0.9940
8 8 0.0032 2493 0.9972
9 2 0.0008 2495 0.9980
10 3 0.0012 2498 0.9992
11 1 0.0004 2499 0.9996
12 1 0.0004 2500 1.0000 
summary(N.gravidanze)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  1.0000  0.9812  1.0000 12.0000
prior_birth_data = data.frame(
  Asimmetria = round(skewness(N.gravidanze), 2), 
  Curtosi = round(kurtosis(N.gravidanze)-3, 2),
  Deviazione_Standard = round(sd(N.gravidanze), 2)
)
knitr::kable((prior_birth_data), "pipe")
Asimmetria Curtosi Deviazione_Standard
2.51 10.99 1.28

Come si evince dall’analisi della variabile N.gravidanze, la maggior parte delle osservazioni (76,56%), per il numero di gravidanze sostenute dalla madre al momento del parto, si concentra entro 1, con osservazioni gradualmente inferiori mentre ci spostiamo verso il maggior numero di gravidanze pregresse (12).

Ci aspettiamo una distribuzione spostata verso sinistra, come indicato dall’asimmetria positiva di 2.51, con valori bassi prevalenti.

ggplot() +
  geom_bar(
    aes(x = as.factor(N.gravidanze)),
    fill = "lightblue", color = "black", alpha = 0.7
  ) +
  labs(
    x = "Numero di gravidanze",
    y = "Frequenza",
    title = "Variabile N.gravidanze",
    subtitle = "Distribuzione delle gravidanze"
  ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

Trattandosi di una variabile discreta utilizziamo un barplot per visualizzare la distribuzione, ottenendo conferma grafica della prevalenza di 0 e 1 gravidanze nelle osservazioni totali e di una forte asimmetria positiva.

Fumatrici:

smokerstable = table(Fumatrici)
knitr::kable((smokerstable), "pipe")
Fumatrici Freq
0 2396
1 104 
smokersperc = round((104*100)/2500, 2)
knitr::kable(smokersperc)
x
4.16 
ggplot() +
  geom_bar(
    aes(x = as.factor(Fumatrici)),
    fill = "lavenderblush4", color = "black", alpha = 0.7
  ) +
  labs(
    x = "Non fumatrici/Fumatrici",
    y = "Osservazioni",
    title = "Variabile Fumatrici",
    subtitle = "Abitudini di tabagismo in gravidanza sulle n osservazioni"
  ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

Per quanto riguarda la variabile Fumatrici, rileviamo come siano osservate 2396 partorienti non fumatrici contro 104 partorienti fumatrici. Come visibile nella seconda tabella, il dato “fumatrici” costituisce un 4,16% delle totali osservazioni.

Lasceremo al personale sanitario valutazioni più approfondite sull’argomento, tuttavia un’informativa del Ministero della Salute datata 29 maggio 2024 indica che i figli di madri fumatrici hanno un eccesso di rischio del 70% di avere malattie alle basse vie respiratorie.

Una rilevazione statistica del 2013 su un campione di n=4953 donne, condotta dall’Istituto Superiore della Sanità, scopre inoltre che su una componente del 23% di fumatrici sul campione, solo il 70% ha smesso di fumare durante la gravidanza, facendo riflettere su una potenziale sottostima delle osservazioni codificate con 1 nel nostro dataset. La ricerca non è abbastanza recente da poter essere valida nel 2024, ma potrebbe comunque indicare una tendenza dei dati a nostra disposizione. Andrebbe effettuato un controllo con serie storica dal 2013 in poi per rilevare se l’andamento generale delle tendenze delle fumatrici oggi è compatibile con i dati a nostra disposizione, rapportandolo alla popolazione di riferimento del contesto dei tre ospedali e verificando se le tendenze locali sono correlate o meno con quelle della popolazione nazionale.

Lo stesso rapporto dell’ISS riporta anche fattori positivi sulla permanenza prolungata delle buone pratiche di salute da parte delle partorienti in seguito a gravidanza. Si consiglia ulteriore studio su questo fattore di rischio.

Gestazione:

pln_distribution = table(Gestazione)
pln_frtable = pln_distribution/sum(pln_distribution)
pl_cumsum = cumsum(pln_distribution)
pl_cumrel = pl_cumsum/sum(pln_distribution)

final_pltable = cbind(ni=pln_distribution, fi=pln_frtable, Ni=pl_cumsum, Fi=pl_cumrel)
colnames(final_pltable) = c("Fr.Assolute", "Fr.Relative", "Fr. Assolute Cumulate", "Fr. relative cumulate")
knitr::kable(final_pltable)
Fr.Assolute Fr.Relative Fr. Assolute Cumulate Fr. relative cumulate
25 1 0.0004 1 0.0004
26 1 0.0004 2 0.0008
27 2 0.0008 4 0.0016
28 4 0.0016 8 0.0032
29 3 0.0012 11 0.0044
30 5 0.0020 16 0.0064
31 8 0.0032 24 0.0096
32 9 0.0036 33 0.0132
33 18 0.0072 51 0.0204
34 16 0.0064 67 0.0268
35 33 0.0132 100 0.0400
36 62 0.0248 162 0.0648
37 192 0.0768 354 0.1416
38 437 0.1748 791 0.3164
39 581 0.2324 1372 0.5488
40 741 0.2964 2113 0.8452
41 329 0.1316 2442 0.9768
42 56 0.0224 2498 0.9992
43 2 0.0008 2500 1.0000 
summary(Gestazione)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   25.00   38.00   39.00   38.98   40.00   43.00
pregnancy_length_data = data.frame(
  Asimmetria = round(skewness(Gestazione), 2), 
  Curtosi = round(kurtosis(Gestazione)-3, 2),
  Deviazione_Standard = round(sd(Gestazione), 2)
)
knitr::kable(pregnancy_length_data)
Asimmetria Curtosi Deviazione_Standard
-2.07 8.26 1.87

Stando alla definizione di “Neonato pretermine” fornita dall’Ospedale Bambin Gesù di Roma, è necessario soffermarsi nel dataset a valutare l’incidenza di nati al di sotto delle 37 settimane per identificare i casi potenzialmente meritevoli di attenzione e intervento sanitario.

Per le osservazioni nel campione a disposizione, il 14,16% delle gravidanze ha una durata pari o inferiore a 37 settimane, mentre la durata media di una gravidanza è di 39 settimane (38,98) e il minimo della distribuzione è di 25. L’84,52% dei casi osservati rientra nelle 40 settimane.

Gli indici di forma restituiscono un’asimmetria negativa con prevalenza di valori alti e una leptocurtosi della curva di densità della distribuzione. Per via della natura della variabile preferiamo mostrarlo con un barplot.

ggplot() +
  geom_bar(
    aes(x = as.factor(Gestazione)),
    fill = "wheat", color = "black", alpha = 0.7
  ) +
  labs(
    x = "Settimane di Gestazione",
    y = "Frequenza",
    title = "Variabile Gestazione",
    subtitle = "Distribuzione delle settimane di gestazione"
  ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

Peso:

weight_class = cut(Peso, breaks = seq(500,5000, by=500), 
                   left=FALSE,
                   labels = paste0("(", seq(500,4500, by=500), "-", seq(1000,5000, by=500), "]")
                   )
weightn_distribution = table(weight_class)
weightn_frtable = weightn_distribution/sum(weightn_distribution)
weight_cumsum = cumsum(weightn_distribution)
weight_cumrel = weight_cumsum/sum(weightn_distribution)

final_weighttable = cbind(ni=weightn_distribution, fi=weightn_frtable, Ni=weight_cumsum, Fi=weight_cumrel)
colnames(final_weighttable) = c("Fr.Assolute", "Fr.Relative", "Fr. Assolute Cumulate", "Fr. relative cumulate")
knitr::kable(final_weighttable)
Fr.Assolute Fr.Relative Fr. Assolute Cumulate Fr. relative cumulate
(500-1000] 6 0.0024 6 0.0024
(1000-1500] 18 0.0072 24 0.0096
(1500-2000] 29 0.0116 53 0.0212
(2000-2500] 94 0.0376 147 0.0588
(2500-3000] 510 0.2040 657 0.2628
(3000-3500] 1016 0.4064 1673 0.6692
(3500-4000] 662 0.2648 2335 0.9340
(4000-4500] 147 0.0588 2482 0.9928
(4500-5000] 18 0.0072 2500 1.0000 
summary(Peso)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     830    2990    3300    3284    3620    4930
weight_data = data.frame(
  Asimmetria = round(skewness(Peso), 2), 
  Curtosi = round(kurtosis(Peso)-3, 2),
  Deviazione_Standard = round(sd(Peso), 2)
)
knitr::kable(weight_data)
Asimmetria Curtosi Deviazione_Standard
-0.65 2.03 525.04

Per quanto riguarda l’analisi preliminare della variabile Peso, considerando che, secondo le linee guida poco sopra citate dell’Ospedale Bambin Gesù di Roma, si considera potenzialmente pretermine (o prematuro) un neonato che abbia un peso pari o inferiore a 2500g, preferiamo effettuare una divisione in classi, chiuse a destra, per costruire una tabella di frequenze che evidenzi come le osservazioni al di sotto dei 2500g risultino essere il 5,8% del campione totale, meritevoli di considerazione da parte del personale sanitario.

Non vi è un grosso scostamento fra media e mediana, indice che non dovremmo trovarci di fronte a valori outlier che influenzano i risultati. L’asimmetria negativa ci indica una tendenza della distribuzione ad assumere valori alti, con una curva leptocurtica.

Costruiamo la curva di densità della distribuzione:

ggplot()+
  geom_density(aes(x=Peso), col="black", fill="thistle")+
  geom_vline(xintercept = mean(Peso), linewidth=1.5, color = "lightcoral")+
  labs(x="Peso del neonato alla nascita (g)", y="Densità", title = "Variabile Peso", subtitle = "Curva di densità della distribuzione")+
  geom_text(aes(x = mean(Peso) + 1000, y = 8e-04, label = sprintf("Media\n %.3f", mean(Peso))), inherit.aes = FALSE) +
  ylim(0, 9e-04)+
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

Le osservazioni precedenti sono confermate dal grafico.

Visto che la Variabile peso sarà la nostra variabile risposta nel modello di regressione che andremo a sviluppare in una fase successiva, effettuiamo subito uno shapiro test per capire se ci troviamo davanti ad una normale.

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

Con un p-value praticamente pari a zero e un risultato W di 0.97, dobbiamo rifiutare l’ipotesi H0 di normalità. Per lo shapiro test Peso non segue una normale, tuttavia il teorema del limite centrale potrebbe mitigare l’effetto della non normalità per dataset con un elevato numero di osservazioni. Valuteremo in seguito gli effetti sul modello.

Lunghezza:

length_class = cut(Lunghezza, breaks = seq(300,570, by=30), 
                   left=FALSE,
                   labels = paste0("(", seq(300,540, by=30), "-", seq(330,570, by=30), "]")
                   )
lengthn_distribution = table(length_class)
lengthn_frtable = lengthn_distribution/sum(lengthn_distribution)
length_cumsum = cumsum(lengthn_distribution)
length_cumrel = length_cumsum/sum(lengthn_distribution)

final_lengthtable = cbind(ni=lengthn_distribution, fi=lengthn_frtable, Ni=length_cumsum, Fi=length_cumrel)
colnames(final_lengthtable) = c("Fr.Assolute", "Fr.Relative", "Fr. Assolute Cumulate", "Fr. relative cumulate")
knitr::kable(final_lengthtable)
Fr.Assolute Fr.Relative Fr. Assolute Cumulate Fr. relative cumulate
(300-330] 4 0.0016 4 0.0016
(330-360] 6 0.0024 10 0.0040
(360-390] 12 0.0048 22 0.0088
(390-420] 23 0.0092 45 0.0180
(420-450] 84 0.0336 129 0.0516
(450-480] 524 0.2096 653 0.2612
(480-510] 1335 0.5340 1988 0.7952
(510-540] 471 0.1884 2459 0.9836
(540-570] 41 0.0164 2500 1.0000 
summary(Lunghezza)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   310.0   480.0   500.0   494.7   510.0   565.0
length_data = data.frame(
  Asimmetria = round(skewness(Lunghezza), 2), 
  Curtosi = round(kurtosis(Lunghezza)-3, 2),
  Deviazione_Standard = round(sd(Lunghezza), 2)
)
knitr::kable(length_data)
Asimmetria Curtosi Deviazione_Standard
-1.51 6.49 26.32

La maggior parte delle osservazioni della variabile Lunghezza si concentrano nella classe 480-510 (53,4%). La lunghezza in mm media di un neonato maschio alla nascita è solitamente 500mm mentre è 495mm per le femmine. Il nostro dataset rispecchia questa tendenza rilevata in altri studi. La media totale della lunghezza dei neonati alla nascita è di 494,7mm, con una mediana di 500mm e uno scostamento minimo per i due indici.

Abbiamo un’asimmetria negativa con una prevalenza di valori alti e una curva leptocurtica.

Costruiamo la curva di densità della distribuzione.

ggplot()+
  geom_density(aes(x=Lunghezza), col="black", fill="springgreen4")+
  geom_vline(xintercept = mean(Lunghezza), linewidth=1.5, color = "lightcoral")+
  labs(x="Lunghezza del neonato alla nascita (mm)", y="Densità", title = "Variabile Lunghezza", subtitle = "Curva di densità della distribuzione")+
  geom_text(aes(x = mean(Lunghezza) -30, y = 0.02, label = sprintf("Media\n %.3f", mean(Lunghezza))), inherit.aes = FALSE) +
  ylim(0, 0.025)+
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

Cranio:

head_class = cut(Cranio, breaks = seq(220,400, by=20), 
                   left=FALSE,
                   labels = paste0("(", seq(220,380, by=20), "-", seq(240,400, by=20), "]")
)
headn_distribution = table(head_class)
headn_frtable = headn_distribution/sum(headn_distribution)
head_cumsum = cumsum(headn_distribution)
head_cumrel = head_cumsum/sum(headn_distribution)

final_headtable = cbind(ni=headn_distribution, fi=headn_frtable, Ni=head_cumsum, Fi=head_cumrel)
colnames(final_headtable) = c("Fr.Assolute", "Fr.Relative", "Fr. Assolute Cumulate", "Fr. relative cumulate")
knitr::kable(final_headtable)
Fr.Assolute Fr.Relative Fr. Assolute Cumulate Fr. relative cumulate
(220-240] 1 0.0004 1 0.0004
(240-260] 3 0.0012 4 0.0016
(260-280] 16 0.0064 20 0.0080
(280-300] 31 0.0124 51 0.0204
(300-320] 201 0.0804 252 0.1008
(320-340] 1057 0.4228 1309 0.5236
(340-360] 1006 0.4024 2315 0.9260
(360-380] 174 0.0696 2489 0.9956
(380-400] 11 0.0044 2500 1.0000 
summary(Cranio)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     235     330     340     340     350     390
head_data = data.frame(
  Asimmetria = round(skewness(Cranio), 2), 
  Curtosi = round(kurtosis(Cranio)-3, 2),
  Deviazione_Standard = round(sd(Cranio), 2)
)
knitr::kable(head_data)
Asimmetria Curtosi Deviazione_Standard
-0.79 2.95 16.43

Per la variabile Cranio, che rileva il diametro della testa del neonato alla nascita (o ancora in gestazione tramite ecografie) in mm, la maggior parte delle osservazioni si concentrano sotto i 360mm (93%). Il minimo 235mm e il massimo 390mm rientrano nella divisione in classi correttamente.

Mediana e Media coincidono su 340mm. La curva ha una asimmetria negativa con valori in prevalenza alti, e una curtosi positiva di 2.95 per una curva leptocurtica.

Costruiamo il grafico di densità della curva di distribuzione.

ggplot()+
  geom_density(aes(x=Cranio), col="black", fill="mintcream")+
  geom_vline(xintercept = mean(Cranio), linewidth=1.5, color = "lightcoral")+
  labs(x="Diametro del cranio del neonato alla nascita (mm)", y="Densità", title = "Variabile Cranio", subtitle = "Curva di densità della distribuzione")+
  geom_text(aes(x = mean(Cranio) -20, y = 0.027, label = sprintf("Media\n %.3f", mean(Cranio))), inherit.aes = FALSE) +
  ylim(0, 0.03)+
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

Tipo.parto:

birth_type = table(Tipo.parto)
knitr::kable(birth_type)
Tipo.parto Freq
Ces 728
Nat 1772 
ggplot() +
  geom_bar(
    aes(x = as.factor(Tipo.parto)),
    fill = "papayawhip", color = "black", alpha = 0.7
  ) +
  labs(
    x = "Tipologia di Parto",
    y = "Osservazioni totali",
    title = "Variabile Tipo.parto",
    subtitle = "Osservazioni sulla tipologia del parto"
  ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

birthtypeperc = round((728*100)/2500, 2)
knitr::kable(birthtypeperc)
x
29.12

Dai dati osservati nei tre ospedali, confermiamo che il 29.12% delle nascite avviene con parto cesareo, mentre le restanti con parto naturale.

Ospedale:

hospital = table(Ospedale)
knitr::kable(hospital)
Ospedale Freq
osp1 816
osp2 849
osp3 835 
ggplot() +
  geom_bar(
    aes(x = as.factor(Ospedale)),
    fill = "red", color = "black", alpha = 0.7
  ) +
  labs(
    x = "Ospedale",
    y = "Osservazioni totali",
    title = "Variabile Ospedale",
    subtitle = "Distribuzione delle nascite negli ospedali"
  ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

La variabile Ospedali presenta un numero di osservazioni molto simile per tutte e tre le strutture, senza significative variazioni sul numero di osservazioni considerato singolarmente.

Sesso:

gender = table(Sesso)
knitr::kable(gender)
Sesso Freq
F 1256
M 1244 
ggplot() +
  geom_bar(
    aes(x = as.factor(Sesso)),
    fill = "mediumslateblue", color = "black", alpha = 0.7
  ) +
  labs(
    x = "Sesso del nascituro",
    y = "Osservazioni totali",
    title = "Variabile Sesso",
    subtitle = "Distribuzione del sesso del nascituro"
  ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

Non si osservano significative differenze fra le nascite di maschi e femmine nel campione di 2500 nati raccolto.

Analisi delle relazioni fra le variabili e il peso del neonato.

Peso dei neonati e relazione con il genere sessuale:

Prima di valutare le relazioni fra le variabili materne e il peso del neonato, vale innanzitutto la pena verificare se vi è una naturale differenza di peso fra nati di genere maschile e femminile. Costruiamo il grafico di seguito

library(dplyr)
male_stats = neonatalcare_data %>%
  filter(Sesso == "M") %>%
  summarise(
    Media = mean(Peso, na.rm = TRUE),
    Mediana = median(Peso, na.rm = TRUE),
    Deviazione_Standard = sd(Peso, na.rm = TRUE),
    Minimo = min(Peso, na.rm = TRUE),
    Massimo = max(Peso, na.rm = TRUE),
    Quartile1 = quantile(Peso, 0.25, na.rm = TRUE),
    Quartile3 = quantile(Peso, 0.75, na.rm = TRUE)
  )

female_stats = neonatalcare_data %>%
  filter(Sesso == "F") %>%
  summarise(
    Media = mean(Peso, na.rm = TRUE),
    Mediana = median(Peso, na.rm = TRUE),
    Deviazione_Standard = sd(Peso, na.rm = TRUE),
    Minimo = min(Peso, na.rm = TRUE),
    Massimo = max(Peso, na.rm = TRUE),
    Quartile1 = quantile(Peso, 0.25, na.rm = TRUE),
    Quartile3 = quantile(Peso, 0.75, na.rm = TRUE)
  )

gender_weight_table = rbind(
  Maschi = male_stats,
  Femmine = female_stats
)

knitr::kable(gender_weight_table)
Media Mediana Deviazione_Standard Minimo Massimo Quartile1 Quartile3
Maschi 3408.215 3430 493.8043 980 4810 3150 3720
Femmine 3161.132 3160 526.3091 830 4930 2900 3470 
ggplot(neonatalcare_data, aes(x = Sesso, y = Peso, fill = Sesso))+
  geom_boxplot(outlier.colour = "lightcoral", outlier.size = 2, alpha = 0.7)+
  labs(
    title = "Distribuzione del peso alla nascita",
    subtitle = "Per genere sessuale (F/M)",
    x = "Genere Sessuale",
    y = "Peso (g)"
  )+
  scale_fill_manual(values = c("M" = "powderblue", "F" = "orchid"))+
  theme_minimal()+
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "none"
  )

Sia i boxplot che il calcolo degli indici di posizione e variabilità per la distribuzione peso divisa per genere sessuale ci restituiscono forti indizi che, alla nascita, il peso dei neonati maschi sia significativamente differente rispetto a quello dei neonati femmina.

Per accertare che la differenza sia statisticamente rilevante, ipotizziamo un “H0 = le medie dei pesi tra neonati maschi e femmine sono uguali” ed effettuiamo un test t di student di seguito.

t.test(Peso ~ Sesso, data = neonatalcare_data)
## 
##  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

Il test restituisce un p-value praticamente pari allo 0, indicando che la differenza media di peso fra neonati maschi e femmine è statisticamente significativa, e quindi ci spinge a rigettare la nostra ipotesi H0 sopra citata.

Ci dice inoltre che, nel dataset a nostra disposizione, il peso medio dei neonati femmina mostra una differenza che va da 207.06 a 287.11 grammi alla nascita rispetto ai maschi.

Relazione fra peso e ospedali:

Prima di procedere a prendere in considerazione le variabili strettamente materne in relazione al peso, vale anche la pena verificare se ci siano forti differenze fra il peso dei neonati nei tre ospedali separati.

ggplot(neonatalcare_data, aes(x = Ospedale, y = Peso, fill = Ospedale))+
  geom_boxplot(outlier.colour = "lightcoral", outlier.size = 2, alpha = 0.7)+
  labs(
    title = "Distribuzione del peso dei neonati",
    subtitle = "Diviso per ospedali",
    x = "Ospedale",
    y = "Peso (g)"
  )+
  scale_fill_brewer(palette = "Pastel1")+
  theme_minimal()+
  theme(
    plot.title = element_text(hjust=0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
  )

All’esame grafico non appaiono significative differenze nella distribuzione del peso divisa per i tre ospedali nel dataset.

Relazione fra peso dei neonati e tipologia di parto:

ggplot(neonatalcare_data, aes(x = Tipo.parto, y = Peso, fill = Tipo.parto))+
  geom_boxplot(outlier.colour = "lightcoral", outlier.size = 2, alpha = 0.7)+
  labs(
    title = "Distribuzione del peso dei neonati per tipo di parto",
    subtitle = "Parto naturale vs cesareo",
    x = "Tipo di parto",
    y = "Peso (g)"
  )+
  scale_fill_manual(values = c("Nat" = "forestgreen", "Ces" = "firebrick"))+
  theme_minimal()+
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "none"
  )

Osserviamo graficamente come non vi sia alcuna influenza sugli indici della distribuzione del peso dei neonati nemmeno discriminando per tipologia di parto.

Relazione fra peso dei neonati, tipologia di parto e ospedale:

Possiamo dettagliare maggiormente l’analisi della relazione fra le variabili verificando anche se, fra i vari ospedali, vi siano differenze fra parto naturale e cesareo limitatamente alla struttura presa in considerazione.

ggplot(data = neonatalcare_data, aes(x = Tipo.parto, y = Peso, fill = Ospedale)) +
  geom_boxplot(outlier.colour = "lightcoral", outlier.size = 2, alpha = 0.7) +
  labs(title = "Peso dei Neonati per Tipo di Parto e Ospedale",
       x = "Tipo di Parto",
       y = "Peso (g)") +
  scale_fill_brewer(palette = "Pastel1")+
  theme_minimal()

Ad un primo occhio, sembrano esserci delle fluttuazioni fra le distribuzioni dei parti di tipo cesareo se divise per ospedale, per cui procederemo ad effettuare un ANOVA test a due vie con le seguenti ipotesi H0:

  • Per il “Tipo di Parto”: Non ci sono differenze significative nei pesi medi dei neonati tra cesareo e naturale;

  • Per “Ospedale”: Non ci sono differenze significative nei pesi medi dei neonati per i tre ospedali;

  • Per Interazione: Non c’è significativa interazione tra il tipo di parto e l’ospedale;

Questo test ci aiuta a verificare statisticamente anche quanto visto nei due grafici precedenti.

anov_result = aov(Peso ~ Tipo.parto * Ospedale, data = neonatalcare_data)
summary(anov_result)
##                       Df    Sum Sq Mean Sq F value Pr(>F)
## Tipo.parto             1      4250    4250   0.015  0.901
## Ospedale               2    934018  467009   1.693  0.184
## Tipo.parto:Ospedale    2     82304   41152   0.149  0.861
## Residuals           2494 687867970  275809

Come si nota dai p-value siamo costretti a rifiutare l’ipotesi H1 del test, in quanto per tutte e tre le ipotesi il test suggerisce che non vi siano differenze statisticamente significative tra le medie nel gruppo e tra i gruppi.

Relazione fra anni della madre e peso del neonato:

ggplot(neonatalcare_data, aes(x = Anni.madre, y = Peso))+
  geom_point(alpha = 0.6, color = "dodgerblue")+
  geom_smooth(method = "lm", color = "darkorange", se = TRUE, linetype = "dashed")+
  labs(
    title = "Relazione tra età della madre e peso del neonato",
    x = "Età della madre (in Anni)",
    y = "Peso del neonato (in grammi)",
    subtitle = "Con linea di regressione lineare e intervallo di confidenza"
  )+
  theme_minimal()+
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

Mettendo graficamente in relazione l’età della madre con il peso del neonato alla nascita, notiamo una certa concentrazione nei valori fra 2500 e 4000 come già verificato precedentemente con la tabella di frequenze per la variabile peso. Si notano anche i valori outlier per le madri con 0 e 1 anni di età segnalati nella sezione apposita.

Non sembra esserci alcuna relazione significativa fra le due variabili, e la linea di regressione quasi orizzontale ne è la conferma, pur considerato l’intervallo di confidenza evidenziato.

Relazione fra numero di gravidanze e peso del neonato:

ggplot(neonatalcare_data, aes(x = factor(N.gravidanze), y = Peso, fill = factor(N.gravidanze))) +
  geom_boxplot(outlier.colour = "lightcoral", outlier.size = 2, alpha = 0.7) +
  labs(
    title = "Distribuzione del peso in relazione al numero di gravidanze",
    x = "Numero di gravidanze",
    y = "Peso del neonato (g)"
  ) +
  scale_fill_viridis_d() +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "none"
  )

Utilizziamo il boxplot sovrastante per mettere in relazione il peso del neonato alla nascita con il numero totale di gravidanze avuto dalla madre nella sua storia clinica. Notiamo immediatamente un grande quantitativo di valori outlier fino alle 3 gravidanze, probabilmente assimilabili a gravidanze difficili o anomale.

Notiamo inoltre come non vi sia una sostanziale differenza negli indici di peso fino ai soggetti con 6 gravidanze totali, dove si spostano il min/max e la differenza interquartile è molto allungata rispetto alle altre condizioni del numero di gravidanze. Lo stesso vale per le persone con 8, 9 o 10 gravidanze totali, che mostrano degli indici molto differenti rispetto alle condizioni iniziali 0, 1, 2, 3.

Ipotizziamo tuttavia che, essendo le gravidanze superiori a 2 progressivamente meno comuni fino a 12, non abbiamo a disposizione abbastanza dati in quest’area del dataset per poter verificare con certezza la relazione peso/numero di gravidanze sopra le 5 gravidanze totali.

Possiamo dividere il boxplot ulteriormente in maschi e femmine per vedere se la tendenza della relazione rilevata fra le variabili è la stessa discriminando per genere sessuale.

ggplot(neonatalcare_data, aes(x = factor(N.gravidanze), y = Peso, fill = factor(N.gravidanze))) +
  geom_boxplot(outlier.colour = "lightcoral", outlier.size = 2, alpha = 0.7) +
  labs(
    title = "Distribuzione del peso in relazione al numero di gravidanze",
    subtitle = "Diviso per F/M",
    x = "Numero di gravidanze",
    y = "Peso del neonato (g)"
  ) +
  facet_wrap(~ Sesso) +
  scale_fill_viridis_d() +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "none"
  )

Quando dividiamo il boxplot per femmine e maschi, notiamo come anche le nascite da madri con 5 gravidanze totali sembrano mostrare graficamente differenze significative negli indici, inoltre nei maschi perde significatività il boxplot per le madri con 9 gravidanze totali, avvalorando nuovamente la tesi che al di sopra delle cinque gravidanze potremmo non avere sufficienti dati fra le osservazioni.

Notiamo inoltre come per le madri con 6 gravidanze totali, nei maschi si ottenga una linea mediana significativamente più alta delle altre, mentre nelle femmine la mediana si attesta intorno alla linea dei 2500grammi che potremmo considerare problematica perchè prima indicatrice di un neonato prematuro.

Si consiglia in questo caso di valutare un estensione del campione, in futuri rilevamenti, per accertarsi che le tendenze indicate nella relazione fra peso e numero di gravidanze nelle madri con più di 5 figli rispecchino queste analisi.

Relazione fra tabagismo e peso del neonato:

ggplot(neonatalcare_data, aes(x= as.factor(Fumatrici), y = Peso, fill = as.factor(Fumatrici)))+
  geom_boxplot(outlier.colour = "lightcoral", outlier.size = 2, alpha = 0.7)+
  labs(
    title = "Distribuzione del peso alla nascita",
    subtitle = "Confronto per madri non fumatrici e fumatrici",
    x = "Tabagismo (0: Non fumatrice, 1: Fumatrice)",
    y = "Peso del neonato alla nascita (g)"
  )+
  scale_fill_manual(values = c("0"="tomato", "1" = "palegreen"),
                    labels = c("Non fumatrici", "Fumatrici"))+
  theme_minimal()+
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "none"
  )

Creiamo un boxplot per confrontare la relazione fra l’eventuale tabagismo delle madri e il peso del neonato. Ad una prima occhiata, sembra esserci una relazione fra il fumo e il peso, in quanto notiamo che la mediana della distribuzione per le madri fumatrici è leggermente inferiore rispetto a quella delle madri non fumatrici.

Va anche segnalato che il campione prevede solo un 4.16% di madri fumatrici rispetto al totale osservato, per cui i dati potrebbero non essere abbastanza da evidenziare una tendenza netta.
Decidiamo comunque di effettuare un t-test con ipotesi H0: “Non esiste differenza tra le medie di peso dei due gruppi”, per verificare che le variazioni osservate nel grafico siano statisticamente rilevanti.

t.test(Peso ~ Fumatrici, data = neonatalcare_data)
## 
##  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

Avendo un p-value superiore a 0.05, e una differenza fra le medie di soli 50 grammi, riteniamo che la differenza di peso fra neonati con madri fumatrici e non fumatrici, in media, non sia statisticamente rilevante, o comunque non in grado di avere una relazione forte con il peso del neonato alla nascita. Necessitiamo quindi di un’analisi più approfondita in merito alle influenze subite dal Peso del neonato alla nascita. Procediamo.

Relazione fra settimane di gestazione, lunghezza, cranio e peso:

Avendo indagato la relazione fra peso e genere, possiamo procedere a indagare anche la relazione fra settimane di gestazione, lunghezza, cranio e peso, nella quale ipotiziamo una quasi ovvia correlazione lineare. Visualizziamo i grafici di seguito

ggplot(neonatalcare_data, aes(x= Gestazione, y = Peso))+
  geom_point(alpha = 0.6, color = "steelblue")+
  geom_smooth(method = "lm", se = TRUE, color = "darkorange", fill = "sandybrown")+
  labs(
    title = "Relazione tra settimane di gestazione e peso del neonato",
    x = "Settimane di gestazione",
    y = "Peso (g)"
  )+
  theme_minimal()+
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

ggplot(neonatalcare_data, aes(x= Lunghezza, y = Peso))+
  geom_point(alpha = 0.6, color = "darkred")+
  geom_smooth(method = "lm", se = TRUE, color = "darkorange", fill = "sandybrown")+
  labs(
    title = "Relazione tra Lunghezza e peso del neonato",
    x = "Lunghezza",
    y = "Peso (g)"
  )+
  theme_minimal()+
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

ggplot(neonatalcare_data, aes(x= Cranio, y = Peso))+
  geom_point(alpha = 0.6, color = "coral4")+
  geom_smooth(method = "lm", se = TRUE, color = "darkorange", fill = "sandybrown")+
  labs(
    title = "Relazione tra Cranio e peso del neonato",
    x = "Cranio",
    y = "Peso (g)"
  )+
  theme_minimal()+
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

Individuiamo una forte correlazione lineare positiva fra peso del neonato e settimane di gestazione. Lo stesso vale anche per peso e cranio, lunghezza, per cui all’aumentare delle tre variabili indicate in media aumenta anche il peso del neonato.

Costruzione di un modello di regressione lineare multipla e selezione del modello ottimale.

Avendo appurato che le variabili del dataset non hanno tutte una forte relazione lineare con il Peso del neonato, e che i dati potrebbero influire differentemente su di esso, costruiamo un modello di regressione lineare multipla per comprendere come tutte le variabili agiscono sul peso del neonato, e quali, nello specifico, non abbiano rilevanza per la nostra indagine statistica.

Innanzitutto costruiamo una matrice di correlazione per condurre un’analisi preliminare sulla correlazione fra le variabili e avere un riepilogo totale.

library(dplyr)

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))
  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)
}

data_select = neonatalcare_data %>%
  select(Peso, Anni.madre, N.gravidanze, Gestazione, Lunghezza, Cranio)

pairs(data_select, lower.panel=panel.cor, upper.panel = panel.smooth)

Dalla matrice di correlazione, notiamo immediatamente che Gestazione e Peso hanno una correlazione di 0.6, proprio come precedentemente verificato, e quindi la variabile gestazione ha un’influenza sul peso del neonato che vale la pena esplorare.

All’aumentare del peso, col passare delle settimane di gestazione, ipotizziamo che anche le dimensioni in termini di lunghezza del feto/neonato aumentino, motivo per cui osserviamo una correlazione fra Gestazione, Lunghezza e diametro del cranio che poi sfociano rispettivamente in uno 0.8 e 0.7 di correlazione con Peso.

Gli scatterplot ci mostrano che tutte e tre le variabili influenzano positivamente Peso, per cui sarà necessario tenerle in considerazione nel nostro modello. Non osserviamo significativa correlazione delle altre variabile su Peso.

Teniamo conto di queste considerazioni, e di quelle precedenti nell’analisi delle relazioni, per costruire il nostro modello.

Costruzione del modello test con tutte le variabili:

Di seguito, sviluppiamo la prima versione del nostro modello di regressione lineare multipla includendo tutte le variabili quantitative del dataset, ed osservando i risultati per poter poi andare a migliorarlo ulteriormente. Useremo la procedura stepwise, facendo un’analisi generale e poi eliminando una alla volta le variabili per raggiungere, nel complesso, un modello parsimonioso e con la maggior capacità predittiva possibile.

model_1 = lm(Peso ~ Anni.madre + N.gravidanze +  Fumatrici + Gestazione + Lunghezza + Cranio + Sesso, data = neonatalcare_data)
summary(model_1)
## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Sesso, data = neonatalcare_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1161.56  -181.19   -15.75   163.70  2630.75 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6714.4109   141.1515 -47.569  < 2e-16 ***
## Anni.madre       0.9585     1.1347   0.845   0.3984    
## N.gravidanze    11.2756     4.6690   2.415   0.0158 *  
## Fumatrici      -30.2959    27.5971  -1.098   0.2724    
## Gestazione      32.9331     3.8267   8.606  < 2e-16 ***
## Lunghezza       10.2342     0.3009  34.009  < 2e-16 ***
## Cranio          10.5177     0.4268  24.642  < 2e-16 ***
## SessoM          78.0845    11.2039   6.969 4.06e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2492 degrees of freedom
## Multiple R-squared:  0.7272, Adjusted R-squared:  0.7264 
## F-statistic:   949 on 7 and 2492 DF,  p-value: < 2.2e-16

Innanzitutto partiamo dal risultato più importante. Il valore del nostro coefficiente di determinazione aggiustato (adjusted R-squared) è fissato a 0.7264, indicandoci che il modello così costruito è già giudicabile come un buon modello, in quanto spiega da solo il 72.64% della variabilità del peso dei neonati.

Con un p-value così piccolo, quasi pari a zero, sappiamo che è altamente significativo nel suo complesso.

Di tutte le variabili considerate, Anni.madre e Fumatrici non influenzano significativamente la risposta, superando 0.05 con un p-value ben più elevato. Nei prossimi modelli test, quindi, elimineremo una alla volta queste due variabili osservando gli effetti sul modello per comprendere come reagisce.

Il numero di gravidanze, invece, con 0.0158 risulta statisticamente significativo, in quanto un aumento del numero di gravidanze è associato, in media, a un aumento di 11.28g di peso del neonato per ogni gravidanza pregressa.

Eliminiamo dal modello la variabile Anni.madre e osserviamo i risultati.

model_2 = update(model_1,~. -Anni.madre)
summary(model_2)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso, data = neonatalcare_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1150.3  -181.3   -15.7   163.0  2636.3 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.6714   135.7178 -49.232  < 2e-16 ***
## N.gravidanze    12.7185     4.3450   2.927  0.00345 ** 
## Fumatrici      -30.4634    27.5948  -1.104  0.26972    
## Gestazione      32.5914     3.8051   8.565  < 2e-16 ***
## Lunghezza       10.2341     0.3009  34.011  < 2e-16 ***
## Cranio          10.5359     0.4262  24.718  < 2e-16 ***
## SessoM          78.1713    11.2028   6.978 3.83e-12 ***
## ---
## 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

Nel secondo test, la variabile N.gravidanze ha un p-value più basso rispetto al primo, acquisendo maggiore significatività e passando da 11.28g a 12.71g di variazione media rispetto al numero di gravidanze pregresse. Il coefficiente di determinazione aggiustato è leggermente migliorato (0.7265), indicando un modello efficiente e più parsimonioso.

Ora procediamo eliminando anche la variabile Fumatrici.

model_3 = update(model_2, ~. -Fumatrici)
summary(model_3)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = neonatalcare_data)
## 
## 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

Il terzo modello test, elimina la variabile Fumatrici e conferma una percentuale del 72.65% aggiustata nella spiegazione della variabilità del peso del neonato.

N.gravidanze ha subito un leggero aumento del p-value ma, considerando che i residui sono ancora 274.6, che R-squared è invariato e che abbiamo un p-value ottimale, si tratta di una variazione che rende il modello 3 più semplice, cioè con meno variabili, e con prestazioni predittive perfino leggermente migliorate rispetto agli altri due.

Ora procediamo valutando cosa accade se eliminiamo proprio N.gravidanze che, fra le variabili statisticamente rilevanti prese in considerazione, è quella con il p-value più alto di tutte e il minor impatto in termini di influenza sulla risposta.

model_4 = update(model_3, ~. -N.gravidanze)
summary(model_4)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso, 
##     data = neonatalcare_data)
## 
## 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

L’eliminazione di N.gravidanze non ha minato la bontà predittiva del modello sul peso del neonato. Notiamo come il coefficiente di determinazione aggiustato è calato leggermente, ma non in maniera statisticamente significativa. Da 72.65% è sceso a 72.57%, mantenendo la percentuale di variabilità complessiva spiegata dal modello quasi intatta.

Il nostro p-value è ancora praticamente pari a zero e la distribuzione dei residui non è variata in maniera tale da destare preoccupazioni, essendo passata da 274,6 a 275.

La statistica F però, cioè l’indice che ci comunica la qualità complessiva del modello costruito, è aumentata a 1654 portandoci a ritenere che l’attuale versione sia la più semplice, e un buon candidato per un modello finale insieme al modello 3.

Per completezza e correttezza nell’attuazione della procedura stepwise, elencheremo comunque di seguito i tentativi rimanenti in cui si andrà a vedere cosa accade se togliamo una o più variabili statisticamente rilevanti dal modello.

model_5 = update(model_4, ~. -Sesso)
summary(model_5)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio, data = neonatalcare_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1105.77  -183.25   -12.83   166.41  2623.80 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6777.1203   135.6417 -49.963   <2e-16 ***
## Gestazione     31.6901     3.8219   8.292   <2e-16 ***
## Lunghezza      10.4252     0.3020  34.522   <2e-16 ***
## Cranio         10.7892     0.4282  25.194   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 277.7 on 2496 degrees of freedom
## Multiple R-squared:  0.7206, Adjusted R-squared:  0.7203 
## F-statistic:  2146 on 3 and 2496 DF,  p-value: < 2.2e-16
model_6 = update(model_5, ~. -Cranio)
summary(model_6)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza, data = neonatalcare_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1152.8  -199.6   -20.4   192.0  3627.4 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5349.8835   138.0018  -38.77   <2e-16 ***
## Gestazione     45.1267     4.2377   10.65   <2e-16 ***
## Lunghezza      13.8973     0.3009   46.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 310.9 on 2497 degrees of freedom
## Multiple R-squared:  0.6496, Adjusted R-squared:  0.6493 
## F-statistic:  2314 on 2 and 2497 DF,  p-value: < 2.2e-16
model_7 = update(model_6, ~. -Lunghezza)
summary(model_7)
## 
## Call:
## lm(formula = Peso ~ Gestazione, data = neonatalcare_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1603.61  -287.34   -11.07   280.46  1897.75 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3197.250    176.851  -18.08   <2e-16 ***
## Gestazione    166.272      4.532   36.69   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 423.3 on 2498 degrees of freedom
## Multiple R-squared:  0.3502, Adjusted R-squared:  0.3499 
## F-statistic:  1346 on 1 and 2498 DF,  p-value: < 2.2e-16

Nonostante abbiamo appurato in precedenza che le differenze di peso fra maschi e femmine sono statisticamente rilevanti quando paragonate in media fra di loro, esse non influiscono a tal punto sul peso del neonato da rendersi altrettanto significative nel nostre modello. Nel test numero 5, infatti, notiamo come eliminando la variabile Sesso otteniamo ancora un 72,03% di variabilità del modello spiegata da Gestazione, Cranio e Lunghezza.

Si potrebbe quindi intendere questa variante come quella più semplice che conserva ancora sufficiente potere esplicativo senza incidere in maniera troppo significativa sui residui.

Notiamo come iniziando ad eliminare dal modello sia Cranio che Lunghezza si perde molta capacità predittiva a scapito anche dell’errore nei residui che sale sensibilmente. Scendiamo rispettivamente a 64% e 35% di variabilità spiegata.

Considerate tutte le variabili test del nostro modello riteniamo che il numero 3, che include Gestazione, Lunghezza, Cranio, Sesso e N.gravidanze sia il più indicato per spiegare la variabilità di Peso in maniera semplice, parsimoniosa e precisa, date le osservazioni a nostra disposizione. Procederemo a lasciare fuori Anni.madre e Fumatrici dal modello.

Prima di procedere alla verifica dei residui del modello però, effettueremo un test BIC (Bayes) e un test StepAIC (Akaike) per accertarci che il modello 3 sia davvero il migliore, e che non sia meglio scegliere il 4.

BIC(model_1, model_2, model_3, model_4, model_5, model_6, model_7)
##         df      BIC
## model_1  9 35233.81
## model_2  8 35226.70
## model_3  7 35220.10
## model_4  6 35220.54
## model_5  5 35262.11
## model_6  4 35820.73
## model_7  3 37356.84
library(MASS)

n = nrow(neonatalcare_data)
stepwise.mod = MASS::stepAIC(model_1, direction = "both", k=log(n))
## Start:  AIC=28131.29
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     53803 187973654 28124
## - Fumatrici     1     90879 188010731 28125
## - N.gravidanze  1    439797 188359648 28129
## <none>                      187919851 28131
## - Sesso         1   3662833 191582684 28172
## - Gestazione    1   5585184 193505035 28197
## - Cranio        1  45791026 233710877 28669
## - Lunghezza     1  87220887 275140738 29077
## 
## Step:  AIC=28124.18
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     91892 188065546 28118
## <none>                      187973654 28124
## - N.gravidanze  1    646039 188619694 28125
## + Anni.madre    1     53803 187919851 28131
## - Sesso         1   3671289 191644943 28165
## - Gestazione    1   5531705 193505359 28189
## - Cranio        1  46066755 234040410 28664
## - Lunghezza     1  87218857 275192512 29069
## 
## 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
## + 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 = neonatalcare_data)
## 
## 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

Come si evince dai risultati ottenuti con i test BIC e StepAIC, in entrambi il modello con il risultato più basso è il numero tre, comprensivo di N.gravidanze, Gestazione, Sesso, Cranio e Lunghezza. Si tratta esattamente di uno dei due potenziali candidati, per cui è confermato che sia possibile utilizzarlo.

Per ulteriore sicurezza, procediamo anche a verificare che nel modello 3 non vi siano problemi di multicollinearità. Eseguiamo un test VIF.

library(car)
vif(model_3)
## N.gravidanze   Gestazione    Lunghezza       Cranio        Sesso 
##     1.023475     1.669189     2.074689     1.624465     1.040054

I risultati sono tutti al di sotto della soglia di 5, per cui non siamo in una situazione di multicollinearità. Questo significa che le variabili indipendenti non contengono informazioni simili perchè non hanno sufficiente correlazione fra loro da diventare problematiche per il modello.

Procediamo con l’analisi dei residui.

Analisi della qualità del modello scelto.

La diagnostica dei residui di un modello deve rispettare le assunzioni di Linearità, normalità, omoschedasticità ed indipendenza.

Andiamo ad effettuare la verifica di queste assunzioni. Innanzitutto creiamo un grafico di analisi.

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

  • Residuals Vs Fitted: I residui nel grafico sono distribuiti casualmente intorno allo zero e non si osservano aperture a ventaglio o pattern anomali o sistematici, per cui dovremmo rispettare le assunzioni di linearità e omoschedasticità dei residui;

  • Q-Q Residuals: I residui dovrebbero essere allineati lungo la diagonale tratteggiata per rispettare la condizione di normalità della distribuzione dei residui. Le code tendono a discostarsi dalla linea per cui effettueremo uno shapiro test per sicurezza;

  • Scale-Location: I punti sono ancora una volta distribuiti casualmente senza aperture e pattern sistematici, per cui rispettiamo la condizione di varianza costante (omoschedasticità) dei residui

  • Residuals vs Leverage: Come osservato nel resto dei grafici, abbiamo potenziali valori outlier, tuttavia di questi solo 1551 sembra poter essere problematico perchè supera la prima soglia di allarme 0.5 e sembra andare molto vicino a 1.

Effettuiamo lo shapiro test sui residui per verificare la normalità della distribuzione.

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

Con un p-value praticamente uguale a zero e un W di 0.97, lo shapiro test ci suggerisce di rifiutare H0 e quindi di non considerare la distribuzione dei residui come una normale.

Tuttavia effettuiamo anche una verifica grafica per essere certi che non sia necessario agire diversamente.

residuals_3 = rstandard(model_3)

plot(
  density(residuals_3),
  main = "Curva di densità dei residui del modello 3 standardizzati",
  xlab = "Residui standardizzati",
  ylab = "Densità",
  col = "olivedrab4",
  lwd = 2
)

curve(
  dnorm(x, mean = mean(residuals_3), sd = sd(residuals_3)),
  col = "lightcoral",
  lwd = 2,
  add = TRUE
  )

legend("topright", 
       inset = c(0, 0),
       legend = c("Curva di densità residui", "Distribuzione normale"),
       col = c("olivedrab4", "lightcoral"),
       lwd = 2,
       bty = "n"
       )

In rosso visualizziamo una distribuzione normale mentre in verde la curva di densità dei residui per il modello 3. Come evidente, la curva dei residui risulta leptocurtica tuttavia non sembra discostarsi in maniera significativa da una normale, per cui procederemo senza effettuare modifiche per ora.

Per verificare l’assunzione di omoschedasticità dei residui, che visivamente non sembrava violata, effettuiamo un test Breusch-Pagan.

lmtest::bptest(model_3)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_3
## BP = 90.253, df = 5, p-value < 2.2e-16

Con un p-value così basso e un BP di 90.25, il test ci spinge a rifiutare l’ipotesi H0 di omoschedasticità e suggerisce eteroschedasticità.

E’ possibile che a influenzare il BP siano singole variabili predittive con pattern specifici di distribuzione dei residui. In base alle analisi precedenti, ipotizziamo che possano essere le variabili Gestazione e N.gravidanze.

Verifichiamo di seguito.

plot(neonatalcare_data$Gestazione, resid(model_3),
     main = "Residuals vs Gestazione",
     xlab = "Gestazione", ylab = "Residuals")
abline(h = 0, col = "red")

plot(neonatalcare_data$N.gravidanze, resid(model_3),
     main = "Residuals vs N.gravidanze",
     xlab = "Gestazione", ylab = "Residuals")
abline(h = 0, col = "red")

Come notiamo nel grafico, le due variabili hanno residui che assumono pattern specifici non casuali intorno allo zero, suggerendo una varianza non costante per via della mancanza di un pattern casuale della distribuzione dei residui come richiesto dal test, inoltre tendono ad aumentare o diminuire con un’apertura a ventaglio. Possiamo ipotizzare quindi che Gestazione e N.gravidanze stiano influenzando i risultati del BP test suggerendo eteroschedasticità.

Effettuiamo ora verifiche sui valori outlier notati in fase di analisi preliminare dei residui. Partiamo con l’analisi dei valori outlier.

plot(rstudent(model_3))
abline(h=c(-2,2), col=2)

car::outlierTest(model_3)
##       rstudent unadjusted p-value Bonferroni p
## 1551 10.051908         2.4906e-23   6.2265e-20
## 155   5.027798         5.3138e-07   1.3285e-03
## 1306  4.827238         1.4681e-06   3.6702e-03

Come si evince dal nostro OutlierTest, abbiamo 1551, 155 e 1306 che sono considerabili come valori outlier. Di questi tuttavia soltanto 1551 risulta estremamente significativo e potenzialmente molto influente. mentre 155 e 1306 sono meno estremi ma comunque rilevanti.

Il grafico mette bene in evidenza 1551, tuttavia non vi sono pattern specifici, la casualità intorno allo zero è rispettata e non osserviamo pattern sistematici nemmeno nella disposizione di cluster di valori outlier fuori dalle soglie, non formano una sottopopolazione, per cui non riteniamo ci siano problemi di predizione nel modello.

Guardiamo ora con maggior precisione alla distanza di Cook.

plot(model_3, which = 4)
points(1551, cooks.distance(model_3)[1551], col = "red", pch = 19)

cook = cooks.distance(model_3)
cookframe = data.frame(Cooks_max_distance = max(cook))
knitr::kable((cookframe), "pipe")
Cooks_max_distance
0.8300965

Grazie a questo secondo grafico confermiamo che 1551 è l’unico valore outlier che supera la soglia di allarme di 0.5, e arriva a un valore di distanza di Cook di 0.83, che lo rende molto vicino alla soglia critica di 1.

Effettuiamo un’analisi visiva delle osservazioni degli outlier per capire se ci possono essere problemi di inserimento nel dataset.

my_outliers = neonatalcare_data[c(1551,155,1306),]
knitr::kable((my_outliers), "pipe")
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
1551 35 1 0 38 4370 315 374 Nat osp3 F
155 30 0 0 36 3610 410 330 Nat osp1 M
1306 23 0 0 41 4900 510 352 Nat osp2 F

Rispetto alle medie tipiche del dataset, 1551 sembra avere una lunghezza del neonato troppo bassa per un peso di 4370. Anche 155 presenta una lunghezza molto al di sotto della media della variabile precedentemente osservata, per cui anomala, per un peso tra l’altro al di sopra della media della variabile. Tutto questo per una gestazione una settimana al di sotto della soglia di prenatalità normalmente considerata.

E’ plausibile pensare che ci troviamo davanti a rilevazioni errate per cui procederemo a tenerle nel dataset perchè frutto di un caso isolato e non in grado di invalidare il lavoro svolto sino ad ora.

Tuttavia questa è una ipotesi partita da una nostra opinione, ma si suggerisce un consulto con il personale medico per verificare effettivamente che casi con questi valori statistici non siano possibili in realtà.

A questo punto ci resta da valutare soltanto l’autocorrelazione dei residui delle variabili del modello per confermare l’assunzione di indipendenza.

lmtest::dwtest(model_3)
## 
##  Durbin-Watson test
## 
## data:  model_3
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0

Il DW test suggerisce una leggera autocorrelazione positiva fra le variabili, ma essendo molto vicino a due non ci preoccupiamo perchè siamo praticamente in procinto di un valore ottimale. Il p-value di 0.12 inoltre ci spinge a non rifiutare l’ipotesi nulla H0 per i residui indipendenti. Non ci sono prove sufficienti per confutare l’assunto e invalidare il modello.

Utilizzo del modello creato per previsioni statistiche.

A questo punto non ci resta che utilizzare il modello 3 per effettuare delle previsioni statistiche sul peso medio dei neonati alla nascita, date come condizione le altre variabili.

Procediamo innanzitutto a verificare quale sarebbe il peso medio alla nascita di una neonata con:

  • Madre non fumatrice;

  • Tre gravidanze;

  • Trentanovesima settimana di gestazione;

forecast = predict(
  model_3, 
  newdata = data.frame(
    Sesso = "F", 
    N.gravidanze = 3, 
    Gestazione = 39, 
    Fumatrici = 0, 
    Lunghezza = round(mean(Lunghezza), 2), 
    Cranio = round(mean(Cranio), 2)
    )
  )

forecast_data = data.frame (
  Previsione = forecast,
  Sesso = "F",
  N.gravidanze = 3,
  Gestazione = 39,
  Fumatrici = 0,
  Lunghezza_media = round(mean(Lunghezza), 2),
  Cranio_medio = round(mean(Cranio), 2)
)

knitr::kable(forecast_data)
Previsione Sesso N.gravidanze Gestazione Fumatrici Lunghezza_media Cranio_medio
3271.078 F 3 39 0 494.69 340.03

Il nostro modello stima che, con le condizioni sopra elencate, la neonata abbia un peso di 3271.08. Per ottenere questo risultato abbiamo impostato come Lunghezza e Cranio i valori della media della distribuzione arrotondati a due cifre decimali.

Effettuiamo altre due prove con dati diversi per osservare come reagisce il modello in fase predittiva.

In questo caso però, calcoleremo il valore atteso delle variabili Cranio e Lunghezza in base alle settimane di gestazione desiderate, così da non invalidare qualsiasi tentativo predittivo utilizzando sempre i valori medi di cui sopra.

Procederemo inoltre a discriminare il calcolo di lunghezza e cranio medi, per le settimane di gestazione indicate, anche per genere sessuale, avendo appurato che le differenze di peso fra maschi e femmine alla nascita hanno una significatività statistica che non può essere ignorata, e questo può riflettersi di rimando sulle dimensioni medie del corpo del neonato nel momento della nascita.

mean_gestazione = mean(Gestazione, na.rm=T)
mean_lunghezza_M = mean(Lunghezza[Sesso == "M"], na.rm = TRUE)
mean_cranio_M = mean(Cranio[Sesso == "M"], na.rm = TRUE)

mean_lunghezza_F = mean(Lunghezza[Sesso == "F"], na.rm = TRUE)
mean_cranio_F = mean(Cranio[Sesso == "F"], na.rm = TRUE)

n_gestazione = 35
lunghezza_forecast_M = round((n_gestazione / mean_gestazione) * mean_lunghezza_M, 2)
cranio_forecast_M = round((n_gestazione / mean_gestazione) * mean_cranio_M, 2)

lunghezza_forecast_F = round((n_gestazione / mean_gestazione) * mean_lunghezza_F, 2)
cranio_forecast_F = round((n_gestazione / mean_gestazione) * mean_cranio_F, 2)

forecast_2 = predict(
  model_3, 
  newdata = data.frame(
    Sesso = "M", 
    N.gravidanze = 0, 
    Gestazione = n_gestazione, 
    Fumatrici = 1, 
    Lunghezza = lunghezza_forecast_M, 
    Cranio = cranio_forecast_M
  )
)

forecast_data2 = data.frame (
  Previsione = forecast_2,
  Sesso = "M",
  N.gravidanze = 0,
  Gestazione = n_gestazione,
  Fumatrici = 1,
  Lunghezza_media = lunghezza_forecast_M,
  Cranio_medio = cranio_forecast_M
)

forecast_3 = predict(
  model_3, 
  newdata = data.frame(
    Sesso = "F", 
    N.gravidanze = 0, 
    Gestazione = n_gestazione, 
    Fumatrici = 1, 
    Lunghezza = lunghezza_forecast_F, 
    Cranio = cranio_forecast_F
  )
)

forecast_data3 = data.frame (
  Previsione = forecast_3,
  Sesso = "F",
  N.gravidanze = 0,
  Gestazione = n_gestazione,
  Fumatrici = 1,
  Lunghezza_media = lunghezza_forecast_F,
  Cranio_medio = cranio_forecast_F
)

knitr::kable(forecast_data2)
Previsione Sesso N.gravidanze Gestazione Fumatrici Lunghezza_media Cranio_medio
2367.287 M 0 35 1 448.64 307.48 
knitr::kable(forecast_data3)
Previsione Sesso N.gravidanze Gestazione Fumatrici Lunghezza_media Cranio_medio
2152.651 F 0 35 1 439.75 303.16

Secondo il modello così utilizzato un neonato pretermine di 35 settimane, nato da una madre fumatrice a zero gravidanze pregresse, avrà un peso stimato alla nascita di 2367.287g per i maschi e 2152.651g per le femmine. Ora rieffettuiamo la stessa previsione di prima, discriminando tra maschi e femmine.

n_gestazione = 39
lunghezza_forecast_M = round((n_gestazione / mean_gestazione) * mean_lunghezza_M, 2)
cranio_forecast_M = round((n_gestazione / mean_gestazione) * mean_cranio_M, 2)

lunghezza_forecast_F = round((n_gestazione / mean_gestazione) * mean_lunghezza_F, 2)
cranio_forecast_F = round((n_gestazione / mean_gestazione) * mean_cranio_F, 2)

forecast_4 = predict(
  model_3, 
  newdata = data.frame(
    Sesso = "M", 
    N.gravidanze = 3, 
    Gestazione = n_gestazione, 
    Fumatrici = 0, 
    Lunghezza = lunghezza_forecast_M, 
    Cranio = cranio_forecast_M
  )
)

forecast_data4 = data.frame (
  Previsione = forecast_4,
  Sesso = "M",
  N.gravidanze = 3,
  Gestazione = n_gestazione,
  Fumatrici = 0,
  Lunghezza_media = lunghezza_forecast_M,
  Cranio_medio = cranio_forecast_M
)

forecast_5 = predict(
  model_3, 
  newdata = data.frame(
    Sesso = "F", 
    N.gravidanze = 3, 
    Gestazione = n_gestazione, 
    Fumatrici = 0, 
    Lunghezza = lunghezza_forecast_F, 
    Cranio = cranio_forecast_F
  )
)

forecast_data5 = data.frame (
  Previsione = forecast_5,
  Sesso = "F",
  N.gravidanze = 3,
  Gestazione = n_gestazione,
  Fumatrici = 0,
  Lunghezza_media = lunghezza_forecast_F,
  Cranio_medio = cranio_forecast_F
)

knitr::kable(forecast_data4)
Previsione Sesso N.gravidanze Gestazione Fumatrici Lunghezza_media Cranio_medio
3429.97 M 3 39 0 499.92 342.62 
knitr::kable(forecast_data5)
Previsione Sesso N.gravidanze Gestazione Fumatrici Lunghezza_media Cranio_medio
3199.61 F 3 39 0 490.01 337.8

Notiamo come, discriminando per genere sessuale e inserendo un valore di Lunghezza e Cranio rapportato alla relativa media per le settimane di Gestazione desiderate, la previsione richiesta in origine dall’azienda varia di quasi 100 grammi scendendo a 3199.61. Un maschio neonato con le stesse condizioni delle variabili materne avrebbe un valore predittivo del peso di 3429.97

Facciamo altre 3 previsioni su 41 settimane di Gestazione.

n_gestazione = 41
lunghezza_forecast_M = round((n_gestazione / mean_gestazione) * mean_lunghezza_M, 2)
cranio_forecast_M = round((n_gestazione / mean_gestazione) * mean_cranio_M, 2)

lunghezza_forecast_F = round((n_gestazione / mean_gestazione) * mean_lunghezza_F, 2)
cranio_forecast_F = round((n_gestazione / mean_gestazione) * mean_cranio_F, 2)

forecast_6 = predict(
  model_3, 
  newdata = data.frame(
    Sesso = "M", 
    N.gravidanze = 1, 
    Gestazione = n_gestazione, 
    Fumatrici = 0, 
    Lunghezza = lunghezza_forecast_M, 
    Cranio = cranio_forecast_M
  )
)

forecast_data6 = data.frame (
  Previsione = forecast_6,
  Sesso = "M",
  N.gravidanze = 1,
  Gestazione = n_gestazione,
  Fumatrici = 0,
  Lunghezza_media = lunghezza_forecast_M,
  Cranio_medio = cranio_forecast_M
)

forecast_7 = predict(
  model_3, 
  newdata = data.frame(
    Sesso = "F", 
    N.gravidanze = 1, 
    Gestazione = n_gestazione, 
    Fumatrici = 0, 
    Lunghezza = lunghezza_forecast_M, 
    Cranio = cranio_forecast_M
  )
)

forecast_data7 = data.frame (
  Previsione = forecast_7,
  Sesso = "F",
  N.gravidanze = 2,
  Gestazione = n_gestazione,
  Fumatrici = 0,
  Lunghezza_media = lunghezza_forecast_M,
  Cranio_medio = cranio_forecast_M
)

forecast_8 = predict(
  model_3, 
  newdata = data.frame(
    Sesso = "M", 
    N.gravidanze = 6, 
    Gestazione = n_gestazione, 
    Fumatrici = 1, 
    Lunghezza = lunghezza_forecast_M, 
    Cranio = cranio_forecast_M
  )
)

forecast_data8 = data.frame (
  Previsione = forecast_8,
  Sesso = "M",
  N.gravidanze = 6,
  Gestazione = n_gestazione,
  Fumatrici = 1,
  Lunghezza_media = lunghezza_forecast_M,
  Cranio_medio = cranio_forecast_M
)

knitr::kable(forecast_data6)
Previsione Sesso N.gravidanze Gestazione Fumatrici Lunghezza_media Cranio_medio
3917.648 M 1 41 0 525.56 360.19 
knitr::kable(forecast_data7)
Previsione Sesso N.gravidanze Gestazione Fumatrici Lunghezza_media Cranio_medio
3839.656 F 2 41 0 525.56 360.19 
knitr::kable(forecast_data8)
Previsione Sesso N.gravidanze Gestazione Fumatrici Lunghezza_media Cranio_medio
3980.023 M 6 41 1 525.56 360.19

Effettuato qualche test predittivo col nostro modello, andiamo a verificare se è possibile evidenziare graficamente la sua efficacia predittiva globale, oltre alla relazione parziale tra Peso e predittori inseriti nel modello.

Grafici di efficienza del modello predittivo.

Prima di confrontare la relazione fra i valori osservati nel dataset e quelli predetti dal nostro modello, possiamo avere un’idea più chiara di come le variabili influenzano la capacità predittiva del modello attraverso questi partial dependence plot.

library("effects")
eff = allEffects(model_3)

plot(eff, "N.gravidanze", main = "Effetto di N.gravidanze sul Peso")

plot(eff, "Gestazione", main = "Effetto di Gestazione sul Peso")

plot(eff, "Cranio", main = "Effetto di Cranio sul Peso")

plot(eff, "Lunghezza", main = "Effetto di Lunghezza sul Peso")

plot(eff, "Sesso", main = "Effetto di Sesso sul Peso")

I grafici sopra hanno lo scopo di mostrare come reagisce la variabile dipendente al cambiamento di una singola variabile indipendente, mantenendo fisse le condizioni delle altre variabili del modello. Prendendo in esame il nostro modello 3, questi grafici ci permettono di comprendere meglio dove ci aspettiamo che il range di previsioni del modello diventi più ampio e impreciso, e dove invece mantenga una certa coerenza e possa ritenersi affidabile.

Innanzitutto tutte le variabili prese in esame per questo modello hanno una relazione lineare positiva con la variabile peso, come già appurato in precedenza. Notiamo come, Gestazione, Lunghezza e Cranio hanno un comportamento simile nella loro relazione con Peso. Esse tendono a concentrarsi intorno alla retta quando si raggiungono le osservazioni intorno alla media della variabile indipendente, risultando quindi in osservazioni più precise e attinenti la realtà per il modello.

L’impatto delle variabili indipendenti tende invece a diventare impreciso quando ci si allontana dai valori medi nel modello. Notiamo in particolare questo effetto in Gestazione, dove allontanandoci dalla media di 40 notiamo come il range di valori possibili del modello si apra considerevolmente intorno alla linea, denotando come Gestazione potrebbe impattare sulle previsioni circa il peso.

La variabile con un effetto simile più preoccupante però è N.gravidanze che, come notiamo nel relativo grafico, superate le 6 gestazioni inizia a dare un range di possibili valori veramente ampio per il peso, andando quindi, a parità di condizioni delle altre variabili, ad influenzare il peso in maniera evidente. Questo rispecchia la relazione che abbiamo studiato precedentemente nell’analisi, dove abbiamo ipotizzato che la penuria di dati a disposizione per le gestazioni superiori a cinque rende i risultati della distribuzione imprecisi e meno affidabili.

Il grafico del genere sessuale invece, ci mostra come la differenza fra maschi e femmine abbia un effetto sulla base del peso nel modello in fase di predizione.

Ora proiettiamo in un grafico i nostri valori osservati contro i valori predetti dal nostro modello.

ggplot(data.frame(Observed = neonatalcare_data$Peso, Predicted = fitted(model_3)), 
       aes(x = Observed, y = Predicted)) +
  geom_point(alpha = 0.6) +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
  labs(title = "Model 3: Valori Predetti vs Osservati", x = "Peso Osservato", y = "Peso Predetto") +
  theme_minimal()

La visualizzazione è popolata da punti che rappresentano sia le nostre osservazioni nel dataset che quelle predette dal modello 3, per le condizioni necessarie delle variabili indipendenti intorno a quel valore di peso. Un modello preciso e affidabile mostrerà una tendenza dei punti ad assieparsi casualmente intorno alla linea rossa, senza mostrare tendenze a curve o pattern specifici non casuali. Eventuali comportamenti simili, potrebbero dipendere dalla violazione dell’omoschedasticità come assunzione dei residui oppure dalla presenza di valori leverage/outliers.

In questo caso, notiamo come il modello risulti preciso e affidabile per il range di peso che va da 2500 a 4000g, iniziando ad allontanarsi dalla linea per valori superiori o inferiori a questo range. In particolare notiamo deviazione soprattutto sotto i 1200g e sopra i 4500g.

Qualsiasi predizione che ottenga risultati inferiori o superiori a questi valori potrebbe quindi essere distorta dal modello e sarà da ritenere meno affidabile del range di peso centrale indicato.
Questo effetto potrebbe dipendere, in parte, anche dall’impatto che N.gravidanze e Gestazione hanno sul peso, e che abbiamo osservato nei partial dependence plot precedenti.

Considerazioni finali

L’esplorazione del dataset fornito dalla Neonatal Health Services ha restituito alcuni risultati interessanti, e alcuni spunti di riflessione per l’azienda che suggeriamo di verificare con ulteriori raccolte dati prolungate nel tempo.

Innanzitutto abbiamo appurato che l’incidenza del tabagismo sul peso finale del neonato alla nascita non è statisticamente rilevante, almeno non con la mole di dati fornita per queste due variabili. Temiamo che in questo frangente della rilevazione ci sia una sottostima rispetto alla percentuale reale per cui suggeriamo di continuare la rilevazione, effettuando ulteriori analisi statistiche in futuro con maggiori dati a disposizione.

Rileviamo inoltre come non vi siano differenze nel peso dei neonati fra i tre ospedali nel dataset, e che l’età della madre e la tipologia di parto non hanno alcuna influenza statisticamente rilevante sull’esito finale del peso.
Al contrario, le analisi confermano che fra maschi e femmine vi sarà necessariamente una differenza di peso statisticamente rilevante anche a parità di condizioni delle altre variabili.

Il numero di gravidanze pregresse della madre potrebbe avere un’influenza sul peso, rilevante al punto da non poterla ignorare nel modello costruito, a partire dalle cinque gravidanze in poi, tuttavia, il sospetto è che semplicemente non vi siano abbastanza dati, per questo numero di parti precedenti, da giustificare l’affermazione con un buon grado di certezza. Invitiamo quindi l’azienda a continuare le rilevazioni statistiche, magari ampliando il campione di osservazioni, per poter effettuare ulteriori verifiche future sulle madri con numerosi parti precedenti.

Facciamo notare come, inoltre, i parti gemellari non vengano assolutamente menzionati nel dataset, per cui l’azienda potrebbe non riuscire a prevedere adeguatamente questo tipo di situazioni basandosi sul modello selezionato di cui parleremo a breve. Sospettiamo che la gestazione simultanea di due feti possa avere impatti sul peso medio finale dei neonati, per cui andrebbero inclusi nel dataset in futuro.

Le relazioni significative nel dataset contro la variabile Peso dipendente sono risultate essere Gestazione, Lunghezza e Cranio. Queste ultime due sono anche particolarmente influenti.
Chiaramente osserviamo come maggiori siano le settimane di gestazione e maggiore sia il peso del neonato per effetto naturale della gestazione stessa. Nello stesso tempo aumenteranno anche la Lunghezza e il diametro craniale del neonato, ma senza alcuna correlazione problematica delle variabili indipendenti.

Considerate queste relazioni, abbiamo costruito un modello di previsione statistica dove queste variabili indipendenti convogliano in una regressione lineare multipla. Il modello n°3 da noi scelto per effettuare delle previsioni, riesce a spiegare il 72,65% della variabilità della dipendente Peso, raggiungendo quindi un buon risultato previsionale.

Ulteriori verifiche successive, inoltre, ci dicono come questo modello risulti affidabile, e quindi utilizzabile, per neonati le cui condizioni ci portano a un peso compreso fra 2500g e 4000g alla nascita, mentre la perdità di affidabilità si ottiene oltre i 4000 e al di sotto dei 1200g. Nel caso di utilizzo per casi limite che possano rientrare in questo range di peso, quindi, consigliamo all’azienda di procedere con cautela sul singolo caso ed effettuare ulteriori accertamenti clinici.

Si consiglia inoltre di controllare due volte la correttezza dei dati appena inseriti nel dataset, in modo da evitare errori di battitura che possano comportare problemi in fase di analisi, avendo rilevato alcuni valori problematici ipoteticamente riconducibili a un errore di inserimento di data entry. Si suggerisce valutazione dei casi segnalati insieme al personale medico.