MODELLO STATISTICO PER LA PREVISIONE DEL PESO NATURALE

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.

1. Raccolta dei Dati & Struttura del Dataset

1.1 Introduzione & descrizione

Il progetto è sviluppato per l’azienda Neonatal Health Solutions e ha come obiettivo la costruzione di un modello statistico predittivo del peso neonatale. Il dataset contiene 2500 osservazioni raccolte tra 3 ospedali, con variabili cliniche, demografiche relative a madre e neonato.

L’analisi si inserisce in un contesto di prevenzione dalle complicazioni neonatali e di ottimizzazione delle risorse ospedaliere.

#Load the data into R (CSV file)
data<-read.csv("neonati.csv",stringsAsFactors = FALSE ) 
  
#First rows of the Dataset
library(knitr)
kable(head(data), caption = "First rows of the Dataset")
First rows of the 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 
attach(data)
#N. observations
N<- dim(data)[1]
cat("Il dataset contiene",N,"osservazioni.\n")
## Il dataset contiene 2500 osservazioni.
#Duplicate check
library(dplyr)
## 
## Caricamento pacchetto: 'dplyr'
## I seguenti oggetti sono mascherati da 'package:stats':
## 
##     filter, lag
## I seguenti oggetti sono mascherati da 'package:base':
## 
##     intersect, setdiff, setequal, union
duplicati<- data%>% group_by_all()%>% filter(n()>1)%>%ungroup()
cat("Il dataset contiene", nrow(duplicati),"duplicati.\n")
## Il dataset contiene 0 duplicati.

Obiettivo: Individuare i fattori che influenzano significativamente il peso neonatale.

In particolare:

  • Valutare l’impatto del fumo materno.

  • Analizzare l’effetto della durata della gestazione (nascite premature vs a termine).

  • Confrontare le differenze tra sessi , tipi di parto, e ospedali.

  • Integrare le variabili antropometriche (lunghezza, cranio) e materne (età, gravidanze) in un modello di regressione.

1.2 Analisi Variabili del Dataset

Di seguito otteniamo le variabili presenti nel Dataset:

#Variables Names
kable(data.frame(Variables=names(data)), caption = "Variables Names\n")
Variables Names
Variables
Anni.madre
N.gravidanze
Fumatrici
Gestazione
Peso
Lunghezza
Cranio
Tipo.parto
Ospedale
Sesso

Organizziamo le variabili estratte dal Dataset all’interno di una tabella, riportando per ogni variabile:

[ Variabile | Tipo di Variabile | Descrizione ]

variabili_df <- data.frame(
  Variabile = c("Anni.madre","N.gravidanze","Fumatrici","Gestazione","Peso",
                "Lunghezza","Cranio","Tipo.parto","Ospedale","Sesso"),
  Tipo = c("Quantitativa continua","Quantitativa discreta","Qualitativa dicotomica",
           "Quantitativa discreta","Quantitativa continua","Quantitativa continua",
           "Quantitativa continua","Qualitativa nominale","Qualitativa nominale",
           "Qualitativa dicotomica"),
  Descrizione = c("Età della madre (anni)",
                  "Numero di gravidanze",
                  "Stato fumo materno (0/1)",
                  "Durata gestazione (settimane)",
                  "Peso alla nascita (grammi)",
                  "Lunghezza neonato (millimetri)",
                  "Diametro craniale (millimetri)",
                  "Tipo di parto (Nat/Ces)",
                  "Ospedale di nascita (osp1/osp2/osp3)",
                  "Sesso neonato (M/F)")
)
kable(variabili_df, caption = "Variable, Type & Description")
Variable, Type & Description
Variabile Tipo Descrizione
Anni.madre Quantitativa continua Età della madre (anni)
N.gravidanze Quantitativa discreta Numero di gravidanze
Fumatrici Qualitativa dicotomica Stato fumo materno (0/1)
Gestazione Quantitativa discreta Durata gestazione (settimane)
Peso Quantitativa continua Peso alla nascita (grammi)
Lunghezza Quantitativa continua Lunghezza neonato (millimetri)
Cranio Quantitativa continua Diametro craniale (millimetri)
Tipo.parto Qualitativa nominale Tipo di parto (Nat/Ces)
Ospedale Qualitativa nominale Ospedale di nascita (osp1/osp2/osp3)
Sesso Qualitativa dicotomica Sesso neonato (M/F)

NB: Le variabili qualitative saranno trattate come fattori nelle analisi inferenziali e nei modelli di regressione.

2. Analisi & Modellazione

2.1 Analisi Preliminare

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

2.1.1 Classificazione Variabili

Nella sezione precedente abbiamo elencato le variabili e le tipologia di variabile in esame. Costruiamo una tabella che consenta di evidenziare la distinzione tra le tipologie di variabili. La tabella distingue chiaramente guiderà la scelta dei test statistici e dei modelli nella fasi successive.

variabili_type_df <- data.frame(
  Variabile_Quantitativa = c(
    "Anni.madre (continua)",
    "N.gravidanze (discreta)",
    "Gestazione (discreta, settimane)",
    "Peso (continua, grammi)",
    "Lunghezza (continua, mm)",
    "Cranio (continua, mm)"
  ),
  Variabile_Qualitativa = c(
    "Fumatrici (dicotomica: 0/1)",
    "Tipo.parto (nominale: Nat/Ces)",
    "Ospedale (nominale: osp1/osp2/osp3)",
    "Sesso (dicotomica: M/F)",
    "",  
    ""   
  )
)
kable(variabili_type_df, caption = "Classification of Variables")
Classification of Variables
Variabile_Quantitativa Variabile_Qualitativa
Anni.madre (continua) Fumatrici (dicotomica: 0/1)
N.gravidanze (discreta) Tipo.parto (nominale: Nat/Ces)
Gestazione (discreta, settimane) Ospedale (nominale: osp1/osp2/osp3)
Peso (continua, grammi) Sesso (dicotomica: M/F)
Lunghezza (continua, mm)
Cranio (continua, mm)

A questo punto, come richiesto nell’enunciato della sezione, esploriamo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione.

2.1.2 Variabili Quantitative

Statistiche Descrittive di Base

Per le variabili quantitative (Anni.madre, N.gravidanze, Gestazione, Peso, Lunghezza, Cranio) calcolo media, mediana, deviazione standard, min/max e quartili.

kable(summary(data[, c("Anni.madre","N.gravidanze","Gestazione","Peso","Lunghezza","Cranio")]))
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

Le medie e mediane risultano vicine per molte variabili, suggerendo distribuzioni abbastanza simmetriche. Tuttavia, il peso e la lunghezza mostrano range ampi, con possibili valori estremi che meritano attenzione.

Indici di Variabilità e Forma

library(moments)
library(dplyr)
#Define the quantitative variables
var_quant <- c("Anni.madre","N.gravidanze","Gestazione","Peso","Lunghezza","Cranio")

#Function to calculate the requested indices
calc_stas<- function(x){
  c(
    Varianza = var(x, na.rm = TRUE),
    Dev_Std = sd(x, na.rm = TRUE),
    Range = diff(range(x, na.rm = TRUE)),
    IQR = IQR(x, na.rm = TRUE),
    Asimmetria = skewness(x, na.rm = TRUE),
    Curtosi = kurtosis(x, na.rm = TRUE)
  )
}
#Apply the function to the qualitative variables
stats_df <- sapply(data[var_quant], calc_stas)%>%
  round(2) %>%
  as.data.frame()
