Modello statistico per la previsione del peso neonatale

Analisi descrittiva, confronti clinici e regressione

1 Obiettivo e dati

1.1 Obiettivo dell’analisi

Il presente lavoro si propone di analizzare il peso alla nascita dei neonati registrati presso tre strutture ospedaliere, con particolare attenzione alle caratteristiche materne, neonatali e relative al parto che possono essere associate a tale misura.

L’analisi si articola in tre fasi principali. In primo luogo, viene effettuata un’esplorazione descrittiva dei dati, finalizzata a esaminarne la distribuzione e a individuare eventuali differenze tra le strutture ospedaliere e i gruppi considerati. Successivamente, viene sviluppata un’analisi delle relazioni tra il peso e le diverse caratteristiche disponibili. Infine, viene costruito un modello di regressione lineare multipla, con l’obiettivo di valutare il contributo simultaneo dei diversi predittori e di descrivere le relazioni tra questi e il peso.

Per facilitare il riconoscimento delle singole osservazioni durante le successive fasi di analisi, a ciascun record è stato assegnato un identificativo progressivo, denominato ID originale. Tale identificativo corrisponde alla posizione della riga nel dataset iniziale e consente di risalire facilmente all’osservazione corrispondente nel dataset di partenza, anche dopo eventuali esclusioni di record dalle analisi.

# Il file CSV deve trovarsi nella stessa cartella del documento Rmd.
dati_originali <- readr::read_csv("neonati.csv", show_col_types = FALSE)

var_numeriche <- c(
  "Anni.madre", "N.gravidanze", "Gestazione",
  "Peso", "Lunghezza", "Cranio"
)

etichette_var <- c(
  Anni.madre = "Età materna (anni)",
  N.gravidanze = "Gravidanze precedenti (numero)",
  Gestazione = "Età gestazionale (settimane)",
  Peso = "Peso (g)",
  Lunghezza = "Lunghezza (mm)",
  Cranio = "Circonferenza cranica (mm)"
)

id_record_escluso <- 1551L

dati <- dati_originali %>%
  mutate(
    ID_originale = row_number(),
    Fumatrici = factor(
      Fumatrici, levels = c(0, 1),
      labels = c("Non fumatrice", "Fumatrice")
    ),
    Tipo.parto = factor(
      Tipo.parto, levels = c("Nat", "Ces"),
      labels = c("Naturale", "Cesareo")
    ),
    Ospedale = factor(
      Ospedale, levels = c("osp1", "osp2", "osp3"),
      labels = c("Ospedale 1", "Ospedale 2", "Ospedale 3")
    ),
    Sesso = factor(
      Sesso, levels = c("F", "M"),
      labels = c("Femmina", "Maschio")
    )
  )

stopifnot(all(c(var_numeriche, "Fumatrici", "Tipo.parto",
                "Ospedale", "Sesso") %in% names(dati)))

1.2 Descrizione del dataset

Il dataset è costituito da 2500 osservazioni individuali relative a neonati appartenenti a tre strutture ospedaliere. Ciascuna osservazione rappresenta un singolo neonato e contiene informazioni fornite da una serie di variabili, che possono essere suddivise in tre gruppi:

• Caratteristiche materne: età della madre, numero di gravidanze precedenti e stato di fumo materno (madre fumatrice o non fumatrice).
• Caratteristiche neonatali: sesso, età gestazionale (in settimane), peso (in grammi), lunghezza e circonferenza cranica (in millimetri).
• Modalità del parto: tipologia di parto (naturale o cesareo) e ospedale di riferimento.

La variabile di interesse principale è il peso, che costituisce la variabile dipendente nell’analisi di regressione.

Per facilitare il controllo delle singole osservazioni e l’eventuale individuazione di valori anomali, a ciascun record viene assegnato un identificativo progressivo. Tale identificativo ha esclusivamente una funzione di tracciamento e non viene utilizzato come variabile esplicativa nei modelli statistici.

L’analisi ha carattere osservazionale: i risultati consentono pertanto di descrivere le associazioni presenti nei dati, senza che queste possano essere interpretate, in assenza di ulteriori elementi, come relazioni causali.

2 Qualità e distribuzione dei dati

Prima di tutto controlliamo se ci sono dati mancanti nel nostro database.

na.flag <- FALSE
for(data in dati_originali){
  if(any(is.na(data))){
    na.flag <- TRUE
    break
  }
}

Da un rapido controllo, segue che non ci sono dati mancanti.

Possiamo quindi proseguire cercando outliers particolarmente estremi. Per un primo controllo individuiamo massimi e minimi delle diverse variabili.

variables <- dati[, c("Anni.madre", "N.gravidanze", "Gestazione", "Peso",
                                         "Lunghezza", "Cranio")]

ind.pos <- data.frame(
  Min = sapply(variables, min),
  Max = sapply(variables, max)
)

rownames(ind.pos) <- c(
  "Età madre",
  "Gravidanze precedenti",
  "Età gestazionale",
  "Peso",
  "Lunghezza",
  "Circonferenza cranica"
)

knitr::kable(
  ind.pos,
  digits = 2,
  col.names = c("", "Minimo", "Massimo"))

	Minimo	Massimo
Età madre	0	46
Gravidanze precedenti	0	12
Età gestazionale	25	43
Peso	830	4930
Lunghezza	310	565
Circonferenza cranica	235	390

Possiamo notare un valore evidentemente anomalo per l’età della madre, ossia 0. Cerchiamo velocemente se ci sono altri valori particolarmente bassi, ad esempio inferiori a 12 anni, consapevoli che si tratti di una soglia estremamente improbabile ma non impossibile.

eta_anomale <- data.frame(
  ID_originale = numeric(),
  Anni = numeric(),
  stringsAsFactors = FALSE
)
for(id in dati$ID_originale){
  if(dati$Anni.madre[dati$ID_originale == id] <= 12){
    eta_anomale <- rbind(
      eta_anomale,
      data.frame(ID_originale = id,
                 Anni = dati$Anni.madre[dati$ID_originale == id])
    )
  }
}

knitr::kable(
  eta_anomale,
  digits = 2,
  col.names = c("ID originale", "Età madre"))

ID originale	Età madre
1152	1
1380	0

n_eta_anomale <- length(eta_anomale$Anni)

# Manteniamo tutte le altre informazioni del record, ma non usiamo questa età.
dati_analisi <- dati %>%
  mutate(Anni.madre = if_else(Anni.madre <= 1, NA_real_, as.numeric(Anni.madre)))

id_eta_ricodificate <- dati$ID_originale[
  !is.na(dati$Anni.madre) & dati$Anni.madre <= 1
]
min_eta <- min(dati_analisi$Anni.madre, na.rm = TRUE)

I 2 valori individuati (0 e 1 anni) non sono plausibili come età materne nel contesto di questo archivio. Li ricodifichiamo quindi come mancanti limitatamente alla variabile età della madre. Gli altri dati dei due record restano disponibili per le analisi che non richiedono l’età materna. Questa scelta non equivale a eliminare le due nascite dall’intero dataset.

Dopo la ricodifica, l’età materna minima osservata tra i valori utilizzabili è 13 anni. Anche questo valore è poco frequente e merita di essere letto come caratteristica del campione, non come riferimento clinico.

La tabella qui sotto riporta alcune statistiche descrittive delle variabili quantitative considerate nell’analisi.

# Per le statistiche descrittive usiamo anche qui le età già ricodificate
# come mancanti, in modo da non reinserire i due valori non plausibili.
variables <- dati_analisi[, c(
  "Anni.madre", "N.gravidanze", "Gestazione", "Peso",
  "Lunghezza", "Cranio"
)]

ind.pos <- data.frame(
    Media = sapply(variables, mean, na.rm = TRUE),
    CV = sapply(variables, sd, na.rm = TRUE)/sapply(variables, mean, na.rm = TRUE),
    Primo_Quart = sapply(variables, quantile, probs = 0.25, na.rm = TRUE),
    Mediana = sapply(variables, median, na.rm = TRUE),
    Terzo_Quart = sapply(variables, quantile, probs = 0.75, na.rm = TRUE)
)

rownames(ind.pos) <- c(
  "Età della madre",
  "Gravidanze precedenti",
  "Età gestazionale",
  "Peso",
  "Lunghezza",
  "Circonferenza cranica"
)

knitr::kable(
  ind.pos,
  digits = 2,
  col.names = c("", "Media", "Coeff. di variazione", "Primo quart.",
                "Mediana", "Terzo quart."),
  align = c("l", "r", "r", "r", "r", "r"))

	Media	Coeff. di variazione	Primo quart.	Mediana	Terzo quart.
Età della madre	28.19	0.19	25	28	32
Gravidanze precedenti	0.98	1.31	0	1	1
Età gestazionale	38.98	0.05	38	39	40
Peso	3284.08	0.16	2990	3300	3620
Lunghezza	494.69	0.05	480	500	510
Circonferenza cranica	340.03	0.05	330	340	350

Oltre alla media, sono indicate la mediana, ovvero il valore che divide le osservazioni ordinate in due metà uguali, e il primo e il terzo quartile, che individuano rispettivamente il valore al di sotto del quale ricade il 25% delle osservazioni e quello al di sotto del quale ricade il 75%. L’intervallo compreso tra questi due quartili contiene quindi il 50% centrale dei dati. È inoltre riportato il coefficiente di variazione, definito come rapporto tra deviazione standard e media, che misura la dispersione relativa dei valori rispetto alla media e permette di confrontare la variabilità di grandezze espresse in unità di misura differenti.
Una piccola precisazione: per il numero di gravidanze precedenti, essendo la media vicina a zero, il coefficiente di variazione è poco informativo e va interpretato con cautela.

Per descrivere la distribuzione delle variabili qualitative considerate, si utilizzano grafici a barre, che rappresentano il numero di osservazioni appartenenti a ciascuna categoria. L’altezza delle barre permette di confrontare immediatamente la frequenza delle diverse modalità e di individuare eventuali squilibri nella composizione del campione.

library(patchwork)
library(scales)

# 1. Età materna
p1 <- ggplot(dati_analisi, aes(x = Anni.madre)) +
  geom_histogram(
    binwidth = 1, boundary = 0.5,
    fill = "#477c8a", color = "white"
  ) +
  scale_x_continuous(breaks = breaks_width(2)) +
  labs(title = "Età materna", x = "Età (anni)", y = "Numero di nascite") +
  theme_minimal()

# 2. Gravidanze precedenti
p2 <- ggplot(dati_analisi, aes(x = N.gravidanze)) +
  geom_histogram(
    binwidth = 1, boundary = -0.5,
    fill = "#477c8a", color = "white"
  ) +
  scale_x_continuous(breaks = breaks_width(1)) +
  labs(title = "Gravidanze precedenti", x = "Numero di gravidanze", y = "Numero di nascite") +
  theme_minimal()

# 3. Età gestazionale
p3 <- ggplot(dati_analisi, aes(x = Gestazione)) +
  geom_histogram(
    binwidth = 1, boundary = 0.5,
    fill = "#477c8a", color = "white"
  ) +
  scale_x_continuous(breaks = breaks_width(2)) +
  labs(title = "Età gestazionale", x = "Settimane", y = "Numero di nascite") +
  theme_minimal()

# 4. peso
p4 <- ggplot(dati_analisi, aes(x = Peso)) +
  geom_histogram(
    binwidth = 200, boundary = 0,
    fill = "#477c8a", color = "white"
  ) +
  scale_x_continuous(breaks = breaks_width(500)) +
  labs(title = "Peso alla nascita", x = "Peso (g)", y = "Numero di nascite") +
  theme_minimal()

# 5. Lunghezza
p5 <- ggplot(dati_analisi, aes(x = Lunghezza)) +
  geom_histogram(
    binwidth = 10, boundary = 0,
    fill = "#477c8a", color = "white"
  ) +
  scale_x_continuous(breaks = breaks_width(40)) +
  labs(title = "Lunghezza", x = "Lunghezza (mm)", y = "Numero di nascite") +
  theme_minimal()

# 6. Circonferenza cranica
p6 <- ggplot(dati_analisi, aes(x = Cranio)) +
  geom_histogram(
    binwidth = 5, boundary = 0,
    fill = "#477c8a", color = "white"
  ) +
  scale_x_continuous(breaks = breaks_width(20)) +
  labs(title = "Circonferenza cranica", x = "Circonferenza cranica (mm)", y = "Numero di nascite") +
  theme_minimal()

# Composizione: 2 colonne e 3 righe
(p1 | p2) /
(p3 | p4) /
(p5 | p6)

