Modello Statistico per la Previsione del Peso Neonatale

Contesto Aziendale

Azienda: Neonatal Health Solutions

Obiettivo: Creare un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali. Il progetto mira a migliorare la gestione delle gravidanze ad alto rischio, ottimizzare le risorse ospedaliere e garantire migliori risultati per la salute neonatale.

Il progetto si inserisce all’interno di un contesto di crescente attenzione verso la prevenzione delle complicazioni neonatali. La possibilità di prevedere il peso alla nascita dei neonati rappresenta un’opportunità fondamentale per migliorare la pianificazione clinica e ridurre i rischi associati a nascite problematiche, come parti prematuri o neonati con basso peso. Di seguito, i principali benefici che questo progetto porterà all’azienda e al settore sanitario:

Miglioramento delle previsioni cliniche: Il peso del neonato è un indicatore chiave della sua salute. Avere un modello predittivo accurato consente al personale medico di intervenire tempestivamente in caso di anomalie, riducendo le complicazioni perinatali come le difficoltà respiratorie o l’ipoglicemia.
Ottimizzazione delle risorse ospedaliere: Sapere in anticipo quali neonati potrebbero avere bisogno di cure intensive aiuta a organizzare le risorse umane e tecnologiche degli ospedali in modo efficiente. Questo si traduce in una riduzione dei costi operativi e una migliore pianificazione dell’utilizzo delle unità di terapia intensiva neonatale (TIN).
Prevenzione e identificazione dei fattori di rischio: Il modello potrà evidenziare i fattori che maggiormente influenzano negativamente il peso del neonato (come il fumo materno, gravidanze multiple o età avanzata della madre). Queste informazioni sono preziose per la prevenzione e la gestione personalizzata delle gravidanze, permettendo interventi proattivi in caso di rischio elevato.
Valutazione delle pratiche ospedaliere: Attraverso un’analisi comparativa tra i tre ospedali coinvolti, l’azienda potrà identificare eventuali differenze nei risultati clinici, come una maggiore incidenza di parti cesarei in una determinata struttura. Ciò consente di monitorare la qualità delle pratiche e armonizzare i protocolli tra i diversi centri ospedalieri, migliorando la coerenza delle cure.
Supporto alla pianificazione strategica: L’analisi dei dati e le previsioni possono essere utilizzate per prendere decisioni informate non solo a livello clinico ma anche strategico. L’azienda potrà sfruttare queste informazioni per implementare nuove politiche di salute pubblica, garantendo un impatto positivo sui tassi di mortalità e morbilità neonatale.

Caricamento librerie necessarie

library(ggplot2)
library(moments)
library(patchwork)
library(car)

## Loading required package: carData

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:patchwork':
## 
##     area

library(pscl)

## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.

library(scatterplot3d)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.3.0
## ✔ purrr     1.1.0     ✔ tidyr     1.3.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::recode() masks car::recode()
## ✖ dplyr::select() masks MASS::select()
## ✖ purrr::some()   masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

1. Raccolta dei Dati e Struttura del Dataset

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

Età della madre: Misura dell’età in anni.
Numero di gravidanze: Quante gravidanze ha avuto la madre.
Fumo materno: Un indicatore binario (0=non fumatrice, 1=fumatrice).
Durata della gravidanza: Numero di settimane di gestazione.
Peso del neonato: Peso alla nascita in grammi.
Lunghezza e diametro del cranio: Lunghezza del neonato e diametro craniale, misurabili anche durante la gravidanza tramite ecografie.
Tipo di parto: Naturale o cesareo.
Ospedale di nascita: Ospedale 1, 2 o 3.
Sesso del neonato: Maschio (M) o femmina (F).

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

dati <- read.csv("neonati.csv",sep = ",",stringsAsFactors = T) 
knitr::kable(head(dati,5))

Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
26	0	42	3380	490	325	Nat	osp3	M
21	2	39	3150	490	345	Nat	osp1	F
34	3	38	3640	500	375	Nat	osp2	M
28	1	41	3690	515	365	Nat	osp2	M
20	0	38	3700	480	335	Nat	osp3	F

N <- dim(dati) [1]
cat("Il dataset contiene",N,"osservazioni\n")

## Il dataset contiene 2500 osservazioni

duplicati <- dati %>% group_by_all() %>% filter(n() > 1) %>% ungroup()
cat("Il dataset contiene",dim(duplicati) [1] ,"duplicati\n") #non ci sono duplicati

## Il dataset contiene 0 duplicati

attach(dati)

2. Analisi e Modellizzazione

2.1 Analisi Preliminare

Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.

str(dati)

## '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  : Factor w/ 2 levels "Ces","Nat": 2 2 2 2 2 2 2 2 1 1 ...
##  $ Ospedale    : Factor w/ 3 levels "osp1","osp2",..: 3 1 2 2 3 2 3 1 2 2 ...
##  $ Sesso       : Factor w/ 2 levels "F","M": 2 1 2 2 1 1 1 2 1 1 ...

summary(dati)

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

Che tipo di variabili presenta il dataset?

Variabili quantitative continua

Peso, Lunghezza, Cranio

Variabili quantitative discrete

Anni.madre, N.gravidanze, Gestazione

Variabili qualitative dummy (o fattori)

Fumatrici

Variabili qualitative nominali

Tipo.parto, Ospedale, Sesso

Ma andiamo ancora più nel dettaglio.

