Modello Statistico per la Previsione del Peso Neonatale

Il presente progetto di statistica inferenziale, realizzato da me Alessio Feudo nell’ambito del Master in Data Science erogato da Profession.AI, ha lo scopo di creare un modello statistico in grado di prevedere il peso dei neonati alla nascita. L’analisi condotta potrà apportare numerosi benefici:

Il dataset

Il dataset è composto da informazioni su 2500 neonati provenienti da 3 ospedali differenti. Le variabili presenti per ciascun neonato sono:

  • 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).

Lo scopo del modello predittivo sarà quello di identificare quali tra queste variabili hanno un impatto maggiore sul peso alla nascita, con speciale attenzione sul fumo materno e sulle settimane di gestazione.

Analisi descrittiva preliminare

In questa prima sezione di analisi, si vogliono descrivere anche visivamente gli andamenti e le relazioni tra le variabili presenti. Si presentano innanzitutto gli indici di posizione delle variabili quantitative.

library(knitr)
library(kableExtra)
library(dplyr)
data=read.csv('neonati.csv', sep = ',')
summary_data=data%>%
  select(-Tipo.parto,-Fumatrici,-Ospedale, -Sesso)%>%
  summary()
kable(summary_data, "html", caption="<b style='color:black;'>Sommario dataset grezzo</b>")%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1:6, extra_css = "padding-right: 2px;")
Sommario dataset grezzo
Anni.madre N.gravidanze Gestazione Peso Lunghezza Cranio
Min. : 0.00 Min. : 0.0000 Min. :25.00 Min. : 830 Min. :310.0 Min. :235
1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:38.00 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330
Median :28.00 Median : 1.0000 Median :39.00 Median :3300 Median :500.0 Median :340
Mean :28.16 Mean : 0.9812 Mean :38.98 Mean :3284 Mean :494.7 Mean :340
3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:40.00 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
Max. :46.00 Max. :12.0000 Max. :43.00 Max. :4930 Max. :565.0 Max. :390 
kable(head(data), "html", caption="<b style='color:black;'>Prime righe dataset</b>")%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1:10, extra_css = "padding-right: 5px;")
Prime righe dataset
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
26 0 0 42 3380 490 325 Nat osp3 M
21 2 0 39 3150 490 345 Nat osp1 F
34 3 0 38 3640 500 375 Nat osp2 M
28 1 0 41 3690 515 365 Nat osp2 M
20 0 0 38 3700 480 335 Nat osp3 F
32 0 0 40 3200 495 340 Nat osp2 F 
kable(tail(data), "html", caption="<b style='color:black;'>Ultime righe dataset</b>")%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1:10, extra_css = "padding-right: 5px;")
Ultime righe dataset
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
2495 32 1 0 38 3750 510 360 Ces osp1 M
2496 25 2 0 37 3250 490 350 Nat osp2 F
2497 26 1 0 40 3140 500 336 Nat osp2 F
2498 21 0 0 37 2760 480 330 Ces osp3 F
2499 31 1 0 39 3000 485 320 Nat osp2 F
2500 29 1 0 40 3410 510 340 Nat osp2 M 
#sostituisco i dati errati con la loro media
media_corretta=mean(data$Anni.madre[data$Anni.madre >= 12])
data$Anni.madre[data$Anni.madre<12]=media_corretta
summary_data_correct=data%>%
  select(-Tipo.parto,-Fumatrici,-Ospedale, -Sesso)%>%
  summary()
kable(summary_data_correct, "html", caption="<b style='color:black;'>Sommario dataset con dati corretti</b>")%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1:6, extra_css = "padding-right: 2px;")
Sommario dataset con dati corretti
Anni.madre N.gravidanze Gestazione Peso Lunghezza Cranio
Min. :13.00 Min. : 0.0000 Min. :25.00 Min. : 830 Min. :310.0 Min. :235
1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:38.00 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330
Median :28.00 Median : 1.0000 Median :39.00 Median :3300 Median :500.0 Median :340
Mean :28.19 Mean : 0.9812 Mean :38.98 Mean :3284 Mean :494.7 Mean :340
3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:40.00 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
Max. :46.00 Max. :12.0000 Max. :43.00 Max. :4930 Max. :565.0 Max. :390 
#conto bimbi per ogni ospedale
data_hosp=data %>%
  group_by(Ospedale) %>%
  summarise(M=sum(Sesso=='M'),
            F=sum(Sesso=='F'),
            tot=M+F,
            med_peso=round(mean(Peso),2),
            med_Anni_madre=round(mean(Anni.madre),2),
            sd_Anni_madre=round(sd(Anni.madre),2),
            med_gravidanze=round(mean(N.gravidanze),2),
            med_Lung=round(mean(Lunghezza),2),
            sd_Lung=round(sd(Lunghezza),2),
            med_Cranio=round(mean(Cranio),2),
            sd_Cranio=round(sd(Cranio),2),
            perc_parto_Nat=round(sum(Tipo.parto=='Nat')/tot,2))

total_M=sum(data_hosp$M)
total_F=sum(data_hosp$F)
total_tot=sum(data_hosp$tot)


data_hosp=data_hosp %>%
  add_row(Ospedale = "<b>Globale</b>", M=total_M, F=total_F, tot=total_tot, med_peso=round(mean(data$Peso),2), med_Anni_madre=round(mean(data$Anni.madre),2), sd_Anni_madre=round(sd(data$Anni.madre),2), med_gravidanze=round(mean(data$N.gravidanze),2), med_Lung=round(mean(data$Lunghezza),2), sd_Lung=round(sd(data$Lunghezza),2), med_Cranio=round(mean(data$Cranio),2), sd_Cranio=round(sd(data$Cranio),2), perc_parto_Nat=round(sum(data$Tipo.parto=='Nat')/nrow(data),2))

kable(data_hosp, "html", escape = FALSE, caption="<b style='color:black;'>Nascite per ogni ospedale</b>")%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1:13, extra_css = "padding-right: 2px;")
Nascite per ogni ospedale
Ospedale M F tot med_peso med_Anni_madre sd_Anni_madre med_gravidanze med_Lung sd_Lung med_Cranio sd_Cranio perc_parto_Nat
osp1 407 409 816 3270.27 28.05 5.27 1.00 494.12 27.22 339.93 17.11 0.70
osp2 414 435 849 3270.48 28.10 5.26 0.92 495.33 25.65 339.90 16.16 0.70
osp3 423 412 835 3311.41 28.41 5.12 1.03 494.60 26.11 340.26 16.03 0.72
Globale 1244 1256 2500 3284.08 28.19 5.22 0.98 494.69 26.32 340.03 16.43 0.71

Si può notare dalle tabelle come il peso medio dei neonati alla nascita si attesti intorno a 3284 grammi mentre la lunghezza in millimetri sia di circa 495 mm. Si nota anche come i tre ospedali presentino valori medi delle variabili molto vicini tra loro, a prova di un adeguato numero di unità statistiche presenti in ognuno. Si nota anche che neonati maschi e femmine siano equamente rappresentati in ciascun ospedale e anche nell’intero dataset. Nel dataset si rileva anche la presenza di due valori di età delle madri forse registrati male. Si sostituiscono con la media dell’età materna stimata dai dati del dataset aventi età maggiore o uguale a 12, che è un’età biologicamente ragionevole. Si trova quindi un’età materna media intorno ai 28 anni e un numero medio di gravidanze già avute per madre pari a circa 1. Dal dataset emerge anche che circa il 70% delle gravidanze termina con un parto naturale. Si studiano adesso le distribuzioni delle varie variabili.

library(ggplot2)
library(moments)
col_cont=c('Anni.madre', 'Peso', 'Lunghezza', 'Cranio')
col_discr=c('N.gravidanze', 'Fumatrici', 'Gestazione')

form_index_df=data.frame(
  Variabile = character(),
  Asimmetria = numeric(),
  Curtosi = numeric(),
  stringsAsFactors = FALSE
)

for(var_name in col_cont){
  skew=0
  kurt=0
  var=data[[var_name]]
  if(var_name=='Anni.madre'){x_name='Anni della madre'
                             tit='Distribuzione Anni della madre'}
  else if(var_name=='Peso'){x_name='Peso neonato in g'
                            tit='Distribuzione peso del neonato'}
  else if(var_name=='Lunghezza'){x_name='Lunghezza neonato in mm'
                                 tit="Distribuzione della lunghezza del neonato"}
  else {x_name='Diametro cranio in mm'
        tit='Distribuzione diametro del cranio'}
  
  skew=skewness(var)
  kurt=kurtosis(var)
  
  form_index_df=form_index_df %>%
    add_row(Variabile=var_name, Asimmetria=round(skew,2), Curtosi=round((kurt-3),2))
  
  
  p=ggplot(data=data,aes(x=var))+
    geom_density(fill = "lightblue", alpha = 0.5)+
    labs(x = x_name, y = "Densità", title = paste(tit)) +
    theme_minimal()
  print(p)
}

kable(form_index_df, "html", escape = FALSE, caption="<b style='color:black;'>Indici di forma variabili continue</b>")%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1:3, extra_css = "padding-right: 2px;")
Indici di forma variabili continue
Variabile Asimmetria Curtosi
Anni.madre 0.15 -0.10
Peso -0.65 2.03
Lunghezza -1.51 6.49
Cranio -0.79 2.95 
for(var_name in col_discr){
  var=data[[var_name]]
  if(var_name=='N.gravidanze'){x_name='Numero gravidanze della madre'
                             tit='Istogramma numero di gravidanze della madre'}
  else if(var_name=='Fumatrici'){
                            var=factor(var)
                            
                            x_name='Madri fumatrici'
                            tit='Istogramma madri fumatrici'}
  else{x_name='Settimane di gestazione del neonato'
                                 tit='Istogramma settimane di gestazione neonato'}
  h=ggplot(data=data)+
    geom_bar(aes(x=var), stat='count', col='blue', fill='lightblue')+
    labs(x = x_name, y = "Occorrenze", title = paste(tit)) +
    theme_minimal()
  print(h)
}