Possiamo notare che le variabili età materna, peso, lunghezza e circonferenza cranica sono distribuite abbastanza simmetricamente, mentre quella delle gravidanze precedenti ha valori via via decrescenti, mentre la distribuzione dell’età gestazionale si concentra intorno alla 40ª settimana per poi decrescere bruscamente.

Per interpretare più correttamente la distribuzione del peso alla nascita, è opportuno tenere conto dell’età gestazionale, poiché il peso atteso varia sensibilmente in funzione delle settimane di gravidanza. In particolare, adottiamo quattro fasce:

• Meno di 32 settimane: raggruppa le nascite molto pretermine, per le quali il peso può essere nettamente inferiore a quello dei neonati più maturi. Manteniamo insieme queste settimane anche per evitare di creare sottogruppi troppo piccoli.
• Da 32 a 36 settimane: comprende le nascite pretermine successive alle 32 settimane e precedenti al termine.
• Da 37 a 41 settimane: comprende le nascite a termine secondo la definizione generale di gravidanza a termine. È una fascia ampia e al suo interno esistono differenze cliniche: la raggruppiamo qui per rendere più stabile il confronto descrittivo, non perché tutte queste settimane siano equivalenti.
• 42 settimane o più: identifica le nascite post-termine, che è utile esaminare separatamente.

classi_gestazione <- c("<32", "32–36", "37–41", "≥42")

dati_per_gestazione <- dati_analisi %>%
  mutate(
    Intervallo_gestazione = case_when(
      Gestazione < 32 ~ "<32",
      Gestazione <= 36 ~ "32–36",
      Gestazione <= 41 ~ "37–41",
      Gestazione >= 42 ~ "≥42"
    ),
    Intervallo_gestazione = factor(
      Intervallo_gestazione,
      levels = classi_gestazione
    )
  )


# Il conteggio per fascia aiuta a capire quanto sono numerosi i gruppi.
conteggi_gestazione <- dati_per_gestazione %>%
  count(Intervallo_gestazione, name = "N") %>%
  mutate(Percentuale = 100 * N / sum(N))

totale_con_gestazione <- sum(!is.na(dati_per_gestazione$Intervallo_gestazione))