stats_df<- tibble::rownames_to_column(stats_df,var= "Measure")

#Final table
kable(stats_df, caption = "Indices by size (rows) and variable (columns)")
Indices by size (rows) and variable (columns)
Measure Anni.madre N.gravidanze Gestazione Peso Lunghezza Cranio
Varianza 27.81 1.64 3.49 275665.68 692.67 269.79
Dev_Std 5.27 1.28 1.87 525.04 26.32 16.43
Range 46.00 12.00 18.00 4100.00 255.00 155.00
IQR 7.00 1.00 2.00 630.00 30.00 20.00
Asimmetria 0.04 2.51 -2.07 -0.65 -1.51 -0.79
Curtosi 3.38 13.99 11.26 5.03 9.49 5.95

  • Il peso presenta asimmetria positiva e curtosi > 3, indicando code più pesanti e valori estremi.

  • L’età materna mostra una distribuzione più concentrata, ma con outlier (madri molto giovani o molto anziane).

  • Le variabili antropometriche (lunghezza, cranio) hanno distribuzioni più regolari, ma con alcuni valori anomali.

Distribuzione Grafiche & Outlier

Riporto una serie di distribuzioni grafiche per rappresentare la tendenza e la variabilità delle variabili:

  • Istogrammi per le var quantitative: evidenziano la forma della distribuzione (asimmetria, code, multimodalità).

  • Boxplot : utili per individuare outlier.

  1. Peso
# Histogram for Peso
hist(Peso, main="Distribuzione Peso Neonati", xlab="Peso (g)", col="lightblue", border="black")
avg_peso<- mean(Peso, na.rm = TRUE)
Q1 <- quantile(Peso, 0.25)
Q3 <- quantile(Peso, 0.75)
abline(v = avg_peso, col = "red", lwd = 2,lty=2)
abline(v = Q1, col = "green", lwd = 1,lty=2)
abline(v = Q3, col = "green", lwd = 1,lty=2)

legend("topleft",legend= c(paste("µ = ",round(avg_peso,2),"g"),
                           paste("Q1 = ", round(Q1,0),"g"),
                           paste("Q3 = ", round(Q3,0),"g")), 
       col= c(2,"green", "green"), lwd= 2, lty=c(2,3,3),bty= "n") 

# Boxplot for Peso
boxplot(Peso, main="Boxplot Peso Neonati", ylab="Peso (g)", col="lightgreen")

# Outlier detection for Peso
IQR <- IQR(Peso, na.rm = TRUE)
outliers <- Peso[Peso < (Q1 - 1.5*IQR) | Peso > (Q3 + 1.5*IQR)]
cat("la variabile Peso contiene", length(outliers), "valori anomali")
## la variabile Peso contiene 69 valori anomali
outlier_matrix <- matrix(outliers, ncol = 10, byrow = TRUE)
## Warning in matrix(outliers, ncol = 10, byrow = TRUE): data length [69] is not a
## sub-multiple or multiple of the number of rows [7]
kable(as.data.frame(t(outlier_matrix)), caption = "elenco outliers")
elenco outliers
V1 V2 V3 V4 V5 V6 V7
1370 1550 830 1840 4760 4650 2000
1340 1410 1615 1620 1800 4720 4690
4680 1900 1960 1280 1430 1780 1580
1500 1720 1770 1280 1950 1180 1170
1850 1980 1750 4810 1970 1890 4720
1560 1390 1170 2040 900 1140 1690
1280 1450 2040 4620 1780 1600 980
1750 1970 1980 1500 4580 1300 930
4600 1190 4600 2000 4930 930 1730
1285 2000 4900 990 4700 1750 1370

Commento: La maggior parte dei neonati pesa tra i 2990 e i 3600 g. Sono presenti outlier sotto i 2000 g (basso peso che è dovuto probabilmente a prematurità) e sopra i 4500 g. Alcuni casi estremi (<1000 g) sono clinicamente implausibili e potrebbero essere errori di registrazione.

  1. Età della Madre
# Histogram for Anni.Madre
hist(Anni.madre, main="Distribuzione Età Madri", xlab="Età (anni)", col="lightblue", border="black")
avg_eta<- mean(Anni.madre, na.rm = TRUE)
Q1 <- quantile(Anni.madre, 0.25)
Q3 <- quantile(Anni.madre, 0.75)
abline(v = avg_eta, col = "red", lwd = 2,lty=2)
abline(v = Q1, col = "green", lwd = 1,lty=2)
abline(v = Q3, col = "green", lwd = 1,lty=2)

legend("topleft",legend= c(paste("µ = ",round(avg_eta,2),"anni"),
                           paste("Q1 = ", round(Q1,0),"anni"),
                           paste("Q3 = ", round(Q3,0),"anni")), 
       col= c(2,"green", "green"), lwd= 2, lty=c(2,3,3),bty= "n") 

# Boxplot for Anni.Madre
boxplot(Anni.madre, main="Boxplot Età Madri", ylab="Età", col="lightgreen")

# Outlier detection for Età Madri
IQR <- IQR(Anni.madre, na.rm = TRUE)
outliers <- Anni.madre[Anni.madre < (Q1 - 1.5*IQR) | Anni.madre > (Q3 + 1.5*IQR)]
cat("la variabile Anni.madre contiene", length(outliers), "valori anomali")
## la variabile Anni.madre contiene 13 valori anomali
outlier_matrix <- matrix(outliers, nrow = 1, byrow = TRUE)
kable(as.data.frame(t(outlier_matrix)), caption = "elenco outliers")
elenco outliers
V1
13
45
43
44
44
43
14
46
1
0
14
44
44

Commento: Distribuzione centrata sui 28-32 anni. Oultier: madri molto giovani (<18 anni) e madri oltre i 40 anni. Tra i 13 valori anomali, sono presenti 2 osservazioni con età <10anni (0,1), i quali potrebbero essere errori di registrazione (è impossibile che l’età della madre sia 0/1 anno).

  1. Lunghezza del Neonato
# Histogram for Lunghezza
hist(Lunghezza, main="Lunghezza neonato", xlab="Lunghezza (cm)", col="lightblue", border="black")
avg_lung<- mean(Lunghezza, na.rm = TRUE)
Q1 <- quantile(Lunghezza, 0.25)
Q3 <- quantile(Lunghezza, 0.75)
abline(v = avg_lung, col = "red", lwd = 2,lty=2)
abline(v = Q1, col = "green", lwd = 1,lty=2)
abline(v = Q3, col = "green", lwd = 1,lty=2)

legend("topleft",legend= c(paste("µ = ",round(avg_lung,2),"cm"),
                           paste("Q1 = ", round(Q1,0),"cm"),
                           paste("Q3 = ", round(Q3,0),"cm")), 
       col= c(2,"green", "green"), lwd= 2, lty=c(2,3,3),bty= "n")

# Boxplot for Lunghezza
boxplot(Lunghezza, main="Boxplot Lunghezza neonato", ylab="Lunghezza", col="lightgreen")

