Modello Statistico per la Previsione del Peso Neonatale

Librerie

Importo le librerie necessarie allo sviluppo corretto dello studio

library(moments)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
library(tidyr)
library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

library(ggpubr)

Importo il dataset cosi da poter iniziare le varie analisi

neonati <- read.csv("neonati.csv", sep = ",")
str(neonati)

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

Analisi preliminare statistiche descrittive

Prima di procedere con l’analisi inferenziale, è fondamentale esplorare le caratteristiche principali delle variabili numeriche attraverso un’analisi descrittiva. Questo passaggio iniziale permette di ottenere una visione d’insieme della distribuzione dei dati, evidenziando:

Valori centrali (media e mediana)
Misure di dispersione (deviazione standard e coefficiente di variazione)
Indicatori di forma (asimmetria e curtosi)

Tali informazioni sono essenziali per comprendere la natura dei dati e individuare eventuali anomalie o outlier. Inoltre, le statistiche descrittive offrono indicazioni utili per orientare la scelta delle tecniche inferenziali più adeguate e per facilitare l’interpretazione dei risultati del successivo modello predittivo.

summary_stats <- function(x) {
  c(
    Media = round(mean(x, na.rm = TRUE), 2),
    Mediana = round(median(x, na.rm = TRUE), 2),
    Varianza = round(var(x, na.rm = TRUE), 2),
    Deviazione_Std = round(sd(x, na.rm = TRUE), 2),
    Coeff_Variazione = round((sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE)) * 100, 2),
    Asimmetria = round(skewness(x, na.rm = TRUE), 2),
    Curtosi = round(kurtosis(x, na.rm = TRUE), 2)
  )
}