#si converte la variabile Fumatrici in stringa di caratteri
data_use=as.data.frame(data)
data_use$Fumatrici=as.numeric(as.character(data_use$Fumatrici))
data_use$Fumatrici=factor(data_use$Fumatrici, levels = c(0, 1), labels = c("Non Fumatrici", "Fumatrici"))

ggplot(data=data_use)+
   geom_boxplot(aes(x=Tipo.parto, y=Gestazione, fill=factor(Fumatrici)))+
   labs(x ='TIpo di parto' , y = "Settimane di gestazione", title = 'Settimane di gestazione in funzione del tipo di parto e fumo materno', fill='Fumatrici') +
  scale_fill_manual(values = c("Non Fumatrici" = "lightblue", "Fumatrici" = "lightcoral"))+
  theme_minimal()

anni_madre_cl=cut(data_use$Anni.madre, breaks=c(13,25,35,47), right=F)
data_use=data_use%>%
  mutate(anni_madre_cl)%>%
  group_by(anni_madre_cl)

ggplot(data=data_use)+
  geom_boxplot(aes(x=anni_madre_cl, y=Gestazione, fill=factor(Fumatrici)))+
  labs(x ='Classi di età madri' , y = "Settimane di gestazione", title = "Boxplot settimane di gestazione in funzione dell'età materna", fill='Fumatrici') +
  scale_fill_manual(values = c("Non Fumatrici" = "lightblue", "Fumatrici" = "lightcoral"))+
  theme_minimal()

ggplot(data=data_use)+
  geom_boxplot(aes(x=factor(Fumatrici), y=Peso), fill='skyblue')+
  labs(x ='Madri fumatrici' , y = "Peso neonato alla nascita", title = "Boxplot peso neonato in funzione del fumo") +
  theme_minimal()

ggplot(data=data_use)+
  geom_boxplot(aes(x=factor(Tipo.parto), y=Peso), fill='skyblue')+
  labs(x ='Tipo di parto' , y = "Peso neonato alla nascita", title = "Boxplot peso neonato in funzione del tipo di parto") +
  theme_minimal()

ggplot(data=data_use)+
  geom_boxplot(aes(x=factor(N.gravidanze), y=Peso), fill='skyblue')+
  labs(x ='Numero di gravidanze della madre' , y = "Peso neonato alla nascita", title = "Boxplot peso neonato in funzione del numero di gravidanze") +
  theme_minimal()

ggplot(data=data_use)+
  geom_point(aes(x=Gestazione, y=Peso), col='skyblue')+
  labs(x ='Settimane di gestazione' , y = "Peso neonato alla nascita", title = "Scatterplot peso neonato in funzione delle settimane di gestazione") +
  theme_minimal()

La distribuzione dell’età materna (dopo correzione di valori errati) si presenta abbastanza simmetrica rispetto al picco. Considerando invece le distribuzioni dei valori antropometrici del bimbo (peso, lunghezza, diametro cranio) si nota come questi siano asimmetrici negativi, come da ‘asimmetria’ in tabella. Questo indica come il neonato che non nasca con parametri medi, tenda ad averli minori di questi ultimi e solo raramente superiori. Per quanto riguarda gli istogrammi, si nota nuovamente come il numero di gravidanze già terminate per le madri sia prevalentemente pari a 0 (moda), ma si nota anche la presenza di donne con 12 figli già nati. Si nota anche come l’incidenza del fumo tra le madri del campione sia limitata, con stragrande prevalenza di madri non fumatrici. Considerando la durata della gravidanza, si evince sia dal relativo istogramma che dalla tabella introduttiva, come la media sia di circa 39 settimane mentre la moda sia di 40, anche qui con molti più valori minori rispetto a quelli maggiori della media, a significare come il fenomeno delle gravidanze brevi sia più incidente rispetto a quello delle gravidanze lunghe. Dal boxplot riassuntivo della relazione tra tipo di parto e durata della gravidanza, si nota come cosiderando le non fumatrici (che sono la grande maggioranza del campione) le mediane siano molto simili in valore e in range interquartile. Si nota però una maggiore incidenza di gravidanze brevi su parti naturali, stando al numero maggiore di outliers. Dal boxplot del peso in funzione del numero di gravidanze, si evince come queste grandezze non abbiano una correlazione stretta. Considerando le fumatrici, è noto dalla letteratura (ref: https://mohre.it/fumo-materno-in-gravidanza-danni-al-cervello-dei-bambini/) che l’effetto del fumo in gravidanza dovrebbe diminuire la durata della gestazione, ma dai boxplot questa tendenza non si evince. Questo potrebbe essere dovuto alla modesta numerosità della categoria delle madri fumatrici nel campione. Passando alla relazione tra l’età delle madri e la durata della gestazione, si può notare dal confronto dei range interquartili come vi sia una lieve tendenza ad avere gravidanze più brevi per madri di età più avanzata. Anche qui le fumatrici non mostrano grandi differenze rispetto alle non fumatrici, sempre a causa della loro bassa nuemerosità probabilmente. Un altro effetto del fumo materno noto in letteratura è quello di portare a neonati mediamente con peso inferiore. Dal relativo boxplot si può notare un’indicazione visiva di tale tendenza: primo, terzo quartile e mediana sono più bassi rispetto agli stessi calcolati per le madri non fumatrici. Questa ultima analisi si può ulteriormente indagare tramite un t-test, che verrà eseguito nella sezione succesiva.

Test di ipotesi

Test di ipotesi: peso e fumo

Si vuole innanzitutto testare l’ipotesi che la differenza nelle medie che si può notare a occhio nel grafico del peso in funzione del fumo sia statisticamente rilevante. Per questo si esegue un test t per confrontare le due medie. Essendo il test a campioni indipendenti, ed essendo le due varianze non significativamente diverse come dimostrato dal test F, si applica la varianza ‘pooled’, cioè pesata. Si fissa il livello di errore del test al 5%.

peso_no_fumo=data$Peso[data$Fumatrici=='0']
peso_fumo=data$Peso[data$Fumatrici=='1']

var.test(peso_no_fumo,peso_fumo)
## 
##  F test to compare two variances
## 
## data:  peso_no_fumo and peso_fumo
## F = 1.2112, num df = 2395, denom df = 103, p-value = 0.2063
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.8987458 1.5741174
## sample estimates:
## ratio of variances 
##           1.211245
data_table=data.frame(
  Variabili=c('Peso no fumo', 'Peso con fumo'),
  Media = c(round(mean(peso_no_fumo),2), round(mean(peso_fumo),2)),
  Varianza = c(round(var(peso_no_fumo),2), round(var(peso_fumo),2))
)
kable(data_table, caption = "<b style='color:black;'>Tabella di Media e Varianza del peso condizionato al fumo</b>")%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1:2, extra_css = "padding-right: 2px;")
Tabella di Media e Varianza del peso condizionato al fumo
Variabili Media Varianza
Peso no fumo 3286.15 277673.8
Peso con fumo 3236.35 229246.7 
t.test(peso_no_fumo,peso_fumo,mu=0, alternative='two.sided', conf.level=0.95, pool.sd=T)
## 
##  Welch Two Sample t-test
## 
## data:  peso_no_fumo and peso_fumo
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -45.61354 145.22674
## sample estimates:
## mean of x mean of y 
##  3286.153  3236.346
plot(density(rt(100000, 114.1)), xlab = "t", ylab = "p(t)", main = "Distribuzione t di Student per 114.1 grandi di libertà")
abline(v=c(qt(0.975,114.1),qt(0.025, 114.1)), col=2)
points(x=1.034, y=0, pch=19, col='blue')

Dal t test si traggono le seguenti conclusioni: il valore della statistica test cade nell’intervallo di accettazione, il p-value associato è maggiore della soglia di errore fissata a priori di 0.05 (le due linee rosse indicano i quantili della distribuzione t di Student associati alle due soglie di errore a destra e sinistra, pari a 2.5% ognuna), inoltre l’intervallo di confidenza (-45.6,145.2) contiene lo zero. Tutto porta a non rifiutare l’ipotesi nulla di medie uguali. Questo risultato è presumibilmente dovuto alla piccola numerosità delle fumatrici nel campione.

Test di ipotesi: cesarei e ospedali

Si vuole ora testare l’ipotesi che il numero di parti cesarei sia statisticamente differente in base all’ospedale. Per fare ciò, si usa un test di indipendenza del Chi-quadro di Pearson per saggiare l’ipotesi che il numero di parti cesarei non dipenda dall’ospedale.

tab_contingenza=table(data$Ospedale, data$Tipo.parto)
kable(tab_contingenza, caption = "<b style='color:black;'>Tabella di contingenza tra ospedali e tipi di parto</b>")%>%
   kable_styling(bootstrap_options = c("striped", "hover"))
Tabella di contingenza tra ospedali e tipi di parto
Ces Nat
osp1 242 574
osp2 254 595
osp3 232 603 
chisq.test(tab_contingenza)
## 
##  Pearson's Chi-squared test
## 
## data:  tab_contingenza
## X-squared = 1.0972, df = 2, p-value = 0.5778
plot(density(rchisq(100000, 2)), xlab = "X-squared", ylab = "p(X-squared))", main = "Distribuzione del Chi quadro per 2 grandi di libertà")
abline(v=qt(0.95,2), col=2)
points(x=1.0972, y=0, pch=19, col='blue')

Essendo il valore della statistica X-quadro contenuto nel range di accettazione, essendo il p-value maggiore del livello di errore classico del 5% (linea rossa), non si rifiuta l’ipotesi nulla di indipendenza tra il tipo di parto e lo specifico ospedale. Cioè nessun ospedale registra un numero di cesarei statisticamnete differente dagli altri e i valori diversi registrati possono essere frutto del caso e non di una tendenza reale.