indici_vis <- function(data, feature){
  
  #Gestione dati e label
  feature_name <- deparse(substitute(feature))  # nome come stringa
  values <- pull(data, {{feature}})  # estrae il vettore
  
  # Calcoli statistici
  mu <- mean(values, na.rm = TRUE)
  sigma <- sd(values, na.rm = TRUE)
  cv <- sigma / mu
  sk <- skewness(values, na.rm = TRUE)
  ku <- kurtosis(values, na.rm = TRUE) - 3  # kurtosi centrata
  
  # Crea df per la label informativa
  info_label <- data.frame(x = mu, y = 0,
                           label = paste0("μ = ", round(mu,2),
                                          "\nσ = ", round(sigma,2),
                                          "\nCV = ", round(cv,2),
                                          "\nSkew = ", round(sk,2),
                                          "\nKurt = ", round(ku,2)))
  
  # Quantili
  q_vals <- quantile(values, probs = c(0.25, 0.5, 0.75), na.rm = TRUE)
  mean_val <- mean(values, na.rm = TRUE)
  
  # Plot dei record ordinati
  record_plot <- ggplot() + 
    geom_point(aes(x = seq_along(values), y = sort(values))) +
    geom_hline(yintercept = mu, color = "blue", linetype = "dashed", linewidth = 0.7) +
    geom_hline(yintercept = q_vals, color = "red", linetype = "dotted" , linewidth = 0.7) +
    geom_label(aes(x = length(values)*0.9, y = mu),
               label = paste0("Media: ", round(mu, 1)),
               color = "blue", size = 3) +
    geom_label(data = data.frame(y = q_vals, label = names(q_vals)),
               aes(x = length(values)*0.9, y = y, label = paste0(label, ": ", round(y, 1))),
               color = "red", size = 3) +
    labs(title = "Distribuzione osservazioni ordinate",
         x = "n. osservazione", y = "Valore osservazione") +
    theme_minimal()
  
  # Boxplot
  box_plot <- ggplot(data, aes(y = {{feature}})) +
    geom_boxplot(fill = "tan1") +
    geom_hline(yintercept = mu, color = "blue", linetype = "dashed", linewidth = 0.7) +
    geom_hline(yintercept = q_vals, color = "red", linetype = "dotted", linewidth = 0.7) +
    labs(title = "Boxplot",
         y = "Valore osservazione") +
    theme_minimal()+
    theme(axis.text.x = element_blank(),  
          axis.ticks.x = element_blank())
  
  # Grafico della densità con box dei valori
  density_plot <- ggplot(data, aes(x = {{feature}})) +
    geom_density(fill = "lightblue", color = "black", alpha = 0.6) +
    geom_vline(xintercept = mu, color = "blue", linetype = "dashed", linewidth = 0.7) +
    geom_label(data = info_label, aes(x = x, y = y, label = label),
               hjust = -0.1, vjust = -0.1, size = 3, fill = "white", label.size = 0.3) +
    labs(title = "Stima della densità di probabilità",
         x = "Valore osservazione", y = "Densità osservazione") +
    theme_minimal() 
  
  print((record_plot + box_plot) +
          plot_annotation(title = paste("Analisi di:", feature_name,"\n"),
                          theme = theme(plot.title = element_text(size = 16))))
  print(density_plot)
} 

indici_vis(dati,Peso) # sotto 1000 gr?

indici_vis(dati,Anni.madre) # Gravidanza a meno di 10 anni da escudere

indici_vis(dati,Gestazione) # Sotto le 30 settimane sono gravidanze premature

indici_vis(dati,Lunghezza)

indici_vis(dati,Cranio)

distr_freq_Ngrav<-as.data.frame(
  cbind(
    ni=table(N.gravidanze),
    fi=round(table(N.gravidanze)/N,2)))
knitr::kable(distr_freq_Ngrav) # 10-12 gravidanze??

	ni	fi
0	1096	0.44
1	818	0.33
2	340	0.14
3	150	0.06
4	48	0.02
5	21	0.01
6	11	0.00
7	1	0.00
8	8	0.00
9	2	0.00
10	3	0.00
11	1	0.00
12	1	0.00

distr_freq_fumo<-as.data.frame(
  cbind(
    ni=table(Fumatrici),
    fi=round(table(Fumatrici)/N,2)))
knitr::kable(distr_freq_fumo) # Non fumatori poco rappresentati. Solo 4%

	ni	fi
0	2396	0.96
1	104	0.04

distr_freq_parto<-as.data.frame(
  cbind(
    ni=table(Tipo.parto),
    fi=round(table(Tipo.parto)/N,2)))
knitr::kable(distr_freq_parto) # Prevalenza di parti narurtali

	ni	fi
Ces	728	0.29
Nat	1772	0.71

distr_freq_osped<-as.data.frame(
  cbind(
    ni=table(Ospedale),
    fi=round(table(Ospedale)/N,2)))
knitr::kable(distr_freq_osped) # Dati perfettamente distribuiti tra i 3 ospedali

	ni	fi
osp1	816	0.33
osp2	849	0.34
osp3	835	0.33

distr_freq_sesso<-as.data.frame(
  cbind(
    ni=table(Sesso),
    fi=round(table(Sesso)/N,2)))
knitr::kable(distr_freq_sesso) # Equa distribuzione tra maschi e femmine

	ni	fi
F	1256	0.5
M	1244	0.5

Considerazione dall’analisi delle variabili

Commentiamo le variabili in analisi in termini di distribuzione, outlier e anomalie

Anni.madre presenta due osservazioni con età inferiore a 10 anni. Considerando che biologicamente è molto improbabile (sviluppo del sistema riproduttivo in età compresa tra i 12-15 anni), possiamo considerare i record associati a questo valore di Anni.madre come degli outlier.
Guardando la variabile Gestazione, si osserva la presenza di gravizande premature con osservazioni il cui valore è inferiore alle 30 settimane. Sotto le 30 settimane sono gravidanze premature.
Il punto precendente si collega al peso dei neonati. Ci sono osservazioni per cui il Peso dei neonati assume valori sotto i 1000 gr

Queste considerazioni son rintracciabili nei seguenti subset:

knitr::kable(subset(dati, Anni.madre < 10))

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

knitr::kable(subset(dati, Gestazione < 30))

	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
310	40	3	28	1560	420	379	Nat	osp3	F
378	27	0	28	1285	400	274	Nat	osp1	F
805	30	2	29	1190	360	272	Nat	osp2	F
928	25	0	28	830	310	254	Nat	osp1	F
1385	33	0	29	1620	410	292	Nat	osp3	F
1429	24	4	29	1280	390	355	Nat	osp1	F
1780	25	2	25	900	325	253	Nat	osp3	F
2120	32	0	27	1140	370	267	Nat	osp3	F
2175	37	8	28	930	355	235	Nat	osp1	F
2437	28	1	27	980	320	265	Nat	osp1	M
2452	28	0	26	930	345	245	Ces	osp3	F

knitr::kable(subset(dati, Peso < 1000))

	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
928	25	0	28	830	310	254	Nat	osp1	F
1619	31	0	31	990	340	278	Ces	osp2	F
1780	25	2	25	900	325	253	Nat	osp3	F
2175	37	8	28	930	355	235	Nat	osp1	F
2437	28	1	27	980	320	265	Nat	osp1	M
2452	28	0	26	930	345	245	Ces	osp3	F