numerical_vars <- neonati[, c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")]
stats <- apply(numerical_vars, 2, summary_stats)
print(stats)

##                  Anni.madre N.gravidanze Gestazione      Peso Lunghezza Cranio
## Media                 28.16         0.98      38.98   3284.08    494.69 340.03
## Mediana               28.00         1.00      39.00   3300.00    500.00 340.00
## Varianza              27.81         1.64       3.49 275665.68    692.67 269.79
## Deviazione_Std         5.27         1.28       1.87    525.04     26.32  16.43
## Coeff_Variazione      18.72       130.51       4.79     15.99      5.32   4.83
## 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

Le statistiche descrittive evidenziano una distribuzione generalmente regolare per la maggior parte delle variabili numeriche, con valori medi e mediani molto simili, variabilità contenuta e modesta asimmetria.

Le variabili Lunghezza e Cranio mostrano una variabilità molto bassa, mentre Gestazione presenta la variabilità maggiore; l’analisi di variabilità è stata condotta considerando il Coefficiente di Variazione per poter eliminare l’influenza dovuta a differenti ordini di grandezza (es: Peso e Lunghezza).

La variabile Gestazione è fortemente asimmetrica a sinistra, e N.gravidanze a destra, richiedendo attenzione nelle fasi inferenziali successive. La curtosi elevata in alcune variabili suggerisce la presenza di code pesanti, ma nel complesso i dati appaiono idonei all’applicazione di test parametrici e modelli lineari, che verranno approfonditi nelle sezioni successive.

freq_table_fun <- function(var) {
  tabella <- table(var)
  totale <- sum(tabella)
  freq_rel <- tabella / totale
  freq_cum <- cumsum(tabella)
  freq_rel_cum <- cumsum(freq_rel)
  
  cbind(
    Frequenza_Assoluta = as.vector(tabella),
    Frequenza_Relativa = round(freq_rel, 4),
    Frequenza_Cumulata = freq_cum,
    Freq_Relativa_Cumulata = round(freq_rel_cum, 4)
  )
}

freq_table_fumo <- freq_table_fun(neonati$Fumatrici)
rownames(freq_table_fumo) <- names(table(neonati$Fumatrici))
print(freq_table_fumo)

##   Frequenza_Assoluta Frequenza_Relativa Frequenza_Cumulata
## 0               2396             0.9584               2396
## 1                104             0.0416               2500
##   Freq_Relativa_Cumulata
## 0                 0.9584
## 1                 1.0000

freq_table_sesso <- freq_table_fun(neonati$Sesso)
rownames(freq_table_sesso) <- names(table(neonati$Sesso))
print(freq_table_sesso)

##   Frequenza_Assoluta Frequenza_Relativa Frequenza_Cumulata
## F               1256             0.5024               1256
## M               1244             0.4976               2500
##   Freq_Relativa_Cumulata
## F                 0.5024
## M                 1.0000

freq_table_parto <- freq_table_fun(neonati$Tipo.parto)
rownames(freq_table_parto) <- names(table(neonati$Tipo.parto))
print(freq_table_parto)

##     Frequenza_Assoluta Frequenza_Relativa Frequenza_Cumulata
## Ces                728             0.2912                728
## Nat               1772             0.7088               2500
##     Freq_Relativa_Cumulata
## Ces                 0.2912
## Nat                 1.0000

freq_table_ospedale <- freq_table_fun(neonati$Ospedale)
rownames(freq_table_ospedale) <- names(table(neonati$Ospedale))
print(freq_table_ospedale)

##      Frequenza_Assoluta Frequenza_Relativa Frequenza_Cumulata
## osp1                816             0.3264                816
## osp2                849             0.3396               1665
## osp3                835             0.3340               2500
##      Freq_Relativa_Cumulata
## osp1                 0.3264
## osp2                 0.6660
## osp3                 1.0000

Analisi delle Variabili Categoriali

L’esplorazione delle variabili categoriali evidenzia alcune caratteristiche rilevanti per le successive analisi inferenziali:

Fumo materno: la stragrande maggioranza delle madri nel campione non fuma (circa il 96%).
Sesso del neonato: la distribuzione tra neonati maschi e femmine appare bilanciata, condizione ottimale per condurre confronti statistici tra i due gruppi. Tale equilibrio consente l’applicazione di test parametrici come il t-test garantendo una buona robustezza nei risultati.
Tipo di parto: circa il 71% dei parti è avvenuto per via naturale, mentre il restante 29% mediante taglio cesareo.
Ospedale di provenienza: il campione risulta equamente distribuito tra i tre ospedali partecipanti allo studio.

Sebbene alcune variabili mostrino una distribuzione sbilanciata, la numerosità dei sottogruppi resta generalmente sufficiente per supportare analisi comparative attendibili nel contesto inferenziale.

Boxplot

ggplot(neonati, aes(x = factor(Fumatrici), y = Peso, fill = factor(Fumatrici))) +
  geom_boxplot() +
  labs(
    title = "Distribuzione del Peso alla Nascita per Fumatrici",
    x = "Fumatrici (0 = No, 1 = Sì)",
    y = "Peso alla Nascita (g)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

Il boxplot mostra una differenza evidente tra i due gruppi:

I neonati di madri fumatrici (1) presentano un peso sistematicamente inferiore rispetto a quelli di madri non fumatrici (0).
La mediana del peso per le i neonati con madri non fumatrici è intorno a 3300 g, con una distribuzione più ampia e una maggiore variabilità.
Per le fumatrici, la mediana si attesta intorno a 3100 g, con una distribuzione più concentrata.
Si osserva un numero maggiore di outlier tra i neonati di madri non fumatrici, soprattutto per pesi molto bassi (sotto i 2000 g).

Questo grafico conferma visivamente l’effetto negativo del fumo materno sul peso neonatale, coerente con le ipotesi e la letteratura esistente.

Istogramma

ggplot(neonati, aes(x = Peso)) +
  geom_histogram(binwidth = 250, fill = "skyblue", color = "black") +
  labs(
    title = "Distribuzione del Peso alla Nascita",
    x = "Peso (g)",
    y = "Frequenza"
  ) +
  theme_minimal()

L’istogramma mostra che la distribuzione del peso neonatale è approssimativamente normale, con una leggera asimmetria negativa, confermata dal coefficiente di asimmetria pari a circa -0.65.

Il picco principale si concentra nella fascia 3000-3500 g, corrispondente al peso considerato fisiologico alla nascita.
La coda sinistra è più estesa, indicando la presenza di un numero limitato di neonati con peso molto basso.
La coda destra è meno pronunciata, con pochi casi estremi oltre i 4500 g.

Scatterplot

ggplot(neonati, aes(x = Gestazione, y = Peso)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE, color = "blue") +
  labs(
    title = "Peso alla Nascita in Funzione delle Settimane di Gestazione",
    x = "Settimane di Gestazione",
    y = "Peso alla Nascita (g)"
  ) +
  theme_minimal()

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

Lo scatterplot mette in evidenza una forte correlazione positiva e lineare tra settimane di gestazione e peso alla nascita

Si osserva una relazione robusta, dove ogni settimana aggiuntiva di gestazione è associata a un incremento del peso di circa 150-200 g.
La maggior parte dei dati si concentra tra le 37 e le 40 settimane, corrispondenti a gravidanze a termine.
I parti prematuri (< 37 settimane) sono chiaramente associati a pesi molto bassi, confermando l’importanza della durata della gestazione per la crescita fetale.
La variabilità del peso aumenta leggermente verso il termine, con una maggiore dispersione dei dati tra la 39ª e la 40ª settimana.

Questi risultati supportano l’ipotesi di una relazione lineare significativa tra durata della gestazione e peso alla nascita.

TEST PARAMETRICI

Test (1): In alcuni ospedali si fanno più parti cesarei

tabella <- table(neonati$Ospedale, neonati$Tipo.parto)


test_chi <- chisq.test(tabella)


test_chi$statistic

## X-squared 
##  1.097202

test_chi$p.value

## [1] 0.5777576

test_chi$expected

##       
##             Ces      Nat
##   osp1 237.6192 578.3808
##   osp2 247.2288 601.7712
##   osp3 243.1520 591.8480

cat(sprintf(
  "Test chi-quadro: p-value = %.4f\n",
  test_chi$p.value
))

## Test chi-quadro: p-value = 0.5778

Non si riscontra una relazione statisticamente significativa tra l’ospedale di appartenenza e il tipo di parto (p = 0.5778). Pertanto, si mantiene l’ipotesi nulla di indipendenza: la distribuzione dei tipi di parto sembra essere simile tra i diversi ospedali del campione analizzato.

Questo risultato è confermato anche tramite analisi grafica, come evidenziato nei grafici sotto riportati.

ggballoonplot(as.data.frame(tabella), x = "Var1", y = "Var2", size = "Freq", fill = "blue") +
  labs(title = "Distribuzione Tipo di Parto per Ospedale",
       x = "Ospedale", y = "Tipo di Parto")

La visualizzazione della tabella di contingenza mostra chiaramente che in ciascun ospedale (osp1, osp2, osp3) il parto naturale (Nat) è largamente prevalente rispetto al cesareo (Ces). Le bolle corrispondenti ai parti naturali sono grandi e simili tra i tre ospedali, indicando frequenze elevate e comparabili. Al contrario, bolle dei cesarei mostrano frequenze basse e relativamente stabili tra gli ospedali.

curve(dchisq(x, df = test_chi$parameter), 
      from = 0, 
      to = max(30, test_chi$statistic + 10), 
      col = "blue", 
      lwd = 2,
      main = "Distribuzione Chi-quadro sotto H0",
      xlab = "Valori della statistica chi-quadro",
      ylab = "Densità")

abline(v = qchisq(0.95, df = test_chi$parameter), col = "red", lwd = 2, lty = 2)
points(test_chi$statistic, 0, col = "darkgreen", pch = 19, cex = 2)

legend("topright",
       legend = c("Distribuzione sotto H0", 
                  "Valore soglia (α=0.05)", 
                  "Statistica osservata"),
       col = c("blue", "red", "darkgreen"),
       lwd = c(2, 2, NA),
       lty = c(1, 2, NA),
       pch = c(NA, NA, 19))

Il test chi-quadro (X² = 1.097, p = 0.578) conferma l’indipendenza tra ospedale e tipo di parto. La statistica osservata (pallino verde) resta ampiamente sotto la soglia critica (linea rossa), non fornendo evidenze di associazione significativa tra le variabili.

Test (2): La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione

Per la costruzione dei test statistici, è particolarmente utile disporre di valori di riferimento noti per la media e la deviazione standard delle variabili peso e lunghezza neonatale. Da una rapida ricerca online emerge che tali parametri assumono i seguenti valori:

Media del peso alla nascita: 3.300 g
Deviazione standard del peso: 500 g
Media della lunghezza alla nascita: 500 mm
Deviazione standard della lunghezza: 20 mm

Essendo noti i valori di media e deviazione standard della popolazione, è possibile applicare lo Z-test per verificare entrambe le ipotesi nulle, sia quella relativa al peso, sia quella relativa alla lunghezza alla nascita.

Peso

mu0 <- 3300
sigma <- 500
peso <- na.omit(neonati$Peso)
n <- length(peso)
media_peso <- mean(peso)

Z <- (media_peso - mu0) / (sigma / sqrt(n))

p_value <- 2 * pnorm(-abs(Z))

IC <- c(
  media_peso - qnorm(0.975) * (sigma / sqrt(n)),
  media_peso + qnorm(0.975) * (sigma / sqrt(n))
)

cat(sprintf(
  "Media campionaria: %.2f | Statistica Z: %.4f | p-value: %s | IC 95%%: [%.2f; %.2f]\n",
  media_peso, Z, format.pval(p_value, digits = 10), IC[1], IC[2]
))

## Media campionaria: 3284.08 | Statistica Z: -1.5919 | p-value: 0.1114026811 | IC 95%: [3264.48; 3303.68]

L’output mostra il risultato di un test Z per valutare se la media del peso neonatale differisce significativamente dal valore noto di riferimento di 3300 grammi. Il p-value ottenuto è superiore ai comuni livelli di significatività (es. 0.05), il che indica che non rifiutiamo l’ipotesi nulla e concludiamo che la media del peso non è significativamente diversa da quella della popolazione.

curve(dnorm(x, 0, 1),
      from = -4, to = 4,
      col = "blue", lwd = 2,
      main = "Distribuzione Normale sotto H0",
      xlab = "Z", ylab = "Densità")
abline(v = c(-qnorm(0.975), qnorm(0.975)), col = "red", lty = 2, lwd = 2)
points(Z, 0, col = "darkgreen", pch = 19, cex = 2)
legend("topright",
       legend = c("N(0,1)", "Soglie critiche", "Z osservato"),
       col = c("blue", "red", "darkgreen"),
       lwd = c(2, 2, NA), lty = c(1, 2, NA), pch = c(NA, NA, 19))

Il grafico mostra la distribuzione normale standard sotto l’ipotesi nulla (Z ~ N(0,1)) per il test Z sul peso neonatale. Le linee rosse tratteggiate rappresentano i valori soglia critici a ±1.96, corrispondenti a un livello di significatività del 5% (α = 0.05) per un test bilaterale. Il punto verde indica il valore Z osservato e si trova all’interno dell’intervallo di accettazione.

Questo conferma visivamente quanto evidenziato dall’output numerico:

Il valore Z cade nell’area di accettazione
L’ipotesi nulla non viene respinta
La media campionaria del peso non è significativamente diversa da quella attesa (3300 grammi)

Lunghezza

mu0 <- 500
sigma <- 20  

lung <- na.omit(neonati$Lunghezza)
n <- length(lung)
media_lung <- mean(lung)


Z_stat <- (media_lung - mu0) / (sigma / sqrt(n))
p_value <- 2 * pnorm(-abs(Z_stat))


IC <- c(
  media_lung - qnorm(0.975) * (sigma / sqrt(n)),
  media_lung + qnorm(0.975) * (sigma / sqrt(n))
)


cat(sprintf(
  "Media campionaria: %.2f | Statistica Z: %.4f | p-value: %s | IC 95%%: [%.2f; %.2f]\n",
  media_lung, Z_stat, format.pval(p_value, digits = 15), IC[1], IC[2]
))

## Media campionaria: 494.69 | Statistica Z: -13.2700 | p-value: < 2.22044604925e-16 | IC 95%: [493.91; 495.48]

Il test Z mostra che la lunghezza neonatale media del campione (494.69) è significativamente diversa da quella della popolazione di riferimento (p-value < 2.22e-16), con un intervallo di confidenza al 95% [493.91; 495.48] che non include la media attesa: la differenza è quindi altamente significativa dal punto di vista statistico.

curve(dnorm(x), from = -15, to = 15, col = "blue", lwd = 2,
      main = "Distribuzione Normale Standard sotto H0 (test Z)",
      xlab = "Valori Z", ylab = "Densità")
abline(v = c(-qnorm(0.975), qnorm(0.975)), col = "red", lty = 2, lwd = 2)
points(Z_stat, 0, col = "darkgreen", pch = 19, cex = 2)

legend("topright",
       legend = c("Distribuzione Normale Standard", 
                  "Soglie critiche ±z(0.975)", 
                  "Z osservato"),
       col = c("blue", "red", "darkgreen"),
       lwd = c(2, 2, NA), lty = c(1, 2, NA), pch = c(NA, NA, 19))

Il grafico mostra come il valore Z osservato per la lunghezza sia molto distante dal centro della distribuzione, indicando una forte significatività statistica. Questo ci porta a rifiutare con sicurezza l’ipotesi nulla. In particolare:

Il valore Z è circa -13, ben oltre la soglia critica;
Il p-value è praticamente zero, confermando la significatività;
La media campionaria è significativamente inferiore al valore ipotizzato.

In conclusione, c’è una chiara evidenza di una differenza reale nella lunghezza media della popolazione rispetto all’ipotesi iniziale.

Le misure antropometriche sono significativamente diverse tra i due sessi

Per confrontare le misure antropometriche tra maschi e femmine è stato utilizzato il test t di Student. Poiché l’analisi riguarda tre diverse variabili (lunghezza, peso e circonferenza cranica), è stata creata una funzione che consente di evitare la ripetizione del codice per ciascun confronto.

t_test_confronto <- function(x1, x2, nome_variabile = "Variabile") {
  x1 <- na.omit(x1)
  x2 <- na.omit(x2)
  
  n1 <- length(x1)
  n2 <- length(x2)
  m1 <- mean(x1)
  m2 <- mean(x2)
  s1 <- sd(x1)
  s2 <- sd(x2)
  

  se <- sqrt(s1^2/n1 + s2^2/n2)
  t_stat <- (m1 - m2) / se
  df <- (s1^2/n1 + s2^2/n2)^2 / ((s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1))
  p_value <- 2 * pt(-abs(t_stat), df)
  
  IC <- c((m1 - m2) - qt(0.975, df)*se, (m1 - m2) + qt(0.975, df)*se)
  
  cat(sprintf(
    "Confronto %s (Maschi - Femmine) | Diff. medie: %.2f | t: %.4f | p-value: %s | IC 95%%: [%.2f; %.2f]\n",
    nome_variabile, m1 - m2, t_stat, format.pval(p_value, digits = 15), IC[1], IC[2]
  ))
  
  curve(dt(x, df), from = -15, to = 15, col = "blue", lwd = 2,
        main = paste("t-test:", nome_variabile, "(Maschi - Femmine)"),
        xlab = "Valori t", ylab = "Densità")
  abline(v = c(-qt(0.975, df), qt(0.975, df)), col = "red", lty = 2, lwd = 2)
  points(t_stat, 0, col = "darkgreen", pch = 19, cex = 2)
  legend("topright",
         legend = c("Distribuzione t", "Soglie critiche", "t osservato"),
         col = c("blue", "red", "darkgreen"),
         lwd = c(2, 2, NA), lty = c(1, 2, NA), pch = c(NA, NA, 19))
}

t_test_confronto(neonati$Lunghezza[neonati$Sesso == "M"],
                 neonati$Lunghezza[neonati$Sesso == "F"],
                 nome_variabile = "Lunghezza")

## Confronto Lunghezza (Maschi - Femmine) | Diff. medie: 9.90 | t: 9.5820 | p-value: < 2.22044604925e-16 | IC 95%: [7.88; 11.93]

t_test_confronto(neonati$Peso[neonati$Sesso == "M"],
                 neonati$Peso[neonati$Sesso == "F"],
                 nome_variabile = "Peso")

## Confronto Peso (Maschi - Femmine) | Diff. medie: 247.08 | t: 12.1061 | p-value: < 2.22044604925e-16 | IC 95%: [207.06; 287.11]

t_test_confronto(neonati$Cranio[neonati$Sesso == "M"],
                 neonati$Cranio[neonati$Sesso == "F"],
                 nome_variabile = "Cranio")

## Confronto Cranio (Maschi - Femmine) | Diff. medie: 4.82 | t: 7.4102 | p-value: 1.71769469990384e-13 | IC 95%: [3.54; 6.09]

Tutti e tre i confronti mostrano differenze significative tra i sessi, con i maschi che mediamente presentano valori più elevati per lunghezza, peso e circonferenza cranica. Queste evidenze statistiche suggeriscono che il sesso del neonato è un fattore associato in modo sistematico alle sue misure antropometriche alla nascita.

SCATTERPLOT E 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, use="complete.obs")
  txt <- format(c(r, 0.123456789), digits=digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8 / strwidth(txt)
  text(0.5, 0.5, txt, cex=cex.cor)
}

pairs(neonati[, c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")],
      lower.panel=panel.cor, upper.panel=panel.smooth)

Le correlazioni più forti si osservano tra le misure antropometriche (peso, lunghezza, cranio), che risultano strettamente interconnesse. In particolare:

Peso e lunghezza mostrano la correlazione più elevata (r = 0.80);
La durata della gestazione è un fattore chiave, correlato positivamente con peso e lunghezza;
Fattori materni come l’età e il numero di gravidanze presentano invece legami deboli con le misure neonatali.

Nel complesso, il pattern emerso suggerisce l’importanza della gestazione come variabile predittiva.

Modello di regressione lineare completo

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

## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione + 
##     Fumatrici + Sesso + Tipo.parto + Ospedale, data = neonati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1488.2  -270.1   -13.6   262.7  1905.9 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -3306.932    187.373 -17.649  < 2e-16 ***
## Anni.madre        3.742      1.706   2.193  0.02838 *  
## N.gravidanze     18.437      7.010   2.630  0.00859 ** 
## Gestazione      163.461      4.520  36.163  < 2e-16 ***
## Fumatrici      -110.330     41.498  -2.659  0.00790 ** 
## SessoM          164.712     16.699   9.863  < 2e-16 ***
## Tipo.partoNat    15.329     18.213   0.842  0.40007    
## Ospedaleosp2     -1.426     20.273  -0.070  0.94395    
## Ospedaleosp3     23.949     20.366   1.176  0.23974    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 413.4 on 2491 degrees of freedom
## Multiple R-squared:  0.3821, Adjusted R-squared:  0.3801 
## F-statistic: 192.6 on 8 and 2491 DF,  p-value: < 2.2e-16

Il modello lineare completo mostra una relazione statisticamente significativa tra diverse variabili cliniche e il peso neonatale.

In particolare:

L’R² è pari a 0.38, indicando che circa il 38% della variabilità del peso alla nascita è spiegata dalle variabili incluse nel modello;
Il test F globale è altamente significativo (p < 2.2e-16), suggerendo che almeno un predittore ha un effetto reale sul peso;
Tra i predittori più influenti, la durata della gestazione si conferma il fattore dominante, con un aumento medio di circa 163 grammi per ogni settimana in più;
Anche il sesso del neonato ha un impatto rilevante: i maschi pesano in media circa 165 grammi in più rispetto alle femmine;
Altri effetti significativi, anche se più contenuti, includono l’età materna, il numero di gravidanze precedenti e il fumo materno.

Non risultano invece significative le variabili relative al tipo di parto e alla struttura ospedaliera, che non sembrano influenzare il peso in modo rilevante nel campione analizzato.

vif(mod_full)

##                  GVIF Df GVIF^(1/(2*Df))
## Anni.madre   1.183690  1        1.087975
## N.gravidanze 1.178562  1        1.085616
## Gestazione   1.043345  1        1.021443
## Fumatrici    1.004497  1        1.002246
## Sesso        1.019952  1        1.009927
## Tipo.parto   1.001659  1        1.000829
## Ospedale     1.003131  2        1.000782

L’analisi dei VIF conferma l’assenza di collinearità tra le variabili indipendenti, con valori ben al di sotto delle soglie critiche. Questo rafforza l’affidabilità delle stime ottenute.

In sintesi, il modello è solido e mostra evidenza chiara di più fattori clinici significativi associati al peso neonatale.

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

shapiro.test(residuals(mod_full))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_full)
## W = 0.9969, p-value = 5.722e-05

Il test di Shapiro-Wilk restituisce W = 0.9969 con un p-value = 5.722e-05, indicando che, sebbene i residui siano visivamente quasi normali (W vicino a 1), la normalità è formalmente rifiutata a livello di significatività del 5%. Dunque, i residui non seguono una distribuzione normale in modo statisticamente significativo, anche se la deviazione appare contenuta.

bptest(mod_full)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod_full
## BP = 13.553, df = 8, p-value = 0.09418

Il test di Breusch-Pagan restituisce BP = 13.553 con p-value = 0.09418, indicando che non si rifiuta l’ipotesi nulla di omoschedasticità. Pertanto, non emergono evidenze statisticamente significative di eteroschedasticità dei residui, suggerendo che la varianza degli errori può essere considerata costante.

dwtest(mod_full)

## 
##  Durbin-Watson test
## 
## data:  mod_full
## DW = 1.8939, p-value = 0.003997
## alternative hypothesis: true autocorrelation is greater than 0

Il test di Durbin-Watson restituisce DW = 1.8939 con un p-value = 0.003997: sebbene il valore DW sia vicino a 2 (indicando apparentemente assenza di autocorrelazione), il p-value < 0.05 porta a rifiutare l’ipotesi nulla di indipendenza dei residui. Si conclude quindi che è presente un’autocorrelazione positiva significativa tra i residui.

SELEZIONE MODELLO RIDOTTO

n <- nrow(neonati)
step_mod <- stepAIC(mod_full, direction = "both", k = log(n))

## Start:  AIC=30183.19
## Peso ~ Anni.madre + N.gravidanze + Gestazione + Fumatrici + Sesso + 
##     Tipo.parto + Ospedale
## 
##                Df Sum of Sq       RSS   AIC
## - Ospedale      2    338648 426000786 30170
## - Tipo.parto    1    121043 425783181 30176
## - Anni.madre    1    821976 426484114 30180
## - N.gravidanze  1   1182030 426844167 30182
## - Fumatrici     1   1207859 426869997 30182
## <none>                      425662138 30183
## - Sesso         1  16624132 442286270 30271
## - Gestazione    1 223470325 649132463 31230
## 
## Step:  AIC=30169.53
## Peso ~ Anni.madre + N.gravidanze + Gestazione + Fumatrici + Sesso + 
##     Tipo.parto
## 
##                Df Sum of Sq       RSS   AIC
## - Tipo.parto    1    130235 426131021 30162
## - Anni.madre    1    845501 426846287 30167
## - N.gravidanze  1   1206111 427206897 30169
## - Fumatrici     1   1231187 427231973 30169
## <none>                      426000786 30170
## + Ospedale      2    338648 425662138 30183
## - Sesso         1  16674239 442675025 30258
## - Gestazione    1 223933725 649934511 31218
## 
## Step:  AIC=30162.47
## Peso ~ Anni.madre + N.gravidanze + Gestazione + Fumatrici + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1    846924 426977945 30160
## - N.gravidanze  1   1188774 427319795 30162
## - Fumatrici     1   1215268 427346289 30162
## <none>                      426131021 30162
## + Tipo.parto    1    130235 426000786 30170
## + Ospedale      2    347841 425783181 30176
## - Sesso         1  16667830 442798852 30251
## - Gestazione    1 223815017 649946038 31210
## 
## Step:  AIC=30159.61
## Peso ~ N.gravidanze + Gestazione + Fumatrici + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1   1233672 428211618 30159
## <none>                      426977945 30160
## + Anni.madre    1    846924 426131021 30162
## - N.gravidanze  1   2379729 429357675 30166
## + Tipo.parto    1    131658 426846287 30167
## + Ospedale      2    371728 426606218 30173
## - Sesso         1  16780605 443758551 30248
## - Gestazione    1 223417523 650395468 31204
## 
## Step:  AIC=30159
## Peso ~ N.gravidanze + Gestazione + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      428211618 30159
## + Fumatrici     1   1233672 426977945 30160
## + Anni.madre    1    865329 427346289 30162
## - N.gravidanze  1   2216633 430428251 30164
## + Tipo.parto    1    115547 428096071 30166
## + Ospedale      2    395946 427815672 30172
## - Sesso         1  16718785 444930403 30247
## - Gestazione    1 222513567 650725185 31197

La selezione del modello mediante procedura stepwise basata sul criterio BIC ha portato a un modello finale composto da tre variabili: Numero di gravidanze, Settimane di gestazione e Sesso del neonato. Le altre variabili presenti nel modello iniziale (come Fumatrici, Anni della madre, Tipo di parto e Ospedale) sono state progressivamente escluse in quanto non apportavano un miglioramento sufficiente al bilanciamento tra bontà di adattamento e complessità del modello secondo il criterio BIC. Tra le variabili selezionate, la Gestazione si conferma come il principale predittore del peso neonatale.

summary(step_mod)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Sesso, data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1494.15  -272.82   -14.43   267.31  1892.39 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3136.835    175.410 -17.883  < 2e-16 ***
## N.gravidanze    23.395      6.509   3.595 0.000331 ***
## Gestazione     162.025      4.499  36.014  < 2e-16 ***
## SessoM         165.144     16.729   9.872  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 414.2 on 2496 degrees of freedom
## Multiple R-squared:  0.3784, Adjusted R-squared:  0.3777 
## F-statistic: 506.5 on 3 and 2496 DF,  p-value: < 2.2e-16

Il modello finale stima il peso neonatale in funzione di:

Numero di gravidanze (coeff. = 23,395 ; p < 0,001)
Settimane di gestazione (coeff. = 162,025; p < 0,001)
Sesso (coeff. = 165,144 per i maschi; p < 0,001)

Tutte le variabili sono statisticamente significative (p < 0,001), con Gestazione come principale predittore: ogni settimana in più di gestazione aumenta il peso di circa 162 grammi, mentre ogni gravidanza in più si associa in media a un aumento di circa 23 grammi.

L’R² pari a 37,8% indica che il modello spiega circa il 38% della variabilità del peso neonatale. I residui mostrano una variabilità ragionevole (errore standard residuo ≈ 414 grammi) e una distribuzione bilanciata intorno allo zero.

Il modello risulta statisticamente significativo nel suo complesso (F = 506,5; p < 0,001) e conferma il ruolo centrale della durata della gravidanza nel determinare il peso alla nascita. Il sesso maschile, infine, si associa coerentemente a un peso medio superiore.

Prova di interazioni

mod_inter<- lm(Peso ~ N.gravidanze + Sesso + Gestazione * Fumatrici, data = neonati)
summary(mod_inter)

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Sesso + Gestazione * Fumatrici, 
##     data = neonati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1498.69  -273.34   -14.34   264.56  1892.01 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -3196.668    177.196 -18.040  < 2e-16 ***
## N.gravidanze            24.367      6.506   3.745 0.000184 ***
## SessoM                 166.750     16.718   9.974  < 2e-16 ***
## Gestazione             163.634      4.544  36.009  < 2e-16 ***
## Fumatrici             1883.300   1140.061   1.652 0.098675 .  
## Gestazione:Fumatrici   -50.808     29.019  -1.751 0.080094 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 413.5 on 2494 degrees of freedom
## Multiple R-squared:  0.381,  Adjusted R-squared:  0.3797 
## F-statistic:   307 on 5 and 2494 DF,  p-value: < 2.2e-16

anova(step_mod, mod_inter)

## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Sesso
## Model 2: Peso ~ N.gravidanze + Sesso + Gestazione * Fumatrici
##   Res.Df       RSS Df Sum of Sq      F   Pr(>F)   
## 1   2496 428211618                                
## 2   2494 426453773  2   1757845 5.1401 0.005919 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Il modello con l’interazione tra Gestazione e Fumatrici migliora significativamente la spiegazione del peso neonatale rispetto al modello senza interazioni (p = 0.0059, test ANOVA).

Nello specifico:

L’interazione è marginalmente significativa (p ≈ 0.08), suggerendo un possibile effetto differenziato della gestazione sul peso a seconda dello stato di fumatore materno.
Le variabili principali (N.gravidanze, Sesso, Gestazione) restano altamente significative e confermano il loro ruolo predittivo.
L’R² aumenta leggermente (da ~37.8% a ~38.1%), indicando un miglioramento modesto ma statisticamente rilevante nella capacità esplicativa del modello.

In conclusione, includere l’interazione tra gestazione e fumo materno risulta corretto e utile per cogliere dinamiche più complesse.

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

L’analisi diagnostica dei residui evidenzia:

Residuals vs Fitted: residui centrati attorno a zero ma con leggera curvatura, suggerendo possibile non-linearità; varianza relativamente costante.
Q-Q Plot: residui generalmente normali al centro, con deviazioni nelle code, specialmente nella coda superiore, indicando outlier o code pesanti.
Scale-Location: lieve tendenza a varianza non costante (eteroschedasticità) lungo i valori fitted.
Residuals vs Leverage: la maggior parte dei punti mostra leverage basso e residui standardizzati compresi tra -4 e +4; tuttavia, alcune osservazioni (1920, 2245, 2089) presentano leverage più elevato e residui relativamente grandi.

L’analisi dei residui mostra lievi non-linearità, variabilità non costante e alcuni outlier influenti. Pur essendo il modello significativo, questi aspetti suggeriscono di considerare ulteriori miglioramenti per aumentare la sua robustezza e precisione predittiva.

shapiro.test(residuals(mod_inter))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_inter)
## W = 0.99683, p-value = 4.488e-05