Test di ipotesi: medie del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione

Per saggiare se le medie di peso e lunghezza estratte dal campione siano o meno uguali a quelle vere della popolazione, che in letteratura sono stimate essere 3300 g e 500 mm, si esegue un test t.

P_med_pop=3300
L_med_pop=500

t.test(x=data$Peso, mu=P_med_pop, alternative='two.sided', conf.level=0.95)
## 
##  One Sample t-test
## 
## data:  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
t.test(x=data$Lunghezza, mu=L_med_pop, alternative='two.sided', conf.level=0.95)
## 
##  One Sample t-test
## 
## data:  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
plot(density(rt(100000, 2499)), xlab = "t", ylab = "p(t)", main = "Distribuzione t di Student per 2499 grandi di libertà-Peso")
abline(v=c(qt(0.975,2499),qt(0.025, 2499)), col=2)
points(x=-1.516, y=0, pch=19, col='blue')

plot(density(rt(100000, 2499)),  xlim = c(-11, 4), xlab = "t", ylab = "p(t)", main = "Distribuzione t di Student per 2499 grandi di libertà-Lunghezza")
abline(v=c(qt(0.975,2499),qt(0.025, 2499)), col=2)
points(x=-10.084, y=0, pch=19, col='blue')

Dal test t sulla media del peso, si nota come non si rifiuta l’ipotesi che la media del campione di statisticamente uguale a quella della popolazione, avendo un p-value superiore alla soglia di errore di 0.05 e ricadendo la statistica test all’interno del range di accettazione. Per la lunghezza invece, il test t mostra come si debba rifiutare ipotesi di uguaglianza tra media del campione e media della popolazione. Infatti la statistica test ricade ben oltre il range di accettazione e anche il p-value è prossimo allo zero, infine l’intervallo di confidenza non contiene il valore medio della popolazione.

Test di ipotesi: misure antropometriche sono significativamente diverse tra i due sessi

Si vuole testare ora se le misure antropometriche, quindi peso, lunghezza e diametro del cranio, siano significativamente diverse tra i due sessi o meno. Per fare ciò si usa un test t sulle medie di ognuna delle variabili tra i due sessi. Si usa la stima ‘pooled’ della varianza sulla variabile peso, in quanto le varianze tra i sessi non sono statisticamente diverse, come da test F. Mentre per lunghezza e cranio il test F indica differenza tra le varianze dei due sessi.

library(ggplot2)
peso_M=data$Peso[data$Sesso=='M']
peso_F=data$Peso[data$Sesso=='F']
var.test(peso_M, peso_F)
## 
##  F test to compare two variances
## 
## data:  peso_M and peso_F
## F = 0.88029, num df = 1243, denom df = 1255, p-value = 0.02436
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.7878460 0.9836202
## sample estimates:
## ratio of variances 
##          0.8802943
cat("var_peso_M= ", round(var(peso_M),0),"\n", "var_peso_F= ", round(var(peso_F),0),"\n")
## var_peso_M=  243843 
##  var_peso_F=  277001
lung_M=data$Lunghezza[data$Sesso=='M']
lung_F=data$Lunghezza[data$Sesso=='F']
var.test(lung_M, lung_F)
## 
##  F test to compare two variances
## 
## data:  lung_M and lung_F
## F = 0.76218, num df = 1243, denom df = 1255, p-value = 1.673e-06
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.6821345 0.8516401
## sample estimates:
## ratio of variances 
##          0.7621782
cat("var_lung_M= ", round(var(lung_M),0),"\n","var_lung_F= ", round(var(lung_F),0),"\n")
## var_lung_M=  578 
##  var_lung_F=  758
cranio_M=data$Cranio[data$Sesso=='M']
cranio_F=data$Cranio[data$Sesso=='F']
var.test(cranio_M, cranio_F)
## 
##  F test to compare two variances
## 
## data:  cranio_M and cranio_F
## F = 0.88484, num df = 1243, denom df = 1255, p-value = 0.03073
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.7919135 0.9886984
## sample estimates:
## ratio of variances 
##           0.884839
cat("var_cranio_M= ", round(var(cranio_M),0),"\n", "var_cranio_F= ", round(var(cranio_F),0),"\n")
## var_cranio_M=  248 
##  var_cranio_F=  280
t.test(peso_M, peso_F, mu=0, conf.level=0.95, alternative='two.sided', pool.sd=T)
## 
##  Welch Two Sample t-test
## 
## data:  peso_M and peso_F
## t = 12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  207.0615 287.1051
## sample estimates:
## mean of x mean of y 
##  3408.215  3161.132
plot(density(rt(100000, 2490.7)), xlab = "t", ylab = "p(t)", xlim=c(-4, 14), main = "Distribuzione t di Student per 2490.7 grandi di libertà - test sul peso")
abline(v=c(qt(0.975,2490.7),qt(0.025, 2490.7)), col=2)
points(x=12.106, y=0, pch=19, col='blue')

ggplot(data=data)+
  geom_boxplot(aes(x=factor(Sesso), y=Peso), fill='skyblue')+
  labs(x ='Sesso neonati' , y = "Peso neonato alla nascita", title = "Boxplot peso neonato in funzione del sesso") +
  theme_minimal()

t.test(lung_M, lung_F, mu=0, conf.level=0.95, alternative='two.sided')
## 
##  Welch Two Sample t-test
## 
## data:  lung_M and lung_F
## t = 9.582, df = 2459.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   7.876273 11.929470
## sample estimates:
## mean of x mean of y 
##  499.6672  489.7643
plot(density(rt(100000, 2459.3)), xlab = "t", ylab = "p(t)", xlim=c(-4, 12), main = "Distribuzione t di Student per 2459.3 grandi di libertà - test lunghezza")
abline(v=c(qt(0.975, 2459.3),qt(0.025, 2459.3)), col=2)
points(x=9.582, y=0, pch=19, col='blue')

ggplot(data=data)+
  geom_boxplot(aes(x=factor(Sesso), y=Lunghezza), fill='skyblue')+
  labs(x ='Sesso neonati' , y = "Lunghezza neonato alla nascita", title = "Boxplot lunghezza neonato in funzione del sesso") +
  theme_minimal()

t.test(cranio_M, cranio_F, mu=0, conf.level=0.95, alternative='two.sided')
## 
##  Welch Two Sample t-test
## 
## data:  cranio_M and cranio_F
## t = 7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  3.541270 6.089912
## sample estimates:
## mean of x mean of y 
##  342.4486  337.6330
plot(density(rt(100000, 2491.4)), xlab = "t", ylab = "p(t)", xlim=c(-4, 10), main = "Distribuzione t di Student per 2491.4 grandi di libertà - test cranio")
abline(v=c(qt(0.975,2491.4),qt(0.025, 2491.4)), col=2)
points(x=7.4102, y=0, pch=19, col='blue')

ggplot(data=data)+
  geom_boxplot(aes(x=factor(Sesso), y=Cranio), fill='skyblue')+
  labs(x ='Sesso neonati' , y = "Diametro cranio del neonato alla nascita", title = "Boxplot diametro cranio neonato in funzione del sesso") +
  theme_minimal()

Dai test risulta come la differenza tra i sessi in tutte e tre le variabili antropometriche sia statisticamente rilevante: infatti per tutti i test condotti il valore della statistica test cade ben oltre il range di accettazione dell’ipotesi nulla (ipotesi nulla = nessuna differenza tra i sessi) racchiuso nei quantili indicati in rosso che rappresentano il 2.5% di errore sulle due code. Inoltre i p-value sono circa zero, ad indicare che il risultato è statisticamente rilevante. Infine gli intervalli di confidenza generati per incapsulare il vero valore della differenza tra le medie di M e F, non contengono lo zero. Si conclude che i maschi presentano tutte è tre le grandezze antropometriche mediamente maggiori delle femmine.

Creazione del Modello di Regressione e selezione del modello ottimale

In questa sezione verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo, si potrà quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato, comprese eventuali interazioni. Ad esempio, si potrà verificare se una maggiore durata della gestazione aumenti in media il peso del neonato, come atteso.

Studio della correlazione tra le variabili e test di normalità del peso

Si inizia indagando la correlazione che sussiste tra le varie variabili del dataset. Nel modello la variabile risposta sarà il peso, mentre tutte le altre agiranno da regressori. Quindi in prima analisi è importante capire quali tra essi sono più legati alla risposta.

panel.cor=function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}
data_numeric=as.data.frame(data)
data_numeric=subset(data, select = c('Anni.madre','Gestazione','Peso','Lunghezza','Cranio'))


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

Il grafico delle correlazioni si può leggere in questo modo: sulla diagonale sono presenti le variabili quantitative del dataset, nella parte triangolare inferiore della matrice (sotto la diagonale) sono proposti i coefficienti di correlazione lineare tra le variabili, mentre nella parte superiore gli scatter plot viasualizzano la relazione tra le stesse. Si possono già fare alcune considerazioni. I regressori più influenti nel determinare il peso del neonato risultano essere:

  • la lunghezza e il diametro del cranio: biologicamente chi è più lungo o ha dimensioni del cranio maggiori alla nascita, pesa di più;
  • la durata della gestazione: più la gravidanza è lunga e più il peso finale cresce;

La variabile inerente agli anni della madre risulta poco influente sul peso alla nascita. Per le variabili qualitative del dataset, si è presentata la loro relazione con il peso tramite boxplot nella sezione iniziale. Si è notato come l’influenza del fumo sul peso sia presente graficamente, anche se non statisticamente rilevante, mentre il numero di gravidanze della madre o il tipo di parto risultano non particolarmente influenti osservando i grafici. Infine il sesso risulta influente sul peso, portando a pesi mediamente differenti tra maschi e femmine, come da test t eseguito nella precedente sezione.

C’è inoltre da segnalare una correlazione non banale tra i regressori stessi, in partciolare tra lunghezza-dimetro cranio e lunghezza-durata gestazione. In seguito verrà analizzata la possibile multicollinearità derivante da questo genere di correlazioni.