#knitr::kable(subset(dati, Fumatrici==TRUE))

Anni.madre pari a 0 e 1 sono chiaramente OUTLIERS.
E l’effetto del fumo? Guardando i subset con Gestazione < 30 e Peso < 1000 si vede che non ci sono Fumatrici. Il fumo non sembra collegato ai parti prematuri.
Si inltre nota poi che fumatori e parti cesarei poco rappresentati con rispettivamene il 4% e il 30%
Mentre c’è una equadistribuzione di maschi e femmine, e di parti tra i 3 ospedali coinvolti nell’analisi.

Inoltre si saggeranno le seguenti ipotesi con i test adatti:

1. In alcuni ospedali si fanno più parti cesarei

Facciamo innanzitutto un controllo sui dati incrociando i parti cesarei agli ospedali coinvolti

mean(Tipo.parto == 'Ces')

## [1] 0.2912

# Filtra solo i casi con parto = "Ces"
dati_ces <- subset(dati, Tipo.parto == "Ces")
N_ces <- nrow(dati_ces)

distr_freq <- as.data.frame(
  cbind(
    ni = table(dati_ces$Ospedale),
    fi = round(table(dati_ces$Ospedale) / N_ces, 2)
  )
)
knitr::kable(distr_freq)

	ni	fi
osp1	242	0.33
osp2	254	0.35
osp3	232	0.32

Si può creare una tabella di contingenza e usare un test del Chiquadroto per saggiare l’ipotesi di indipendenza tra le variabili Ospedale e Tipo.parto

H0: indipendenza
H1: associazione

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

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

# Per la visualizzazione: install.packages("ggpubr")
ggpubr::ggballoonplot(data=as.data.frame(tab_cont), fill="blue")

# Calcolo delle attese
attese <- tab_cont

for(i in 1:nrow(tab_cont)){
  for(j in 1:ncol(tab_cont)){
    attese[i,j] <- margin.table(tab_cont,1)[i] * margin.table(tab_cont,2)[j] / margin.table(tab_cont)
  }
}
# Mentre le osservate sono
osservate<-tab_cont

# Test fatto a mano
chiquadro.manuale <- sum((osservate-attese)^2/attese)
chiquadro.manuale

## [1] 1.097202

plot(density(rchisq(1000000,df = 2)), xlim=c(0,40))
abline(v=qchisq(0.95,2),col=2)
points(chiquadro.manuale,0,cex=3,col=4)

# Usando la funzione chisq.test
test.indipendenza <- chisq.test(tab_cont)
test.indipendenza

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

Con un p-value di 0.5778 e un alpha di 0.05, accetto l’ipotesi nulla di indipendenza: i dati raccolti ci dicono che non si apprezzano differenze di tipo di parto tra i diversi ospedali.

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