Il test di Shapiro-Wilk indica una lieve deviazione dalla normalità (p-value molto basso), suggerendo che i residui non seguono perfettamente una distribuzione normale, probabilmente a causa di outlier o code pesanti. Tuttavia, la deviazione non è drammatica.

bptest(mod_inter)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod_inter
## BP = 6.9002, df = 5, p-value = 0.2282

Il test di Breusch-Pagan non rileva evidenze di varianza non costante dei residui (p-value > 0.2), quindi l’ipotesi di omoschedasticità è ragionevole. Ciò significa che la dispersione degli errori è stabile lungo i valori predetti.

dwtest(mod_inter)

## 
##  Durbin-Watson test
## 
## data:  mod_inter
## DW = 1.8957, p-value = 0.004543
## alternative hypothesis: true autocorrelation is greater than 0

Il test di Durbin-Watson evidenzia una autocorrelazione positiva significativa tra i residui (p-value < 0.01), indicando che gli errori non sono indipendenti, il che può influenzare l’affidabilità delle stime e richiede attenzione.

In conclusione, il modello mostra una buona adeguatezza generale, ma presenta alcune criticità: la non normalità lieve dei residui e la presenza di autocorrelazione possono influenzare l’affidabilità delle inferenze. La stabilità della varianza degli errori è invece confermata.