Si procede con la verifica della gaussianità della variabile Peso. Si attua un test di Shapiro-Wilk, come completamento alla visualizzazione della distribuzione del peso e degli indici di forma già ricavati precedentemente.

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

Dal test si evince come il valore della statistica W sia vicino a uno (1=gaussianità), ma il p-value piccolo suggerisca di rifiutare l’ipotesi nulla di normalità, essendo minore del livello di significatività al 5%. Cioè pur essendo la distribuzione del peso quasi normale stando al valore di W, ciò che la differenzia da una perfetta normale è statisticamente rilevante tanto da concludere che la variabile non può considerarsi come distribuita normalente. Anche se la variabile risposta non è gaussiana, ciò che conta per la validità del modello di regressione sono le ipotesi sui residui, che saranno opportunamente verificate in seguito.

Formulazione del modello

Si passa ora alla formulazione di un adeguato modello di regressione lineare multipla. Si inizia considerando un modello con tutte le variabili del dataset come regressori, con la variabile sesso come variabile di controllo. L’unica variabile esclusa è quella riguardante l’ospedale specifico da cui derivano i dati. Questo perchè si cerca un modello generale e anche perchè il peso medio dei nenonati è risultato essere simile tra i tre ospedali esaminati. Le variabili qualitative aventi le rispettive categorie, sono implementate nel modello sotto forma di variabili ‘dummy’ con tanti regressori x_i nella formuula del modello quante sono le loro categorie meno una, tranne nel caso del numero di gravidanze, dove si usano classi singole fino a 4 gravidanze, mentre le restanti possibilità vengono inglobate in un’unica classe ‘>4’ gravidanze, sulla base di considerazioni probabilistiche. Il contributo di questo tipo di variabile dovrà poi essere interpretato in maniera differente da una variabile quantitativa, riportando il suo effetto rispetto ad una delle categorie della variable stessa settata come livello di paragone base nel modello. Inoltre la variabile ‘Tipo.parto’ verrà tenuta nel modello scelto alla fine, perchè in alcuni casi si può stabilire in anticipo se il bimbo dovrà nascere cesareo o naturale, e quindi tale variabile potrà avere ruolo predittivo sul peso finale in tal senso.

Nel modello iniziale vengono considerate anche variabili non lineari, derivanti da un’ispezione visiva degli scatterplot nella matrice delle correlazioni. Infine si aggungono anche termini di interazione tra quelle variabili che si suppongono essere collegate e vicendevolmente influenzabili. Il modello migliore sarà quello con un valore di R^2 aggiustato maggiore, valori per i cirteri di Akaike (AIC) e Bayes (BIC) minimi, coefficienti stimati diversi da zero e significatività dei coefficienti elevata.

library(MASS)
data$Tipo.parto=as.factor(data$Tipo.parto)

data=data %>%
  mutate(gravidanze_raggruppate = case_when(
    as.numeric(as.character(N.gravidanze)) == 0 ~ "0",
    as.numeric(as.character(N.gravidanze)) == 1 ~ "1",
    as.numeric(as.character(N.gravidanze)) == 2 ~ "2",
    as.numeric(as.character(N.gravidanze)) == 3 ~ "3",
    as.numeric(as.character(N.gravidanze)) == 4 ~ "4",
    as.numeric(as.character(N.gravidanze)) > 4 ~ ">4"
  ))
data$gravidanze_raggruppate=as.factor(data$gravidanze_raggruppate)
data$gravidanze_raggruppate=relevel(data$gravidanze_raggruppate, ref = "0")


data$Fumatrici=as.factor(data$Fumatrici)
data$Sesso=as.factor(data$Sesso)