# Outlier detection for Lunghezza
IQR <- Q3 - Q1
outliers <- Lunghezza[Lunghezza < (Q1 - 1.5*IQR) | Lunghezza > (Q3 + 1.5*IQR)]
cat("la variabile Lunghezza contiene", length(outliers), "valori anomali")
## la variabile Lunghezza contiene 59 valori anomali
outlier_matrix <- matrix(outliers, ncol = 10, byrow = FALSE)
## Warning in matrix(outliers, ncol = 10, byrow = FALSE): data length [59] is not
## a sub-multiple or multiple of the number of rows [6]
kable(as.data.frame(t(outlier_matrix)), caption = "elenco outliers")
elenco outliers
V1 V2 V3 V4 V5 V6
390 400 410 405 420 420
360 430 400 410 560 380
420 390 405 360 430 310
390 430 390 420 410 420
370 430 430 410 385 390
560 315 420 340 410 430
430 380 325 420 420 410
400 430 355 370 410 380
355 410 425 400 370 430
565 405 320 345 430 390

Commento: La maggior parte dei neonati è lunga tra i 480 e i 510 cm. Sono presenti outlier sotto i 420 cm e sopra i 560 cm. Tutti i casi elencati (ben 59) non vanno eliminati ma interpretati con cautela.

  1. Diametro del Cranio
# Histogram for Cranio
hist(Cranio, main="Diametro Cranio", xlab="Diametro (cm)", col="lightblue", border="black")
avg_cranio<- mean(Cranio, na.rm = TRUE)
abline(v = avg_cranio, col = "red", lwd = 2,lty=2)
Q1 <- quantile(Cranio, 0.25)
Q3 <- quantile(Cranio, 0.75)

abline(v = Q1, col = "green", lwd = 1,lty=2)
abline(v = Q3, col = "green", lwd = 1,lty=2)

legend("topleft",legend= c(paste("µ = ",round(avg_cranio,2),"cm"),
                           paste("Q1 = ", round(Q1,0),"cm"),
                           paste("Q3 = ", round(Q3,0),"cm")), 
       col= c(2,"green", "green"), lwd= 2, lty=c(2,3,3),bty= "n")

# Boxplot for Cranio
boxplot(Cranio, main="Boxplot Diametro Cranio", ylab="Diametro (cm)", col="lightgreen")

# Outlier detection for Lunghezza
IQR <- IQR(Cranio, na.rm = TRUE)
outliers <- Cranio[Cranio < (Q1 - 1.5*IQR) | Cranio > (Q3 + 1.5*IQR)]
cat("la variabile Cranio contiene", length(outliers), "valori anomali")
## la variabile Cranio contiene 48 valori anomali
outlier_matrix <- matrix(outliers, ncol = 10, byrow = FALSE)
## Warning in matrix(outliers, ncol = 10, byrow = FALSE): data length [48] is not
## a sub-multiple or multiple of the number of rows [5]
kable(as.data.frame(t(outlier_matrix)), caption = "elenco outliers")
elenco outliers
V1 V2 V3 V4 V5
298 382 287 273 285
280 390 384 295 276
382 274 289 295 386
390 277 280 272 254
297 385 295 275 266
383 390 292 292 293
278 253 277 390 381
270 267 290 276 235
294 299 298 290 273
290 265 245 298 382

Commento: La maggior parte dei neonati è ha il diametro del cranio tra i 330 e i 350 cm. Sono presenti outlier sotto i 290 cm e sopra i 560 cm. Tutti i casi elencati (ben 48) non vanno eliminati ma interpretati con cautela.

  1. Settimane di Gestazione
# Histogram for Gestazione
hist(Gestazione, main="Settimane di Gestazione", xlab="Weeks", col="lightblue", border="black")
avg_gest<- mean(Gestazione, na.rm = TRUE)
abline(v = avg_gest, col = "red", lwd = 2,lty=2)
Q1 <- quantile(Gestazione, 0.25)
Q3 <- quantile(Gestazione, 0.75)

abline(v = Q1, col = "green", lwd = 1,lty=2)
abline(v = Q3, col = "green", lwd = 1,lty=2)

legend("topleft",legend= c(paste("µ = ",round(avg_gest,0),"weeks"),
                           paste("Q1 = ", round(Q1,0),"weeks"),
                           paste("Q3 = ", round(Q3,0),"weeks")), 
       col= c(2,"green", "green"), lwd= 2, lty=c(2,3,3),bty= "n")

# Boxplot for Gestazione
boxplot(Gestazione, main="Boxplot Settimane di Gestazione", ylab="Weeks", col="lightgreen")

# Outlier detection for Gestazione
IQR <- IQR(Gestazione, na.rm = TRUE)
outliers <- Gestazione[Gestazione < (Q1 - 1.5*IQR) | Gestazione > (Q3 + 1.5*IQR)]
cat("la variabile Gestazione contiene", length(outliers), "valori anomali")
## la variabile Gestazione contiene 67 valori anomali
outlier_matrix <- matrix(outliers, ncol = 10, byrow = FALSE)
## Warning in matrix(outliers, ncol = 10, byrow = FALSE): data length [67] is not
## a sub-multiple or multiple of the number of rows [7]
kable(as.data.frame(t(outlier_matrix)), caption = "elenco outliers")
elenco outliers
V1 V2 V3 V4 V5 V6 V7
34 33 34 30 34 34 31
34 33 28 32 28 34 34
32 33 33 31 32 34 33
33 33 33 29 34 28 32
31 33 34 34 33 30 33
34 32 29 34 29 33 31
31 32 25 32 33 34 34
33 31 32 27 33 30 28
30 32 33 34 30 33 31
27 26 31 33 34 33 34

Commento: Distribuzione centrata sui 38-40 settimane. Outlier: tutti i valori anomali sono sotto le 35 settimane. Tra i 67 valori anomali sono presenti casi di nascite premature (< 30 settimane): potrebbero non risultare errori, ma casi clinici reali che spiegano i valori bassi di peso osservati.

2.1.3 Variabili Qualitative - Distribuzione di Frequenza

# Frequency Function
freq_table <- function(x, N = length(x)) {
  tb <- table(x)
  data.frame(
    cat = names(tb),
    ni = as.integer(tb),
    fi = round(as.numeric(tb) / N, 2),
    check.names = FALSE
  )
}

# Frequency Distribution
kable(freq_table(Fumatrici,length(Fumatrici)), caption = "Distribuzione Fumatrici")
Distribuzione Fumatrici
cat ni fi
0 2396 0.96
1 104 0.04 
kable(freq_table(Tipo.parto,length(Tipo.parto)), caption = "Distribuzione Tipo Parto")
Distribuzione Tipo Parto
cat ni fi
Ces 728 0.29
Nat 1772 0.71 
kable(freq_table(Ospedale,length(Ospedale)), caption = "Distribuzione Ospedali")
Distribuzione Ospedali
cat ni fi
osp1 816 0.33
osp2 849 0.34
osp3 835 0.33 
kable(freq_table(Sesso,length(Sesso)), caption = "Distribuzione Sesso")
Distribuzione Sesso
cat ni fi
F 1256 0.5
M 1244 0.5

Commento:

  • Solo il 4% delle madri è fumatrice: campione totalmente sbilanciato.

  • I parti naturali sono più frequenti (71%) rispetto a quelli cesarei (29%)

  • La distribuzione è equilibrata per sesso e ospedali.

2.1.4 Formulazione delle ipotesi da testare