Validazione del modello

train_indices <- sample(seq_len(nrow(neonati)), size = 0.7 * nrow(neonati))
train_data <- neonati[train_indices, ]
test_data <- neonati[-train_indices, ]


mod_train <- lm(Peso ~ Anni.madre + N.gravidanze + Gestazione * Fumatrici + Sesso + Tipo.parto + Ospedale,
                data = train_data)


pred_train <- predict(mod_train, newdata = train_data)
pred_test <- predict(mod_train, newdata = test_data)


MSE <- function(y_obs, y_pred) mean((y_obs - y_pred)^2)
mse_train <- MSE(train_data$Peso, pred_train)
rmse_train <- sqrt(mse_train)
mse_test <- MSE(test_data$Peso, pred_test)
rmse_test <- sqrt(mse_test)

print(paste("MSE su training:", round(mse_train, 2)))

## [1] "MSE su training: 163932.2"

print(paste("RMSE su training:", round(rmse_train, 2)))

## [1] "RMSE su training: 404.89"

print(paste("MSE su test:", round(mse_test, 2)))

## [1] "MSE su test: 185213.13"

print(paste("RMSE su test:", round(rmse_test, 2)))

## [1] "RMSE su test: 430.36"

La validazione del modello è stata effettuata dividendo il dataset in training (70%) e test (30%).