Per il peso, in accordo a quanto dichiarato dal WHO (fonte: https://www.who.int/tools/child-growth-standards/standards/weight-for-age) la media dei neonati risulta:

neonati maschi: 3300 gr
neanate femmine: 3200 gr

Visto che il campione presenta una equidistribuzione di maschi-femmine, è possibile considerare un valore medio come media della popolazione (a valle faremo ulteriore considerazioni). Le ipotesi del test quindi diventano:

H0: mu popolazione = mu campione
H1: mu popolazione != mu campione

Usiamo un t-test a causa della dimensione ridotta del campione, della distribuzione non normale della variabile e della varianza della popolazione non nota a priori (varia tanto geograficamente).

Effetuiamo il test di normalità della variabile peso per confermare la nostra affermazione di non normalità e poi passiamo al t test per confronto campione-popolazione.

shapiro.test(Peso)

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

# Peso ha una distribuzione non normale
mu_peso_popol <- 3250 #media delle medie della popolazione tra maschi e femmine
risu_ttest <- t.test(Peso,mu=mu_peso_popol) #di default two.sided e 0.95 di conf.level
risu_ttest

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

plot(density(rt(100000,2499)),xlim=c(-4,4))
abline(v=qt(0.025,2499))
qt(0.025,2499)

## [1] -1.960914

abline(v=qt(0.975,2499))
qt(0.975,2499)

## [1] 1.960914

# abline(v=risu_ttest$statistic, col=2)
points(risu_ttest$statistic,0,cex=3,col=2)

Con un p-value di 0.0012 e un alpha di 0.05, rifiuto l’ipotesi nulla: la media della popolazione differisce significativamente dalla media del campione. NOTA si potrebbe approfondire differenziando il campione in base al sesso, considerando però che si avrebbero sempre meno dati a disposizione.

cat("Media popolazione neonati maschi: 3300 gr")

## Media popolazione neonati maschi: 3300 gr

mean(subset(Peso, Sesso == 'M'))

## [1] 3408.215

cat("Media popolazione neonati maschi: 3200 gr")

## Media popolazione neonati maschi: 3200 gr

mean(subset(Peso, Sesso == 'F'))

## [1] 3161.132

Passando alla lunghezza, in accordo a quanto dichiarato dal WHO (fonte: https://www.who.int/tools/child-growth-standards/standards/length-height-for-age), la media dei neonati risulta di:

neonati maschi: 500 mm
neanate femmine: 490 mm

Visto che il campione presenta una equidistribuzione di maschi-femmine, è possibile considerare un valore medio come media della popolazione. Analogamente a quanto visto per il peso:

shapiro.test(Lunghezza)

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

Le ipotesi del test quindi diventano:

H0: mu popolazione = mu campione
H1: mu popolazione != mu campione

mu_lungh_popol <- 495 #media delle medie della popolazione tra maschi e femmine
risu_ttest <- t.test(Lunghezza,mu=mu_lungh_popol)
risu_ttest

## 
##  One Sample t-test
## 
## data:  Lunghezza
## t = -0.58514, df = 2499, p-value = 0.5585
## alternative hypothesis: true mean is not equal to 495
## 95 percent confidence interval:
##  493.6598 495.7242
## sample estimates:
## mean of x 
##   494.692

plot(density(rt(100000,2499)),xlim=c(-4,4))
abline(v=qt(0.025,2499))
qt(0.025,2499)

## [1] -1.960914

abline(v=qt(0.975,2499))
qt(0.975,2499)

## [1] 1.960914

# abline(v=risu_ttest$statistic, col=2)
points(risu_ttest$statistic,0,cex=3,col=2)

Con un p-value di 0.5585 e un alpha di 0.05, non posso rifiutare l’ipotesi nulla: la media campionaria è signitificativamente uguale a quella della popolazione

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

Prepariamo i dati necessari al test:

Lunghezza.M <- subset(Lunghezza, Sesso == 'M')
Lunghezza.F <- subset(Lunghezza, Sesso == 'F')
Cranio.M <- subset(Cranio, Sesso == 'M')
Cranio.F <- subset(Cranio, Sesso == 'F')

Usiamo il t test per il confronto delle medie campionarie. Occorre verificare che:

Le variabili misurate devono essere distribuite normalmente (test per verificare la normalità)
Le varianze tra gruppi devono essere almeno simili (test per verificare l’omogeneità tra varianze)

Le ipotesi del test quindi diventano:

H0: mu_campione1 - mu_campione2 = 0
H1: mu_campione1 - mu_campione2 != 0

Sappiamo che Lunghezza non è una variabile con una distribuzione normale. Lo stesso vale per Cranio, di cui a seguire il test di normalità. Date queste premesse, usiamo un test non parametrico: usiamo pairwise.wilcox.test.

summary(Lunghezza[Sesso=='M'])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   320.0   490.0   500.0   499.7   515.0   560.0

summary(Lunghezza[Sesso=='F'])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   310.0   480.0   490.0   489.8   505.0   565.0

boxplot(Lunghezza~Sesso, ylab = 'Lunghezza [mm]')

pairwise.wilcox.test(Lunghezza, Sesso,
                     paired = F, #non appaiati
                     pool.sd = T,
                     p.adjust.method = "none") # Son significativamente diversi

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  Lunghezza and Sesso 
## 
##   F     
## M <2e-16
## 
## P value adjustment method: none

shapiro.test(Cranio)

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

summary(Cranio[Sesso=='M'])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   265.0   334.0   343.0   342.4   352.0   390.0

summary(Cranio[Sesso=='F'])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   235.0   330.0   340.0   337.6   348.2   390.0

boxplot(Cranio~Sesso, ylab = 'Cranio [mm]')

pairwise.wilcox.test(Cranio, Sesso) # Son significativamente diversi

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  Cranio and Sesso 
## 
##   F      
## M 9.6e-15
## 
## P value adjustment method: holm

In entrambi casi, non si può accettare l’ipotesi nulla (in entrambi i casi con un p-value prossimo allo zero). Possiamo affermare che esistono differenze significative nelle misure antropometriche tra i due sessi: sia sulla lunghezza che sulle misure del cranio.

2.2 Creazione del Modello di Regressione

Verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo, potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato ed eventuali interazioni. Ad esempio, ci aspettiamo che una maggiore durata della gestazione aumenterebbe in media il peso del neonato.

Dalle precedenti analisi emerge che la variabile Peso non presenta una distribuzione normale (con un p-value prossimo allo zero).

Potremmo quindi vagliare l’ipotesi di fare una correzione sulla variabile risposta oppure usare i GLM, ma prima di tutto andiamo a studiare la linearità del modello e se le ipotesi alla base vengono rispettate.

Iniziamo con lo studio della correlazione delle variabili in gioco.

# Studiamo le correlazioni in anticipo
# round(cor(dati),2) #non è possobile per via delle variabili qualitative
# Mentre utilizzando la funzione pair
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr")
  on.exit(par(usr = 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)
}
#correlazioni
pairs(dati,lower.panel=panel.cor, upper.panel=panel.smooth)

Si può osservare una alta correlazione con Gestazione, Lunghezza, Cranio (e in parte Sesso).

Verifichiamo a parte le variabili qualitative. Lo faremo sia graficamente (boxplot) che non un ttest.

par(mfrow=c(1,2))
boxplot(Peso, ylab='Peso [gr]', xlab='Campione')
boxplot(Peso~Sesso, ylab = 'Peso [gr]')

t.test(Peso~Sesso) #differenze significative

## 
##  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

par(mfrow=c(1,2))
boxplot(Peso, ylab='Peso [gr]', xlab='Campione')
boxplot(Peso~Ospedale, ylab = 'Peso [gr]')

pairwise.wilcox.test(Peso, Ospedale) #differenze non significative

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  Peso and Ospedale 
## 
##      osp1 osp2
## osp2 0.64 -   
## osp3 0.49 0.26
## 
## P value adjustment method: holm

par(mfrow=c(1,2))
boxplot(Peso, ylab='Peso [gr]', xlab='Campione')
boxplot(Peso~Tipo.parto, ylab = 'Peso [gr]')

t.test(Peso~Tipo.parto) #differenze non significative

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

Viene confermata la correlazione tra Peso e Sesso e dai test possiamo affermare che esistono differenze significative di peso nei neonati di diverso sesso. Al contrario, sulla base dei dati a disposizione, non si osservano differenze significative di peso nei diversi ospedali e in base al tipo di parto.

Siamo pronti quindi a costruire passo passo un primo modello lineare che comprenda tutte le variabili. Negli step successivi andremo a migliorare il modello.

mod<-lm(Peso ~ . ,
        data=dati)
summary(mod)

## 
## Call:
## lm(formula = Peso ~ ., data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1124.40  -181.66   -14.42   160.91  2611.89 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6738.4762   141.3087 -47.686  < 2e-16 ***
## Anni.madre        0.8921     1.1323   0.788   0.4308    
## N.gravidanze     11.2665     4.6608   2.417   0.0157 *  
## Fumatrici       -30.1631    27.5386  -1.095   0.2735    
## Gestazione       32.5696     3.8187   8.529  < 2e-16 ***
## Lunghezza        10.2945     0.3007  34.236  < 2e-16 ***
## Cranio           10.4707     0.4260  24.578  < 2e-16 ***
## Tipo.partoNat    29.5254    12.0844   2.443   0.0146 *  
## Ospedaleosp2    -11.2095    13.4379  -0.834   0.4043    
## Ospedaleosp3     28.0958    13.4957   2.082   0.0375 *  
## SessoM           77.5409    11.1776   6.937 5.08e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared:  0.7289, Adjusted R-squared:  0.7278 
## F-statistic: 669.2 on 10 and 2489 DF,  p-value: < 2.2e-16

# modelsummary(mod, output = "markdown")
# stargazer(mod, type = "html", title = "Modello di regressione completo")
AIC(mod)

## [1] 35171.95

BIC(mod)

## [1] 35241.84

Leggendo il summary si nota che le variabili più significative sono le settimane di gestazione, seguite da la variabili antropometriche Lunghezza e Cranio. I regressori non significativi sono Anni.madre, Fumatrici, Ospedale e, in seconda fascia, N.gravidanze e Tipo.parto.

Per far fronte alla non normalità della variabile peso (con possibili impatti sulla distribuzione dei residui) potremmo provare altre strategie: trasformazione di boxcox oppure modelli GLM. Per effettuare un confronto diretto, ci soffermeremo sui GLM e usero il BIC come parametro di confronto.

# boxcox(lm(Peso ~ . , data = dati))
# lambda <- 0.5
# shapiro.test(Peso^0.5)
# model_boxcox <- lm(Peso^0.5 ~ . , data = dati)
# summary(model_boxcox)

model_gamma <- glm(Peso ~ . , family = Gamma(link="log"), data = dati)
summary(model_gamma)

## 
## Call:
## glm(formula = Peso ~ ., family = Gamma(link = "log"), data = dati)
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.417e+00  4.580e-02  96.441  < 2e-16 ***
## Anni.madre     2.799e-04  3.670e-04   0.763   0.4457    
## N.gravidanze   3.581e-03  1.511e-03   2.370   0.0178 *  
## Fumatrici     -7.025e-03  8.925e-03  -0.787   0.4313    
## Gestazione     1.883e-02  1.238e-03  15.214  < 2e-16 ***
## Lunghezza      3.401e-03  9.746e-05  34.897  < 2e-16 ***
## Cranio         3.604e-03  1.381e-04  26.098  < 2e-16 ***
## Tipo.partoNat  6.746e-03  3.917e-03   1.722   0.0851 .  
## Ospedaleosp2  -2.857e-03  4.355e-03  -0.656   0.5119    
## Ospedaleosp3   9.416e-03  4.374e-03   2.153   0.0314 *  
## SessoM         1.793e-02  3.623e-03   4.949 7.96e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.0078821)
## 
##     Null deviance: 75.182  on 2499  degrees of freedom
## Residual deviance: 18.828  on 2489  degrees of freedom
## AIC: 35309
## 
## Number of Fisher Scoring iterations: 5

# Confrontiamo i due modelli
BIC(mod,model_gamma)

##             df      BIC
## mod         12 35241.84
## model_gamma 12 35378.94

Dal confronto emerge che il modello di regressione GLM presenta un BIC maggiore del primo modello base sviluppato. Nel paragrafo successivo, approfondiremo la questione.

2.3 Selezione del Modello Ottimale

Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC), selezioneremo il modello più parsimonioso, eliminando le variabili non significative. Verranno considerati anche modelli con interazioni tra le variabili e possibili effetti non lineari.

Iniziamo la procedura stepwise eliminando i regrossori non significativi individuati in precedenza partendo da Anni.madre e Fumatrici, per poi continuare con Ospedale, Tipo.gravidanza e N.gravidanza.

# Procedura stepwise
mod1<-lm(Peso ~ . -Anni.madre -Fumatrici,
        data=dati)
summary(mod1)

## 
## Call:
## lm(formula = Peso ~ . - Anni.madre - Fumatrici, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1113.18  -181.16   -16.58   161.01  2620.19 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6707.4293   135.9438 -49.340  < 2e-16 ***
## N.gravidanze     12.3619     4.3325   2.853  0.00436 ** 
## Gestazione       31.9909     3.7896   8.442  < 2e-16 ***
## Lunghezza        10.3086     0.3004  34.316  < 2e-16 ***
## Cranio           10.4922     0.4254  24.661  < 2e-16 ***
## Tipo.partoNat    29.2803    12.0817   2.424  0.01544 *  
## Ospedaleosp2    -11.0227    13.4363  -0.820  0.41209    
## Ospedaleosp3     28.6408    13.4886   2.123  0.03382 *  
## SessoM           77.4412    11.1756   6.930 5.36e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared:  0.7287, Adjusted R-squared:  0.7278 
## F-statistic: 836.3 on 8 and 2491 DF,  p-value: < 2.2e-16

# l'R2 adj rimane più o meno lo stesso. Eliminiamo anche Ospedale
mod2<-update(mod1, ~ . -Ospedale,
         data=dati)
summary(mod2)

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

# Stesso discorso. L'R2 adj rimane più o meno lo stesso. Eliminiamo anche Tipo.parto e N.gravidanze
mod3<-update(mod2, ~ . -Tipo.parto -N.gravidanze,
             data=dati)
summary(mod3)

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

#AIC(mod3,mod2,mod1,mod,model_gamma)
#anova(mod3,mod2,mod1,mod,model_gamma)
BIC(mod3,mod2,mod1,mod,model_gamma)

##             df      BIC
## mod3         6 35220.54
## mod2         8 35221.75
## mod1        10 35228.03
## mod         12 35241.84
## model_gamma 12 35378.94

Usando come criterio la minimizzazione del BIC, il modello migliore sembra essere il modello 3 in cui il Peso è linearmente dipendente da Gestazione, Lunghezza, Cranio e Sesso. A questo punto possiamo prevare la procedura automatica e comparare i risultati. Come ultimo test prima della procedura automatica, proviamo a reintrodurre il fumo ma con delle relazioni non lineari.

mod4<-update(mod3, ~ . + Fumatrici*Anni.madre,
             data=dati)
summary(mod4)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     Fumatrici + Anni.madre + Fumatrici:Anni.madre, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1164.20  -184.25   -14.04   162.73  2615.88 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -6724.5661   141.3994 -47.557  < 2e-16 ***
## Gestazione              32.4306     3.8264   8.475  < 2e-16 ***
## Lunghezza               10.2031     0.3012  33.878  < 2e-16 ***
## Cranio                  10.6000     0.4260  24.883  < 2e-16 ***
## SessoM                  78.8492    11.2127   7.032 2.62e-12 ***
## Fumatrici              -39.9458   146.9672  -0.272   0.7858    
## Anni.madre               1.9410     1.0806   1.796   0.0726 .  
## Fumatrici:Anni.madre     0.4584     5.1005   0.090   0.9284    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.9 on 2492 degrees of freedom
## Multiple R-squared:  0.7266, Adjusted R-squared:  0.7258 
## F-statistic:   946 on 7 and 2492 DF,  p-value: < 2.2e-16

BIC(mod4,mod3,mod2,mod1,mod,model_gamma)

##             df      BIC
## mod4         9 35239.65
## mod3         6 35220.54
## mod2         8 35221.75
## mod1        10 35228.03
## mod         12 35241.84
## model_gamma 12 35378.94

La variabile Fumatrici si conferma essere non significativa. Il mod 3 rimane il migliore.

Procedura automatica

# Per ripetere tutta la procedura in automatico
stepwise.mod <- MASS::stepAIC(mod,
              direction = "both",
              k=log(N)) #AIC con k=2, log(n) per BIC

## Start:  AIC=28139.32
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     46578 186809099 28132
## - Fumatrici     1     90019 186852540 28133
## - Ospedale      2    685979 187448501 28133
## - N.gravidanze  1    438452 187200974 28137
## - Tipo.parto    1    447929 187210450 28138
## <none>                      186762521 28139
## - Sesso         1   3611021 190373542 28179
## - Gestazione    1   5458403 192220925 28204
## - Cranio        1  45326172 232088693 28675
## - Lunghezza     1  87951062 274713583 29096
## 
## Step:  AIC=28132.12
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     90897 186899996 28126
## - Ospedale      2    692738 187501837 28126
## - Tipo.parto    1    448222 187257321 28130
## <none>                      186809099 28132
## - N.gravidanze  1    633756 187442855 28133
## + Anni.madre    1     46578 186762521 28139
## - Sesso         1   3618736 190427835 28172
## - Gestazione    1   5412879 192221978 28196
## - Cranio        1  45588236 232397335 28670
## - Lunghezza     1  87950050 274759149 29089
## 
## Step:  AIC=28125.51
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Ospedale      2    701680 187601677 28119
## - Tipo.parto    1    440684 187340680 28124
## <none>                      186899996 28126
## - N.gravidanze  1    610840 187510837 28126
## + Fumatrici     1     90897 186809099 28132
## + Anni.madre    1     47456 186852540 28133
## - Sesso         1   3602797 190502794 28165
## - Gestazione    1   5346781 192246777 28188
## - Cranio        1  45632149 232532146 28664
## - Lunghezza     1  88355030 275255027 29086
## 
## Step:  AIC=28119.23
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Tipo.parto    1    463870 188065546 28118
## <none>                      187601677 28119
## - N.gravidanze  1    651066 188252743 28120
## + Ospedale      2    701680 186899996 28126
## + Fumatrici     1     99840 187501837 28126
## + Anni.madre    1     54392 187547285 28126
## - Sesso         1   3649259 191250936 28160
## - Gestazione    1   5444109 193045786 28183
## - Cranio        1  45758101 233359778 28657
## - Lunghezza     1  88054432 275656108 29074
## 
## Step:  AIC=28117.58
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      188065546 28118
## - N.gravidanze  1    623141 188688687 28118
## + Tipo.parto    1    463870 187601677 28119
## + Ospedale      2    724866 187340680 28124
## + Fumatrici     1     91892 187973654 28124
## + Anni.madre    1     54816 188010731 28125
## - Sesso         1   3655292 191720838 28158
## - Gestazione    1   5464853 193530399 28181
## - Cranio        1  46108583 234174130 28658
## - Lunghezza     1  87632762 275698308 29066

summary(stepwise.mod)

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

summary(mod3)

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

La procedura automatica permette di giungere ad un modello molto simile al miglior modello trovato manualmente con al sola aggiunta della variabile numero di gravidanze. A questo punto proviamo a costruire un glm usando le stesse variabili identificare dalla procedura stepwise.

model_gamma_2 <- glm(Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso, family = Gamma(link="log"), data = dati)
model_gamma_2

## 
## Call:  glm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + 
##     Cranio + Sesso, family = Gamma(link = "log"), data = dati)
## 
## Coefficients:
##  (Intercept)  N.gravidanze    Gestazione     Lunghezza        Cranio  
##     4.432756      0.004017      0.018799      0.003386      0.003623  
##       SessoM  
##     0.018088  
## 
## Degrees of Freedom: 2499 Total (i.e. Null);  2494 Residual
## Null Deviance:       75.18 
## Residual Deviance: 18.93     AIC: 35310

BIC(mod3,model_gamma,stepwise.mod,model_gamma_2)

##               df      BIC
## mod3           6 35220.54
## model_gamma   12 35378.94
## stepwise.mod   7 35220.10
## model_gamma_2  7 35353.69

Il miglior modello, utilizzando come parametro di scelta il BIC, risulta essere il modello scelto con la procedura stepwise. Il modello presenta un R² aggiustato di circa 0.727.

Tuttavia, consideranzo l’importanza scientifica della variabile Fumatrici, la aggiungeremo al modello per futuri sviluppi del modello. Si fa notare che la sua introduzione fa aumentare leggermente il BIC ma non altera l’R2.

stepwise.mod.final<-update(stepwise.mod, ~ . + Fumatrici,
             data=dati)
summary(stepwise.mod.final)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Fumatrici, data = dati)
## 
## 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 ** 
## 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 ***
## Fumatrici      -30.4634    27.5948  -1.104  0.26972    
## ---
## 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