tabella_per_gestazione <- dati_per_gestazione %>%
  filter(!is.na(Intervallo_gestazione)) %>%
  group_by(Intervallo_gestazione) %>%
  summarise(
    N = n(),
    Peso_Q1 = quantile(Peso, 0.25, na.rm = TRUE),
    Peso_mediana = median(Peso, na.rm = TRUE),
    Peso_Q3 = quantile(Peso, 0.75, na.rm = TRUE),
    Lunghezza_Q1 = quantile(Lunghezza, 0.25, na.rm = TRUE),
    Lunghezza_mediana = median(Lunghezza, na.rm = TRUE),
    Lunghezza_Q3 = quantile(Lunghezza, 0.75, na.rm = TRUE),
    Cranio_Q1 = quantile(Cranio, 0.25, na.rm = TRUE),
    Cranio_mediana = median(Cranio, na.rm = TRUE),
    Cranio_Q3 = quantile(Cranio, 0.75, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    `N (%)` = paste0(N, " (", formatta_numero(100 * N / totale_con_gestazione, 1), "%)"),
    `Peso mediano (Q1–Q3), g` = paste0(
      formatta_numero(Peso_mediana, 0), " (",
      formatta_numero(Peso_Q1, 0), "–", formatta_numero(Peso_Q3, 0), ")"
    ),
    `Lunghezza mediana (Q1–Q3), mm registrati` = paste0(
      formatta_numero(Lunghezza_mediana, 0), " (",
      formatta_numero(Lunghezza_Q1, 0), "–", formatta_numero(Lunghezza_Q3, 0), ")"
    ),
    `Circonferenza cranica mediana (Q1–Q3)` = paste0(
      formatta_numero(Cranio_mediana, 0), " (",
      formatta_numero(Cranio_Q1, 0), "–", formatta_numero(Cranio_Q3, 0), ")"
    )
  ) %>%
  select(
    `Fascia di età gestazionale (settimane)` = Intervallo_gestazione,
    `N (%)`, `Peso mediano (Q1–Q3), g`,
    `Lunghezza mediana (Q1–Q3), mm registrati`,
    `Circonferenza cranica mediana (Q1–Q3)`
  )

kable(
  tabella_per_gestazione,
  caption = "Misure neonatali descritte separatamente per fascia di età gestazionale"
)

Misure neonatali descritte separatamente per fascia di età gestazionale

Fascia di età gestazionale (settimane)	N (%)	Peso mediano (Q1–Q3), g	Lunghezza mediana (Q1–Q3), mm registrati	Circonferenza cranica mediana (Q1–Q3)
<32	24 (1,0%)	1.282 (1.102–1.638)	385 (355–410)	275 (267–296)
32–36	138 (5,5%)	2.560 (2.152–2.860)	460 (445–479)	322 (310–334)
37–41	2280 (91,2%)	3.310 (3.050–3.630)	500 (485–510)	340 (332–350)
≥42	58 (2,3%)	3.575 (3.310–3.892)	510 (500–524)	346 (336–360)

Notiamo che il gruppo 37–41 settimane comprende la maggior parte del campione, mentre i gruppi sotto le 32 settimane e da 42 settimane in poi sono molto più piccoli. In particolare, i gruppi meno numerosi offrono una base più limitata per confronti dettagliati.

Allo scopo di visualizzare meglio come le variabili sono distribuite, utilizziamo dei boxplot. La scatola si estende dal primo al terzo quartile e la sua ampiezza rappresenta l’intervallo interquartile (IQR), che contiene il 50% centrale delle osservazioni e misura la dispersione dei valori. La linea all’interno della scatola indica la mediana. I baffi rappresentano l’estensione dei valori non anomali secondo il criterio adottato, mentre gli eventuali punti esterni individuano osservazioni potenzialmente anomale, gli outlier. Questi ultimi sono precisamente i valori che si trovano più una volta e mezza la distanza interquartile al di sotto del primo quartile o al di sopra del terzo quartile.

dati_per_gestazione %>%
  filter(!is.na(Intervallo_gestazione)) %>%
  select(Intervallo_gestazione, Peso, Lunghezza, Cranio) %>%
  pivot_longer(
    cols = c(Peso, Lunghezza, Cranio),
    names_to = "Misura",
    values_to = "Valore"
  ) %>%
  mutate(
    Misura = recode(
      Misura,
      Peso = "Peso (g)",
      Lunghezza = "Lunghezza (mm registrati)",
      Cranio = "Circonferenza cranica"
    )
  ) %>%
  ggplot(aes(x = Intervallo_gestazione, y = Valore)) +
  geom_boxplot(fill = "#d8e8eb", color = "#477c8a", outlier.alpha = 0.45) +
  facet_wrap(~ Misura, scales = "free_y", ncol = 1) +
  labs(
    title = "Distribuzione delle misure neonatali nelle diverse fasce di età gestazionale",
    x = "Età gestazionale (settimane)",
    y = "Valore registrato"
  )

Mostra le tabelle degli outlier

Caratteristiche dei neonati con outlier — intervallo 37–41

ID record	Età gestazionale (settimane)	Circonferenza cranica	Peso (g)	Lunghezza (mm)
34	39	382	3150	480
119	40	372	3410	550
126	41	370	4680	545
151	41	280	2280	450
165	40	360	4350	550
190	39	390	4050	525
205	38	384	3850	505
220	40	363	3520	445
241	41	355	4070	550
271	41	350	4140	550
288	38	332	2480	446
295	40	305	1850	460
304	37	326	2060	420
317	40	352	4000	550
322	37	305	1750	430
328	41	340	3900	550
329	40	340	4560	540
340	40	382	4150	520
353	39	300	2680	460
368	41	370	4600	540
383	39	346	3600	550
428	40	360	4400	550
458	38	315	2150	450
471	40	346	3680	560
532	41	380	3770	520
565	40	386	4250	520
600	41	370	4120	550
666	40	345	4240	555
684	39	390	3000	475
697	39	300	2820	510
712	38	300	2430	473
728	37	330	2340	440
729	39	355	3900	555
954	38	323	2430	445
956	41	310	2210	430
958	40	354	4240	550
964	40	335	4010	550
988	40	365	4150	550
991	40	385	4050	550
1014	37	295	2050	390
1113	40	380	4420	520
1157	40	380	4060	505
1188	40	320	4140	550
1194	39	325	2160	450
1215	40	355	3630	550
1222	40	372	4330	550
1275	38	305	1980	430
1283	40	376	4250	550
1293	38	380	4600	485
1306	41	352	4900	510
1321	40	383	4340	515
1323	39	355	3660	555
1356	40	390	4090	525
1395	39	304	3790	505
1420	38	380	3770	530
1433	41	364	4810	530
1510	37	305	2040	465
1513	40	367	4620	560
1541	38	343	4540	530
1551	38	374	4370	315
1556	41	300	2420	490
1577	37	310	2060	445
1634	38	328	2500	430
1635	39	322	3430	445
1639	40	365	4760	550
1644	41	360	3850	550
1712	39	300	3800	520
1735	37	345	2750	440
1743	38	305	1950	445
1745	40	378	4310	525
1775	40	380	4480	530
1826	41	360	4270	550
1838	40	360	4580	515
1844	37	310	2090	440
1868	40	390	3470	525
1920	39	350	4930	550
1960	37	345	2990	445
2006	38	309	2180	445
2023	39	354	4650	510
2040	38	359	3240	410
2076	40	360	4720	540
2213	37	299	2610	450
2215	40	298	2440	465
2235	41	373	4690	540
2273	38	300	2450	470
2306	37	315	2500	440
2391	41	366	4400	565
2392	40	355	4720	540
2417	40	378	4000	520
2438	39	315	2200	440

Caratteristiche dei neonati con outlier — intervallo 32–36

ID record	Età gestazionale (settimane)	Circonferenza cranica	Peso (g)	Lunghezza (mm)
101	34	287	1370	390
312	32	276	1280	360
492	33	295	1410	380
748	33	277	1390	390
947	32	297	1615	390
1181	36	373	4070	500
1428	36	292	1280	385
1553	35	360	4520	520
1701	32	301	1430	380

Caratteristiche dei neonati con outlier — intervallo <32

ID record	Età gestazionale (settimane)	Circonferenza cranica	Peso (g)	Lunghezza (mm)
310	28	379	1560	420
1429	29	355	1280	390

Tra questi selezioniamo poi solo i valori estremi che distano ben tre volte la distanza interquartile dal primo e dal terzo quartile.

Caratteristiche dei neonati con valori estremi — intervallo 37–41

ID record	Età gestazionale (settimane)	Circonferenza cranica	Peso (g)	Lunghezza (mm)
1014	37	295	2050	390
1551	38	374	4370	315

Notiamo subito che la combinazione relativa al record 1551 è fortemente discordante: la lunghezza è insolitamente bassa rispetto alle altre caratteristiche antropometriche registrate. Il dato è quindi da considerare sospetto.
Per evitare che un record potenzialmente compromesso condizioni le analisi successive, decidiamo pertanto di escludere l’intera osservazione 1551 dal dataset analitico.

stopifnot(sum(dati_analisi$ID_originale == id_record_escluso) == 1)

dati_analisi <- dati_analisi %>%
  filter(ID_originale != id_record_escluso)

stopifnot(!any(dati_analisi$ID_originale == id_record_escluso))

3 Tre domande descrittive

In questa sezione vengono affrontate alcune domande di interesse clinico attraverso opportuni strumenti di statistica inferenziale. In particolare, si vuole verificare se la frequenza dei parti cesarei differisca tra i tre ospedali, se il peso medio e la lunghezza media dei neonati osservati siano compatibili con i rispettivi valori di riferimento della popolazione italiana (3300 g e 50 cm) e se le misure antropometriche differiscano tra neonati di sesso maschile e femminile.

Per rispondere a queste domande, ciascuna analisi sarà preceduta dalla formulazione di un sistema di ipotesi, costituito da un’ipotesi nulla \(\mathrm{H}_0\) e da un’ipotesi alternativa \(\mathrm{H}_1\). La scelta del test statistico dipende dalla natura delle variabili e dal confronto che si intende effettuare: si potranno quindi confrontare una media con un valore di riferimento, le medie di due gruppi oppure verificare l’indipendenza tra due variabili categoriali.

Un test di ipotesi consente di valutare quanto i dati osservati siano compatibili con l’ipotesi nulla. A questo scopo si utilizza il p-value, che rappresenta la probabilità, supponendo vera \(\mathrm{H}_0\), di ottenere un risultato almeno altrettanto estremo di quello osservato. Un p-value piccolo indica che i dati sarebbero poco compatibili con l’ipotesi nulla e costituiscono quindi un’evidenza a suo sfavore. Convenzionalmente, si adotta una soglia di significatività pari a 0,05: quando il p-value è inferiore a tale valore, si rifiuta l’ipotesi nulla. Un’altra informazione riportata nelle tabelle è l’intervallo di confidenza (IC), generalmente calcolato al 95%. Esso fornisce un intervallo di valori compatibili con i dati osservati per il parametro che si intende stimare, tenendo conto della variabilità campionaria.

È importante ricordare che il p-value non esprime la probabilità che l’ipotesi nulla sia vera e non misura l’entità o la rilevanza clinica di una differenza. Un risultato statisticamente significativo, infatti, non è necessariamente rilevante dal punto di vista clinico; analogamente, un risultato non significativo non dimostra che una differenza o un’associazione siano assenti.

3.1 Il ricorso al cesareo varia tra ospedali?

Confrontiamo la percentuale di cesarei nei tre ospedali. Poiché entrambe le variabili sono qualitative e vogliamo verificare se la distribuzione del tipo di parto cambia tra tre strutture, utilizziamo un test chi-quadrato di indipendenza su una tabella di frequenze.

Le ipotesi sono:

• \(\mathrm{H}_0\): ospedale e tipo di parto sono indipendenti; la distribuzione dei parti naturali e cesarei è la stessa nei tre ospedali.
• \(\mathrm{H}_1\): ospedale e tipo di parto non sono indipendenti; almeno un ospedale presenta una distribuzione diversa.

Il test valuta quindi un’associazione complessiva: da solo non identifica quali ospedali differiscano tra loro.

tabella_parto <- with(dati_analisi, table(Ospedale, Tipo.parto))
test_chi <- chisq.test(tabella_parto, correct = FALSE)

v_cramer <- sqrt(
  as.numeric(test_chi$statistic) /
    (sum(tabella_parto) * (min(dim(tabella_parto)) - 1))
)

riepilogo_cesarei <- dati_analisi %>%
  group_by(Ospedale) %>%
  summarise(
    N = n(),
    Cesarei = sum(Tipo.parto == "Cesareo", na.rm = TRUE),
    Percentuale_cesarei = 100 * mean(Tipo.parto == "Cesareo", na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(Percentuale_cesarei = round(Percentuale_cesarei, 1))

knitr::kable(
  riepilogo_cesarei,
  col.names = c("Ospedale", "Parti", "Cesarei", "Cesarei (%)"),
  caption = "Numero e percentuale di cesarei per ospedale"
)

Numero e percentuale di cesarei per ospedale

Ospedale	Parti	Cesarei	Cesarei (%)
Ospedale 1	816	242	29.7
Ospedale 2	849	254	29.9
Ospedale 3	834	232	27.8

Il test chi-quadrato restituisce un p-value pari a 0,588: non emerge evidenza statistica di un’associazione tra ospedale e tipo di parto. La V di Cramer, pari a 0,02, quantifica l’intensità dell’associazione osservata. Il valore, prossimo a zero, indica che l’associazione tra le due variabili è molto debole. In ogni caso, l’assenza di significatività non dimostra che i tassi siano equivalenti, né il confronto consente di valutare l’appropriatezza dei singoli cesarei o la qualità dei reparti.

3.2 Il peso medio è compatibile con 3.300 g? E la lunghezza con 50 cm?

Confrontiamo le medie campionarie con i valori di riferimento indicati, utilizzando un test t di Student a un campione. Questo test è appropriato quando non si conosce la deviazione standard della popolazione, che viene pertanto stimata a partire dal campione. La statistica test, ottenuta rapportando la differenza tra la media osservata e il valore di riferimento al relativo errore standard, segue sotto l’ipotesi nulla una distribuzione t di Student.

Il p-value rappresenta la probabilità, supponendo vera l’ipotesi nulla, di ottenere una statistica test almeno altrettanto distante da zero di quella osservata. Trattandosi di un test bilaterale, tale probabilità corrisponde all’area complessiva delle due code della distribuzione, oltre i valori della statistica test osservata e del suo opposto.

Per il peso le ipotesi sono:

• \(\mathrm{H}_0\): la media del peso alla nascita nel campione di riferimento è 3.300 g (\(\mu_{peso}=3300\)).
• \(\mathrm{H}_1\): la media del peso alla nascita è diversa da 3.300 g (\(\mu_{peso}\neq3300\)).

Per la lunghezza le ipotesi sono:

• \(\mathrm{H}_0\): la media della lunghezza è 50 cm (\(\mu_{lunghezza}=50\)).
• \(\mathrm{H}_1\): la media della lunghezza è diversa da 50 cm (\(\mu_{lunghezza}\neq50\)).

test_peso <- t.test(dati_analisi$Peso, mu = 3300)
test_lunghezza <- t.test(dati_analisi$Lunghezza / 10, mu = 50)

tab_confronto <- tibble(
  Misura = c("Peso (g)", "Lunghezza (cm)"),
  Media_campione = c(
    mean(dati_analisi$Peso, na.rm = TRUE),
    mean(dati_analisi$Lunghezza / 10, na.rm = TRUE)
  ),
  Riferimento = c(3300, 50),
  Differenza = c(
    mean(dati_analisi$Peso, na.rm = TRUE) - 3300,
    mean(dati_analisi$Lunghezza / 10, na.rm = TRUE) - 50
  ),
  IC95_inferiore = c(
    test_peso$conf.int[1] - 3300,
    test_lunghezza$conf.int[1] - 50
  ),
  IC95_superiore = c(
    test_peso$conf.int[2] - 3300,
    test_lunghezza$conf.int[2] - 50
  ),
  p = c(test_peso$p.value, test_lunghezza$p.value)
) %>%
  mutate(
    across(
      c(Media_campione, Riferimento, Differenza,
        IC95_inferiore, IC95_superiore),
      ~ round(.x, 2)
    ),
    `p-value` = vapply(p, formatta_p, character(1))
  ) %>%
  select(-p)

knitr::kable(
  tab_confronto,
  caption = "Confronto con i valori di riferimento",
  col.names = c("Misura", "Media campione", "Riferimento",
                "Differenza", "IC95% inferiore", "IC95% superiore",
                "p-value")
)

Confronto con i valori di riferimento

Misura	Media campione	Riferimento	Differenza	IC95% inferiore	IC95% superiore	p-value
Peso (g)	3283.65	3300	-16.35	-36.94	4.23	0,119
Lunghezza (cm)	49.48	50	-0.52	-0.63	-0.42	< 0,001

Per il peso, essendo il p-value maggiore di 0,05, non emerge evidenza statistica di una differenza rispetto a 3.300 g.

Per la lunghezza, invece, lo scostamento è statisticamente rilevabile, ma la sua rilevanza clinica va valutata alla luce del contesto e del valore di riferimento adottato.

3.3 Le misure differiscono tra maschi e femmine?

Confrontiamo due gruppi indipendenti (maschi e femmine) per ciascuna delle tre misure. Lo strumento appropriato è il test t di Welch, che non assume l’uguaglianza delle varianze. Definiamo la differenza come media nei maschi meno media nelle femmine.

Per ciascuna misura (peso, lunghezza e circonferenza cranica) formuliamo le stesse ipotesi, applicate alla rispettiva media:

• \(\mathrm{H}_0\): la media nei maschi è uguale alla media nelle femmine (\(\mu_M-\mu_F=0\)).
• \(\mathrm{H}_1\): le due medie sono diverse (\(\mu_M-\mu_F\neq0\)).

dati_F <- filter(dati_analisi, Sesso == "Femmina")
dati_M <- filter(dati_analisi, Sesso == "Maschio")

test_peso_mf <- t.test(dati_M$Peso, dati_F$Peso)
test_lunghezza_mf <- t.test(dati_M$Lunghezza, dati_F$Lunghezza)
test_cranio_mf <- t.test(dati_M$Cranio, dati_F$Cranio)

medie_M <- c(
  mean(dati_M$Peso, na.rm = TRUE),
  mean(dati_M$Lunghezza / 10, na.rm = TRUE),
  mean(dati_M$Cranio / 10, na.rm = TRUE)
)
medie_F <- c(
  mean(dati_F$Peso, na.rm = TRUE),
  mean(dati_F$Lunghezza / 10, na.rm = TRUE),
  mean(dati_F$Cranio / 10, na.rm = TRUE)
)

tab_mf <- tibble(
  Misura = c("Peso (g)", "Lunghezza (cm)",
             "Circonferenza cranica (cm)"),
  Media_M = medie_M,
  Media_F = medie_F,
  Differenza_M_F = medie_M - medie_F,
  IC95_inferiore = c(
    test_peso_mf$conf.int[1],
    test_lunghezza_mf$conf.int[1] / 10,
    test_cranio_mf$conf.int[1] / 10
  ),
  IC95_superiore = c(
    test_peso_mf$conf.int[2],
    test_lunghezza_mf$conf.int[2] / 10,
    test_cranio_mf$conf.int[2] / 10
  ),
  p = c(
    test_peso_mf$p.value,
    test_lunghezza_mf$p.value,
    test_cranio_mf$p.value
  )
) %>%
  mutate(
    across(
      c(Media_M, Media_F, Differenza_M_F,
        IC95_inferiore, IC95_superiore),
      ~ round(.x, 2)
    ),
    `p-value` = vapply(p, formatta_p, character(1))
  ) %>%
  select(-p)

knitr::kable(
  tab_mf,
  caption = "Confronto tra le misure di maschi e femmine",
  col.names = c("Misura", "Media maschi", "Media femmine",
                "Differenza M−F", "IC95% inferiore",
                "IC95% superiore", "p-value")
)

Confronto tra le misure di maschi e femmine

Misura	Media maschi	Media femmine	Differenza M−F	IC95% inferiore	IC95% superiore	p-value
Peso (g)	3408.22	3160.17	248.05	208.05	288.04	< 0,001
Lunghezza (cm)	49.97	48.99	0.98	0.78	1.18	< 0,001
Circonferenza cranica (cm)	34.24	33.76	0.48	0.36	0.61	< 0,001

Tutte e tre le differenze risultano statisticamente significative nel confronto non aggiustato. Questi confronti non correggono per età gestazionale o per le altre caratteristiche: l’eventuale coefficiente del sesso nel modello multivariato risponde a una domanda diversa, confrontando neonati a parità degli altri predittori.

4 Modello di regressione del peso

La regressione lineare multipla è una tecnica statistica che permette di studiare la relazione tra una variabile di interesse, detta variabile dipendente, e un insieme di altre variabili, dette variabili indipendenti o predittori. Nel nostro caso, l’obiettivo è analizzare come le caratteristiche materne e neonatali, insieme ad alcune informazioni relative al parto, siano associate al peso alla nascita.

Il modello esprime il peso come combinazione lineare dei diversi predittori:

\[ Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \dots + \beta_k X_{ik} + \varepsilon_i, \]

dove \(Y_i\) rappresenta il peso del neonato \(i\)-esimo, \(X_{ij}\) il valore assunto dal predittore \(j\)-esimo per la stessa osservazione, \(\beta_0\) l’intercetta, \(\beta_j\) i coefficienti di regressione ed \(\varepsilon_i\) la componente di errore, che raccoglie la parte della variabilità non spiegata dal modello.

I coefficienti di regressione descrivono come varia mediamente il peso previsto dal modello al variare di ciascun predittore, mantenendo costanti tutti gli altri. Per una variabile quantitativa, il coefficiente indica la variazione attesa del peso associata a un incremento unitario di quella variabile. Per una variabile qualitativa, invece, il coefficiente esprime la differenza rispetto a una categoria di riferimento, a parità degli altri predittori. Il segno del coefficiente indica la direzione della relazione: positivo se all’aumentare del predittore aumenta il peso previsto, negativo nel caso opposto.

Per valutare la rilevanza statistica dei singoli coefficienti si considera il p-value, calcolato nell’ambito di un test d’ipotesi che verifica se il coefficiente della variabile considerata sia significativamente diverso da zero. In particolare, l’ipotesi nulla afferma che, mantenendo costanti gli altri predittori, non vi sia un’associazione lineare tra quella variabile e il peso.

Nelle tabelle dei risultati compaiono anche l’errore standard e la statistica t. L’errore standard misura l’incertezza associata alla stima di ciascun coefficiente: valori più piccoli indicano stime generalmente più precise. La statistica t si ottiene dividendo il coefficiente stimato per il suo errore standard e misura quanto la stima si discosta da zero in rapporto alla propria incertezza. Il p-value viene quindi ricavato dalla distribuzione t di Student e permette di valutare la compatibilità dei risultati con l’ipotesi nulla.

È importante sottolineare che la significatività statistica non coincide necessariamente con la rilevanza clinica: un’associazione statisticamente significativa non è automaticamente importante dal punto di vista medico, così come un risultato non significativo non dimostra l’assenza di una relazione.

Per valutare la capacità del modello di descrivere i dati osservati si considera innanzitutto il coefficiente di determinazione \(R^2\). Esso rappresenta la proporzione della variabilità della variabile dipendente che viene spiegata dal modello. Assume valori compresi tra 0 e 1: valori più elevati indicano che il modello riesce a spiegare una quota maggiore della variabilità osservata. Tuttavia, l’aggiunta di nuovi predittori non può diminuire il valore di \(R^2\), anche quando le variabili aggiunte apportano un contributo limitato. Per questo motivo, nella valutazione di modelli con un numero diverso di predittori è utile considerare anche il coefficiente di determinazione corretto (\(R^2\) aggiustato), che tiene conto del numero di variabili incluse.

Per confrontare modelli alternativi si possono inoltre utilizzare due criteri di selezione: l’Akaike Information Criterion (AIC) e il Bayesian Information Criterion (BIC). Entrambi valutano il compromesso tra la capacità del modello di adattarsi ai dati e la sua complessità, introducendo una penalizzazione per il numero di parametri stimati. In questo modo, un modello più complesso viene preferito soltanto se il miglioramento dell’adattamento è sufficientemente consistente da compensare la penalizzazione.

I due criteri differiscono principalmente per l’entità della penalizzazione: il BIC penalizza maggiormente i modelli più complessi, soprattutto quando la numerosità campionaria è elevata. In entrambi i casi, valori più bassi indicano modelli preferibili rispetto a quelli con cui vengono confrontati.

Nel presente lavoro si è scelto di utilizzare il BIC come criterio principale per la selezione del modello, privilegiando una soluzione relativamente parsimoniosa, capace di descrivere adeguatamente i dati senza introdurre un numero eccessivo di predittori e termini aggiuntivi. Questa scelta mira a limitare il rischio di sovradattamento (overfitting) e a favorire un modello più semplice e interpretabile. Il valore di \(R^2\) e quello di \(R^2\) aggiustato vengono invece utilizzati come indicatori complementari della capacità esplicativa del modello.

Infine, un ulteriore aspetto da considerare nella costruzione del modello è la multicollinearità, che si verifica quando alcuni predittori sono fortemente correlati tra loro. In questa situazione può diventare difficile distinguere con precisione il contributo dei singoli termini e gli errori standard dei relativi coefficienti possono aumentare. Per valutarla utilizzeremo il Variance Inflation Factor (VIF), definito per ciascun termine come

\[ \mathrm{VIF}_j = \frac{1}{1-R_j^2}, \]

dove \(R_j^2\) è il coefficiente di determinazione ottenuto regredendo il termine considerato su tutti gli altri predittori del modello. Un VIF pari a 1 indica assenza di collinearità lineare con gli altri termini. Come regole orientative, valori inferiori a 5 sono spesso considerati contenuti, valori superiori a 5 suggeriscono di approfondire la situazione e valori superiori a 10 indicano una criticità più marcata. Queste soglie non sono criteri automatici di esclusione: il VIF va interpretato insieme alla motivazione teorica del termine, alla precisione delle stime e alle prestazioni complessive del modello.

Per i termini che richiedono più coefficienti, invece, viene utilizzato il Generalized Variance Inflation Factor (GVIF), che valuta l’inflazione dell’incertezza associata congiuntamente all’intero gruppo di coefficienti, tenendo conto anche delle loro correlazioni reciproche. Per questi ultimi termini viene inoltre riportata la quantità \(GVIF^{1/(2Df)}\), che consente di rendere i valori confrontabili tra termini con un diverso numero di gradi di libertà.

4.1 Costruzione di un campione comune

Dopo l’esclusione del record 1551, i record con età materna ricodificata come mancante, ossia i record 1152 e 1380 (vedi Sezione 2), vengono esclusi dal campione utilizzato per la regressione. In questo modo tutti i modelli candidati sono stimati sulle medesime osservazioni, rendendo confrontabili i rispettivi criteri di selezione. Sul campione così ottenuto, riassegniamo a ID record una numerazione progressiva, coincidente con l’indice delle osservazioni nel modello. La colonna ID originale conserva invece l’identificativo del record dataset iniziale, consentendo di risalire ai dati di partenza.

Le variabili originali vengono mantenute per le analisi descrittive e per rappresentare i grafici nelle loro unità di misura. Nella regressione utilizzeremo anche le versioni centrate di età gestazionale, lunghezza e circonferenza cranica, ottenute sottraendo a ciascuna variabile la propria media campionaria. Questa trasformazione rende più agevole l’interpretazione dei termini di ordine superiore e può ridurre la collinearità non essenziale tra effetti lineari, termini quadratici e interazioni.

predittori_candidati <- c(
  "Anni.madre", "N.gravidanze", "Fumatrici", "Gestazione",
  "Lunghezza", "Cranio", "Tipo.parto", "Ospedale", "Sesso"
)

dati_modello <- dati_analisi %>%
  # Conserviamo l'identificativo del CSV per poter risalire al record originale.
  select(ID_originale, Peso, all_of(predittori_candidati)) %>%
  drop_na() %>%
  droplevels() %>%
  # Nel campione completo del modello, ID_record coincide con l'indice di riga.
  mutate(ID_record = row_number()) %>%
  relocate(ID_record, ID_originale)

# Centriamo le tre misure continue usate nei termini quadratici e nelle interazioni.
# I centri sono calcolati sul campione completo comune.
centri <- c(
  Gestazione = mean(dati_modello$Gestazione),
  Lunghezza = mean(dati_modello$Lunghezza),
  Cranio = mean(dati_modello$Cranio)
)

dati_modello <- dati_modello %>%
  mutate(
    Gestazione_c = Gestazione - centri[["Gestazione"]],
    Lunghezza_c = Lunghezza - centri[["Lunghezza"]],
    Cranio_c = Cranio - centri[["Cranio"]]
  )

# Applica ai nuovi profili le stesse centrature usate nella stima.
centra_predittori <- function(nuovi_dati) {
  nuovi_dati %>%
    mutate(
      Gestazione_c = Gestazione - centri[["Gestazione"]],
      Lunghezza_c = Lunghezza - centri[["Lunghezza"]],
      Cranio_c = Cranio - centri[["Cranio"]]
    )
}

# Verifica: l'ID del campione coincide esattamente con l'indice di riga.
stopifnot(identical(dati_modello$ID_record, seq_len(nrow(dati_modello))))

n_esclusi <- nrow(dati_analisi) - nrow(dati_modello)
id_esclusi_modello <- setdiff(
  dati_analisi$ID_originale,
  dati_modello$ID_originale
)

stopifnot(!any(dati_modello$ID_originale == id_record_escluso))

stopifnot(
  nrow(dati_modello) > 0,
  nlevels(dati_modello$Fumatrici) >= 2,
  nlevels(dati_modello$Sesso) >= 2,
  nlevels(dati_modello$Tipo.parto) >= 2,
  nlevels(dati_modello$Ospedale) >= 2
)

4.2 Esplorazione delle relazioni

Prima di procedere alla regressione, esaminiamo le relazioni tra il peso alla nascita e le altre variabili quantitative, nonché le relazioni tra i predittori stessi. A questo scopo, utilizziamo due rappresentazioni complementari.

In primo luogo, attraverso diagrammi di dispersione e curve di smoothing, analizziamo la relazione tra il peso e ciascun predittore quantitativo, così da osservarne l’andamento e individuare eventuali tendenze non lineari. Successivamente, esaminiamo la matrice delle correlazioni di Pearson tra le variabili quantitative, che consente di valutare la forza e la direzione delle relazioni lineari, mettendo in evidenza anche eventuali sovrapposizioni informative tra predittori (per esempio tra lunghezza e circonferenza cranica).

Le correlazioni non implicano un rapporto causale e sono calcolate sul campione comune della regressione, così da riferirsi alle stesse osservazioni utilizzate nella successiva stima dei modelli.

variabili_correlazione <- c(
  "Peso", "Anni.madre", "N.gravidanze", "Gestazione",
  "Lunghezza", "Cranio"
)

# Scatterplot del peso rispetto a ciascun predittore quantitativo.
dati_modello %>%
  select(all_of(variabili_correlazione)) %>%
  pivot_longer(-Peso, names_to = "Predittore", values_to = "Valore") %>%
  mutate(Predittore = unname(etichette_var[Predittore])) %>%
  ggplot(aes(x = Valore, y = Peso)) +
  geom_point(alpha = 0.25, color = "#477c8a") +
  geom_smooth(method = "loess", se = TRUE,
              color = "#bd6b42", linewidth = 0.9) +
  facet_wrap(~ Predittore, scales = "free_x", ncol = 2) +
  labs(
    title = "Peso e predittori quantitativi",
    x = "Valore del predittore", y = "Peso (g)"
  )

# Creazione delle variabili dummy.
dati_cor <- dati_modello %>%
  select(
    -ID_originale,
    -ID_record,
    -Cranio_c,
    -Lunghezza_c,
    -Gestazione_c,
    -Sesso,
    -Fumatrici,
    -Tipo.parto,
    -Ospedale
  ) %>%
  mutate(
    Maschio = as.integer(as.character(dati_modello$Sesso) == "Maschio"),
    Fumatrice = as.integer(as.character(dati_modello$Fumatrici) == "Fumatrice"),
    Cesareo = as.integer(as.character(dati_modello$Tipo.parto) == "Cesareo"),
    Ospedale_2 = as.integer(as.character(dati_modello$Ospedale) == "Ospedale 2"),
    Ospedale_3 = as.integer(as.character(dati_modello$Ospedale) == "Ospedale 3")
  )

# Matrice delle correlazioni.
matrice_cor <- cor(
  dati_cor,
  use = "pairwise.complete.obs",
  method = "pearson"
)

# Recupero delle etichette, usando il nome originale se manca la descrizione.
etichette_cor <- etichette_var[rownames(matrice_cor)]

etichette_cor[is.na(etichette_cor)] <-
  rownames(matrice_cor)[is.na(etichette_cor)]

etichette_cor <- unname(etichette_cor)

rownames(matrice_cor) <- etichette_cor
colnames(matrice_cor) <- etichette_cor

# Solo triangolo inferiore, diagonale esclusa.
matrice_cor_visibile <- matrice_cor
matrice_cor_visibile[
  upper.tri(matrice_cor_visibile, diag = TRUE)
] <- NA_real_

# Formato long.
cor_long <- as.data.frame(as.table(matrice_cor_visibile)) %>%
  rename(
    Variabile_1 = Var1,
    Variabile_2 = Var2,
    Correlazione = Freq
  ) %>%
  filter(!is.na(Correlazione))

# Grafico.
grafico_cor <- ggplot(
  cor_long,
  aes(x = Variabile_2, y = Variabile_1, fill = Correlazione)
) +
  geom_tile(color = "white", linewidth = 0.5) +
  geom_text(
    aes(label = sprintf("%.3f", Correlazione)),
    size = 3.2
  ) +
  scale_fill_gradient2(
    low = "#B35858",
    mid = "white",
    high = "#756BB1",
    midpoint = 0,
    limits = c(-1, 1)
  ) +
  scale_x_discrete(position = "bottom") +
  scale_y_discrete(limits = rev(etichette_cor[-1])) +
  coord_equal() +
  labs(
    title = "Correlazioni tra le variabili",
    x = NULL,
    y = NULL,
    fill = "Pearson r"
  ) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    axis.text.x = element_text(
      angle = 40,
      hjust = 1,
      size = 9
    ),
    axis.text.y = element_text(size = 9),
    plot.title = element_text(
      hjust = 0.5,
      face = "bold"
    ),
    legend.position = "right"
  )

grafico_cor

Dall’analisi congiunta dei grafici emerge una forte associazione positiva tra il peso alla nascita e tre variabili: età gestazionale, circonferenza cranica e lunghezza, in particolare quest’ultima. Le curve LOESS suggeriscono tuttavia possibili deviazioni da un andamento lineare, soprattutto nel caso dell’età gestazionale e, in misura minore, della circonferenza cranica per i valori più elevati. Anche la curva relativa alla lunghezza sembra discostarsi dall’andamento lineare per valori molto piccoli, ma il numero limitato di osservazioni in questa regione non consente di trarre conclusioni affidabili.

Si osserva inoltre una modesta correlazione positiva tra il peso alla nascita e il sesso maschile, mentre la correlazione con il fumo materno è negativa, ma di intensità molto debole.

Per quanto riguarda le relazioni tra i predittori, si osservano forti associazioni positive tra età gestazionale, lunghezza e circonferenza cranica. Sebbene queste variabili condividano quindi parte dell’informazione, le correlazioni osservate non raggiungono valori particolarmente elevati, tali da far ipotizzare, sulla sola base dell’analisi bivariata, una forte multicollinearità. È comunque opportuno verificare questo aspetto anche nella successiva fase di costruzione e diagnostica del modello, poiché le correlazioni a coppie non sono sufficienti a escludere eventuali problemi di collinearità.

Si può infine notare una dispersione piuttosto ampia dei valori osservati attorno alle tendenze individuate. Ciò suggerisce che, pur essendo alcune variabili fortemente associate al peso, nessuna di esse, considerata isolatamente, ne descriva completamente la variabilità. Risulta pertanto opportuno considerare congiuntamente i diversi predittori attraverso un modello di regressione multipla.

4.3 Modello iniziale e selezione dei predittori

Stimiamo innanzitutto un modello che include tutti i predittori candidati. Ricordiamo che le variabili età gestazionale, lunghezza e circonferenza cranica sono state centrate per ridurre la collinearità.

Come già accennato all’inizio della sezione, le variabili qualitative vengono rappresentate nel modello mediante variabili indicatrici (dummy). In particolare, a una variabile qualitativa con (n) modalità corrispondono (n-1) variabili indicatrici, ciascuna delle quali assume valore 1 quando l’osservazione appartiene a una determinata modalità e 0 altrimenti. La modalità non rappresentata costituisce la categoria di riferimento. Nel nostro caso, otteniamo le variabili «fumatrice», «maschio», «cesareo», «ospedale 2» e «ospedale 3».

formula_completa <- Peso ~ Anni.madre + N.gravidanze + Fumatrici +
  Gestazione_c + Lunghezza_c + Cranio_c + Tipo.parto + Ospedale + Sesso

mod1 <- lm(formula_completa, data = dati_modello)

mod2 <- step(
  mod1,
  direction = "backward",
  scope = list(
    lower = Peso ~ Fumatrici,
    upper = formula(mod1)
  ),
  k = log(nobs(mod1)),
  trace = 0
)

vif_valori <- function(modello) {
  valori <- car::vif(modello)

  if (is.matrix(valori)) {
    tibble::tibble(
      termine_vif = rownames(valori),
      VIF = ifelse(
        valori[, "Df"] == 1,
        valori[, "GVIF"],
        valori[, "GVIF"]^(1 / (2 * valori[, "Df"]))
      )
    )
  } else {
    tibble::tibble(
      termine_vif = names(valori),
      VIF = as.numeric(valori)
    )
  }
}

vif_valore <- function(modello, termine) {
  vif_valori(modello) %>%
    filter(termine_vif == termine) %>%
    pull(VIF)
}

tabella_coefficienti <- function(modello) {
  vif_modello <- vif_valori(modello)

  broom::tidy(modello) %>%
    filter(term != "(Intercept)") %>%
    mutate(
      termine_vif = case_when(
        startsWith(term, "Fumatrici") ~ "Fumatrici",
        startsWith(term, "Tipo.parto") ~ "Tipo.parto",
        startsWith(term, "Ospedale") ~ "Ospedale",
        startsWith(term, "Sesso") ~ "Sesso",
        TRUE ~ term
      )
    ) %>%
    left_join(vif_modello, by = "termine_vif") %>%
    mutate(
      p.value = vapply(p.value, formatta_p, character(1)),
      Termine = dplyr::recode(
        term,
        "Gestazione_c" = "Età gestazionale",
        "Lunghezza_c" = "Lunghezza",
        "Cranio_c" = "Circonferenza cranica",
        "I(Gestazione_c^2)" = "(Età gestazionale)<sup>2</sup>",
        "I(Gestazione_c * Cranio_c)" = "Età gestazionale × Circonferenza cranica",
        "I(Gestazione_c * Lunghezza_c)" = "Età gestazionale × Lunghezza",
        "I(Cranio_c * Lunghezza_c)" = "Circonferenza cranica × Lunghezza",
        "Anni.madre" = "Età materna",
        "N.gravidanze" = "Gravidanze precedenti",
        "FumatriciFumatrice" = "Fumatrice",
        "Tipo.partoCesareo" = "Cesareo",
        "OspedaleOspedale 2" = "Ospedale 2",
        "OspedaleOspedale 3" = "Ospedale 3",
        "SessoMaschio" = "Maschio",
        .default = term
      )
    ) %>%
    transmute(
      Termine,
      Stima = estimate,
      `Errore standard` = std.error,
      `Statistica t` = statistic,
      `p-value` = p.value,
      `VIF (o GVIF corretto)` = VIF
    )
}

knitr::kable(
  tabella_coefficienti(mod1),
  digits = 4,
  caption = "Coefficienti del modello 1",
  escape = FALSE
)

Coefficienti del modello 1

Termine	Stima	Errore standard	Statistica t	p-value	VIF (o GVIF corretto)
Età materna	0.5356	1.1250	0.4761	0,634	1.1901
Gravidanze precedenti	12.4605	4.5801	2.7206	0,007	1.1899
Fumatrice	-27.3556	27.0213	-1.0124	0,311	1.0075
Età gestazionale	29.7564	3.7580	7.9182	< 0,001	1.7052
Lunghezza	10.9308	0.3020	36.1987	< 0,001	2.1444
Circonferenza cranica	9.8623	0.4225	23.3413	< 0,001	1.6631
Cesareo	-29.5171	11.8581	-2.4892	0,013	1.0041
Ospedale 2	-11.7248	13.1888	-0.8890	0,374	1.0010
Ospedale 3	25.1656	13.2495	1.8994	0,058	1.0010
Maschio	77.7578	10.9715	7.0872	< 0,001	1.0403

Ricordiamo che i coefficienti di regressione descrivono come varia il peso atteso al variare di un predittore, mantenendo costanti tutte le altre variabili indipendenti. I p-value elevati associati ai coefficienti dell’età materna e della variabile fumatrice indicano che non disponiamo di evidenze statistiche sufficienti per considerarli significativamente diversi da zero. Questo risultato suggerisce che, nel modello considerato, non emerga un’associazione statisticamente significativa tra queste variabili e il peso alla nascita.

Per individuare un modello più parsimonioso, ci baseremo principalmente sul valore del BIC e, in secondo luogo, sull’\(R^2\). In particolare, procederemo mediante la cosiddetta backward selection, eliminando a ogni passaggio il predittore la cui esclusione determina la maggiore riduzione del BIC. Manterremo tuttavia sempre nel modello il predittore relativo al fumo materno, per ragioni sostantive, dato il potenziale interesse del suo effetto sul peso alla nascita.

knitr::kable(
  tabella_coefficienti(mod2),
  digits = 4,
  caption = "Coefficienti del modello 2",
  escape = FALSE
)

Coefficienti del modello 2

Termine	Stima	Errore standard	Statistica t	p-value	VIF (o GVIF corretto)
Gravidanze precedenti	13.3633	4.2630	3.1347	0,002	1.0263
Fumatrice	-27.4540	27.0703	-1.0142	0,311	1.0067
Età gestazionale	29.8717	3.7442	7.9782	< 0,001	1.6852
Lunghezza	10.8760	0.3022	35.9942	< 0,001	2.1377
Circonferenza cranica	9.9160	0.4227	23.4570	< 0,001	1.6574
Maschio	78.2981	10.9939	7.1219	< 0,001	1.0400

riepilogo_modelli <- function(modello) {
  tibble(
    R_quadro = summary(modello)$r.squared,
    R_quadro_corretto = summary(modello)$adj.r.squared,
    BIC = BIC(modello)
  )
}

knitr::kable(
  bind_rows(
    `Modello 1` = riepilogo_modelli(mod1),
    `Modello 2` = riepilogo_modelli(mod2),
    .id = "Modello"
  ),
  digits = 4,
  caption = "Confronto tra modello 1 e modello 2",
  col.names = c("Modello", "R<sup>2</sup>", "R<sup>2</sup> corretto", "BIC")
)

Confronto tra modello 1 e modello 2

Modello	R2	R2 corretto	BIC
Modello 1	0.7389	0.7378	35104.46
Modello 2	0.7373	0.7367	35088.20

Dalla tabella di confronto, possiamo vedere che sebbene \(R^2\) sia leggermente diminuito, il BIC si è decisamente ridotto, facendoci favorire il secondo modello sul primo.

4.4 Termini quadratici e interazioni

Dai grafici esplorativi non possiamo desumere nessun andamento non lineare particolarmente marcato per le curve LOESS della lunghezza e del cranio nelle zone con più punti (e quindi più affidabili). Invece, nella curva dell’età gestazionale notiamo che la pendenza diminuisce gradualmente, motivo per cui introdurremo una sua potenza quadra.
Inoltre, testeremo l’interazione tra queste tre variabili, molto legate tra loro, aggiungendo anche delle variabili prodotto. I termini quadratici e le interazioni saranno costruiti sulle variabili centrate definite in precedenza. Sarebbe interessante testare anche altre interazioni, ma eviteremo per attenerci al principio di essenzialità imposto per questa analisi.

mod3 <- update(
  mod2,
  . ~ . +
    Gestazione_c + Lunghezza_c + Cranio_c + N.gravidanze + Fumatrici +
    I(Gestazione_c^2) +
    I(Gestazione_c * Cranio_c) +
    I(Gestazione_c * Lunghezza_c) +
    I(Cranio_c * Lunghezza_c)
)

knitr::kable(
  tabella_coefficienti(mod3),
  digits = 4,
  caption = "Coefficienti del modello 3",
  escape = FALSE
)

Coefficienti del modello 3

Termine	Stima	Errore standard	Statistica t	p-value	VIF (o GVIF corretto)
Gravidanze precedenti	14.5844	4.2164	3.4589	< 0,001	1.0280
Fumatrice	-24.6225	26.7592	-0.9201	0,358	1.0072
Età gestazionale	33.5770	4.1880	8.0174	< 0,001	2.1587
Lunghezza	11.1774	0.3053	36.6107	< 0,001	2.2345
Circonferenza cranica	10.3498	0.4240	24.4111	< 0,001	1.7070
Maschio	71.7269	10.9203	6.5682	< 0,001	1.0506
(Età gestazionale)2	-5.7829	1.5117	-3.8255	< 0,001	10.0738
Età gestazionale × Circonferenza cranica	0.5117	0.2084	2.4554	0,014	8.6321
Età gestazionale × Lunghezza	0.3369	0.1560	2.1598	0,031	15.7813
Circonferenza cranica × Lunghezza	0.0149	0.0156	0.9548	0,340	8.4731

I coefficienti positivi dei termini di interazione tra circonferenza cranica e lunghezza indicano che l’effetto marginale di ciascuna di queste variabili sul peso aumenta al crescere dell’altra. Il coefficiente negativo associato al termine quadratico dell’età gestazionale suggerisce invece una relazione concava tra peso ed età gestazionale: il peso atteso tende ad aumentare a un ritmo progressivamente minore al crescere dell’età gestazionale.

Osservando la tabella, si notano gli elevati valori di VIF associati ad alcuni termini, in particolare a quelli quadratici, che segnalano la presenza di multicollinearità. Prima di valutare l’eventuale eliminazione di questi termini, esamineremo come cambiano i VIF a seguito di una nuova procedura di selezione backward.

Si osserva inoltre un p-value relativamente elevato per il termine di interazione tra lunghezza e circonferenza cranica. Questo suggerisce che non vi sia evidenza statistica sufficiente a supportare la presenza di un effetto di interazione tra le due variabili

knitr::kable(
  bind_rows(
    `Modello 2` = riepilogo_modelli(mod2),
    `Modello 3` = riepilogo_modelli(mod3),
    .id = "Modello"
  ),
  digits = 4,
  caption = "Confronto tra modello 2 e modello 3",
  col.names = c("Modello", "R<sup>2</sup>", "R<sup>2</sup> corretto", "BIC")
)

Confronto tra modello 2 e modello 3

Modello	R2	R2 corretto	BIC
Modello 2	0.7373	0.7367	35088.20
Modello 3	0.7438	0.7428	35056.54

Osserviamo che, dopo l’introduzione dei termini quadratici, l’\(R^2\) è aumentato e il BIC è diminuito. Quest’ultimo risultato suggerisce che il miglioramento dell’adattamento del modello compensa ampiamente l’aumento della sua complessità, fornendo un primo riscontro positivo a favore dell’introduzione di tali termini.

4.5 Il modello finale

Come anticipato, applichiamo nuovamente una backward selection basata sul BIC. La formula inferiore mantiene gli effetti principali coinvolti nelle interazioni che abbiamo deciso di valutare: età gestazionale, lunghezza, circonferenza cranica e fumo materno. In questo modo, preserviamo una struttura interpretabile e riduciamo il rischio di imporre al modello andamenti difficili da giustificare. I termini quadratici e le interazioni possono invece essere eliminati se non migliorano sufficientemente il criterio.

mod4 <- step(
  mod3,
  direction = "backward",
  scope = list(
    lower = Peso ~ Gestazione_c + Lunghezza_c + Cranio_c + Fumatrici,
    upper = formula(mod3)
  ),
  k = log(nobs(mod3)),
  trace = 0
)

knitr::kable(
  tabella_coefficienti(mod4),
  digits = 4,
  caption = "Coefficienti del modello 4",
  escape = FALSE
)

Coefficienti del modello 4

Termine	Stima	Errore standard	Statistica t	p-value	VIF (o GVIF corretto)
Gravidanze precedenti	14.6225	4.2162	3.4682	< 0,001	1.0279
Fumatrice	-24.7005	26.7586	-0.9231	0,356	1.0072
Età gestazionale	33.6496	4.1872	8.0362	< 0,001	2.1580
Lunghezza	11.1672	0.3051	36.6002	< 0,001	2.2318
Circonferenza cranica	10.3429	0.4239	24.3989	< 0,001	1.7065
Maschio	71.7527	10.9201	6.5707	< 0,001	1.0506
(Età gestazionale)2	-6.3794	1.3765	-4.6344	< 0,001	8.3534
Età gestazionale × Circonferenza cranica	0.6355	0.1632	3.8946	< 0,001	5.2920
Età gestazionale × Lunghezza	0.4172	0.1314	3.1756	0,002	11.1948

knitr::kable(
  bind_rows(
    `Modello 3` = riepilogo_modelli(mod3),
    `Modello 4` = riepilogo_modelli(mod4),
    .id = "Modello"
  ),
  digits = 4,
  caption = "Confronto tra modello 3 e modello 4",
  col.names = c("Modello", "R<sup>2</sup>", "R<sup>2</sup> corretto", "BIC")
)

Confronto tra modello 3 e modello 4

Modello	R2	R2 corretto	BIC
Modello 3	0.7438	0.7428	35056.54
Modello 4	0.7437	0.7428	35049.64

In questo caso, è stato escluso soltanto il termine di interazione tra circonferenza cranica e lunghezza, il cui p-value (p = 0,340) era elevato.
Tra le specificazioni confrontate, quest’ultimo modello semplificato mediante BIC presenta il valore più basso di questo criterio. Tuttavia, il VIF dell’interazione tra età gestazionale e lunghezza è pari a 11,195, un valore superiore alla soglia orientativa di 10. Per questo motivo, il modello non viene ancora assunto come specificazione di riferimento: di seguito valutiamo una versione priva di tale interazione.

mod5 <- update(
  mod4,
  . ~ . - I(Gestazione_c * Lunghezza_c)
)

# Il modello 5 diventa la specificazione di riferimento per le analisi successive.
modello_finale <- mod5

vif_modello_finale <- vif_valori(modello_finale)

knitr::kable(
  tabella_coefficienti(modello_finale),
  digits = 4,
  caption = "Coefficienti del modello 5",
  escape = FALSE
)

Coefficienti del modello 5

Termine	Stima	Errore standard	Statistica t	p-value	VIF (o GVIF corretto)
Gravidanze precedenti	14.3075	4.2227	3.3883	< 0,001	1.0273
Fumatrice	-25.4167	26.8065	-0.9482	0,343	1.0071
Età gestazionale	33.9147	4.1940	8.0864	< 0,001	2.1572
Lunghezza	11.0366	0.3029	36.4388	< 0,001	2.1913
Circonferenza cranica	10.4058	0.4242	24.5292	< 0,001	1.7028
Maschio	72.9420	10.9336	6.6714	< 0,001	1.0494
(Età gestazionale)2	-3.3646	0.9986	-3.3694	< 0,001	4.3802
Età gestazionale × Circonferenza cranica	0.9181	0.1370	6.7003	< 0,001	3.7180

Dopo l’eliminazione dell’interazione, il VIF massimo scende a 4,380: tutti i valori risultano quindi inferiori a 5 e non abbiamo più problemi di collinearità.

confronto_bic <- tibble(
  Modello = c("Modello 1", "Modello 2", "Modello 3", "Modello 4", "Modello 5"),
  BIC = c(BIC(mod1), BIC(mod2), BIC(mod3), BIC(mod4), BIC(modello_finale)),
  R_quadro_corretto = c(
    summary(mod1)$adj.r.squared,
    summary(mod2)$adj.r.squared,
    summary(mod3)$adj.r.squared,
    summary(mod4)$adj.r.squared,
    summary(modello_finale)$adj.r.squared
  )
)

knitr::kable(
  confronto_bic,
  digits = 4,
  caption = "Confronto tra i modelli",
  col.names = c("Modello", "BIC", "R<sup>2</sup> corretto")
)

Confronto tra i modelli

Modello	BIC	R2 corretto
Modello 1	35104.46	0.7378
Modello 2	35088.20	0.7367
Modello 3	35056.54	0.7428
Modello 4	35049.64	0.7428
Modello 5	35051.92	0.7419

L’eliminazione dell’interazione comporta un lieve aumento del BIC nel modello 5, pari a 2,28, mentre \(R^2\) corretto rimane quasi identico. Il modello 4 conserva quindi un adattamento leggermente migliore secondo il BIC; tuttavia, il peggioramento associato alla sua esclusione è contenuto, mentre la riduzione della multicollinearità è rilevante. Considerando congiuntamente questi elementi, adottiamo il modello 5 come specificazione di riferimento per le analisi successive.

5 Diagnostica del modello

Dopo aver selezionato il modello finale, si procede a verificarne la compatibilità con le principali ipotesi della regressione lineare.

5.1 Analisi del leverage

Il leverage misura quanto una singola osservazione sia distante, nello spazio delle variabili esplicative, dalla configurazione complessiva del campione.

Per individuare valori particolarmente elevati si utilizza come riferimento la soglia convenzionale 2p/n, dove p è il numero di parametri stimati e n la numerosità campionaria.

lev <- hatvalues(modello_finale)
n <- nobs(modello_finale)
p <- sum(lev)

soglia <- 2*p/n
lev_elevati <- lev[lev > soglia]

plot(lev, pch=19, col="#4472C4",
     main="Valori di leverage",
     xlab="ID_record (indice del modello)",
     ylab="Leverage")
abline(h=soglia, col="#D9534F", lty=2, lwd=2)

Nel campione si individuano 204 osservazioni che superano la soglia convenzionale. Un leverage elevato non implica necessariamente che un’osservazione sia influente: segnala soltanto una configurazione insolita delle variabili esplicative.

lev <- hatvalues(modello_finale)
n <- nobs(modello_finale)
p <- sum(lev)

soglia <- 2 * p / n

# Indici delle cinque osservazioni con leverage più elevato
indici_top_lev <- order(lev, decreasing = TRUE)[
  seq_len(min(3, length(lev)))
]

# Tabella delle osservazioni
tabella_leverage <- data.frame(
  ID_record = dati_modello$ID_record[indici_top_lev],
  Leverage = lev[indici_top_lev]
)

knitr::kable(
  tabella_leverage,
  digits = 4,
  caption = "Le tre osservazioni con leverage più elevato"
)

Le tre osservazioni con leverage più elevato

	ID_record	Leverage
310	310	0.3657
1427	1427	0.1646
1777	1777	0.1315

La tabella riporta a titolo di esempio le tre osservazioni caratterizzate dai valori di leverage più elevati. Il loro superamento della soglia convenzionale non implica necessariamente un’influenza significativa sulle stime del modello, ma segnala una posizione relativamente insolita nello spazio dei regressori. Per valutare l’effettiva influenza di queste osservazioni, è necessario considerare anche la distanza di Cook.

5.2 Individuazione degli outlier

Per individuare osservazioni caratterizzate da valori della risposta particolarmente distanti da quelli previsti dal modello, si considerano i residui studentizzati. Questi misurano la distanza tra il peso effettivamente osservato e quello previsto dal modello, espressa in unità di deviazione standard. A differenza dei residui ordinari, tengono conto anche della variabilità della stima e della posizione dell’osservazione rispetto alle altre.

Come riferimento, si considerano i valori −2 e 2: osservazioni con residui studentizzati esterni a questo intervallo presentano uno scostamento relativamente elevato rispetto alle previsioni del modello e meritano pertanto un approfondimento.

res_stud <- rstudent(modello_finale)

plot(res_stud, pch=19, col="#4472C4",
     main="Residui studentizzati",
     xlab="ID_record (indice del modello)",
     ylab="Residuo studentizzato")
abline(h=c(-2,2), col="#D9534F", lty=2, lwd=2)

# Test bilaterale sui residui studentizzati con correzione di Bonferroni
p_non_corretto <- 2 * pt(
  -abs(res_stud),
  df = df.residual(modello_finale) - 1
)
p_bonferroni <- p.adjust(p_non_corretto, method = "bonferroni")

Il grafico permette di individuare le osservazioni caratterizzate da residui studentizzati particolarmente elevati. Per valutare la significatività statistica degli scostamenti, si effettuano test t bilaterali sui residui studentizzati, verificando per ciascuna osservazione l’ipotesi nulla di assenza di uno scostamento anomalo rispetto al modello. Ai p-value così ottenuti si applica la correzione di Bonferroni, che tiene conto dei confronti multipli e riduce il rischio di individuare falsi positivi.

indici_outlier <- which(p_bonferroni < 0.05)

if (length(indici_outlier) == 0) {
  cat("Il test di Bonferroni non individua residui studentizzati significativi al livello 0,05.\n")
} else {
  tabella_outlier_bonferroni <- data.frame(
    ID_originale = dati_modello$ID_originale[indici_outlier],
    Residuo_studentizzato = res_stud[indici_outlier],
    p_Bonferroni = formatta_p(p_bonferroni[indici_outlier])
  )

  knitr::kable(
    tabella_outlier_bonferroni,
    digits = 4,
    caption = "Osservazioni individuate dal test sui residui studentizzati con correzione di Bonferroni",
    col.names = c("ID originale", "Residuo studentizzato", "p Bonferroni")
  )
}

Osservazioni individuate dal test sui residui studentizzati con correzione di Bonferroni

	ID originale	Residuo studentizzato	p Bonferroni
155	155	5.4813	< 0,001
1305	1306	4.8913	0,003
1397	1399	-4.3613	0,034
1691	1694	4.3833	0,030

La procedura individua 4 osservazioni con residui studentizzati significativamente elevati anche dopo la correzione di Bonferroni (p-value aggiustati inferiori a 0,05). I residui individuati non sono tutti positivi e gli scostamenti rispetto alle previsioni possono quindi presentarsi in entrambe le direzioni. Queste osservazioni presentano pertanto uno scostamento particolarmente marcato rispetto alle previsioni, ma questo non significa che siano implausibili.

5.3 Distanza di Cook

La distanza di Cook misura l’influenza esercitata da ciascuna osservazione sulle stime del modello in base alla grandezza dei residui e del leverage. Per valutarne l’entità si utilizzano spesso due soglie indicative: 0,5, oltre la quale un’osservazione può meritare un approfondimento, e 1, che segnala una possibile influenza particolarmente elevata sulle stime del modello. Tali soglie non costituiscono criteri rigidi, ma riferimenti da interpretare nel contesto dell’analisi.

#{r residui-leverage, echo=TRUE, results='hide', message=FALSE, warning=FALSE, fig.width=8, fig.height=4.5, fig.align='center'}

# Residui studentizzati in funzione del leverage e curve della distanza di Cook
plot(modello_finale, which = 5,
     caption = "",
     sub.caption = "",
     main = "Residui vs Leverage")

Il grafico Residui vs Leverage mostra congiuntamente leverage e residui studentizzati. Le curve della distanza di Cook aiutano a riconoscere i casi che possono incidere maggiormente sulle stime.

cook <- cooks.distance(modello_finale)

plot(cook, pch=19, col="#4472C4",
     main="Distanza di Cook",
     xlab="ID record",
     ylab="Distanza di Cook")

# Etichettiamo le cinque osservazioni con la distanza di Cook più elevata.
indici_top_cook <- order(cook, decreasing = TRUE)[seq_len(min(5, length(cook)))]
text(
  x = indici_top_cook,
  y = cook[indici_top_cook],
  labels = dati_modello$ID_record[indici_top_cook],
  pos = 3, cex = 0.75, col = "#17324d"
)

tabella_cook <- tibble(
  ID_record = dati_modello$ID_record,
  ID_originale = dati_modello$ID_originale,
  Distanza_Cook = as.numeric(cook),
  Leverage = as.numeric(hatvalues(modello_finale)),
  Residuo_studentizzato = as.numeric(rstudent(modello_finale))
) %>%
  arrange(desc(Distanza_Cook)) %>%
  slice_head(n = 3) %>%
  mutate(across(c(Distanza_Cook, Leverage, Residuo_studentizzato), ~ round(.x, 4)))

id_cook_max_leverage <- tabella_cook$ID_originale[which.max(tabella_cook$Leverage)]
id_cook_max_residuo <- tabella_cook$ID_originale[which.max(abs(tabella_cook$Residuo_studentizzato))]

knitr::kable(
  tabella_cook,
  caption = paste(nrow(tabella_cook), "osservazioni con la maggiore distanza di Cook"),
  col.names = c("ID record", "ID originale", "Distanza di Cook", "Leverage",
                "Residuo studentizzato")
)

3 osservazioni con la maggiore distanza di Cook

ID record	ID originale	Distanza di Cook	Leverage	Residuo studentizzato
1777	1780	0.0509	0.1315	1.7396
155	155	0.0277	0.0083	5.4813
1917	1920	0.0238	0.0146	3.8177

È interessante osservare che le 3 osservazioni riportate presentano caratteristiche differenti: il record 1780 si distingue soprattutto per l’elevato leverage, mentre il record 155 presenta un residuo studentizzato più elevato in valore assoluto. La distanza di Cook permette di considerare congiuntamente questi aspetti e di valutare la loro influenza complessiva sulle stime del modello.

Nel complesso, i risultati non evidenziano osservazioni con un’influenza particolarmente marcata sulle stime. Non emerge pertanto, sulla base di questo indicatore, una ragione sufficiente per escludere ulteriori record dal campione.

5.4 Test di Breusch-Pagan: omoschedasticità

L’ipotesi di omoschedasticità richiede che la varianza degli errori sia costante al variare delle variabili esplicative. Prima del test osserviamo i grafici che aiutano a valutare la forma funzionale e la dispersione dei residui.

# Residui in funzione dei valori stimati
plot(modello_finale, which=1,
     sub.caption = "")

# Ampiezza dei residui in funzione dei valori stimati
plot(modello_finale, which=3,
     sub.caption = "")

Nel grafico Residuals vs Fitted si cercano eventuali andamenti sistematici, che possono suggerire una forma funzionale non adeguata. Il grafico Scale-Location consente invece di valutare se la dispersione dei residui cambia al variare dei valori stimati: un andamento crescente o decrescente può essere indicativo di eteroschedasticità.

Il controllo grafico viene affiancato dal test di Breusch-Pagan, che consente di verificare formalmente l’ipotesi di omoschedasticità degli errori, ossia la costanza della loro varianza al variare dei valori delle variabili esplicative. L’ipotesi nulla del test è quella di omoschedasticità, mentre l’ipotesi alternativa è la presenza di eteroschedasticità. Sotto l’ipotesi nulla, la statistica test segue approssimativamente una distribuzione chi-quadrato, con gradi di libertà pari al numero di variabili esplicative considerate nella regressione ausiliaria.

library(lmtest)

bp <- bptest(modello_finale)

Il test restituisce una statistica pari a BP = 7.094, con 8 gradi di libertà e un p-value pari a 0.5265. Non si dispone, dunque, di evidenze statistiche sufficienti per rifiutare l’ipotesi nulla di omoschedasticità.

5.5 Test di Durbin-Watson: autocorrelazione

L’indipendenza degli errori viene esaminata attraverso il test di Durbin-Watson, che verifica la presenza di autocorrelazione di primo ordine tra i residui. Il test si basa sulle differenze tra residui consecutivi nell’ordine delle osservazioni: se residui vicini tendono ad assumere valori simili, si manifesta una possibile autocorrelazione positiva, mentre un’alternanza sistematica di valori positivi e negativi può indicare autocorrelazione negativa.

La statistica del test assume valori compresi tra 0 e 4. Un valore prossimo a 2 è generalmente compatibile con l’assenza di autocorrelazione di primo ordine; valori inferiori a 2 suggeriscono autocorrelazione positiva, mentre valori superiori a 2 suggeriscono autocorrelazione negativa.

dw <- dwtest(modello_finale)

La statistica di Durbin-Watson è pari a DW = 1.9634, con un p-value di 0.1804.

Anche in questo caso, non disponiamo di evidenze statistiche sufficienti per rifiutare l’ipotesi nulla di assenza di autocorrelazione positiva di primo ordine. Va tuttavia precisato che il test è particolarmente appropriato quando le osservazioni presentano un ordinamento temporale o comunque significativo. Se l’ordine delle osservazioni nel campione è arbitrario, come nel nostro caso, la sua interpretazione risulta limitata.

5.6 Test di Shapiro-Wilk: normalità

Infine, si verifica l’ipotesi di normalità dei residui mediante il test di Shapiro-Wilk. Il test confronta la distribuzione dei residui osservati con quella attesa nel caso in cui essi provengano da una distribuzione normale. L’ipotesi nulla è pertanto che i residui siano normalmente distribuiti, mentre l’ipotesi alternativa è che si discostino da tale distribuzione. La statistica W, compresa tra 0 e 1, misura il grado di compatibilità dei dati con la normalità: valori prossimi a 1 indicano una maggiore compatibilità, mentre valori più bassi suggeriscono scostamenti più pronunciati.

sh <- shapiro.test(residuals(modello_finale))

Il test restituisce una statistica pari a W = 0.98849 e un p-value pari a 2.611 × 10^-13. Sebbene il valore di W suggerisca una buona compatibilità dei residui con la distribuzione normale, il p-value molto piccolo porta a rifiutare l’ipotesi nulla di normalità. Questa apparente discrepanza è riconducibile anche alla numerosità del campione: quando il numero di osservazioni è elevato, il test può rilevare come statisticamente significativi anche scostamenti relativamente contenuti dalla normalità. Il rifiuto dell’ipotesi nulla non implica quindi necessariamente che la deviazione dalla normalità sia rilevante dal punto di vista pratico.

Il risultato del test deve essere interpretato congiuntamente all’analisi grafica, in quanto la significatività statistica del test non fornisce da sola una misura dell’entità dello scostamento dalla distribuzione normale. Per questo motivo affianchiamo al test il Q-Q plot, che confronta i quantili empirici dei residui con quelli di una distribuzione normale, e il grafico di densità.

# Q-Q plot dei residui
plot(modello_finale, which=2,
     sub.caption = "")

# Distribuzione empirica dei residui
plot(density(residuals(modello_finale)),
     main="Distribuzione dei residui",
     xlab="Residui",
     ylab="Densità",
     col="#4472C4",
     lwd=2)

Nel Q-Q plot, i punti seguono approssimativamente la retta teorica nella parte centrale della distribuzione, mentre si osservano deviazioni più marcate in corrispondenza di entrambe le code, particolarmente evidenti in quella superiore. Questo andamento suggerisce la presenza di code più pesanti rispetto a quelle attese sotto l’ipotesi di normalità.

La distribuzione dei residui presenta invece una forma unimodale e approssimativamente a campana, con una concentrazione delle osservazioni intorno allo zero. Si osservano tuttavia code estese e una moderata asimmetria positiva, coerentemente con le deviazioni evidenziate dal Q-Q plot.

6 Esempio di previsione applicata

Il brief richiede il peso previsto per una neonata la cui madre è alla terza gravidanza e partorisce alla 39ª settimana. Quindi il valore della variabile corrispondente alle gravidanze precedenti dovrà essere fissato a 2.

Non sono specificati età materna, abitudine al fumo materno, ospedale, tipo di parto o misure neonatali. Per illustrare il procedimento fissiamo tali caratteristiche a valori espliciti: mediana per le variabili numeriche e livelli di riferimento per quelle categoriche. Queste scelte sono convenzioni illustrative, non informazioni effettivamente note sulla madre.

profilo <- tibble(
  Anni.madre = median(dati_modello$Anni.madre),
  N.gravidanze = 2,
  Fumatrici = factor(
    "Non fumatrice", levels = levels(dati_modello$Fumatrici)
  ),
  Gestazione = 39,
  Lunghezza = median(dati_modello$Lunghezza),
  Cranio = median(dati_modello$Cranio),
  Tipo.parto = factor(
    "Naturale", levels = levels(dati_modello$Tipo.parto)
  ),
  Ospedale = factor(
    "Ospedale 1", levels = levels(dati_modello$Ospedale)
  ),
  Sesso = factor(
    "Femmina", levels = levels(dati_modello$Sesso)
  )
) %>%
  centra_predittori()

previsione_media <- predict(
  modello_finale, newdata = profilo,
  interval = "confidence", level = 0.95
)

previsione_individuale <- predict(
  modello_finale, newdata = profilo,
  interval = "prediction", level = 0.95
)

tabella_previsione <- tibble(
  Quantità = c(
    "Peso medio atteso per profili simili",
    "Peso previsto per una singola neonata"
  ),
  Stima_g = c(
    previsione_media[1, "fit"],
    previsione_individuale[1, "fit"]
  ),
  Limite_inferiore_g = c(
    previsione_media[1, "lwr"],
    previsione_individuale[1, "lwr"]
  ),
  Limite_superiore_g = c(
    previsione_media[1, "upr"],
    previsione_individuale[1, "upr"]
  )
) %>%
  mutate(across(where(is.numeric), ~ round(.x, 0)))

knitr::kable(
  tabella_previsione,
  caption = "Previsione illustrativa e intervalli al 95% (grammi)"
)

Previsione illustrativa e intervalli al 95% (grammi)

Quantità	Stima_g	Limite_inferiore_g	Limite_superiore_g
Peso medio atteso per profili simili	3320	3302	3338
Peso previsto per una singola neonata	3320	2797	3843

L’intervallo di confidenza riguarda la media attesa per profili simili; l’intervallo di previsione è più ampio perché riguarda una singola neonata e comprende anche la variabilità individuale residua. Per una previsione individuale è quindi più pertinente il secondo.

Limite per l’uso clinico. Il brief indica che lunghezza e misura cranica sono disponibili prima del parto tramite ecografia, ma i valori nel dataset potrebbero essere misure rilevate sul neonato alla nascita. Finché non è confermato che i predittori siano disponibili al momento della previsione e misurati nello stesso modo, il calcolo è un esempio statistico e non una previsione prenatale utilizzabile per decisioni cliniche.

7 Grafici delle relazioni principali

Costruiamo una griglia di valori per ciascuna delle variabili indipendenti di interesse (età gestazionale, lunghezza e circonferenza cranica) e utilizziamo il modello per ottenere il corrispondente peso previsto, distinguendo le previsioni in base allo stato di fumo materno e al sesso del neonato. Le altre variabili considerate nel modello vengono mantenute costanti.

Le curve descrivono quindi il peso medio previsto per specifici profili, e non una media marginale sull’intera popolazione. Questa scelta consente di esaminare l’associazione tra il peso previsto e ciascuna delle variabili indipendenti considerate, evitando che le differenze tra i profili dipendano anche da variazioni degli altri predittori. Le bande rappresentano intervalli di confidenza al 95%.

Nel primo grafico, il numero di gravidanze precedenti, la lunghezza e la circonferenza cranica vengono fissati alle rispettive mediane calcolate sull’intero campione.

# Definiamo l’intervallo di età gestazionale comune ai due sessi
limiti_gestazione <- dati_modello %>%
  group_by(Sesso) %>%
  summarise(
    minimo = min(Gestazione, na.rm = TRUE),
    massimo = max(Gestazione, na.rm = TRUE),
    .groups = "drop"
  )

gest_min <- max(limiti_gestazione$minimo)
gest_max <- min(limiti_gestazione$massimo)


# Costruiamo la griglia delle predizioni
griglia <- tidyr::expand_grid(
  Gestazione = seq(gest_min, gest_max, length.out = 80),
  Fumatrici = factor(
    c("Non fumatrice", "Fumatrice"),
    levels = levels(dati_modello$Fumatrici)
  ),
  Sesso = factor(
    c("Femmina", "Maschio"),
    levels = levels(dati_modello$Sesso)
  )
) %>%
  mutate(
    N.gravidanze = median(dati_modello$N.gravidanze),
    Lunghezza = median(dati_modello$Lunghezza),
    Cranio = median(dati_modello$Cranio)
  ) %>%
  centra_predittori()


predizioni <- predict(
  modello_finale, newdata = griglia,
  interval = "confidence", level = 0.95
) %>%
  as_tibble() %>%
  bind_cols(
    griglia %>% select(Gestazione, Fumatrici, Sesso)
  )


ggplot(
  predizioni,
  aes(
    x = Gestazione,
    y = fit,
    color = Fumatrici,
    fill = Fumatrici
  )
) +
  geom_ribbon(
    aes(ymin = lwr, ymax = upr),
    alpha = 0.16,
    color = NA
  ) +
  geom_line(linewidth = 1) +
  facet_wrap(~ Sesso, nrow = 1) +
  scale_color_manual(
    values = c(
      "Non fumatrice" = "#477c8a",
      "Fumatrice" = "#bd6b42"
    )
  ) +
  scale_fill_manual(
    values = c(
      "Non fumatrice" = "#477c8a",
      "Fumatrice" = "#bd6b42"
    )
  ) +
  scale_x_continuous(
    breaks = scales::breaks_width(2)
  ) +
  scale_y_continuous(
    breaks = scales::breaks_width(200)
  ) +
  labs(
    title = "peso previsto in funzione dell’età gestazionale",
    subtitle = "Confronto per sesso e stato di fumo materno",
    x = "Età gestazionale (settimane)",
    y = "Peso medio previsto (g)",
    color = "Fumo materno",
    fill = "Fumo materno"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 13),
    plot.subtitle = element_text(size = 10),
    strip.text = element_text(face = "bold", size = 11),
    axis.title = element_text(size = 10),
    legend.position = "bottom",
    panel.spacing = unit(1.5, "lines")
  )

In questo primo grafico, il modello restituisce una curva con un andamento crescente, ma con un incremento che tende progressivamente a ridursi nelle settimane più avanzate.
Le bande di confidenza si allargano sensibilmente alle età gestazionali più basse e si restringono nella parte centrale e finale dell’intervallo. Questo indica una maggiore incertezza nella stima del peso medio previsto alle età gestazionali più basse, potenzialmente dovuta anche alla minore numerosità delle osservazioni disponibili in tale regione.

Nel secondo e nel terzo grafico, l’età gestazionale è stata fissata a 40 settimane, un valore comunemente assunto come riferimento per una gravidanza a termine.

limiti_lunghezza <- dati_modello %>%
  group_by(Sesso) %>%
  summarise(
    minimo = min(Lunghezza, na.rm = TRUE),
    massimo = max(Lunghezza, na.rm = TRUE),
    .groups = "drop"
  )

lun_min <- max(limiti_lunghezza$minimo)
lun_max <- min(limiti_lunghezza$massimo)


griglia <- tidyr::expand_grid(
  Lunghezza = seq(lun_min, lun_max, length.out = 80),
  Fumatrici = factor(
    c("Non fumatrice", "Fumatrice"),
    levels = levels(dati_modello$Fumatrici)
  ),
  Sesso = factor(
    c("Femmina", "Maschio"),
    levels = levels(dati_modello$Sesso)
  )
) %>%
  mutate(
    Anni.madre = median(dati_modello$Anni.madre),
    N.gravidanze = median(dati_modello$N.gravidanze),
    Gestazione = 40,
    Cranio = median(dati_modello$Cranio),
    Tipo.parto = factor(
      "Naturale",
      levels = levels(dati_modello$Tipo.parto)
    ),
    Ospedale = factor(
      "Ospedale 1",
      levels = levels(dati_modello$Ospedale)
    )
  ) %>%
  centra_predittori()


predizioni <- predict(
  modello_finale, newdata = griglia,
  interval = "confidence", level = 0.95
) %>%
  as_tibble() %>%
  bind_cols(
    griglia %>% select(Lunghezza, Fumatrici, Sesso)
  )


ggplot(
  predizioni,
  aes(
    x = Lunghezza,
    y = fit,
    color = Fumatrici,
    fill = Fumatrici
  )
) +
  geom_ribbon(
    aes(ymin = lwr, ymax = upr),
    alpha = 0.16,
    color = NA
  ) +
  geom_line(linewidth = 1) +
  facet_wrap(~ Sesso, nrow = 1) +
  scale_color_manual(
    values = c(
      "Non fumatrice" = "#477c8a",
      "Fumatrice" = "#bd6b42"
    )
  ) +
  scale_fill_manual(
    values = c(
      "Non fumatrice" = "#477c8a",
      "Fumatrice" = "#bd6b42"
    )
  ) +
  scale_x_continuous(
    breaks = scales::breaks_width(20)
  ) +
  scale_y_continuous(
    breaks = scales::breaks_width(200)
  ) +
  labs(
    title = "peso previsto in funzione della lunghezza",
    subtitle = "Confronto per sesso e stato di fumo materno",
    x = "Lunghezza (mm)",
    y = "Peso medio previsto (g)",
    color = "Fumo materno",
    fill = "Fumo materno"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 13),
    plot.subtitle = element_text(size = 10),
    strip.text = element_text(face = "bold", size = 11),
    axis.title = element_text(size = 10),
    legend.position = "bottom",
    panel.spacing = unit(1.5, "lines")
  )