Il modello è stato costruito sui dati di training e poi utilizzato per prevedere il peso neonatale sia nel training che nel test.

L’errore quadratico medio (MSE) e la radice dell’errore quadratico medio (RMSE) risultano leggermente più bassi nel training.
La differenza contenuta tra gli errori in training e test indica che il modello generalizza abbastanza bene, senza segni evidenti di overfitting.
Gli RMSE intorno a 400 grammi rappresentano la deviazione media tra i valori osservati e quelli previsti, suggerendo una buona capacità predittiva per questo tipo di dati.

In sintesi, il modello mostra una performance stabile su dati nuovi, confermando la sua affidabilità per la previsione del peso neonatale.

Previsione

Supponiamo di applicare il modello finale mod_inter per stimare il peso alla nascita di un neonato in un caso specifico:

Madre alla terza gravidanza
39 settimane di gestazione
Fumatrice
Sesso maschile
Parto naturale
Ricovero presso l’ospedale 1

new_neonate <- data.frame(
  Anni.madre = 30,          
  N.gravidanze = 3,
  Gestazione = 39,
  Fumatrici = 1,       
  Sesso = factor("M", levels = c("F", "M")),
  Tipo.parto = factor("Nat", levels = c("Nat", "Ces")),
  Ospedale = factor("osp1", levels = c("osp1", "osp2", "osp3"))
)