BIC(stepwise.mod,stepwise.mod.final)

##                    df     BIC
## stepwise.mod        7 35220.1
## stepwise.mod.final  8 35226.7

2.4 Analisi della Qualità del Modello

Una volta ottenuto il modello finale, valuteremo la sua capacità predittiva utilizzando metriche come R² e il Root Mean Squared Error (RMSE). Un’attenzione particolare sarà rivolta all’analisi dei residui e alla presenza di valori influenti, che potrebbero distorcere le previsioni, indagando su di essi.

L’analisi del modello parte dal controllo dei problemi di multicollinearità.

# Gli indicatori devono stare sotto alla soglia di 5
knitr::kable(car::vif(stepwise.mod.final), col.names = c('multi_collinearità_indice'))

	multi_collinearità_indice
N.gravidanze	1.026120
Gestazione	1.675575
Lunghezza	2.078644
Cranio	1.624603
Sesso	1.040271
Fumatrici	1.006607

Non sussistono problemi di multicollinearità. Passiamo all’analisi dei residui

# Aanlisi dei residui
par(mfrow=c(2,2))
plot(stepwise.mod.final)

# attenzione all'osservazione 1551
knitr::kable(dati[1551,])

	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
1551	35	1	0	38	4370	315	374	Nat	osp3	F