attach(data)
general_mod=lm(Peso~Lunghezza+I(Lunghezza^2)+Cranio+I(Cranio^0.5)+Gestazione+I(Gestazione^0.5)+Anni.madre+Sesso+Fumatrici+gravidanze_raggruppate+Tipo.parto+Lunghezza:Cranio+Lunghezza:Gestazione+Lunghezza:Fumatrici+Gestazione:Fumatrici, data=data)
summary_gen_mod=summary(general_mod)
kable(summary(general_mod)$coefficients, caption = "<b style='color:black;'>Coefficienti stimati per il modello 'general_mod'</b>")
Coefficienti stimati per il modello ‘general_mod’
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13236.8134506 7982.5439735 1.6582199 0.0973995
Lunghezza -16.0584434 5.2700520 -3.0471129 0.0023349
I(Lunghezza^2) 0.0521790 0.0064779 8.0549865 0.0000000
Cranio 125.6987135 29.6329506 4.2418561 0.0000230
I(Cranio^0.5) -3371.4369136 896.7186584 -3.7597488 0.0001740
Gestazione -150.8018413 308.3209184 -0.4891067 0.6248094
I(Gestazione^0.5) 3610.0357243 2760.1577281 1.3079092 0.1910253
Anni.madre 0.1316789 1.1295888 0.1165724 0.9072083
SessoM 73.1472409 10.9586236 6.6748566 0.0000000
Fumatrici1 -23.5450336 794.9370176 -0.0296187 0.9763735
gravidanze_raggruppate>4 54.5792021 40.5679360 1.3453778 0.1786262
gravidanze_raggruppate1 12.8503049 12.7337615 1.0091523 0.3130000
gravidanze_raggruppate2 58.6646634 17.3060331 3.3898389 0.0007103
gravidanze_raggruppate3 50.8577630 24.3573036 2.0879882 0.0369007
gravidanze_raggruppate4 96.8871887 40.2196358 2.4089524 0.0160709
Tipo.partoNat 29.7697586 11.7854759 2.5259700 0.0115997
Lunghezza:Cranio -0.0485999 0.0151138 -3.2156034 0.0013184
Lunghezza:Gestazione -0.2044895 0.1985639 -1.0298421 0.3031846
Lunghezza:Fumatrici1 4.3439569 1.4345533 3.0280903 0.0024864
Gestazione:Fumatrici1 -54.5228094 21.7133836 -2.5110232 0.0121013 
n=nrow(data)
step_mod=stepAIC(general_mod, direction='both', k=log(n))
## Start:  AIC=28069.94
## Peso ~ Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^0.5) + 
##     Gestazione + I(Gestazione^0.5) + Anni.madre + Sesso + Fumatrici + 
##     gravidanze_raggruppate + Tipo.parto + Lunghezza:Cranio + 
##     Lunghezza:Gestazione + Lunghezza:Fumatrici + Gestazione:Fumatrici
## 
##                          Df Sum of Sq       RSS   AIC
## - gravidanze_raggruppate  5   1245969 177851820 28048
## - Anni.madre              1       968 176606819 28062
## - Lunghezza:Gestazione    1     75526 176681377 28063
## - I(Gestazione^0.5)       1    121817 176727669 28064
## - Gestazione:Fumatrici    1    449009 177054860 28068
## - Tipo.parto              1    454370 177060222 28068
## <none>                                176605851 28070
## - Lunghezza:Fumatrici     1    652967 177258818 28071
## - Lunghezza:Cranio        1    736340 177342191 28072
## - I(Cranio^0.5)           1   1006633 177612484 28076
## - Sesso                   1   3172760 179778612 28107
## - I(Lunghezza^2)          1   4620437 181226288 28127
## 
## Step:  AIC=28048.4
## Peso ~ Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^0.5) + 
##     Gestazione + I(Gestazione^0.5) + Anni.madre + Sesso + Fumatrici + 
##     Tipo.parto + Lunghezza:Cranio + Lunghezza:Gestazione + Lunghezza:Fumatrici + 
##     Gestazione:Fumatrici
## 
##                          Df Sum of Sq       RSS   AIC
## - Lunghezza:Gestazione    1     69097 177920918 28042
## - I(Gestazione^0.5)       1    117462 177969282 28042
## - Anni.madre              1    211018 178062839 28044
## - Tipo.parto              1    401793 178253613 28046
## - Gestazione:Fumatrici    1    416629 178268449 28046
## <none>                                177851820 28048
## - Lunghezza:Fumatrici     1    612175 178463995 28049
## - Lunghezza:Cranio        1    741117 178592937 28051
## - I(Cranio^0.5)           1    987301 178839121 28054
## + gravidanze_raggruppate  5   1245969 176605851 28070
## - Sesso                   1   3289674 181141494 28086
## - I(Lunghezza^2)          1   4482788 182334608 28103
## 
## Step:  AIC=28041.55
## Peso ~ Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^0.5) + 
##     Gestazione + I(Gestazione^0.5) + Anni.madre + Sesso + Fumatrici + 
##     Tipo.parto + Lunghezza:Cranio + Lunghezza:Fumatrici + Gestazione:Fumatrici
## 
##                          Df Sum of Sq       RSS   AIC
## - Anni.madre              1    214085 178135003 28037
## - Gestazione:Fumatrici    1    401851 178322769 28039
## - Tipo.parto              1    408323 178329241 28040
## <none>                                177920918 28042
## - Lunghezza:Fumatrici     1    612046 178532963 28042
## - I(Gestazione^0.5)       1   1022277 178943195 28048
## + Lunghezza:Gestazione    1     69097 177851820 28048
## - I(Cranio^0.5)           1   1114074 179034991 28049
## - Lunghezza:Cranio        1   1216938 179137856 28051
## + gravidanze_raggruppate  5   1239540 176681377 28063
## - Sesso                   1   3255335 181176252 28079
## - I(Lunghezza^2)          1   6247813 184168730 28120
## 
## Step:  AIC=28036.73
## Peso ~ Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^0.5) + 
##     Gestazione + I(Gestazione^0.5) + Sesso + Fumatrici + Tipo.parto + 
##     Lunghezza:Cranio + Lunghezza:Fumatrici + Gestazione:Fumatrici
## 
##                          Df Sum of Sq       RSS   AIC
## - Gestazione:Fumatrici    1    395282 178530284 28034
## - Tipo.parto              1    404547 178539550 28035
## <none>                                178135003 28037
## - Lunghezza:Fumatrici     1    585009 178720012 28037
## + Anni.madre              1    214085 177920918 28042
## + Lunghezza:Gestazione    1     72164 178062839 28044
## - I(Gestazione^0.5)       1   1067009 179202012 28044
## - I(Cranio^0.5)           1   1138082 179273085 28045
## - Lunghezza:Cranio        1   1226442 179361445 28046
## + gravidanze_raggruppate  5   1452428 176682574 28055
## - Sesso                   1   3298515 181433518 28075
## - I(Lunghezza^2)          1   6284986 184419988 28116
## 
## Step:  AIC=28034.45
## Peso ~ Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^0.5) + 
##     Gestazione + I(Gestazione^0.5) + Sesso + Fumatrici + Tipo.parto + 
##     Lunghezza:Cranio + Lunghezza:Fumatrici
## 
##                          Df Sum of Sq       RSS   AIC
## - Lunghezza:Fumatrici     1    273576 178803860 28030
## - Tipo.parto              1    388513 178918797 28032
## <none>                                178530284 28034
## + Gestazione:Fumatrici    1    395282 178135003 28037
## + Anni.madre              1    207515 178322769 28039
## - Gestazione              1    998657 179528941 28041
## + Lunghezza:Gestazione    1     57113 178473172 28042
## - I(Cranio^0.5)           1   1107271 179637555 28042
## - I(Gestazione^0.5)       1   1178983 179709267 28043
## - Lunghezza:Cranio        1   1212329 179742614 28044
## + gravidanze_raggruppate  5   1414182 177116102 28054
## - Sesso                   1   3244777 181775061 28072
## - I(Lunghezza^2)          1   6452528 184982812 28115
## 
## Step:  AIC=28030.45
## Peso ~ Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^0.5) + 
##     Gestazione + I(Gestazione^0.5) + Sesso + Fumatrici + Tipo.parto + 
##     Lunghezza:Cranio
## 
##                          Df Sum of Sq       RSS   AIC
## - Fumatrici               1     34882 178838742 28023
## - Tipo.parto              1    382092 179185952 28028
## <none>                                178803860 28030
## + Lunghezza:Fumatrici     1    273576 178530284 28034
## + Anni.madre              1    189312 178614548 28036
## - Gestazione              1    969223 179773084 28036
## + Gestazione:Fumatrici    1     83849 178720012 28037
## + Lunghezza:Gestazione    1     63625 178740235 28037
## - I(Cranio^0.5)           1   1125057 179928918 28038
## - I(Gestazione^0.5)       1   1146547 179950407 28039
## - Lunghezza:Cranio        1   1218024 180021884 28040
## + gravidanze_raggruppate  5   1379794 177424066 28050
## - Sesso                   1   3324925 182128785 28069
## - I(Lunghezza^2)          1   6478506 185282366 28112
## 
## Step:  AIC=28023.11
## Peso ~ Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^0.5) + 
##     Gestazione + I(Gestazione^0.5) + Sesso + Tipo.parto + Lunghezza:Cranio
## 
##                          Df Sum of Sq       RSS   AIC
## - Tipo.parto              1    378070 179216813 28021
## <none>                                178838742 28023
## + Anni.madre              1    187482 178651260 28028
## - Gestazione              1    969233 179807976 28029
## + Lunghezza:Gestazione    1     63380 178775363 28030
## + Fumatrici               1     34882 178803860 28030
## - I(Cranio^0.5)           1   1131963 179970705 28031
## - I(Gestazione^0.5)       1   1145976 179984719 28031
## - Lunghezza:Cranio        1   1224426 180063168 28032
## + gravidanze_raggruppate  5   1349428 177489314 28043
## - Sesso                   1   3313742 182152484 28061
## - I(Lunghezza^2)          1   6498868 185337610 28104
## 
## Step:  AIC=28020.57
## Peso ~ Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^0.5) + 
##     Gestazione + I(Gestazione^0.5) + Sesso + Lunghezza:Cranio
## 
##                          Df Sum of Sq       RSS   AIC
## <none>                                179216813 28021
## + Tipo.parto              1    378070 178838742 28023
## + Anni.madre              1    184500 179032313 28026
## - Gestazione              1    970444 180187257 28026
## + Lunghezza:Gestazione    1     69607 179147206 28027
## + Fumatrici               1     30861 179185952 28028
## - I(Cranio^0.5)           1   1136401 180353214 28028
## - I(Gestazione^0.5)       1   1147988 180364800 28029
## - Lunghezza:Cranio        1   1212160 180428972 28030
## + gravidanze_raggruppate  5   1299827 177916986 28042
## - Sesso                   1   3311786 182528599 28058
## - I(Lunghezza^2)          1   6508021 185724833 28102
summary_step_mod=summary(step_mod)
kable(summary(step_mod)$coefficients, caption = "<b style='color:black;'>Coefficienti stimati per il modello 'step_mod'</b>")
Coefficienti stimati per il modello ‘step_mod’
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8722.6014861 6887.6411612 1.266414 0.2054835
Lunghezza -17.7277601 5.0959488 -3.478795 0.0005123
I(Lunghezza^2) 0.0482445 0.0050725 9.510907 0.0000000
Cranio 133.4274239 28.7436958 4.641972 0.0000036
I(Cranio^0.5) -3523.2632727 886.5062006 -3.974324 0.0000726
Gestazione -456.1653910 124.2051406 -3.672677 0.0002451
I(Gestazione^0.5) 6113.1423271 1530.3768798 3.994534 0.0000667
SessoM 74.5174532 10.9832123 6.784668 0.0000000
Lunghezza:Cranio -0.0554373 0.0135059 -4.104662 0.0000418 
refined_mod1=update(step_mod, ~.+Fumatrici)
summary_ref_mod1=summary(refined_mod1)
kable(summary(refined_mod1)$coefficients, caption = "<b style='color:black;'>Coefficienti stimati per il modello 'refined_mod1'</b>")
Coefficienti stimati per il modello ‘refined_mod1’
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8638.1331981 6889.638458 1.2537861 0.2100374
Lunghezza -17.7213834 5.096542 -3.4771384 0.0005155
I(Lunghezza^2) 0.0481817 0.005074 9.4957462 0.0000000
Cranio 133.0999863 28.751340 4.6293490 0.0000039
I(Cranio^0.5) -3513.6516946 886.729329 -3.9624850 0.0000763
Gestazione -456.1645959 124.219382 -3.6722497 0.0002455
I(Gestazione^0.5) 6114.6053757 1530.553989 3.9950276 0.0000666
SessoM 74.6469188 10.986251 6.7945764 0.0000000
Fumatrici1 -17.6403729 26.937547 -0.6548619 0.5126172
Lunghezza:Cranio -0.0553044 0.013509 -4.0938940 0.0000438 
refined_mod2=update(refined_mod1, ~.-Lunghezza:Cranio-I(Cranio^0.5))
summary_ref_mod2=summary(refined_mod2)
kable(summary(refined_mod2)$coefficients, caption = "<b style='color:black;'>Coefficienti stimati per il modello 'refined_mod2'</b>")
Coefficienti stimati per il modello ‘refined_mod2’
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.569374e+04 3710.9474288 -4.2290374 0.0000243
Lunghezza -3.106130e+01 4.0722599 -7.6275337 0.0000000
I(Lunghezza^2) 4.254030e-02 0.0041765 10.1856768 0.0000000
Cranio 1.060933e+01 0.4179068 25.3868431 0.0000000
Gestazione -4.522638e+02 113.5082217 -3.9844148 0.0000696
I(Gestazione^0.5) 6.070579e+03 1395.2703684 4.3508266 0.0000141
SessoM 7.380731e+01 11.0198713 6.6976565 0.0000000
Fumatrici1 -1.972720e+01 27.0297609 -0.7298327 0.4655611 
refined_mod3=update(general_mod, ~.-I(Lunghezza^2)-I(Gestazione^0.5)-I(Cranio^0.5)-Lunghezza:Cranio-Lunghezza:Gestazione-Lunghezza:Fumatrici-Gestazione:Fumatrici)
summary_ref_mod3=summary(refined_mod3)
kable(summary(refined_mod3)$coefficients, caption = "<b style='color:black;'>Coefficienti stimati per il modello 'refined_mod3', solo regressori lineari</b>")
Coefficienti stimati per il modello ‘refined_mod3’, solo regressori lineari
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6726.0447023 141.5742751 -47.5089468 0.0000000
Lunghezza 10.2640665 0.3008198 34.1203147 0.0000000
Cranio 10.4711214 0.4275774 24.4894159 0.0000000
Gestazione 32.8913279 3.8323528 8.5825417 0.0000000
Anni.madre 0.5650175 1.1585062 0.4877121 0.6257968
SessoM 77.7302673 11.1883167 6.9474497 0.0000000
Fumatrici1 -33.8115155 27.6113811 -1.2245500 0.2208608
gravidanze_raggruppate>4 40.4463272 41.6389922 0.9713570 0.3314650
gravidanze_raggruppate1 9.7104811 13.0379510 0.7447858 0.4564716
gravidanze_raggruppate2 51.2106970 17.7567411 2.8840144 0.0039602
gravidanze_raggruppate3 37.2761498 24.9551839 1.4937237 0.1353747
gravidanze_raggruppate4 87.6252911 41.2539219 2.1240475 0.0337646
Tipo.partoNat 30.9460934 12.1003979 2.5574443 0.0106033 
refined_mod4=update(refined_mod3, ~.-Anni.madre-gravidanze_raggruppate)
summary_ref_mod4=summary(refined_mod4)
kable(summary(refined_mod4)$coefficients, caption = "<b style='color:black;'>Coefficienti stimati per il modello 'refined_mod4', solo regressori lineari significativi</b>")
Coefficienti stimati per il modello ‘refined_mod4’, solo regressori lineari significativi
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6675.93348 135.7770904 -49.1683351 0.0000000
Lunghezza 10.22758 0.3010839 33.9691964 0.0000000
Cranio 10.63790 0.4242445 25.0749386 0.0000000
Gestazione 31.40674 3.7882430 8.2905817 0.0000000
SessoM 79.24775 11.2024687 7.0741325 0.0000000
Fumatrici1 -27.52451 27.5783292 -0.9980485 0.3183527
Tipo.partoNat 29.31670 12.1131895 2.4202296 0.0155818 
refined_mod5=update(refined_mod2, ~.+I(Cranio^0.5))
summary_ref_mod5=summary(refined_mod5)
kable(summary(refined_mod5)$coefficients, caption = "<b style='color:black;'>Coefficienti stimati per il modello 'refined_mod5'</b>")
Coefficienti stimati per il modello ‘refined_mod5’
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8660.2350378 5458.975098 -1.5864214 0.1127707
Lunghezza -27.7762071 4.479812 -6.2003063 0.0000000
I(Lunghezza^2) 0.0391943 0.004589 8.5409579 0.0000000
Cranio 43.3938241 18.673187 2.3238575 0.0202132
Gestazione -535.9758703 123.067512 -4.3551370 0.0000138
I(Gestazione^0.5) 7113.0849873 1515.768955 4.6927238 0.0000028
SessoM 73.4558466 11.017084 6.6674489 0.0000000
Fumatrici1 -19.2965789 27.019578 -0.7141702 0.4751888
I(Cranio^0.5) -1205.8276024 686.635920 -1.7561382 0.0791877 
#testo rilevanza di aggiunta variabili al modello