limiti_cranio <- dati_modello %>%
  group_by(Sesso) %>%
  summarise(
    minimo = min(Cranio, na.rm = TRUE),
    massimo = max(Cranio, na.rm = TRUE),
    .groups = "drop"
  )

cran_min <- max(limiti_cranio$minimo)
cran_max <- min(limiti_cranio$massimo)


griglia <- tidyr::expand_grid(
  Cranio = seq(cran_min, cran_max, length.out = 80),
  Fumatrici = factor(
    c("Non fumatrice", "Fumatrice"),
    levels = levels(dati_modello$Fumatrici)
  ),
  Sesso = factor(
    c("Femmina", "Maschio"),
    levels = levels(dati_modello$Sesso)
  )
) %>%
  mutate(
    Anni.madre = median(dati_modello$Anni.madre),
    N.gravidanze = median(dati_modello$N.gravidanze),
    Gestazione = 40,
    Lunghezza = median(dati_modello$Lunghezza),
    Tipo.parto = factor(
      "Naturale",
      levels = levels(dati_modello$Tipo.parto)
    ),
    Ospedale = factor(
      "Ospedale 1",
      levels = levels(dati_modello$Ospedale)
    )
  ) %>%
  centra_predittori()


predizioni <- predict(
  modello_finale, newdata = griglia,
  interval = "confidence", level = 0.95
) %>%
  as_tibble() %>%
  bind_cols(
    griglia %>% select(Cranio, Fumatrici, Sesso)
  )