peso_previsto <- predict(mod_inter, newdata = new_neonate)
print(paste("Peso previsto (grammi):", round(peso_previsto, 1)))

## [1] "Peso previsto (grammi): 3326.7"

La previsione ottenuta indica un peso stimato di circa 3327 grammi. Questo valore rappresenta la stima del modello basata sulle variabili incluse, fornendo un esempio pratico di come il modello possa essere applicato per valutare il peso atteso di un neonato con determinate condizioni materne e cliniche

Effetto combinato tra Gestazione e Fumo

ggplot(neonati, aes(x = Gestazione, y = Peso, color = factor(Fumatrici))) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE) +
  facet_wrap(~Fumatrici) +
  scale_color_manual(
    name = "Fumatrici",
    values = c("0" = "blue", "1" = "red"),
    labels = c("Non Fumatrici", "Fumatrici")
  ) +
  labs(
    title = "Peso in funzione della Gestazione e Fumo Materno",
    x = "Settimane di Gestazione",
    y = "Peso del neonato (g)"
  ) +
  theme_minimal()

ggplot(neonati, aes(x = Gestazione, y = Peso, color = factor(Fumatrici), group = Fumatrici)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE) +
  labs(title = "Interazione: Fumo materno e Settimane di Gestazione",
       x = "Gestazione (settimane)", y = "Peso (g)", color = "Fumatrici") +
  theme_light()