anova(step_mod, refined_mod1) #prima si inserisce il modello con la variabile, poi quello senza
## Analysis of Variance Table
## 
## Model 1: Peso ~ Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^0.5) + 
##     Gestazione + I(Gestazione^0.5) + Sesso + Lunghezza:Cranio
## Model 2: Peso ~ Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^0.5) + 
##     Gestazione + I(Gestazione^0.5) + Sesso + Fumatrici + Lunghezza:Cranio
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1   2491 179216813                           
## 2   2490 179185952  1     30861 0.4288 0.5126
#criteri di confronto tra modelli

AIC(general_mod, step_mod, refined_mod1, refined_mod2, refined_mod3, refined_mod4, refined_mod5)
##              df      AIC
## general_mod  21 35050.16
## step_mod     10 35064.84
## refined_mod1 11 35066.41
## refined_mod2  9 35082.28
## refined_mod3 14 35179.16
## refined_mod4  8 35182.82
## refined_mod5 10 35081.19
BIC(general_mod, step_mod, refined_mod1, refined_mod2, refined_mod3, refined_mod4, refined_mod5)
##              df      BIC
## general_mod  21 35172.46
## step_mod     10 35123.09
## refined_mod1 11 35130.48
## refined_mod2  9 35134.69
## refined_mod3 14 35260.70
## refined_mod4  8 35229.41
## refined_mod5 10 35139.43
adj_r_squared=c(
  round(summary_gen_mod$adj.r.squared, 3),
  round(summary_step_mod$adj.r.squared, 3),
  round(summary_ref_mod1$adj.r.squared, 3),
  round(summary_ref_mod2$adj.r.squared, 3),
  round(summary_ref_mod3$adj.r.squared, 3),
  round(summary_ref_mod4$adj.r.squared, 3),
  round(summary_ref_mod5$adj.r.squared, 3)
)


rmse_function=function(resid){
  mse_vals=sqrt(mean(resid^2))
  return(mse_vals)
}


rmse=c(
    round(rmse_function(residuals(general_mod)), 0),
    round(rmse_function(residuals(step_mod)), 0),
    round(rmse_function(residuals(refined_mod1)), 0),
    round(rmse_function(residuals(refined_mod2)), 0),
    round(rmse_function(residuals(refined_mod3)), 0),
    round(rmse_function(residuals(refined_mod4)), 0),
    round(rmse_function(residuals(refined_mod5)), 0)
)

models_table=data.frame(
  Modello = c("general_mod", "step_mod", "refined_mod1", "refined_mod2", "refined_mod3", "refined_mod4", "refined_mod5"),
  R_Squared_Adj = adj_r_squared,
  RMSE=rmse
)

kable(models_table)
Modello R_Squared_Adj RMSE
general_mod 0.742 266
step_mod 0.739 268
refined_mod1 0.739 268
refined_mod2 0.737 269
refined_mod3 0.727 273
refined_mod4 0.726 274
refined_mod5 0.737 269

Nell’analisi per selezionare il miglior modello partendo dal modello assunto come il più completo e conservativo, ‘general_mod’, si segue la procedura ‘Stepwise’ che consiste nell’aggiungere o togliere variabili dal modello inziale e verificarne la bontà di volta in volta usando indicatori come R^2 aggiustato e criteri AIC e BIC. Tale procedura per il caso esaminato conduce al modello migliore tra quelli estraibili partendo dal modello iniziale, denominato ‘step_mod’. Tale modello pur conservando le variabili maggiormente correlate con il peso e conservando alcuni andamenti non lineari e di interazione, scarta la variabile inerente al fumo materno. Dato che però nel particolare caso oggetto di studio si vuole considerare anche il possibile effetto del fumo delle madri sul peso alla nascita, si sceglie di tenere comunque tale variabile nel modello, pur avendo un test ANOVA che valuta il migliorameto in varianza spiegata tra modello con e senza variabile fumo che indica come la significatività statistica di tale aggiunta non sia rilevante. Il modello risultante è denominato ‘refined_mod1’. Si implementano anche altre varianti di modelli, a parte dalla procedura ‘Stepwise’, per capire se modelli con dipendenze più semplici siano comunque adeguati, nel tentativo di trovare non solo il modello più adatto, ma anche il modello che renda la formulazione stessa più immediata e semplice, a parità di bontà rispetto ai dati. Per questo si testano modelli come ‘refined_mod2’ che rispetto al modello precedente elimina alcune dipendenze non lineari e di interazione. Nel modello ‘refined_mod3’ si testa anche un tipo di approccio completamente lineare con tutte le variabili del dataset presenti. La stessa idea di linearità è mantenuta anche nel modello ‘refined_mod4’, pur conservando solo le variabili davvero significative. Infine il modello ‘refined_mod5’ si basa sul modello ‘refined_mod1’ eliminando solo il terminie di interazione.

Per stabilire il modello migliore si confrontano le stime dei valori di R^2 aggiustato, Root mean square error (RMSE), AIC e BIC. Il modello migliore a livello teorico è quello avente un R^2 aggiustato più elevato (spiega meglio i dati), un RMSE basso (differenza media tra valori previsti e effettivi è piccola) ed infine valori di AIC e BIC minimi. Da tali analisi il modello migliore risulterebbe essere indicativamente il ‘general_mod’ o al massimo il modello ‘step_mod’. Il primo però, è sovrabbondante in regressori e non si presta ad una interpretazione intuitiva. Il secondo pur avendo meno regressori, manca del fumo. Si potrebbe quindi concludere prendendo il modello ‘refined_mod1’ come modello migliore, ma ciò non è ancora detto, perchè ci sono ancora da analizzare i residui.

Studio dei residui

Si procede a studiare i residui del modello ‘refined_mod1’ per capire se rispettano le ipotesi di gaussianità, media zero, omoschedasticità (varianza costante) e indipendenza. Si studiano inoltre eventuali ‘outliers’ (cioè valori dei residui particolarmente grandi) e valori ‘leverages’ (cioè valori nello spazio dei regresori che pesano di più sul modello di regressione).

library(lmtest)

models=list(refined_mod1, refined_mod5, refined_mod4)
i=1
for(mod in models){
if(i==1){name='refined_mod1'
         n_par=10}
  else if(i==2){name='refined_mod5' 
                n_par=9}
else{name='refined_mod4'  
     n_par=7}  
  
par(mfrow=c(2,2))
plot(mod, main=paste("Diagnostica residui ", name))

#media 0
cat("TEST MEDIA 0 RESIDUI ", name)
t=t.test(residuals(mod))
print(t)

#gaussianità
cat("TEST GAUSSIANITÀ RESIDUI ", name)
sh=shapiro.test(residuals(mod))
print(sh)

#omoschedasticità
cat("TEST OMOSCHEDASTICITÀ RESIDUI ", name)
bp=bptest(mod)
print(bp)

#indipendenza
cat("TEST INDIPENDENZA RESIDUI ", name)
dw=dwtest(mod)
print(dw)

#outliers e leverages
par(mfrow=c(1,1))
cook=cooks.distance(mod)
plot(cook, main = paste("Distanza di Cook ", name))
i=i+1

#residui studentizzati
res_stud=rstudent(mod)
plot(res_stud, main=paste("Residui studentizzati modello ", name))
res_max=res_stud[res_stud==max(res_stud)]
tab_res=kable(res_max, caption=paste("Valore outlier massimo modello ", name))
print(tab_res)

#leverages
leverages=hatvalues(mod)
limit=2*n_par/n


lev_max=leverages[leverages==max(leverages)]


tab_lev=kable(lev_max, caption=paste("Valore leverage massimo modello ", name))
print(tab_lev)
}

## TEST MEDIA 0 RESIDUI  refined_mod1
##  One Sample t-test
## 
## data:  residuals(mod)
## t = -2.4694e-15, df = 2499, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -10.50164  10.50164
## sample estimates:
##     mean of x 
## -1.322481e-14 
## 
## TEST GAUSSIANITÀ RESIDUI  refined_mod1
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod)
## W = 0.99084, p-value = 1.595e-11
## 
## TEST OMOSCHEDASTICITÀ RESIDUI  refined_mod1
##  studentized Breusch-Pagan test
## 
## data:  mod
## BP = 108.66, df = 9, p-value < 2.2e-16
## 
## TEST INDIPENDENZA RESIDUI  refined_mod1
##  Durbin-Watson test
## 
## data:  mod
## DW = 1.9511, p-value = 0.1106
## alternative hypothesis: true autocorrelation is greater than 0