ggplot(
  predizioni,
  aes(
    x = Cranio,
    y = fit,
    color = Fumatrici,
    fill = Fumatrici
  )
) +
  geom_ribbon(
    aes(ymin = lwr, ymax = upr),
    alpha = 0.16,
    color = NA
  ) +
  geom_line(linewidth = 1) +
  facet_wrap(~ Sesso, nrow = 1) +
  scale_color_manual(
    values = c(
      "Non fumatrice" = "#477c8a",
      "Fumatrice" = "#bd6b42"
    )
  ) +
  scale_fill_manual(
    values = c(
      "Non fumatrice" = "#477c8a",
      "Fumatrice" = "#bd6b42"
    )
  ) +
  scale_x_continuous(
    breaks = scales::breaks_width(20)
  ) +
  scale_y_continuous(
    breaks = scales::breaks_width(200)
  ) +
  labs(
    title = "Peso previsto in funzione della circonferenza cranica",
    subtitle = "Confronto per sesso e stato di fumo materno",
    x = "Circonferenza cranica",
    y = "Peso medio previsto (g)",
    color = "Fumo materno",
    fill = "Fumo materno"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", size = 13),
    plot.subtitle = element_text(size = 10),
    strip.text = element_text(face = "bold", size = 11),
    axis.title = element_text(size = 10),
    legend.position = "bottom",
    panel.spacing = unit(1.5, "lines")
  )