I risultati mostrano che il fumo materno ha un effetto costante e persistente sul peso alla nascita durante tutta la gravidanza. La differenza di peso tra neonati di madri fumatrici e non fumatrici si mantiene relativamente stabile dall’inizio alla fine della gestazione.

Diagnostica residui

ggplot(data = data.frame(Fitted = fitted(mod_inter), Residuals = resid(mod_inter)), 
       aes(x = Fitted, y = Residuals)) +
  geom_point(alpha = 0.5) +
  geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
  labs(title = "Residui vs Valori Predetti", x = "Valori Predetti", y = "Residui") +
  theme_minimal()

Il grafico mostra in modo chiaro che l’assunzione di omoschedasticità non è rispettata. I punti non si distribuiscono in modo uniforme intorno alla linea orizzontale, ma formano un evidente pattern a imbuto: la varianza dei residui aumenta all’aumentare dei valori predetti. Questo indica che il modello commette errori crescenti per i neonati con peso maggiore, segnalando la possibilità di una varianza non costante che dovrebbe essere corretta. Inoltre, si osserva una forte concentrazione di dati nella fascia 3000-3500 grammi, mentre alcuni residui molto distanti, in particolare due sopra i 1500 e uno oltre i 1800, suggeriscono la presenza di outlier potenzialmente influenti.