## 
## 
## Table: Valore outlier massimo modello  refined_mod1
## 
## |     |        x|
## |:----|--------:|
## |1551 | 5.195386|
## 
## 
## Table: Valore leverage massimo modello  refined_mod1
## 
## |     |        x|
## |:----|--------:|
## |1551 | 0.614033|

## TEST MEDIA 0 RESIDUI  refined_mod5
##  One Sample t-test
## 
## data:  residuals(mod)
## t = -5.5519e-15, df = 2499, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -10.53693  10.53693
## sample estimates:
##     mean of x 
## -2.983289e-14 
## 
## TEST GAUSSIANITÀ RESIDUI  refined_mod5
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod)
## W = 0.98835, p-value = 2.052e-13
## 
## TEST OMOSCHEDASTICITÀ RESIDUI  refined_mod5
##  studentized Breusch-Pagan test
## 
## data:  mod
## BP = 125.15, df = 8, p-value < 2.2e-16
## 
## TEST INDIPENDENZA RESIDUI  refined_mod5
##  Durbin-Watson test
## 
## data:  mod
## DW = 1.9534, p-value = 0.1221
## alternative hypothesis: true autocorrelation is greater than 0

## 
## 
## Table: Valore outlier massimo modello  refined_mod5
## 
## |     |        x|
## |:----|--------:|
## |1551 | 6.606117|
## 
## 
## Table: Valore leverage massimo modello  refined_mod5
## 
## |     |         x|
## |:----|---------:|
## |1551 | 0.2757956|

## TEST MEDIA 0 RESIDUI  refined_mod4
##  One Sample t-test
## 
## data:  residuals(mod)
## t = 2.4982e-15, df = 2499, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -10.76191  10.76191
## sample estimates:
##   mean of x 
## 1.37109e-14 
## 
## TEST GAUSSIANITÀ RESIDUI  refined_mod4
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod)
## W = 0.97428, p-value < 2.2e-16
## 
## TEST OMOSCHEDASTICITÀ RESIDUI  refined_mod4
##  studentized Breusch-Pagan test
## 
## data:  mod
## BP = 89.498, df = 6, p-value < 2.2e-16
## 
## TEST INDIPENDENZA RESIDUI  refined_mod4
##  Durbin-Watson test
## 
## data:  mod
## DW = 1.9562, p-value = 0.1364
## alternative hypothesis: true autocorrelation is greater than 0

## 
## 
## Table: Valore outlier massimo modello  refined_mod4
## 
## |     |       x|
## |:----|-------:|
## |1551 | 9.97964|
## 
## 
## Table: Valore leverage massimo modello  refined_mod4
## 
## |     |         x|
## |:----|---------:|
## |1551 | 0.0486682|

Per studiare i residui del modello si usano diversi grafici e test. Nei grafici sono presenti i 4 plot diagnostici dei residui, il grafico sulla distanza di Cook e il grafico dei residui studentizzati. I grafici diagnostici sui residui riportano 4 tipi di analisi:

  • il primo in alto a sinistra analizza la media dei residui e l’eventuale presenza di pattern mostra che il modello non cattura bene tutta l’informazione
  • il secondo in alto a destra analizza la gaussianità
  • in basso a sinistra si ha lo studio dell’omoschedasticità
  • in basso a destra outliers e leverages.

Dallo studio sui residui si può notare come il ‘refined_mod1’ abbia residui a media 0 (t-test), pur mostrando un accenno di curvatura su bassi valori delle x nel primo grafico diagnostico dei residui, indice di informazione non filtrata dal modello. Tale modello risulta però mancante delle ipotesi di gaussianità (Shapiro-Wilk test) e omoschedasticità (Breush-Pagan test), pur avendo residui indipendenti (Durbin-Watson test). Lo studio della distanza di Cook, che aiuta ad individuare i dati che hanno un contributo elevato nel modello in termini di outlier dei residui o di leverage, suggerisce la presenza di un dato ben oltre la soglia di distanza pari a 1, che è la soglia di allarme di influenza. Tale dato corrisponde all’osservazione 1551 del dataset e questa descrive un neonato con un peso ben oltre la media per il valore di lunghezza associato. Per questo motivo si valuta anche il ‘refined_mod5’ per osservare se la rimozione del termine che li distingue, quello di interazione tra lunghezza del neonato e del cranio, sia causa di miglioramento, e infine il ‘refined_mod4’, che implementa solo regressori lineari significativi. Si può vedere come pur conservando gli stessi risultati sulle ipotesi sui residui già visti per il ‘refined_mod1’, gli ultimi due modelli abbiano una distanza di Cook per il punto precedentemente menzionato molto minore, specie nel modello lineare dove tale punto non è neanche più oltre la soglia di allarme. C’è però da notare come considerando i residui studentizzati, lo stesso punto rappresenti un significativo outlier. Questo porta a concludere che tale punto pur essendo mal catturato dal modello (alto residuo) risulta poco influente sulla stima del modello stesso. Come ulteriore analisi, si può notare anche come tale punto problematico nei vari modelli abbia valori di residuo e leverage in dipendenza inversa, al crescere di uno l’altro decresce, indicando come le richieste di avere un modello ben fittato sui punti del dataset(basso residuo) e allo stesso tempo generalizzabile(basso leverage), e quindi poco legato al campione specifico, siano in competizione. Tornando alla scelta del modello, dato che il modello è fatto con obiettivi di predizione, la presenza di punti leverage con distanza di Cook elevata non giova in tal senso, perchè il modello è molto dipendente dal campione specifico esaminato e si presta meno ad essere generalizzato su nuovi dati. Per tale motivo, si sceglie di considerare come miglior modello il ‘refined_mod4’, che pur rinunciando ad integrare termini non lineari e di interazione, ben si presta per scopi predittivi. Inoltre la differenza in R^2 aggiustato e RMSE per tale modello e i modelli favoriti da tali criteri non è eccessiva, seppur non sia trascurabile. In aggiunta, tutti i coefficienti del ‘refined_mod4’ risultano avere valore stimato decisamente diverso da zero e una significatività massima, con uniche eccezioni per il tipo di parto, che rasenta la soglia di significatività e per il fumo che è stato aggiunto a mano per scopi di indagine. Proprio in riferimento a tale modello, si può notare come l’effetto del fumo sul peso, già osservato graficamente in precedenza, sia negativo, portando a neonati con peso minore in concomitanza a madri fumatrici. Inoltre la variabile qualitativa dicotomica Sesso, indica come neonati maschi pesino mediamente 80 grammi in più rispetto alle femmine, assunte come riferimento. Infine, l’effetto della duarata della gestazione è positivo sul peso, mostrando come una gestazione più lunga di una settimana possa portare ad un peso maggiore di circa 31 grammi, mentre il tipo di parto indica come un bimbo nato naturale, pesi in media 30 grammi in più di un bimbo nato con il cesareo. In riferimento al tipo di parto, il modello ha mostrato un legame tra questa variabile ed il peso che dai grafici precedenti non era ovvio.

Come nota a margine del presente modello è importante dire che, date le ipotesi disattese sui residui, sarebbe utile valutare alternative ad una regressione lineare, ad esempio si potrebbe testare una modellizzazione derivante dai modelli lineari generalizzati o GLM, più adatti nei contesti in cui le ipotesi della regressione lineare vengono meno. Tale ulteriore indagine però non rientra negli scopi del presente progetto.

Studio della multicollinearità

Si può infine indagare la presenza di multicollinearità usando gli indicatori VIF (Variance Inflation Factor) che devono mostrare valori inferiori a 5 per dimostrare assenza di multicollinearità.

car::vif(refined_mod4)
##  Lunghezza     Cranio Gestazione      Sesso  Fumatrici Tipo.parto 
##   2.078848   1.607612   1.659003   1.039087   1.004316   1.003059

I VIF per il modello ‘refined_mod4’ sono tutti sotto la soglia di 5, ad indicare come non vi sia multicollinearità tra i regressori del mdoello.

Predizioni usando il modello

Una volta selezionato il modello più adatto per rappresentare la relazione tra il peso dei neonati e le altre varibili del dataset, si vuole testare il modello per fare previsioni utili. C’è comunque da premettere che tale modello lineare è un’approssimazione della relazione vera sul range di valori dei regressori che si hanno disponibili e nulla si sa sull’evoluzione di tale relazione fuori da tali ranges, inoltre fare previsioni biologicamente insensate non sarebbe realistico.

lun=c(470,490,500,500,355)
cra=c(298,325,340,344,270)
gest=c(34,42,41,37,31)
sex=c('M','M','F','M','F')   # o 0 o 1
fum=c('0','0','0','1','0')  # o 0 o 1
parto=c('Ces','Nat','Nat','Ces','Nat')
real=c(2400,3380,3300,3280,1180)

par_list=list()
estim_list=list()

par_list[[1]]=c('N.','Lunghezza','Cranio','Gestazione','Sesso','Fumatrice','Tipo.parto')
estim_list[[1]]=c('N.','Valore vero','Stima','Lower bd','Upper bd')

for(i in c(1,2,3,4,5)){

par_list[[i+1]]=c(i,lun[i],cra[i],gest[i],sex[i],fum[i],parto[i])

new_point=data.frame(Lunghezza=lun[i], Cranio=cra[i], Gestazione=gest[i], Sesso=factor(sex[i], levels = c('F', 'M')), Fumatrici=factor(fum[i], levels = c('0', '1')), Tipo.parto=factor(parto[i], levels=c('Nat','Ces')))

pred=predict(refined_mod4, newdata = new_point, interval = "prediction", level = 0.95)
estim_list[[i+1]]=c(i,real[i],round(pred[1],0),round(pred[2],0),round(pred[3],0))
}