Di seguito formuliamo le ipotesi in modo operativo e indichiamo i test statistici appropriati:

  1. In alcuni ospedali si fanno più parti cesarei

  • Domanda statistica: Esiste una associazione tra ospedale e tipologia di parto?

  • Variabili: Ospedale (osp1, osp2, osp3), Tipo.parto (Nat/Ces)

Dato che le variabili da sottoporre a test statistico sono categoriche uso Test chi-quadro su tabella di contingenza dove l’interpretazione attesa è :

  • H0 : Ospedale e tipo di parto sono indipendenti (proporzione di cesarei uguale nei 3 ospedali)

  • H1 : Ospedale e tipo di parto sono associati (almeno un ospedale ha proporzione diversa)

#Tab COntingenza ospedale x tipo di parto
tab_parto_osp<- table( Tipo.parto, Ospedale)
knitr::kable(round(tab_parto_osp, 3), caption = "Tab. di contingenza per ospedale")
Tab. di contingenza per ospedale
osp1 osp2 osp3
Ces 242 254 232
Nat 574 595 603 
## Proporzioni per ospedale (per riga)
prop_per_osp <- prop.table(tab_parto_osp, margin = 1)
knitr::kable(round(prop_per_osp, 3), caption = "Proporzioni per ospedale")
Proporzioni per ospedale
osp1 osp2 osp3
Ces 0.332 0.349 0.319
Nat 0.324 0.336 0.340 
#Visualizzazione Baloon Plot
ggpubr::ggballoonplot(data=as.data.frame(tab_parto_osp), fill = "blue")

Commento: Ho generato un Baloon Plot (=grafico i cui punti d’intersezione tra 2 variabili corrisponde a un palloncino, la cui grandezza dipende dalla frequenza assoluta osservata).

Secondo le regole relative all’Indipendenza e Associazione, Il balloon plot mostra differenze visive nelle frequenze (palloncini con diverse dimensioni e senza un pattern preciso), ma la conferma statistica richiede il test del chi‑quadro.

#Calcolo delle Frequenze attese
attese<- tab_parto_osp

for(i in 1:nrow(tab_parto_osp)){
  for(j in 1:ncol(tab_parto_osp)){
    attese[i,j] <- margin.table(tab_parto_osp,1)[i] * margin.table(tab_parto_osp,2)[j] / margin.table(tab_parto_osp)
  }
}
#Frequence osservate
osservate<- tab_parto_osp
#Calcolo manuale del Chi-Quadro
x_quadro_manuale<- sum((osservate-attese)^2/attese)
x_quadro_manuale
## [1] 1.097202

Vado a simulare graficamente la Distribuzione del Chi-Quadro

plot(density(rchisq(1000000,df = 2)), xlim=c(0,30))
abline(v=qchisq(0.95,2),col=2)
points(x_quadro_manuale,0,cex=3,col=4)

Commento: Il calcolo manuale e la simulazione grafica confermano che la statistica osservata rientra nella zona di accettazione e quindi per ora H0 non può venir rifiutata. Per avvalorare le nostre ipotesi usiamo il Test d’indipendenza con funzione integrata con lo scopo di calcolare statistica, p-value, G.d.L.

test_indipendenza<- chisq.test(tab_parto_osp)
test_indipendenza
## 
##  Pearson's Chi-squared test
## 
## data:  tab_parto_osp
## X-squared = 1.0972, df = 2, p-value = 0.5778

Conclusione: Dal calcolo si denota un p-value = 0.578 > alpha (liv. di significatività) = 0.05 -> non rifiutiamo H0. Non emergono differenze significative nelle proporzioni di cesarei tra ospedali.

  1. La media del peso e della lunghezza del campione = popolazione

  • Domanda Statistica: Le medie di Peso e Lunghezza osservate nel campione sono coerenti con i valori medi della popolazione?

  • Variabili: Peso (g), Lunghezza (mm)