In tutti e tre i grafici si osserva un peso previsto leggermente inferiore per i neonati di madri fumatrici, mentre il sesso maschile risulta associato a un peso previsto maggiore rispetto a quello femminile.

Le distanze tra le curve possono essere interpretate quantitativamente attraverso i coefficienti del modello di riferimento. In particolare, a parità degli altri predittori, il fumo materno è associato a una riduzione del peso previsto di circa 25 grammi, ma tale differenza non risulta statisticamente significativa (p = 0,343). Il sesso maschile è invece associato a un peso previsto superiore di circa 73 grammi rispetto a quello femminile, con un’associazione statisticamente significativa (p = < 0,001).

8 Conclusioni

L’analisi ha permesso di esaminare le caratteristiche dei neonati osservati nei tre ospedali, approfondendo alcune relazioni di interesse clinico e costruendo un modello statistico per descrivere il peso alla nascita in funzione delle caratteristiche materne e neonatali.

Per quanto riguarda i confronti descrittivi, non emerge un’associazione statisticamente significativa tra ospedale e tipo di parto: le percentuali di cesareo risultano infatti piuttosto simili nelle tre strutture (p = 0,588). Anche il confronto del peso medio con il valore di riferimento di 3.300 g non evidenzia una differenza statisticamente significativa (p = 0,119). La lunghezza media risulta invece inferiore al riferimento di 50 cm, con uno scostamento di circa 0,52 cm, statisticamente significativo (p < 0,001). Nel confronto tra i sessi, tutte e tre le misure antropometriche considerate mostrano valori medi superiori nei maschi rispetto alle femmine. Si tratta tuttavia di differenze non aggiustate per le altre caratteristiche dei neonati e della gravidanza, che vanno quindi distinte da quelle stimate attraverso il modello di regressione.