df_par=do.call(rbind, par_list)
df_est=do.call(rbind, estim_list)

kable(df_par, caption = "<b style='color:black;'>Tabella valori regressori</b>")%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1:7, extra_css = "padding-right: 2px;")
Tabella valori regressori
N. Lunghezza Cranio Gestazione Sesso Fumatrice Tipo.parto
1 470 298 34 M 0 Ces
2 490 325 42 M 0 Nat
3 500 340 41 F 0 Nat
4 500 344 37 M 1 Ces
5 355 270 31 F 0 Nat 
kable(df_est, caption = "<b style='color:black;'>Tabella valori stimati</b>")%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1:5, extra_css = "padding-right: 2px;")
Tabella valori stimati
N. Valore vero Stima Lower bd Upper bd
1 2400 2448 1908 2989
2 3380 3221 2681 3760
3 3300 3372 2833 3911
4 3280 3311 2769 3853
5 1180 830 288 1372

Le stime di predizione sono accompagnate dall’intervallo di predizione che rappresenta l’incertezza associata alla stima previsionale del modello. Questi intervalli considerano sia l’incertezza associata alle stime dei parametri del modello, sia quanto il modello sia adatto a rappresentare i dati, tramite i suoi residui. Nelle due tabelle precedenti si testa il modello su 5 unità statistiche casuali selezionate dal dataset, per sondare l’efficienza del modello nel riproporre valori del peso già noti. Si nota come le stime siano tendenzialmente vicine ai valori veri noti, con unica eccezione rilevante per la quinta stima, pur conservando un certo errore tipicamente dell’ordine delle decine/centinaia di grammi. Nelle due tabelle si adottano le stesse numerazioni per la stessa unità statistica. Di seguito si testa il modello su nuovi dati.

lun=c(472,350,500,570,570,570,570)
cra=c(295,289,350,400,400,400,400)
gest=c(34,36,39,40,40,40,36)
sex=c('M','M','F','F','M','M','M')   # o 0 o 1
fum=c('0','0','1','0','0','1','0')  # o 0 o 1
parto=c('Nat','Ces','Ces','Nat','Nat','Nat','Nat')

par_list=list()
estim_list=list()

par_list[[1]]=c('N.','Lunghezza','Cranio','Gestazione','Sesso','Fumatrice', 'Tipo.parto')
estim_list[[1]]=c('N.','Stima','Lower bd','Upper bd')

for(i in c(1,2,3,4,5,6,7)){

par_list[[i+1]]=c(i,lun[i],cra[i],gest[i],sex[i],fum[i], parto[i])

new_point=data.frame(Lunghezza=lun[i], Cranio=cra[i], Gestazione=gest[i], Sesso=factor(sex[i], levels = c('F', 'M')), Fumatrici=factor(fum[i], levels = c('0', '1')), Tipo.parto=factor(parto[i], levels=c('Nat','Ces')))

pred=predict(refined_mod4, newdata = new_point, interval = "prediction", level = 0.95)
estim_list[[i+1]]=c(i,round(pred[1],0),round(pred[2],0),round(pred[3],0))
}

df_par=do.call(rbind, par_list)
df_est=do.call(rbind, estim_list)

kable(df_par, caption = "<b style='color:black;'>Tabella valori regressori</b>")%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1:7, extra_css = "padding-right: 2px;")
Tabella valori regressori
N. Lunghezza Cranio Gestazione Sesso Fumatrice Tipo.parto
1 472 295 34 M 0 Nat
2 350 289 36 M 0 Ces
3 500 350 39 F 1 Ces
4 570 400 40 F 0 Nat
5 570 400 40 M 0 Nat
6 570 400 40 M 1 Nat
7 570 400 36 M 0 Nat 
kable(df_est, caption = "<b style='color:black;'>Tabella valori stimati</b>")%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(1:4, extra_css = "padding-right: 2px;")
Tabella valori stimati
N. Stima Lower bd Upper bd
1 2466 1925 3007
2 1188 644 1731
3 3358 2817 3900
4 4695 4153 5236
5 4774 4233 5315
6 4746 4203 5290
7 4648 4105 5191

Nelle due tabelle mostrate si è testato il modello su valori dei parametri all’interno dei rispettivi ranges derivati dal dataset. Su tali dati i valori del peso stimati risultano vicini a valori simili presenti nel dataset, mostrando come il modello dia risultati verosimili. Negli ultimi quattro set di parametri test, si è voluta testare a parità delle altre variabili, la differenza tra maschi e femmine, tra madri fumatrici e non e tra gestazioni di durata diversa. Correttamente il modello riporta un peso maggiore per i maschi e tra maschi, il neonato con madre fumatrice riporta un peso inferiore. Anche la gestazione ha un effetto ben modellizzato, portando ad un peso minore in corrispondenza di gestazioni più brevi. Anche se non riportati in tabella, si è testa anche l’abilità del modello di prevedere pesi per valori delle variabili estremi, cioè fuori dai ranges indicati nel dataset. Si è notato come il modello in tali contesti non sia più adeguato portando a risultati non relastici come valori negativi, a ribadire l’importanza di utilizzare il modello solo su range biologicamente validi dei regressori.

Visualizzazioni

Si vuole concludere questa analisi con delle visualizzazioni del modello di regressione in funzione di alcuni dei regressori che lo compongono.

ggplot(data, aes(x = Lunghezza, y = Peso)) + 
  geom_point(col='skyblue') + 
  geom_smooth(method = "lm", formula = y ~ x, se=FALSE, aes(col=factor(Sesso))) +
  labs(title = "Retta di regressione Peso~Lunghezza+Sesso", x = "Lunghezza", y = "Peso",  color = "Sesso")+
  theme_classic()

ggplot(data, aes(x = Gestazione, y = Peso)) + 
  geom_point(col='skyblue') + 
  geom_smooth(method = "lm", formula = y ~ x, se=FALSE, aes(col=factor(Sesso))) +
  labs(title = "Retta di regressione Peso~Gestazione+Sesso", x = "Gestazione", y = "Peso",  color = "Sesso")+
  theme_classic()

ggplot(data, aes(x = Gestazione, y = Peso)) + 
  geom_point(col='skyblue') + 
  geom_smooth(method = "lm", formula = y ~ x, se=FALSE, aes(col=factor(Fumatrici))) +
  labs(title = "Retta di regressione Peso~Gestazione+Fumatrici", x = "Gestazione", y = "Peso",  color = "Fumatrici")+
  theme_classic()

ggplot(data, aes(x = Gestazione, y = Peso)) + 
  geom_point(col='skyblue') + 
  geom_smooth(method = "lm", formula = y ~ x, se=FALSE, aes(col=factor(Tipo.parto))) +
  labs(title = "Retta di regressione Peso~Gestazione+Tipo di parto", x = "Gestazione", y = "Peso",  color = "Tipo di parto")+
  theme_classic()

ggplot(data, aes(x = Cranio, y = Peso)) + 
  geom_point(col='skyblue') + 
  geom_smooth(method = "lm", formula = y ~ x, se=FALSE, aes(col=factor(Sesso))) +
  labs(title = "Retta di regressione Peso~Cranio+Sesso", x = "Diametro cranio", y = "Peso",  color = "Sesso")+
  theme_classic()

Nei grafici è possibile visualizzare nuovamente l’importanza della relazione tra peso e gestazione, mentre per la relazione con il fumo, si nota visivamente la penuria di dati relativi a madri fumatrici per poter fare previsoni accurate e statisticamente rilevanti, seppur si possa già notare come il fumo riduca il peso a parità di gestazione.

Conclusioni

Nel presente documento si è studiato il peso neonatale per comprenderne le correlazioni con altre variabili di interesse, come le variabili antropometriche del neonato (lunghezza, diametro del cranio), la durata della gravidanza e variabili legate alla madre, prime fra tutte l’età e l’incidenza del fumo. Dall’analisi esplorativa iniziale e dai test statistici successivi è emerso come il peso del neonato sia dipendnente dalla gestazione, con pesi minori associati a gravidanze più brevi. Inoltre il peso risulta statisticamente differente tra maschi e femmine. Infine il fumo materno sembra diminuire il peso alla nascita, ma non sembra alterare significativamente la durata della gestazione. Questa caratteristica seppur prevista dalla letteratura, non si manifesta bel dataset studiato presumibilmente per le poche occorrenze di madri fumatrici. Dai test statistici è inoltre risultato che le medie del peso e della lunghezza stimate dal dataset sono rappresentative dell’intera popolazione, il che fa assumere allo studio eseguito sul campione specifico un carattere generalizzabile. Successivamente si è messo a punto un modello di regressione lineare multipla che permettesse di rappresentare efficacemente il legame tra peso e altre variabili e che fosse allo stesso tempo il più intuitivo possibile. Su questo schema, considerando anche la significatività dei coefficienti e il loro valore, si è scelto il modello ‘refined_mod4’ come il più rappresentativo, pur non essendo il miglior modello secondo criteri dell’ R^2 aggiustato o AIC/BIC, ed essendo un modello puramente lineare e pur avendo alcune pecche relative alla normalità ed omoschedasticità dei residui, anomalie queste ultime che in relatà sono comuni a tutti i modelli di regressione lineare indagati nello studio. Un punto di forza di tale modello è l’assenza di valori leverages problematici, anzi la sua scelta è legata proprio a tale caratteristica che lo rende meno dipendente dal campione e quindi più gerenalizzabile. Si è testato tale modello sia su dati veri che su valori delle variabili arbitrari per valutarne la bontà, ottenendo risultati verosmili. Tale modello potrà essere utile per prevedere il peso dei nascituri stanti alcuni valori delle variabili implementate come regressori, per aiutare le strutture ospedaliere a garantire ai neonati trattamenti adeguati, dando particolare attenzione all’influenza della durata della gestazione e approfondendo ulteriormente l’analisi dell’impatto del fumo, collezionando campioni più numerosi con tale caratteristica.