Dato che vado a fare un confronto tra medie di gruppi utilizzo il test T, dove ho come ipotesi:

  • H0 : La media campionaria è uguale alla media della popolazione (µ campione = µ popolazione)

  • H1 : La media campionaria è diversa da quella della popolazione ((µ campione ≠ µ popolazione)

  1. Peso

Andiamo a verificare l’ipotesi nulla prima per la variabile Peso e incomincio effettuando il test di normalità della variabile peso per verificare la sua non normalità e in fase successiva passiamo al t test per confronto popolazione campione.

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

Commento: p-value = 0 -> i residui non seguono una distribuzione normale. Tuttavia dato che abbiamo campioni molto grandi anche anche piccole deviazioni diventano significative. Per questo motivo, grazie al Teorema del Limite Centrale, il t‑test rimane applicabile.

NB: Come valori di riferimento prendo come media del peso 3,3 kg (3300 g). Fonte : Manuali MSD – Sezione Pediatria, Parametri di accrescimento nei neonati (MSD Manuals, 2025).

# Valori di riferimento della popolazione
mu_peso_pop <- 3300

# Statistiche campionarie
avg_peso <- mean(Peso, na.rm = TRUE)
sd_peso <- sd(Peso, na.rm = TRUE)
n_peso <- sum(!is.na(Peso))

# Statistica t test
t_peso <- t.test(Peso, mu = mu_peso_pop)
t_peso
## 
##  One Sample t-test
## 
## data:  Peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.490 3304.672
## sample estimates:
## mean of x 
##  3284.081
#Quantili
cat("I valori di sogllia sono", qt(0.025,2499),", ",qt(0.975,2499))
## I valori di sogllia sono -1.960914 ,  1.960914

Commento: Eseguo il t test a un campione della variabile Peso, confrontando la media osservata µ campione con µ popolazione (3250 g). I risultati ottenuti sono:

  • t value= -1.5 rientra nei valori di soglia = [ -1.96 ; +1.94 ]

  • p-value = 0.1296 > alpha = 0.05 : l’ipotesi nulla è non statisticamente rilevante, non rifiuto H0.

  • IC = [ 3263.5 3304 ] : sapendo che H0 = 3300 rientra per poco in IC, non rifiuto H0.

Fatte queste considerazioni non posso rifiutare H0e quindi la variabile Peso ha µ campione = µ popolazione. Volendo possiamo visualizza il grafico della distribuzione del peso e verifichiamo se il punto del t statistico si trova all’interno della zona d’accettazione.

plot(density(rt(100000,2499)),xlim=c(-4,4))
abline(v=qt(0.025,2499), col=2)
abline(v=qt(0.975,2499),col=2)
points(t_peso$statistic,0,cex=3,col=2)

Risultato : Come Volevasi Dimostrare, il punto si trova all’interno della zona d’accettazione-> H0 non si rifiuta.

  1. Lunghezza

Ora verifico l’ipotesi nulla prima per la variabile Lunghezza e per prima azione effettuo il test di normalità della variabile peso per verificare la sua non normalità e in fase successiva passp al t test per confronto popolazione campione.

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

Commento: p-value = 0 -> i residui non seguono una distribuzione normale. Anche qui grazie alla numerosità del campione e al teorema del limite centrale il test-t rimane applicabile.

NB: In questo caso prendo come valore di riferimento la media della Lunghezza 50 cm (500 mm). Fonte : Manuali MSD – Sezione Pediatria, Parametri di accrescimento nei neonati (MSD Manuals, 2025).

# Valori di riferimento della popolazione
mu_peso_pop <- 500

# Statistiche campionarie
avg_lung <- mean(Lunghezza, na.rm = TRUE)
sd_lung <- sd(Lunghezza, na.rm = TRUE)
n_lung<- sum(!is.na(Lunghezza))

# Statistica t test
t_Lunghezza<- t.test(Lunghezza, mu = mu_peso_pop)
t_Lunghezza
## 
##  One Sample t-test
## 
## data:  Lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6598 495.7242
## sample estimates:
## mean of x 
##   494.692
#Quantili
cat("I valori di sogllia sono", qt(0.025,2499),", ",qt(0.975,2499))
## I valori di sogllia sono -1.960914 ,  1.960914

Commento: Nel caso della variabile Lunghezza i risultati ottenuti sono:

  • t value= -10.084 non rientra nei valori di soglia = [ -1.96 ; +1.94 ]

  • p-value = 0 < alpha = 0.05 : l’ipotesi nulla è statisticamente non rilevante, rifiuto H0.

  • IC = [ 493.66 495.72 ] : sapendo che H0 = 500, esso non rientra per poco in IC, quindi rifiuto H0.

Fatte queste considerazioni posso rifiutare H0 e quindi la variabile Lunghezza ha µ campione ≠ µ popolazione. Anche qui posso visualizzare e accertarmi tramite grafico se il punto si trova all’interno o meno della zona d’accettazione.

plot(density(rt(100000,2499)),xlim=c(-11,11))
abline(v=qt(0.025,2499), col=2)
abline(v=qt(0.975,2499),col=2)
points(t_Lunghezza$statistic,0,cex=3,col=2)

Risultato : Come Volevasi Dimostrare, il punto si totalmente al di fuori della zona d’accettazione-> H0 si rifiuta (La media campionaria è diversa da quella della popolazione).

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

  • Domanda Statistica: Le misure antropometriche (Peso, Lunghezza, Cranio) differiscono significativamente tra maschi e femmine?

  • Variabili: Quantitative : Peso (g), Lunghezza (mm), Cranio (mm) , Qualitativa : Sesso (M/F)

Dato che vado a fare un confronto tra multipli utilizzo Wilcoxon test (non parametrico), dove le ipotesi sono:

  • H0: Le medie delle misure antropometriche sono uguali tra maschi e femmine.

  • H1: Le medie delle misure antropometriche sono diverse tra maschi e femmine.

  1. Lunghezza

Per la variabile Lunghezza sappiamo già che non ha una distribuzione normale. Proprio per questa assunzione viene usato il test di Wilcoxon .

summary(Lunghezza[Sesso=='M'], caption = "Lenght per Men")
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   320.0   490.0   500.0   499.7   515.0   560.0
summary(Lunghezza[Sesso=='F'], caption = "Lenght per Women")
##    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]')

Commento: Dal boxplot si osserva che i neonati maschi tendono ad avere una lunghezza leggermente superiore rispetto alle femmine. La varianza è minore nei maschi, con una distribuzione più concentrata.

Appunto facciamo il wilcox test,

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

Commento: Nel caso della variabile Lunghezza i risultati ottenuti sono:

  • p-value = 0 < alpha = 0.05 : l’ipotesi nulla è statisticamente non rilevante, rifiuto H0.
  • Risultato: la lunghezza neonatale differisce significativamente tra i due sessi.
  1. Peso

Anche per la variabile Peso possiamo utilizzare il test di wilcoxon perché abbiamo dimostrato in precedenza la mancanza di una distribuzione normale.

summary(Peso[Sesso=='M'], caption = "Weight per Men")
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     980    3150    3430    3408    3720    4810
summary(Peso[Sesso=='F'], caption = "Weight per Women")
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     830    2900    3160    3161    3470    4930
boxplot(Peso~Sesso, ylab = 'Peso [g]')

Commento: Il peso medio dei neonati maschi è superiore a quello delle femmine (≈3500 g vs ≈3200 g). Le distribuzioni sono simili in varianza.

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

Commento: Nel caso della variabile Peso, i risultati ottenuti sono:

  • p-value = 0 < alpha = 0.05 : l’ipotesi nulla è statisticamente non rilevante, rifiuto H0.
  • Risultato: il peso neonatale differisce significativamente tra i due sessi.
  1. Cranio

Invece per la variabile Cranio non abbiamo idea se ha una distribuzione normale. Quindi applichiamo il test di shapiro.

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

Commento: p-value=0 -> Anche la variabile peso non ha distribuzione normale.

summary(Cranio[Sesso=='M'], caption = "Weight per Men")
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   265.0   334.0   343.0   342.4   352.0   390.0
summary(Cranio[Sesso=='F'], caption = "Weight per Women")
##    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]')

Commento: Dal grafico dei boxplot si desume che i maschi presentano una media leggermente superiore (≈342 mm vs ≈337 mm), con varianze simili.

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

Commento: Nel caso della variabile Cranio, i risultati ottenuti sono:

  • p-value = 0 < alpha = 0.05 : l’ipotesi nulla è statisticamente non rilevante, rifiuto H0.
  • Risultato: il diametro cranico differisce significativamente tra i due sessi.

Conclusioni: Per tutte e tre le variabili antropometriche (Peso, Lunghezza, Cranio), il test di Wilcoxon ha portato al rifiuto dell’ipotesi nulla. Possiamo quindi affermare che esistono differenze statisticamente significative tra maschi e femmine nelle misure antropometriche alla nascita.

2.2 Creazione del Modello di Regressione

Obiettivo: Costruire un modello di regressione lineare multipla per spiegare la variabilità del Peso neonatale in funzione delle variabili disponibili (antropometriche, materne e contestuali).

2.2.1 Scelta variabili del modello

  • Variabile Risposta: Peso (g) -> variabile quantitativa continua che rappresenta l’outcome di principale interesse.

  • Variabili regressori: Quantitative ( Anni.madre, N.gravidanze, Gestazione, Lunghezza, Cranio ) + Qualitative ( Fumatrici, Tipo.parto, Ospedale, Sesso )

Ipotesi:

  • H0: tutti i coefficienti = 0 (nessuna influenza sul peso).

  • H1: almeno un coefficiente è diverso da 0 (almeno una variabile influenza il peso).

2.2.2 Correlazione tra variabili

Prima di costruire il modello di regressione lineare multipla, con il peso come variabile risposta, andiamo a studiare la linearità del modello. Per fare ciò mostro una matrice di scatterplot per correlazioni.

Ciò mi permette di verificare:

  • Multicollinearità: se le variabili indipendenti sono molto correlate.

  • Scelta Variabili: posso decidere quali mantenere o scartare.

  • Diagnosi Preliminare: Individuo gli outlier o pattern strani che potrebbero influenzare il modello.

NB: Nel grafico analizzo solo le variabili quantitative dell’analisi delle correlazioni.

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(data[sapply(data, is.numeric)], lower.panel=panel.cor, upper.panel=panel.smooth)

Commento: Nella matrtice di scatterplot si può osservare un’alta correlazione con la variabile Peso per le variabili Lunghezza, Cranio, Gestazione

Per verificare le variabili qualitative, utilizzo dei boxplot.

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

Commento: Si denota un peso medio maggiore per i maschi rispetto al peso medio femminile ma mantengono la stessa variabilità. Inoltre risultano un numero maggiori di outliers verso il basso rispetto all’alto.

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

Risultati:

  • t = -12,1 : t risulta più lontano da 0 e quindi rifiuto H0.

  • p-value = 0 : rifiuto H0 -> rifiuto H0

  • IC = [ -287,1 ; -207,06 ] -> la Var stimata è [-287,1 ; -207,06] più pesante nei maschi.

par(mfrow=c(1,2))
boxplot(Peso, ylab='Fumatori', xlab='Campione')
boxplot(Peso~Fumatrici, ylab = 'Fumatori')

Commento: Si denota un leggero numero medio di non fumatori superiore rispetto ai fumatori con a sua volta una variabilità superiore dei non fumatori. Inoltre risultano un numero maggiori di outliers da parte dei non fumatori.

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

Risultati:

  • t = 1,03 : t risulta più vicino da 0 e quindi non rifiuto H0.

  • p-value = 0,3 > alpha = 0.05 : non rifiuto H0

  • IC = [ -45,61 ; 145,23 ] -> la Var stimata è [-287,1 ; -207,06] più pesante nei maschi.

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

Commento: Si denota tutte gli ospedali presentano +o - lo stesso numero medio di nascite (3200). Inoltre risultano per tutti gli ospedali un numero elevato di outliers verso il basso.

2.2.3 Modello di regressione lineare multipla

Costruiamo il modello di regressione lineare multipla che comprenda tutte le variabili.

mod_full <- lm(Peso ~ Anni.madre + N.gravidanze + Gestazione +
                 Lunghezza + Cranio + Fumatrici + Tipo.parto +
                 Ospedale + Sesso, data = data)

summary(mod_full)
## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione + 
##     Lunghezza + Cranio + Fumatrici + Tipo.parto + Ospedale + 
##     Sesso, data = data)
## 
## 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 *  
## 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 ***
## Fumatrici       -30.1631    27.5386  -1.095   0.2735    
## 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

Commento:

  • RMSE ≈ 274 g → errore medio di previsione del peso neonatale

  • R^2 = 0.73 : Modello spiega il 73% della variabilità.

  • p-value = 0 : il modello è significativo.

Dal summary si nota che le variabili più significative sono Gestazione, Lunghezza e Cranio, mentre i regressori non significativi sono Anni.madre, Fumatrici, Ospedale e in seconda fascia N.gravidanze e Tipo.parto. A questo punto eliminiamo tutte le variabili non significative in prima istanza per verificare se il nuovo modello riesce a spiegare meglio la variabilità.

mod2<- lm(Peso ~ N.gravidanze + Gestazione +
                 Lunghezza + Cranio + Tipo.parto + Sesso , data = data)

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

Commento:

  • RMSE ≈ 274 g → errore medio di previsione del peso neonatale

  • R^2 = 0.73 : Modello spiega il 73% della variabilità.

  • p-value = 0 : il modello è significativo.

Nonostante eliminazione delle variabili in prima istanza, la variabilità del modello non è migliorata ma il modello è più semplice rispetto al modello precedente. A pari variabilità si sceglie sempre il modello più semplice. A questo punto valutiamo ancora di migliorare il modello eliminando le variabili non significative di seconda fascia come Tipo.parto e N.gravidanze.

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

Commento:

  • RMSE ≈ 274 g → errore medio di previsione del peso neonatale

  • R^2 = 0.73 : Modello spiega il 73% della variabilità.

  • p-value = 0 : il modello è significativo.

Anche nonostante l’eliminazione delle variabili in seconda fascia, la variabilità del modello non è migliorata ma anche in questo caso il modello è migliore del precedente. Quindi come visto precedentemente scegliamo questo modello perché è più semplice.

Conclusione: ll modello finale (mod3) con Gestazione, Lunghezza, Cranio e Sesso spiega circa il 73% della variabilità del peso neonatale. Tutti i regressori inclusi risultano significativi, confermando che fattori antropometrici e biologici sono determinanti per il peso alla nascita.

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.

Nella sezione precedente abbiamo già applicato la procedura di Stepwise, eliminando progressivamente i regressori non significativi (Anni.madre, Fumatrici, Ospedale, Tipo.parto e N.gravidanze). Ora confrontiamo i modelli annidati ottenuti con metodi formali.

Adesso uso i seguenti metodi per confrontare i modelli annidati:

  • Anova : Metodo di confronto di modelli annidati, valutando la riduzione della devianza residua (RSS).

  • BIC: Metodo di confronto modelli, che permette di confrontare la bontà e l’adattamento dei modelli.

anova(mod3,mod2,mod_full)
## Analysis of Variance Table
## 
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Sesso
## Model 3: Peso ~ Anni.madre + N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Fumatrici + Tipo.parto + Ospedale + Sesso
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1   2495 188688687                                  
## 2   2493 187601677  2   1087011 7.2433 0.0007301 ***
## 3   2489 186762521  4    839155 2.7959 0.0247927 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Commento: Il mod3 ha RSS più basso e quindi spiega meglio la variabilità dei dati rispetto agli altri. Inoltre p-value = 0.025 < alpha = 0.05 -> il p-value è significativo.

Risultato: Il mod3 è il modello più semplice e parsimonioso.

BIC(mod_full,mod2,mod3)
##          df      BIC
## mod_full 12 35241.84
## mod2      8 35221.75
## mod3      6 35220.54

Commento: Il modello con BIC minore è il mod3, quindi confermiamo la scelta del mod3 come quello che spiega meglio la variabilità del dataset.

Conclusione:

Sia l’ANOVA che il BIC convergono sulla scelta di mod3 come modello ottimale. Questo modello, che include solo Gestazione, Lunghezza, Cranio e Sesso, spiega circa il 73% della variabilità del peso neonatale, mantenendo un numero ridotto di regressori e garantendo interpretabilità clinica e statistica.

2.4 Analisi della Qualità del Modello

2.4.1 Valutazione R², RMSE, significatività coefficienti

Una volta ottenuto il modello finale (mod3), valutiamo la sua capacità predittiva e la bontà dell’adattamento. Abbiamo già visto che:

  • R^2 = 0.73 : Il modello spiega il 73% della variabilità del peso neonatale.

  • RMSE = 274 g : l’errore medio nella previsione del Peso (g).

Inoltre l’analisi del modello parte dal controllo dei problemi di multicollinearità:

knitr::kable(car::vif(mod3), col.names = c('indice multicollinearità'))
indice multicollinearità
Gestazione 1.653502
Lunghezza 2.069517
Cranio 1.606131
Sesso 1.038814

Commento: Tutti le variabili hanno VIF basso (<5) e quindi non sussistono problemi di multicollinearità.

Ora approfondiamo con l’analisi dei residui e la ricerca dei valori influenti.

2.4.2 Analisi dei residui e valori influenti

Obiettivo: Verificare che i residui rispettino le assunzioni del modello lineare:

  • Media dei residui = 0

  • Varianza costante (omoschedaticità)

  • Distribuzione dei residui relativamente normale.

  • Assenza di Pattern sistematici.

Per prima cosa facciamo un’analisi grafica dei residui.

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

Risultati del Grafico:

  1. Residuals vs Fitted : il grafico Residuals vs Fitted mostra un leggero pattern curvilineo, segnale di possibile non linearità o eteroschedasticità.

  2. Q-Q Residuals : il Q-Q plot sembra vicino alla normale, ma il test di Shapiro rifiuta l’ipotesi di normalità.

  3. Scale-Location : i punti seguono una distribuzione uniforme lungo un pattern predefinito. La linea rossa che definisce questo pattern è leggermente curvilineo e ciò dimostra la non eteroschedacità e omoschedasticità del modello.

  4. Residuals-Leverage : La maggior parte dei punti si trova all’interno delle curve di Cook. Due osservazioni (etichettate “1780” e “1551”) hanno un po’ più di leva e residui più alti con “1551” dentro la curva di Cook’s DIstance. Questo è un’osservazione pericolosa cha va oltre la distanza di cook di 0.5 avvicinandosi molto a 1.

knitr::kable(data[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

Effettuiamo quindi i singoli test per studiare come si comportano le osservazioni e i residui e valutare la qualità del modello.

shapiro.test(mod3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod3$residuals
## W = 0.9742, p-value < 2.2e-16
lmtest::bptest(mod3)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod3
## BP = 89.148, df = 4, p-value < 2.2e-16
lmtest::dwtest(mod3)
## 
##  Durbin-Watson test
## 
## data:  mod3
## DW = 1.9557, p-value = 0.1337
## alternative hypothesis: true autocorrelation is greater than 0
plot(density(residuals(mod3)))

Risultati Test Statistici:

  • Shapiro Test | W = 0.9742, p-value < 2.2e-16 : Il p-value è basso quindi si rifiuta H0 di normalità. I residui non sono compatibili con una distribuzione normale.

  • Breusch–Pagan | BP = 89.148, df = 4, p-value < 2.2e-16 : anche qui p-value basso quindi non ci sono evidenze di omoschedaticità. La varianza dei residui non può essere considerata costante.

  • Durbin-Watson | DW = 1.9557, p-value = 0.1337 : p-value> alpha =0.05 quindi in questo caso il valore è molto vicino a 2, e il p-value alto indica che non c’è autocorrelazione significativa. I residui sono indipendenti, come richiesto dal modello lineare classico.

Conclusione: Il modello mod3 supera solo 1 dei 3 test diagnostici, infatti l’indipendenza è OK mentre non vengono rispettate le assunzioni di Normalità e Varianza Costante. Possiamo utilizzare il modello perché le violazioni sono lievi ma bisogna fare particolare attenzione in fase di previsione.

2.4.2.1 Analisi Leverage e Outliers
#Leverage Analysis
par(mfrow=c(1,1))
lev <- hatvalues(mod3)
plot(lev)
p <- sum(lev)
n <- length(lev)
soglia = 2*p/n
abline(h=soglia, col=2)

lev[lev > soglia]
##          15          34          42          61          67          96 
## 0.006144457 0.006237758 0.004252678 0.004587185 0.005345322 0.004801899 
##         101         106         117         131         151         155 
## 0.007185969 0.012812305 0.004746840 0.006964953 0.010847336 0.006670119 
##         190         205         206         220         249         295 
## 0.005288500 0.005297016 0.009402903 0.007376867 0.004655679 0.004008437 
##         304         305         310         312         315         348 
## 0.004419003 0.005382162 0.028757757 0.013126063 0.005342858 0.004187962 
##         378         383         445         471         486         492 
## 0.015366923 0.004305772 0.007094099 0.004289364 0.004740595 0.008175223 
##         565         587         592         615         629         638 
## 0.004635950 0.008375235 0.006313937 0.004575867 0.004041054 0.006216787 
##         656         666         684         697         702         726 
## 0.004735898 0.004345904 0.008750225 0.005809934 0.004719868 0.004038656 
##         748         750         765         805         821         895 
## 0.008236046 0.006712718 0.006049647 0.014303716 0.004039952 0.005203974 
##         928         947         956         964         968         991 
## 0.022095348 0.007803378 0.007670034 0.004618460 0.004075295 0.004259145 
##        1014        1067        1091        1096        1130        1157 
## 0.008221844 0.007868942 0.008927908 0.004268008 0.006561457 0.004080725 
##        1166        1181        1188        1200        1238        1248 
## 0.004020427 0.005586402 0.006404596 0.005445729 0.005374725 0.014573080 
##        1273        1283        1293        1294        1356        1357 
## 0.007080183 0.004055883 0.005555616 0.004735838 0.005285646 0.006526609 
##        1361        1385        1395        1400        1402        1420 
## 0.004080786 0.012017154 0.004646054 0.005559147 0.004788164 0.005130876 
##        1428        1429        1551        1553        1556        1560 
## 0.008191659 0.021037040 0.048521821 0.006810651 0.005868384 0.004588261 
##        1593        1606        1610        1619        1628        1634 
## 0.004855489 0.004746852 0.007761683 0.014504316 0.004663751 0.004527422 
##        1686        1692        1693        1701        1712        1735 
## 0.008715871 0.004260736 0.005041490 0.010197488 0.006945920 0.004370535 
##        1780        1802        1806        1809        1827        1858 
## 0.025528862 0.004005956 0.004090967 0.008388773 0.006008547 0.004328741 
##        1868        1893        1911        1920        1977        2037 
## 0.004746360 0.004245302 0.004481835 0.004131544 0.005384125 0.004234559 
##        2040        2066        2089        2114        2115        2120 
## 0.011429685 0.004013637 0.006252340 0.012850646 0.011763363 0.018002540 
##        2140        2149        2175        2200        2215        2216 
## 0.006090376 0.013033782 0.021167345 0.011030260 0.004833388 0.007548562 
##        2220        2224        2225        2257        2307        2318 
## 0.004023907 0.005780067 0.005408611 0.005767952 0.013898144 0.004770755 
##        2337        2359        2391        2408        2437        2452 
## 0.004700477 0.004750437 0.004386621 0.009632048 0.023757012 0.023227398 
##        2458        2476        2478 
## 0.008002531 0.004070348 0.005745434
cat("ci sono", length(lev[lev>soglia]),"osservazioni lontane dallo spazio dei regressori") 
## ci sono 135 osservazioni lontane dallo spazio dei regressori

Commento: Il mod3 ha ben 135 osservazioni strutturalmente importanti ma non pericolosi (punti critici < 0.2, 0.3)

#Oultiers
plot(rstudent(mod3), xlab = "Osservazioni", ylab = "Valore rstudent outlier", main = "Osservazioni Outliers")
abline(h=c(-2,2)) 

car::outlierTest(mod3)
##      rstudent unadjusted p-value Bonferroni p
## 1551 9.986149         4.7193e-23   1.1798e-19
## 155  4.951276         7.8654e-07   1.9663e-03
## 1306 4.781188         1.8440e-06   4.6100e-03

Risultati:

Queste osservazioni (1551, 155, 1306) hanno un peso/residuo che non è compatibile con il resto del campione. Tutti ben oltre la soglia ±3: sono osservazioni molto anomale rispetto al modello.

  • P-value estremamente basso per ciascun outlier

  • Bonferroni: anche dopo correzione i valori del p-value risultano estremamente significativi e questo significa che le osservazioni 1551, 155, 1306 sono outlier estremamente significativi.

2.4.2.2 Distanza di Cook
cook<- cooks.distance(mod3)
plot(cook,ylim= c(0,1))

Commento: 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.

Conclusione:

  • Esistono osservazioni (1551,155,1306) 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.

  • Il modello rispetta l’indipendenza dei residui, ma mostra deviazioni da normalità e omoschedasticità. Tuttavia, con N=2500 il modello rimane robusto e utilizzabile, pur con cautela in fase predittiva.

  • Posso utilizzare il modello creato mod3 perché le violazioni sono lievi ma bisogna fare particolare attenzione in fase di previsione.

Sintesi Finale : Dopo aver costruito e confrontato diversi modelli, il modello finale (mod3) con Gestazione, Lunghezza, Cranio e Sesso risulta il più parsimonioso e spiega circa il 73% della variabilità del peso neonatale. L’analisi diagnostica evidenzia alcune violazioni (normalità ed eteroschedasticità dei residui) e la presenza di osservazioni influenti (es. ID 1551), che meritano ulteriori indagini. Nonostante ciò, il modello mantiene una buona capacità predittiva (RMSE ≈ 274 g) ed è interpretabile clinicamente: i fattori antropometrici e biologici sono determinanti per il peso alla nascita.

3. Previsioni e Risultati

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

data_prediction<- data.frame(
  N.gravidanze = 3,
  Gestazione   = 39,
  Lunghezza    = 500, #avg campione
  Cranio       = 340, #avg campione
  Sesso        = "F",
  Fumatrici    = 0
) 
cat("Il neonato peserà gr", round(predict(mod3, newdata = data_prediction),2))
## Il neonato peserà gr 3299.19

Commento: Il modello prevede che una neonata con queste caratteristiche abbia un peso medio atteso di circa 3299 g.

4. Visualizzazioni

Obiettivo: Utilizziamo i grafici e rappresentazioni visive per comunicare i risultati e mostrare i risultati del modello e mostrare le relazioni più significative tra le variabili.

In particolare, analizzo 2 aspetti rilevanti d’interesse clinico:

  • L’impatto del numero di settimane di gestazione sul peso neonatale,

  • L’effetto del fumo materno sul peso previsto.

4.1.1 Relazione tra Gestazione e Peso

library(ggplot2)
## Warning: il pacchetto 'ggplot2' è stato creato con R versione 4.5.2
ggplot(data= data)+
  geom_point(aes(x = Gestazione, y = Peso, col= Sesso)) +
  geom_smooth(aes(x = Gestazione, y = Peso, col= Sesso),method="lm", se=FALSE) +
  geom_smooth(aes(x = Gestazione,y = Peso),
            col = "black", se = FALSE, method = "lm")+
  labs(title="Relazione tra settimane di gestazione e peso neonatale",
       x="Gestazione [weeks]",
       y="Peso [g]")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Commento: Il grafico mostra chiaramente una relazione positiva tra le settimane di gestazione e il peso neonatale: all’aumentare della gravidanza, il peso della nascita tende ad aumentare. Inoltre distinguendo per sesso:

  • La retta di regressione dei maschi si colloca leggermente più in alto rispetto a quella delle femmine.

  • Questo suggerisce che a parità di settimane di gestazione, i neonati maschi tendono ad avere un peso medio superiore rispetto alle femmine.

  • La distribuzione delle osservazioni non è uniforme: ci sono poche nascite premature (<35 settimane), quindi in quell’intervallo el stime sono meno robuste.

Risultato: La durata della gestazione è un fattore significativo per il peso della nascita e il sesso del neonato genera una variazione sistematica (i maschi risultano mediamente più pesanti delle femmine).

4.1.2 Effetto del fumo materno sul Peso

Smoker<- factor(Fumatrici,levels=c(0,1), labels = c("NO","SI"))
ggplot(data= data)+
  geom_point(aes(x= Gestazione ,
                 y= Peso,
                 col= Smoker))+
  geom_smooth(aes(x= Gestazione,
                  y= Peso, col=Smoker), se= FALSE, method = "lm")+
  
labs(title = "Effetto del fumo materno sul peso neonatale",x= "Gestazione [weeks]", y= "Peso [g]")
## `geom_smooth()` using formula = 'y ~ x'

Commento: Il grafico mostra la relazione tra settimane di gestazione e peso neonatale, distinguendo tra madri fumatrici e non fumatrici. In entrambi i gruppi si osserva un incremento del peso con l’aumentare della durata della gravidanza; tuttavia la retta di regressione relativa alle madri fumatrici si colloca sistematicamente al di sotto di quella delle non fumatrici.

  • Solo un numero ridotto di madri nel dataset è fumatrice q euindi la retta di regressione per questo gruppo è meno stabile.

Questo risultato evidenzia come il fumo materno sia associato a un peso neonatale mediamente inferiore, indipendente dall’età gestazionale.

4.1.3 Relazione tra Peso e Misure Antropometriche

Infine rappresentiamo la relazione tra Peso e le altre misure antropometriche (Lunghezza e Cranio), che in definitiva sono le variabili quantitative del modello finale ottenuto.

#car::scatter3d(Peso ~ Lunghezza + Cranio, data=data)
library(scatterplot3d)
## Warning: il pacchetto 'scatterplot3d' è stato creato con R versione 4.5.2
scatterplot3d(Lunghezza, Cranio, Peso, pch=20, color="blue",
              xlab="Lunghezza [mm]", ylab="Cranio [mm]", zlab="Peso [gr]", main = "Relazione tra Peso e Misure Antropometriche", cex.symbols = 0.8)

Commento: Anche nel grafico 3D si nota l’osservazione anomale 1551, già vista nell’analisi dei residui.

Conclusione: Le visualizzazioni confermano i risultati del modello: la durata della gestazione e le misure antropometriche sono i principali predittori del peso neonatale, mentre il fumo materno mostra un effetto negativo, seppur meno marcato. Le rappresentazioni grafiche rendono immediata la comprensione di tali relazioni.

5. Conclusioni

Il progetto di previsione del peso neonatale rappresenta un’iniziativa fondamentale per Neonatal Health Solutions, in quanto dimostra come l’uso di dati clinici dettagliati e strumenti di analisi statistica avanzati possa tradursi in conoscenze utili e applicabili.

Dall’analisi è emerso che:

  • Le variabili antropometriche ( Lunghezza, Cranio ) e la durata della gestazione (Gestazione) sono i predittori più forti e significativi del peso dei neonati.

  • Il sesso del neonato incide sistematicamente, con i maschi mediamente più pesanti delle femmine.

  • Variabili come età materna, tipo di parto, ospedale e fumo materno non hanno mostrato un impatto statisticamente rilevante nel dataset, pur mantenendo un interesse clinico che merita un ulteriore approfondimenti.

Il modello finale spiega 73% della variabilità del peso neonatale, garantendo un buon equilibrio tra accuratezza e semplicità. Le analisi diagnostiche hanno evidenziato alcune violazioni delle assunzioni classiche come:

  • La non normalità e eteroschedasticità dei residui

  • La presenza di outlier e osservazioni influenti (es. 1551)

Dal punto di vista applicativo, il modello consente di applicare stime predittive utili in ambito clinico e gestionale, supportando i professionisti sanitari nella valutazione del rischio neonatale e nella pianificazione delle risorse ospedaliere.