Abbiamo individuato una osservazione pericolosa che va oltre la distanza di cook di 0.5, avvicinandosi molto a 1. Effettuiamo quindi i singoli test per studiare passo passo come si comportano le osservazioni e i residui e valutare, quindi, la qualità del modello.

# Indagine solle osservazioni leverage
par(mfrow=c(1,1))
lev<-hatvalues(stepwise.mod.final)
plot(lev, xlab = "Osservazioni", ylab = "Valore di leverage", main = "Osservazioni Leverage")
p<-sum(lev)
n<-length(lev)
soglia=2*p/n
abline(h=soglia,col=2) # oss 1551 ben visibile

#lev[lev>soglia]
length(lev[lev>soglia]) # ci 152 osservazione che sono lontane dallo spazio dei regressori, delle variabili esplicative

## [1] 211

# Studio degli outliers
plot(rstudent(stepwise.mod.final), xlab = "Osservazioni", ylab = "Valore rstudent outlier", main = "Osservazioni Outliers")
abline(h=c(-2,2))  # oss 1551 ben visibile

# facciamo un ttest sugli outliers
car::outlierTest(stepwise.mod.final) # osservazione 1511, 155 e 1306 da capire come gestire

##       rstudent unadjusted p-value Bonferroni p
## 1551 10.039719         2.8060e-23   7.0149e-20
## 155   5.022108         5.4723e-07   1.3681e-03
## 1306  4.823102         1.4986e-06   3.7465e-03