La regressione multipla ha permesso di approfondire queste relazioni considerando simultaneamente diversi predittori. Il modello di riferimento include il numero di gravidanze precedenti, il fumo materno, il sesso, l’età gestazionale, la lunghezza e la circonferenza cranica, oltre a termini che consentono di rappresentare un andamento non lineare rispetto all’età gestazionale e un’interazione con la circonferenza cranica. Nel complesso, il modello spiega una quota consistente della variabilità osservata nel peso, con un \(R^2\) corretto pari a 0,742.

La scelta della specificazione di riferimento ha richiesto un compromesso tra capacità esplicativa e stabilità delle stime. In particolare, l’esclusione dell’interazione tra età gestazionale e lunghezza ha permesso di ridurre la multicollinearità, portando tutti i VIF al di sotto della soglia convenzionale di 5, senza modificare sostanzialmente la capacità esplicativa del modello. Questa scelta non elimina ogni possibile incertezza nelle stime, ma consente di adottare una specificazione più equilibrata e più agevole da interpretare.

Le verifiche diagnostiche hanno fornito indicazioni complessivamente articolate. Non emergono evidenze statistiche di eteroschedasticità, mentre il test di Shapiro-Wilk segnala uno scostamento dalla normalità dei residui, da valutare anche alla luce della numerosità campionaria e dell’analisi grafica. Il test di Durbin-Watson, invece, deve essere interpretato con cautela, poiché l’ordine delle osservazioni non ha un significato temporale noto. Infine, le distanze di Cook risultano ampiamente inferiori alle soglie indicative comunemente adottate, suggerendo che nessuna singola osservazione eserciti un’influenza rilevante sul risultato complessivo. Eventuali valori relativamente elevati di leverage o residui meritano comunque di essere considerati nell’interpretazione diagnostica, ma non costituiscono, di per sé, una ragione sufficiente per escludere le corrispondenti osservazioni dall’analisi.

È importante sottolineare che i risultati ottenuti descrivono associazioni osservate nel campione e non consentono di stabilire relazioni causali. Inoltre, le previsioni prodotte dal modello per profili con caratteristiche fissate rappresentano esempi del suo funzionamento, non una dimostrazione della sua efficacia come strumento di previsione clinica. Un’eventuale applicazione in questo ambito richiederebbe una validazione su dati indipendenti, accompagnata da una valutazione specifica dell’accuratezza e dell’affidabilità delle previsioni.

Nel complesso, l’analisi mostra come la considerazione congiunta delle caratteristiche materne e neonatali consenta di descrivere il peso alla nascita in modo più articolato rispetto all’esame delle singole variabili. Il modello costituisce quindi uno strumento utile per sintetizzare le relazioni presenti nel campione e formulare previsioni condizionate a specifici profili, purché i risultati siano interpretati tenendo conto delle assunzioni adottate e dei limiti dei dati disponibili.