ggplot(data = data.frame(Residuals = resid(mod_inter)), aes(sample = Residuals)) +
  stat_qq() +
  stat_qq_line(col = "blue") +
  labs(title = "Q-Q Plot dei residui") +
  theme_minimal()

Nel Q-Q plot i residui non seguono perfettamente la linea teorica, soprattutto agli estremi, rivelando code più pesanti di quanto previsto da una distribuzione normale. Questo implica che gli errori non si distribuiscono in modo perfettamente simmetrico e regolare, ma presentano deviazioni significative nelle code. La parte centrale dei punti aderisce bene alla linea, suggerendo che per la maggior parte dei casi la normalità è rispettata.

ggplot(data = data.frame(Residuals = resid(mod_inter)), aes(x = Residuals)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  labs(title = "Distribuzione dei Residui", x = "Residui", y = "Frequenza") +
  theme_minimal()

L’istogramma dei residui conferma quanto osservato negli altri grafici: la distribuzione appare approssimativamente normale, ma presenta irregolarità. Si nota una leggera asimmetria verso destra, con una coda destra più lunga che segnala la presenza di errori positivi particolarmente grandi. Nonostante la maggior parte dei residui sia concentrata intorno a zero, come auspicabile, alcuni valori estremi distorcono la forma ideale.

CONCLUSIONE

L’analisi condotta ha evidenziato come la durata della gestazione, il sesso del neonato e il numero di gravidanze siano i principali determinanti del peso alla nascita nel campione esaminato. Inoltre, l’interazione con la variabile relativa alle madri fumatrici suggerisce un effetto ulteriore sul peso neonatale. Il modello finale, pur spiegando una quota parziale della variabilità complessiva, fornisce indicazioni utili per comprendere le dinamiche cliniche alla base del peso neonatale.

Modello Statistico per la Previsione del Peso Neonatale

Matteo Zuliani

2025-07-01