# Ed infine, distanza di cook per valutare entrambe le cose contemporaneamente
cook<-cooks.distance(stepwise.mod.final)
plot(cook,ylim = c(0,1), xlab = "Osservazioni", ylab = "Distanza di cooks", main = "Osservazioni secondo Cooks")

Viene confermata la ‘pericolosità’ dell’osservazione 1551. La rimozione di questa osservazione potrebbe cambiare sensibilmente le metriche del modello in analisi. Si può pensare di rimuovere il record e rianalizzare il modello. Oppure si possono fare indagini più approfondite sul perché si tale osservazione.

Passando ai residui.

# Test di normalità dei residui
shapiro.test(residuals(stepwise.mod.final))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(stepwise.mod.final)
## W = 0.9741, p-value < 2.2e-16

# rifiuto l'ipotesi nulla, i residui non sono normali
plot(density(residuals(stepwise.mod.final)))

# Test di omeoschedasticità
lmtest::bptest(stepwise.mod.final)

## 
##  studentized Breusch-Pagan test
## 
## data:  stepwise.mod.final
## BP = 89.798, df = 6, p-value < 2.2e-16

# rifiuto l'ipotesi nulla, la varianza non è costante

# test di autocorrelazione
lmtest::dwtest(stepwise.mod.final)

## 
##  Durbin-Watson test
## 
## data:  stepwise.mod.final
## DW = 1.9542, p-value = 0.126
## alternative hypothesis: true autocorrelation is greater than 0

# Non posso rifiutare l'ipotesi nulla, i residui non sono autocorrelati

Valutiamo infine altre metriche come: MSE e RMSE. Tali parametri sono utili come parametri di confronto. Proviamo a confrontare il modello migliore con il glm del paragrafo precedente.

#Valutiamo anche l'errore quadratico medio
MSE<-function(y_oss,y_prev){
  return(sum((y_oss-y_prev)^2)/length(y_prev))
}

RMSE<-function(y_oss,y_prev){
  return(sqrt(MSE(y_oss, y_prev)))
}

mse_train_stepwise<-MSE(dati$Peso, fitted(stepwise.mod.final))
mse_train_stepwise

## [1] 75189.46

mse_train_gamma<-MSE(dati$Peso, fitted(model_gamma_2))
mse_train_gamma

## [1] 76332.72

rmse_train_stepwise<-RMSE(dati$Peso, fitted(stepwise.mod.final))
rmse_train_stepwise

## [1] 274.207

rmse_train_gamma<-RMSE(dati$Peso, fitted(model_gamma_2))
rmse_train_gamma

## [1] 276.2838

Ricapitolando, dall’analisi della qualità del modello emerge che:

ci sono delle osservazioni che hanno alti valori di leva e outlier. In particolare, considerando la distanza di cooks, l’osservazione 1551 risulta particolarmente critica. Servono ulteriori indagini su queste osservazioni. Si può provare a rimuoverle dal set e ricalcolare le metriche.
NON VENGONO RISPETTATE le ipotesi alla base del modello lineare: i residui sono eteroschedastici e con distribuzione non lineare. Questo incide particolmente sulla qualità del modello. Servono ulteriori dati e ulteriori indagini per trovare un modello migliore
Se si guarda all’errore quadratico medio e si confrontano il miglior modello della procedura stepwise con il miglior modello glm, si conferma ancora una volta la bontà del modello stepwise. Si può utilizzare questa metriche per analizzare le performance del modello sul set di dati di training e su quello di test.

3. Previsioni e Risultati

Una volta validato il modello, lo useremo per fare previsioni pratiche. Ad esempio, potremo stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.

# Previsione 
dato.esempio = data.frame(
  N.gravidanze = 3,
  Gestazione = 39,
  Lunghezza = 450,
  Cranio = 370,
  Sesso = "M",
  Fumatrici = 0
)

cat("Il neonato peserà gr", round(predict(stepwise.mod, newdata = dato.esempio),2))

## Il neonato peserà gr 3206.95

4. Visualizzazioni

Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili. Ad esempio, potremmo visualizzare l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.

Fumatrici.factor <- factor(Fumatrici, levels = c(0,1), labels = c("No","Si"))

ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col=Fumatrici.factor))+ #position = "jitter"
  geom_smooth(aes(x=Gestazione,
                  y=Peso,
                  col=Fumatrici.factor),se=F,method = "lm")+
  labs(title = "Relazione tra peso e settimane di gestazione",
                 x = "Gestazione [sett]", y = "Peso [gr]")

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

Se si mette in relazione il peso previsto con le settimane di gestazione e la variabile Fumatrici, si possono osservare diverse cose:

C’è una relazione lineare tra Peso e Gestazione. La qualità della previsione del peso in funzione delle settimane di gestazione potrebbe essere migliorata a causa della distribuzione della variabile Gestazione: ci sono poche osservazioni relative ai parti prematuri (settimane inferiori a 30-35).
Ci son pocchissime osservazioni di donne fumatrici: solo il 4%. Con questa premessa, e considerando anche che dall’analisi Fumatrici risulta essere un regressore non significativo, nel grafico mostrato si possono osservare due differenti rette di regressione in base all’essere fumatrice. Sembrerebbe che fumare faccia sì che ci sia un minore incremento di Peso all’aumentare delle settimane di gestazione. E’ infatti l’unico regressore ad avere un coefficiente negativo.

ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col=Sesso))+ #position = "jitter"
  geom_smooth(aes(x=Gestazione,
                  y=Peso,
                  col=Sesso),se=F,method = "lm")+
  labs(title = "Relazione tra peso e settimane di gestazione",
                 x = "Gestazione [sett]", y = "Peso [gr]")

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

Analizzando invece la relazione peso previsto con le settimane di gestazione in funzione del sesso, si possono osservare invece:

Equa distribuzione delle osservazioni tra i sessi
Medesimo incremento di Peso all’aumentare delle settimane di gestazione. Comporamento molto simile tra i due sessi: pendenza della retta molto simile, differente intercetta a causa delle diverse medie di peso tra i due sessi.

Infine, grafichiamo la relazione che vige tra il peso e le altre misure antropometriche Lunghezza e Cranio

#car::scatter3d(Peso~Lunghezza+Cranio)
scatterplot3d(Lunghezza, Cranio, Peso, pch=16, color="blue",
              xlab="Lunghezza [mm]", ylab="Cranio [mm]", zlab="Peso [gr]", 
              main="Relazione tra peso e misure antropometriche")

Da notare la presenza della osservazione 1551 emersa dall’analisi del modello. Per il resto, si può osservare una certa relazione lineare tra il peso del neonato e le lunghezze antropometriche.

Conclusioni

Il progetto di previsione del peso neonatale è un’iniziativa fondamentale per Neonatal Health Solutions. Attraverso l’uso di dati clinici dettagliati e strumenti di analisi statistica avanzati, possiamo contribuire a migliorare la qualità della cura prenatale, ridurre i rischi per i neonati e ottimizzare l’efficienza delle risorse ospedaliere. Questo progetto rappresenta un punto di svolta per l’azienda, consentendo non solo un miglioramento della pratica clinica ma anche l’implementazione di politiche sanitarie più informate e proattive.

In conclusione, lo studio preliminare dei dati ci ha permesso di analizzare le variabili a disposizione osservando in particolare la tipologia, come si distribuiscono e come si relazionano. A seguire le principali osservazioni raccolte dal paragrafo 2.1:

Anni.madre presenta due osservazioni con età inferiore a 10 anni. Considerando che biologicamente è molto improbabile (sviluppo del sistema riproduttivo in età compresa tra i 12-15 anni), possiamo considerare i record associati a questo valore di Anni.madre come degli outlier.
Guardando la variabile Gestazione, si osserva la presenza di gravizande premature con osservazioni il cui valore è inferiore alle 30 settimane. Sotto le 30 settimane sono gravidanze premature.
Il punto precendente si collega al peso dei neonati. Ci sono osservazioni per cui il Peso dei neonati assume valori sotto i 1000 gr
Anni.madre pari a 0 e 1 sono chiaramente OUTLIERS.
E l’effetto del fumo? Guardando i subset con Gestazione < 30 e Peso < 1000 si vede che non ci sono Fumatrici. Il fumo non sembra collegato ai parti prematuri.
Si inltre nota poi che fumatori e parti cesarei poco rappresentati con rispettivamene il 4% e il 30%
Mentre c’è una equadistribuzione di maschi e femmine, e di parti tra i 3 ospedali coinvolti nell’analisi.
Dai test effettuati risulta che non si apprezzano differenze di tipo di parto tra i diversi ospedali; la media del peso della popolazione differisce significativamente dalla media del campione; mentre la media della lunghezza del campione è signitificativamente uguale a quella della popolazione; esistono differenze significative nelle misure antropometriche tra i due sessi: sia sulla lunghezza che sulle misure del cranio.

Per quanto riguarda il modello lineare sviluppato dall’analisi ci sono diverse considerazioni da fare:

ci sono delle osservazioni che hanno alti valori di leva e outlier. In particolare, considerando la distanza di cooks, l’osservazione 1551 risulta particolarmente critica. Servono ulteriori indagini su queste osservazioni. Si può provare a rimuoverle dal set e ricalcolare le metriche.
NON VENGONO RISPETTATE le ipotesi alla base del modello lineare: i residui sono eteroschedastici e con distribuzione non lineare. Questo incide particolmente sulla qualità del modello. Servono ulteriori dati e ulteriori indagini per trovare un modello migliore

Nonostante queste premesse abbiamo utilizzato il modello per visualizzare le relazioni che sussistono tra le variabili e, ad esempio, se si mette in relazione il peso previsto con le settimane di gestazione e le variabili Fumatrici e Sesso, si possono osservare emerge che:

C’è una relazione lineare tra Peso e Gestazione. La qualità della previsione del peso in funzione delle settimane di gestazione potrebbe essere migliorata a causa della distribuzione della variabile Gestazione: ci sono poche osservazioni relative ai parti prematuri (settimane inferiori a 30-35).
Ci son pocchissime osservazioni di donne fumatrici: solo il 4%. Con questa premessa, e considerando anche che dall’analisi Fumatrici risulta essere un regressore non significativo, nel grafico mostrato si possono osservare due differenti rette di regressione in base all’essere fumatrice. Sembrerebbe che fumare faccia sì che ci sia un minore incremento di Peso all’aumentare delle settimane di gestazione. E’ infatti l’unico regressore ad avere un coefficiente negativo.
Equa distribuzione delle osservazioni tra i sessi e medesimo incremento di Peso all’aumentare delle settimane di gestazione. Comporamento molto simile tra i due sessi: pendenza della retta molto simile, differente intercetta a causa delle diverse medie di peso tra i due sessi.
Viene confermata una relazione lineare tra il peso del neonato e le lunghezze antropometriche, lunghezza e cranio.

Progetto Stat inferenziale Meringolo

Meringolo

2025-09-14