Modello statistico per la previsione del peso neonatale
Thomas De Lazzari
2026-10-09
Introduzione
Questo progetto analizza i dati di 2.500 nascite raccolte da tre ospedali, con l’obiettivo di studiare i fattori associati al peso alla nascita e costruire un modello per stimarlo.
L’analisi comprende una prima esplorazione dei dati, alcuni test statistici, un modello di regressione lineare multipla e una previsione applicata a un caso specifico.
I risultati saranno interpretati in modo semplice, tenendo conto anche dei limiti dell’analisi.
1. Pulizia ed esplorazione dei dati
1.1 Importazione e controllo iniziale
Per iniziare, carico il dataset e controllo le dimensioni, le variabili e la presenza di valori mancanti.
# Carico i dati
dati <- read.csv("neonati.csv")
# Controllo le dimensioni
dimensioni <- data.frame(
Osservazioni = nrow(dati),
Variabili = ncol(dati)
)
knitr::kable(dimensioni, row.names = FALSE,
caption = "Dimensioni del dataset")| Osservazioni | Variabili |
|---|---|
| 2500 | 10 |
# Controllo la struttura
struttura <- data.frame(
Variabile = names(dati),
Tipo = sapply(dati, class)
)
knitr::kable(struttura, row.names = FALSE,
caption = "Struttura del dataset")| Variabile | Tipo |
|---|---|
| Anni.madre | integer |
| N.gravidanze | integer |
| Fumatrici | integer |
| Gestazione | integer |
| Peso | integer |
| Lunghezza | integer |
| Cranio | integer |
| Tipo.parto | character |
| Ospedale | character |
| Sesso | character |
# Statistiche descrittive iniziali
numeriche <- dati[sapply(dati, is.numeric)]
statistiche <- data.frame(
Variabile = names(numeriche),
Minimo = sapply(numeriche, min, na.rm = TRUE),
Mediana = sapply(numeriche, median, na.rm = TRUE),
Media = sapply(numeriche, mean, na.rm = TRUE),
Massimo = sapply(numeriche, max, na.rm = TRUE)
)
knitr::kable(statistiche, digits = 2, row.names = FALSE,
caption = "Statistiche descrittive iniziali")| Variabile | Minimo | Mediana | Media | Massimo |
|---|---|---|---|---|
| Anni.madre | 0 | 28 | 28.16 | 46 |
| N.gravidanze | 0 | 1 | 0.98 | 12 |
| Fumatrici | 0 | 0 | 0.04 | 1 |
| Gestazione | 25 | 39 | 38.98 | 43 |
| Peso | 830 | 3300 | 3284.08 | 4930 |
| Lunghezza | 310 | 500 | 494.69 | 565 |
| Cranio | 235 | 340 | 340.03 | 390 |
# Creo una tabella dei valori mancanti
mancanti <- data.frame(
Variabile = names(dati),
Valori_mancanti = colSums(is.na(dati))
)
knitr::kable(mancanti, row.names = FALSE,
caption = "Valori mancanti per variabile")| Variabile | Valori_mancanti |
|---|---|
| Anni.madre | 0 |
| N.gravidanze | 0 |
| Fumatrici | 0 |
| Gestazione | 0 |
| Peso | 0 |
| Lunghezza | 0 |
| Cranio | 0 |
| Tipo.parto | 0 |
| Ospedale | 0 |
| Sesso | 0 |
1.2 Controllo dei valori anomali
Il dataset non presenta valori mancanti, ma alcune variabili hanno valori che potrebbero essere errati. Controllo quindi le osservazioni sospette e le categorie presenti.
# Controllo i valori sospetti
anomali <- dati[
dati$Anni.madre < 12 | dati$Anni.madre > 55 |
dati$Gestazione < 22 | dati$Gestazione > 44 |
dati$Peso < 500 | dati$Peso > 6000 |
dati$Lunghezza < 350 | dati$Lunghezza > 600 |
dati$Cranio < 250 | dati$Cranio > 450,
]
# Mostro le osservazioni da controllare
knitr::kable(anomali, row.names = TRUE,
caption = "Osservazioni con valori sospetti")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| 928 | 25 | 0 | 0 | 28 | 830 | 310 | 254 | Nat | osp1 | F |
| 1152 | 1 | 1 | 0 | 41 | 3250 | 490 | 350 | Nat | osp2 | F |
| 1380 | 0 | 0 | 0 | 39 | 3060 | 490 | 330 | Nat | osp3 | M |
| 1551 | 35 | 1 | 0 | 38 | 4370 | 315 | 374 | Nat | osp3 | F |
| 1619 | 31 | 0 | 0 | 31 | 990 | 340 | 278 | Ces | osp2 | F |
| 1780 | 25 | 2 | 0 | 25 | 900 | 325 | 253 | Nat | osp3 | F |
| 2175 | 37 | 8 | 0 | 28 | 930 | 355 | 235 | Nat | osp1 | F |
| 2437 | 28 | 1 | 0 | 27 | 980 | 320 | 265 | Nat | osp1 | M |
| 2452 | 28 | 0 | 0 | 26 | 930 | 345 | 245 | Ces | osp3 | F |
# Controllo le variabili categoriche
categorie <- c("Fumatrici", "Tipo.parto", "Ospedale", "Sesso")
# Creo una tabella unica con tutte le frequenze
tabella_categorie <- do.call(rbind, lapply(categorie, function(x) {
frequenze <- table(dati[[x]])
data.frame(
Variabile = x,
Categoria = names(frequenze),
Frequenza = as.vector(frequenze),
Percentuale = round(as.vector(prop.table(frequenze)) * 100, 2)
)
}))
# Mostro la tabella
knitr::kable(tabella_categorie, row.names = FALSE,
caption = "Frequenze delle variabili categoriche")| Variabile | Categoria | Frequenza | Percentuale |
|---|---|---|---|
| Fumatrici | 0 | 2396 | 95.84 |
| Fumatrici | 1 | 104 | 4.16 |
| Tipo.parto | Ces | 728 | 29.12 |
| Tipo.parto | Nat | 1772 | 70.88 |
| Ospedale | osp1 | 816 | 32.64 |
| Ospedale | osp2 | 849 | 33.96 |
| Ospedale | osp3 | 835 | 33.40 |
| Sesso | F | 1256 | 50.24 |
| Sesso | M | 1244 | 49.76 |
Il dataset contiene 2.500 osservazioni e non presenta valori mancanti. Sono stati però individuati alcuni valori sospetti, che verranno controllati prima delle analisi.
1.3 Pulizia e preparazione dei dati
Sono state trovate due età materne impossibili e una lunghezza molto sospetta rispetto al peso e alle settimane di gestazione. Questi tre valori vengono considerati mancanti, mentre le misure basse dei neonati prematuri vengono mantenute.
Converto anche lunghezza e cranio in centimetri e preparo le variabili categoriche per le analisi.
# Creo una copia del dataset originale
dati_puliti <- dati
# Correggo le età materne impossibili
dati_puliti$Anni.madre[dati_puliti$Anni.madre < 12] <- NA
# Considero mancante la lunghezza sospetta
dati_puliti$Lunghezza[1551] <- NA
# Converto le misure in centimetri
dati_puliti$Lunghezza <- dati_puliti$Lunghezza / 10
dati_puliti$Cranio <- dati_puliti$Cranio / 10
# Converto le variabili categoriche in fattori
# Converto la variabile in fattore
dati_puliti$Fumatrici <- factor(dati_puliti$Fumatrici,
levels = c(0, 1),
labels = c("No", "Si"))
dati_puliti$Tipo.parto <- factor(dati_puliti$Tipo.parto,
levels = c("Nat", "Ces"),
labels = c("Naturale", "Cesareo"))
dati_puliti$Ospedale <- factor(dati_puliti$Ospedale,
levels = c("osp1", "osp2", "osp3"))
dati_puliti$Sesso <- factor(dati_puliti$Sesso,
levels = c("F", "M"))
# Controllo i valori mancanti dopo la pulizia
controllo <- data.frame(
Variabile = names(dati_puliti),
Valori_mancanti = colSums(is.na(dati_puliti))
)
knitr::kable(controllo, row.names = FALSE,
caption = "Valori mancanti dopo la pulizia")| Variabile | Valori_mancanti |
|---|---|
| Anni.madre | 2 |
| N.gravidanze | 0 |
| Fumatrici | 0 |
| Gestazione | 0 |
| Peso | 0 |
| Lunghezza | 1 |
| Cranio | 0 |
| Tipo.parto | 0 |
| Ospedale | 0 |
| Sesso | 0 |
1.4 Distribuzione delle variabili quantitative
Per capire meglio come sono distribuiti i dati, utilizzo istogrammi e boxplot. Questi grafici permettono di individuare i valori più comuni e quelli più lontani dagli altri.
# Seleziono le variabili quantitative
variabili <- c("Anni.madre", "N.gravidanze", "Gestazione",
"Peso", "Lunghezza", "Cranio")
# Titoli più chiari per i grafici
titoli <- c("Età materna (anni)", "Gravidanze precedenti",
"Gestazione (settimane)", "Peso (g)",
"Lunghezza (cm)", "Cranio (cm)")
# Istogrammi
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
for (i in seq_along(variabili)) {
hist(dati_puliti[[variabili[i]]],
main = titoli[i],
xlab = titoli[i],
col = "lightblue",
border = "white")
}
# Boxplot
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
for (i in seq_along(variabili)) {
boxplot(dati_puliti[[variabili[i]]],
main = titoli[i],
ylab = titoli[i],
col = "lightgreen")
}
# Ripristino la finestra grafica
par(mfrow = c(1, 1))Interpretazione dei risultati
La maggior parte delle madri ha un’età compresa tra 25 e 32 anni e ha avuto poche gravidanze precedenti.Le nascite si concentrano soprattutto tra la 38ª e la 40ª settimana di gestazione.
Il peso dei neonati si concentra intorno ai 3.300 g, la lunghezza intorno ai 50 cm e il diametro craniale ai 34 cm.
Sono presenti alcuni valori molto bassi, che potrebbero essere legati alle nascite premature.
1.5 Distribuzione delle variabili categoriche
Utilizzo i grafici a barre per osservare come sono distribuite le nascite nelle diverse categorie.
# Seleziono le variabili categoriche
categoriche <- c("Fumatrici", "Tipo.parto", "Ospedale", "Sesso")
# Creo i grafici a barre
titoli <- c("Fumo materno", "Tipo di parto",
"Ospedale", "Sesso del neonato")
par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
for (i in seq_along(categoriche)) {
barplot(table(dati_puliti[[categoriche[i]]]),
main = titoli[i],
xlab = titoli[i],
ylab = "Numero di nascite",
col = "lightblue")
}
par(mfrow = c(1, 1))Interpretazione dei risultati
La maggior parte delle madri non fuma (95,84%) e i parti naturali sono più frequenti dei cesarei (70,88%).
Le nascite sono distribuite abbastanza equamente tra i tre ospedali e tra maschi e femmine. Le madri fumatrici sono invece molto meno numerose.
2. Verifica delle ipotesi cliniche
2.1 Differenze nel ricorso al cesareo tra ospedali
Voglio verificare se la percentuale di parti cesarei cambia tra i tre ospedali.
Le ipotesi sono:
- H0: il tipo di parto è indipendente dall’ospedale.
- H1: il tipo di parto è associato all’ospedale.
Utilizzo il test chi-quadrato di Pearson perché entrambe le variabili sono categoriche. Considero un livello di significatività del 5%.
# Creo la tabella dei parti per ospedale
tabella_cesarei <- table(dati_puliti$Ospedale,
dati_puliti$Tipo.parto)
# Calcolo le frequenze e le percentuali
risultati_cesarei <- data.frame(
Ospedale = rownames(tabella_cesarei),
Naturali = tabella_cesarei[, "Naturale"],
Cesarei = tabella_cesarei[, "Cesareo"],
Percentuale_cesarei = round(
prop.table(tabella_cesarei, 1)[, "Cesareo"] * 100, 2)
)
# Mostro la tabella
knitr::kable(risultati_cesarei, row.names = FALSE,
caption = "Parti naturali e cesarei per ospedale")| Ospedale | Naturali | Cesarei | Percentuale_cesarei |
|---|---|---|---|
| osp1 | 574 | 242 | 29.66 |
| osp2 | 595 | 254 | 29.92 |
| osp3 | 603 | 232 | 27.78 |
# Eseguo il test chi-quadrato
test_cesarei <- chisq.test(tabella_cesarei)
# Riassumo i risultati del test
risultati_test <- data.frame(
Chi_quadrato = round(unname(test_cesarei$statistic), 3),
Gradi_liberta = unname(test_cesarei$parameter),
P_value = round(test_cesarei$p.value, 3)
)
knitr::kable(risultati_test, row.names = FALSE,
caption = "Risultati del test chi-quadrato")| Chi_quadrato | Gradi_liberta | P_value |
|---|---|---|
| 1.097 | 2 | 0.578 |
# Controllo la frequenza attesa minima
cat("Frequenza attesa minima:", round(min(test_cesarei$expected), 2))## Frequenza attesa minima: 237.62Interpretazione dei risultati
I cesarei sono il 29,66% nel primo ospedale, il 29,92% nel secondo e il 27,78% nel terzo.
Il p-value è 0,578, quindi superiore a 0,05. Non risultano differenze statisticamente significative tra gli ospedali. Anche le percentuali sono abbastanza simili, con una differenza massima di 2,14 punti percentuali.
2.2 Confronto con i valori di riferimento italiani
Voglio verificare se il peso e la lunghezza medi dei neonati sono diversi dai valori di riferimento di 3.300 g e 50 cm.
Questi valori sono riportati anche dall’Ospedale Pediatrico Bambino Gesù, che indica un peso medio alla nascita di circa 3.300 g e una lunghezza media di 50 cm.
Per il peso: - H0: il peso medio è uguale a 3.300 g. - H1: il peso medio è diverso da 3.300 g.
Per la lunghezza: - H0: la lunghezza media è uguale a 50 cm. - H1: la lunghezza media è diversa da 50 cm.
Utilizzo due test t a un campione perché confronto le medie del campione con valori di riferimento e non conosco la deviazione standard della popolazione. Considero un livello di significatività del 5%.
Fonte: Ospedale Pediatrico Bambino Gesù, Da 0 a 30 giorni: come si presenta e come cresce, https://www.ospedalebambinogesu.it/da-0-a-30-giorni-come-si-presenta-e-come-cresce-80012/
# Test t per il peso e la lunghezza
test_peso <- t.test(dati_puliti$Peso, mu = 3300)
test_lunghezza <- t.test(dati_puliti$Lunghezza, mu = 50)
# Creo una tabella con i risultati
risultati_riferimento <- data.frame(
Variabile = c("Peso (g)", "Lunghezza (cm)"),
Media = c(mean(dati_puliti$Peso, na.rm = TRUE),
mean(dati_puliti$Lunghezza, na.rm = TRUE)),
Riferimento = c(3300, 50),
Differenza = c(mean(dati_puliti$Peso, na.rm = TRUE) - 3300,
mean(dati_puliti$Lunghezza, na.rm = TRUE) - 50),
IC_inferiore = c(test_peso$conf.int[1],
test_lunghezza$conf.int[1]),
IC_superiore = c(test_peso$conf.int[2],
test_lunghezza$conf.int[2]),
P_value = c(test_peso$p.value,
test_lunghezza$p.value)
)
# Arrotondo e mostro la tabella
risultati_riferimento$P_value <- format.pval(
risultati_riferimento$P_value,
digits = 3,
eps = 0.001
)
knitr::kable(risultati_riferimento,
digits = 3,
row.names = FALSE,
caption = "Confronto con i valori di riferimento")| Variabile | Media | Riferimento | Differenza | IC_inferiore | IC_superiore | P_value |
|---|---|---|---|---|---|---|
| Peso (g) | 3284.081 | 3300 | -15.919 | 3263.490 | 3304.672 | 0.13 |
| Lunghezza (cm) | 49.476 | 50 | -0.524 | 49.374 | 49.579 | <0.001 |
Interpretazione dei risultati
Il peso medio dei neonati è di circa 3.284 g, quindi molto vicino al riferimento di 3.300 g. Il p-value è 0,130, quindi la differenza non è significativa.
La lunghezza media è di 49,47 cm, circa mezzo centimetro in meno rispetto ai 50 cm di riferimento. In questo caso la differenza è significativa (p-value < 0,001), ma è piccola e non indica necessariamente un problema clinico.
2.3 Differenze delle misure antropometriche tra maschi e femmine
Voglio verificare se peso, lunghezza e diametro craniale sono diversi tra neonati maschi e femmine.
Per ciascuna misura considero:
- H0: la media è uguale nei due gruppi.
- H1: la media è diversa tra i due gruppi.
Utilizzo il test t di Welch perché confronto due gruppi indipendenti, senza supporre che abbiano la stessa varianza. Considero un livello di significatività del 5%.
# Seleziono le tre misure da confrontare
misure <- c("Peso", "Lunghezza", "Cranio")
nomi <- c("Peso (g)", "Lunghezza (cm)", "Cranio (cm)")
# Eseguo i test t e preparo la tabella
risultati_sesso <- do.call(rbind, lapply(seq_along(misure), function(i) {
variabile <- misure[i]
test <- t.test(dati_puliti[[variabile]] ~ dati_puliti$Sesso)
media_f <- mean(dati_puliti[[variabile]][dati_puliti$Sesso == "F"],
na.rm = TRUE)
media_m <- mean(dati_puliti[[variabile]][dati_puliti$Sesso == "M"],
na.rm = TRUE)
data.frame(
Variabile = nomi[i],
Media_F = media_f,
Media_M = media_m,
Differenza_M_F = media_m - media_f,
IC_inferiore = -test$conf.int[2],
IC_superiore = -test$conf.int[1],
P_value = test$p.value
)
}))
# Formatto i p-value
risultati_sesso$P_value <- format.pval(
risultati_sesso$P_value,
digits = 3,
eps = 0.001
)
# Mostro i risultati
knitr::kable(risultati_sesso, digits = 2,
row.names = FALSE,
caption = "Confronto delle misure tra maschi e femmine")| Variabile | Media_F | Media_M | Differenza_M_F | IC_inferiore | IC_superiore | P_value |
|---|---|---|---|---|---|---|
| Peso (g) | 3161.13 | 3408.22 | 247.08 | 207.06 | 287.11 | <0.001 |
| Lunghezza (cm) | 48.99 | 49.97 | 0.98 | 0.78 | 1.18 | <0.001 |
| Cranio (cm) | 33.76 | 34.24 | 0.48 | 0.35 | 0.61 | <0.001 |
# Creo i boxplot
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
boxplot(Peso ~ Sesso, data = dati_puliti,
main = "Peso", xlab = "Sesso",
ylab = "Peso (g)", col = "lightblue")
boxplot(Lunghezza ~ Sesso, data = dati_puliti,
main = "Lunghezza", xlab = "Sesso",
ylab = "Lunghezza (cm)", col = "lightblue")
boxplot(Cranio ~ Sesso, data = dati_puliti,
main = "Diametro craniale", xlab = "Sesso",
ylab = "Diametro (cm)", col = "lightblue")
par(mfrow = c(1, 1))Interpretazione dei risultati
I maschi hanno in media un peso maggiore di 247 g, una lunghezza superiore di circa 1 cm e una misura craniale maggiore di 0,48 cm rispetto alle femmine.
Per tutte e tre le misure il p-value è inferiore a 0,001, quindi rifiutiamo l’ipotesi nulla. Le differenze sono statisticamente significative, anche se i valori dei due gruppi si sovrappongono in parte.
3. Correlazioni e regressione lineare multipla
3.1 Relazioni tra variabili quantitative
Prima di costruire il modello, controllo le correlazioni tra il peso alla nascita e le altre variabili quantitative. Utilizzo la correlazione di Pearson e gli scatterplot per capire quali variabili sono maggiormente legate al peso e se ci sono relazioni tra i predittori.
# Seleziono le variabili quantitative
variabili_num <- dati_puliti[, c("Anni.madre", "N.gravidanze",
"Gestazione", "Peso",
"Lunghezza", "Cranio")]
# Calcolo la matrice di correlazione
matrice_cor <- round(cor(variabili_num,
use = "complete.obs"), 2)
# Mostro la tabella
knitr::kable(matrice_cor,
caption = "Correlazioni tra le variabili numeriche")| Anni.madre | N.gravidanze | Gestazione | Peso | Lunghezza | Cranio | |
|---|---|---|---|---|---|---|
| Anni.madre | 1.00 | 0.38 | -0.13 | -0.02 | -0.06 | 0.02 |
| N.gravidanze | 0.38 | 1.00 | -0.10 | 0.00 | -0.06 | 0.04 |
| Gestazione | -0.13 | -0.10 | 1.00 | 0.59 | 0.62 | 0.46 |
| Peso | -0.02 | 0.00 | 0.59 | 1.00 | 0.81 | 0.70 |
| Lunghezza | -0.06 | -0.06 | 0.62 | 0.81 | 1.00 | 0.62 |
| Cranio | 0.02 | 0.04 | 0.46 | 0.70 | 0.62 | 1.00 |
# Preparo gli scatterplot
variabili <- c("Anni.madre", "N.gravidanze", "Gestazione",
"Lunghezza", "Cranio")
titoli <- c("Età materna (anni)", "Gravidanze precedenti",
"Gestazione (settimane)", "Lunghezza (cm)",
"Diametro craniale (cm)")
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
# Mostro la relazione tra peso e predittori
for (i in seq_along(variabili)) {
plot(dati_puliti[[variabili[i]]], dati_puliti$Peso,
main = titoli[i],
xlab = titoli[i],
ylab = "Peso (g)",
pch = 20,
col = "steelblue")
}
par(mfrow = c(1, 1))
Interpretazione dei risultati
Il peso è maggiormente legato alla lunghezza (0,81), al cranio (0,70) e alle settimane di gestazione (0,59). In generale, all’aumentare di queste misure aumenta anche il peso.
L’età materna e il numero di gravidanze precedenti mostrano invece correlazioni molto basse.
Lunghezza e cranio sono anche correlati tra loro (0,62), quindi sarà utile controllare questa relazione nella regressione.
3.2 Modello di regressione multipla iniziale
Costruisco un primo modello di regressione lineare multipla per studiare quali variabili sono legate al peso alla nascita. Non inserisco ospedale e tipo di parto, perché non sono necessari all’obiettivo di previsione prima della nascita.
Parto dalle altre variabili e successivamente controllerò quali sono davvero utili per ottenere un modello più semplice.
# Creo il modello iniziale
mod1 <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici +
Gestazione + Lunghezza + Cranio + Sesso,
data = dati_puliti)
# Riassumo i coefficienti del modello
coefficienti <- data.frame(
Variabile = rownames(summary(mod1)$coefficients),
summary(mod1)$coefficients,
row.names = NULL
)
# Rinomino le colonne
names(coefficienti) <- c("Variabile", "Stima", "Errore_standard",
"Valore_t", "P_value")
# Formatto i p-value
coefficienti$P_value <- format.pval(
coefficienti$P_value, digits = 3, eps = 0.001
)
# Mostro la tabella
knitr::kable(coefficienti, digits = 3, row.names = FALSE,
caption = "Coefficienti del modello iniziale")| Variabile | Stima | Errore_standard | Valore_t | P_value |
|---|---|---|---|---|
| (Intercept) | -6704.910 | 138.598 | -48.377 | < 0.001 |
| Anni.madre | 0.605 | 1.127 | 0.537 | 0.59153 |
| N.gravidanze | 12.455 | 4.587 | 2.715 | 0.00667 |
| FumatriciSi | -27.336 | 27.075 | -1.010 | 0.31276 |
| Gestazione | 30.084 | 3.766 | 7.989 | < 0.001 |
| Lunghezza | 108.757 | 3.022 | 35.988 | < 0.001 |
| Cranio | 99.053 | 4.233 | 23.402 | < 0.001 |
| SessoM | 78.242 | 10.996 | 7.116 | < 0.001 |
# Riassumo le caratteristiche del modello
riepilogo_mod1 <- data.frame(
Osservazioni = nobs(mod1),
R2 = summary(mod1)$r.squared,
R2_aggiustato = summary(mod1)$adj.r.squared,
Errore_residuo = summary(mod1)$sigma,
AIC = AIC(mod1),
BIC = BIC(mod1)
)
knitr::kable(riepilogo_mod1, digits = 3, row.names = FALSE,
caption = "Risultati del modello iniziale")| Osservazioni | R2 | R2_aggiustato | Errore_residuo | AIC | BIC |
|---|---|---|---|---|---|
| 2497 | 0.737 | 0.737 | 269.386 | 35043.32 | 35095.73 |
Interpretazione dei risultati
Il modello iniziale spiega circa il 74% delle differenze di peso tra i neonati.
Le variabili più importanti nel modello sono gestazione, lunghezza, cranio e sesso. Anche il numero di gravidanze precedenti risulta significativo, mentre età materna e fumo non mostrano una relazione significativa con il peso.
Nel prossimo passaggio proverò a eliminare le variabili meno utili, controllando che il modello mantenga risultati simili.
3.3 Riduzione del modello
3.3.1 Valutazione dell’età materna
Provo a semplificare il modello eliminando prima l’età materna, che non risulta significativa.
Confronto il nuovo modello con quello iniziale usando R² aggiustato, AIC, BIC e il test ANOVA, per verificare se questa modifica peggiora i risultati.
# Uso le stesse osservazioni del modello iniziale
dati_confronto <- model.frame(mod1)
# Creo il modello senza età materna
mod2 <- update(mod1, ~ . - Anni.madre,
data = dati_confronto)
# Confronto le prestazioni dei due modelli
confronto <- data.frame(
Modello = c("Iniziale", "Senza età materna"),
R2_aggiustato = c(summary(mod1)$adj.r.squared,
summary(mod2)$adj.r.squared),
AIC = c(AIC(mod1), AIC(mod2)),
BIC = c(BIC(mod1), BIC(mod2))
)
knitr::kable(confronto, digits = 3, row.names = FALSE,
caption = "Confronto dei due modelli")| Modello | R2_aggiustato | AIC | BIC |
|---|---|---|---|
| Iniziale | 0.737 | 35043.32 | 35095.73 |
| Senza età materna | 0.737 | 35041.61 | 35088.20 |
# Verifico se l'eliminazione peggiora il modello
test_confronto <- anova(mod2, mod1)
risultato_anova <- data.frame(
Confronto = "Eliminazione età materna",
P_value = format.pval(test_confronto$`Pr(>F)`[2],
digits = 3, eps = 0.001)
)
knitr::kable(risultato_anova, row.names = FALSE,
caption = "Risultato del test ANOVA")| Confronto | P_value |
|---|---|
| Eliminazione età materna | 0.592 |
Interpretazione dei risultati
Eliminando l’età materna, il modello mantiene praticamente gli stessi risultati. L’R² aggiustato resta a 0,737, mentre AIC e BIC diminuiscono leggermente. Anche il test ANOVA (p-value = 0,592) non mostra un peggioramento significativo. Mantengo quindi la versione più semplice.
3.3.2 Valutazione del fumo materno
Dopo aver escluso l’età materna, verifico se posso eliminare anche il fumo materno, che nel modello iniziale non risultava significativo.
Confronto i due modelli per capire se togliere questa variabile peggiora i risultati.
# Creo un modello senza il fumo materno
mod3 <- update(mod2, ~ . - Fumatrici)
# Confronto le prestazioni dei due modelli
confronto_fumo <- data.frame(
Modello = c("Con fumo", "Senza fumo"),
R2_aggiustato = c(summary(mod2)$adj.r.squared,
summary(mod3)$adj.r.squared),
AIC = c(AIC(mod2), AIC(mod3)),
BIC = c(BIC(mod2), BIC(mod3))
)
knitr::kable(confronto_fumo, digits = 3, row.names = FALSE,
caption = "Confronto dei modelli con e senza fumo")| Modello | R2_aggiustato | AIC | BIC |
|---|---|---|---|
| Con fumo | 0.737 | 35041.61 | 35088.2 |
| Senza fumo | 0.737 | 35040.64 | 35081.4 |
# Test ANOVA
test_fumo <- anova(mod3, mod2)
risultato_fumo <- data.frame(
Confronto = "Eliminazione fumo materno",
P_value = format.pval(test_fumo$`Pr(>F)`[2],
digits = 3, eps = 0.001)
)
knitr::kable(risultato_fumo, row.names = FALSE,
caption = "Risultato del test ANOVA")| Confronto | P_value |
|---|---|
| Eliminazione fumo materno | 0.311 |
Interpretazione dei risultati
Eliminando anche il fumo materno, l’R² aggiustato resta a 0,737 e AIC e BIC diminuiscono. Il test ANOVA (p-value = 0,311) non mostra un peggioramento significativo. Mantengo quindi il modello senza fumo, perché è più semplice e offre risultati simili.
Questo non significa che il fumo non influisca sul peso, ma che nel modello non emerge un contributo significativo.
3.3.3 Controllo delle variabili rimaste
Dopo aver eliminato età materna e fumo, controllo i coefficienti delle variabili rimaste per capire se il modello può essere semplificato ancora.
# Controllo i coefficienti del modello ridotto
coefficienti_mod3 <- data.frame(
Variabile = rownames(summary(mod3)$coefficients),
summary(mod3)$coefficients,
row.names = NULL
)
# Rinomino le colonne
names(coefficienti_mod3) <- c(
"Variabile", "Stima", "Errore_standard",
"Valore_t", "P_value"
)
# Formatto i p-value
coefficienti_mod3$P_value <- format.pval(
coefficienti_mod3$P_value,
digits = 3, eps = 0.001
)
# Mostro la tabella
knitr::kable(coefficienti_mod3,
digits = 3,
row.names = FALSE,
caption = "Coefficienti del modello ridotto")| Variabile | Stima | Errore_standard | Valore_t | P_value |
|---|---|---|---|---|
| (Intercept) | -6683.833 | 133.160 | -50.194 | < 0.001 |
| N.gravidanze | 13.145 | 4.258 | 3.087 | 0.00204 |
| Gestazione | 29.634 | 3.737 | 7.930 | < 0.001 |
| Lunghezza | 108.899 | 3.019 | 36.077 | < 0.001 |
| Cranio | 99.192 | 4.227 | 23.465 | < 0.001 |
| SessoM | 78.138 | 10.993 | 7.108 | < 0.001 |
Interpretazione dei risultati
Nel modello ridotto tutte le variabili rimaste sono significative (p-value < 0,05).
Non elimino quindi altre variabili e passo a verificare se il modello può essere migliorato aggiungendo effetti non lineari.
3.4 Valutazione delle relazioni non lineari e delle interazioni
3.4.1 Relazioni non lineari
Dal grafico tra gestazione e peso si nota che la relazione potrebbe non essere perfettamente lineare.
Provo quindi ad aggiungere il quadrato delle settimane di gestazione e confronto il nuovo modello con quello ridotto, per vedere se migliora i risultati.
# Aggiungo il quadrato della gestazione
mod_quadratico <- update(mod3, ~ . + I(Gestazione^2))
# Confronto i due modelli
confronto_quadratico <- data.frame(
Modello = c("Lineare", "Quadratico"),
R2_aggiustato = c(summary(mod3)$adj.r.squared,
summary(mod_quadratico)$adj.r.squared),
AIC = c(AIC(mod3), AIC(mod_quadratico)),
BIC = c(BIC(mod3), BIC(mod_quadratico))
)
knitr::kable(confronto_quadratico, digits = 3,
row.names = FALSE,
caption = "Confronto tra modello lineare e quadratico")| Modello | R2_aggiustato | AIC | BIC |
|---|---|---|---|
| Lineare | 0.737 | 35040.64 | 35081.40 |
| Quadratico | 0.737 | 35035.27 | 35081.86 |
# Test ANOVA
test_quadratico <- anova(mod3, mod_quadratico)
risultato_quadratico <- data.frame(
Confronto = "Aggiunta del quadrato della gestazione",
P_value = format.pval(test_quadratico$`Pr(>F)`[2],
digits = 3, eps = 0.001)
)
knitr::kable(risultato_quadratico, row.names = FALSE,
caption = "Risultato del test ANOVA")| Confronto | P_value |
|---|---|
| Aggiunta del quadrato della gestazione | 0.00671 |
Aggiungendo il quadrato della gestazione, l’R² aggiustato rimane circa 0,737. Il test ANOVA è significativo (p-value = 0,007), ma il miglioramento è piccolo.
L’AIC favorisce il modello quadratico, mentre il BIC preferisce leggermente quello lineare. Per ora mantengo entrambi come possibili modelli.
3.4.2 Interazioni
Voglio verificare se la relazione tra settimane di gestazione e peso cambia tra madri fumatrici e non fumatrici.
Per farlo riparto dal modello che contiene il fumo, aggiungo l’interazione e confronto i risultati.
# Aggiungo l'interazione tra gestazione e fumo
mod_interazione <- update(mod2, ~ . + Gestazione:Fumatrici)
# Confronto i due modelli
confronto_interazione <- data.frame(
Modello = c("Senza interazione", "Con interazione"),
R2_aggiustato = c(summary(mod2)$adj.r.squared,
summary(mod_interazione)$adj.r.squared),
AIC = c(AIC(mod2), AIC(mod_interazione)),
BIC = c(BIC(mod2), BIC(mod_interazione))
)
knitr::kable(confronto_interazione, digits = 3,
row.names = FALSE,
caption = "Confronto dei modelli con e senza interazione")| Modello | R2_aggiustato | AIC | BIC |
|---|---|---|---|
| Senza interazione | 0.737 | 35041.61 | 35088.20 |
| Con interazione | 0.737 | 35042.37 | 35094.78 |
# Test ANOVA
test_interazione <- anova(mod2, mod_interazione)
risultato_interazione <- data.frame(
Confronto = "Interazione tra gestazione e fumo",
P_value = format.pval(test_interazione$`Pr(>F)`[2],
digits = 3, eps = 0.001)
)
knitr::kable(risultato_interazione, row.names = FALSE,
caption = "Risultato del test ANOVA")| Confronto | P_value |
|---|---|
| Interazione tra gestazione e fumo | 0.266 |
Interpretazione dei risultati
Aggiungendo l’interazione tra gestazione e fumo, l’R² aggiustato rimane a 0,737, mentre AIC e BIC aumentano. Il test ANOVA non è significativo (p-value = 0,266). L’interazione quindi non migliora in modo significativo il modello e decido di non mantenerla.
3.5 Scelta del modello finale
Dopo aver confrontato i vari modelli, scelgo quello più semplice che mantiene buoni risultati.
Confronto le prestazioni dei modelli e controllo anche che le variabili non siano troppo legate tra loro.
# Confronto i modelli principali
confronto_finale <- data.frame(
Modello = c("Iniziale", "Senza età", "Ridotto",
"Quadratico", "Con interazione"),
R2_aggiustato = c(summary(mod1)$adj.r.squared,
summary(mod2)$adj.r.squared,
summary(mod3)$adj.r.squared,
summary(mod_quadratico)$adj.r.squared,
summary(mod_interazione)$adj.r.squared),
AIC = c(AIC(mod1), AIC(mod2), AIC(mod3),
AIC(mod_quadratico), AIC(mod_interazione)),
BIC = c(BIC(mod1), BIC(mod2), BIC(mod3),
BIC(mod_quadratico), BIC(mod_interazione))
)
knitr::kable(confronto_finale, digits = 3,
row.names = FALSE,
caption = "Confronto finale dei modelli")| Modello | R2_aggiustato | AIC | BIC |
|---|---|---|---|
| Iniziale | 0.737 | 35043.32 | 35095.73 |
| Senza età | 0.737 | 35041.61 | 35088.20 |
| Ridotto | 0.737 | 35040.64 | 35081.40 |
| Quadratico | 0.737 | 35035.27 | 35081.86 |
| Con interazione | 0.737 | 35042.37 | 35094.78 |
# Scelgo il modello ridotto come candidato finale
mod_finale <- mod3
# Controllo la multicollinearità
valori_vif <- car::vif(mod_finale)
tabella_vif <- data.frame(
Variabile = names(valori_vif),
VIF = as.numeric(valori_vif)
)
knitr::kable(tabella_vif, digits = 3,
row.names = FALSE,
caption = "Controllo della multicollinearità")| Variabile | VIF |
|---|---|
| N.gravidanze | 1.024 |
| Gestazione | 1.679 |
| Lunghezza | 2.133 |
| Cranio | 1.657 |
| Sesso | 1.040 |
Interpretazione dei risultati
Il modello ridotto mantiene un R² aggiustato di 0,737, quindi spiega circa il 74% delle differenze di peso tra i neonati.
Il modello quadratico ottiene un AIC migliore, ma richiede un termine in più. Scelgo il modello ridotto perché è più semplice e ha il BIC leggermente più basso.
I valori VIF sono tutti inferiori a 5, quindi non risultano problemi importanti di multicollinearità.
3.6 Interpretazione dei coefficienti del modello finale
Il modello finale utilizza cinque variabili per stimare il peso alla nascita. Tutte risultano significative (p-value < 0,05).
I coefficienti indicano quanto cambia il peso previsto, mantenendo fisse le altre variabili:
- Gravidanze precedenti: ogni gravidanza in più è associata a circa 13,14 g in più.
- Gestazione: ogni settimana in più è associata a circa 29,63 g in più.
- Lunghezza: ogni centimetro in più è associato a circa 108,90 g in più.
- Cranio: ogni centimetro in più è associato a circa 99,19 g in più.
- Sesso: i maschi hanno un peso previsto maggiore di circa 78,14 g rispetto alle femmine, mantenendo fisse le altre caratteristiche.
L’intercetta è negativa (−6683,83 g), ma non ha un significato clinico utile perché si riferisce a valori pari a zero per le variabili numeriche.
Il modello spiega circa il 74% delle differenze di peso osservate. Le relazioni trovate non dimostrano necessariamente rapporti di causa ed effetto.
4. Diagnostica del modello finale
4.1 Analisi grafica dei residui
Controllo i residui del modello per verificare se presentano andamenti particolari, se la loro variabilità è costante e se alcune osservazioni sono molto diverse dalle altre.
# Visualizzo i quattro grafici diagnostici
par(mfrow = c(2, 2))
plot(mod_finale)
par(mfrow = c(1, 1))I grafici mostrano alcune irregolarità nei residui. Si nota una leggera curvatura e la dispersione non sembra sempre costante.
Anche la distribuzione dei residui non sembra perfettamente normale. Sono presenti alcune osservazioni lontane dalle altre, che verranno controllate nei passaggi successivi.
4.2 Test sui residui
Dai grafici si notano alcune irregolarità nei residui. Utilizzo alcuni test per controllare meglio la loro distribuzione e verificare se la variabilità è costante.
Controllo anche se i residui sono correlati tra loro, anche se questo test è meno utile perché le nascite non hanno un ordine temporale noto.
# Eseguo i test sui residui
test_bp <- lmtest::bptest(mod_finale)
test_dw <- lmtest::dwtest(mod_finale)
test_shapiro <- shapiro.test(residuals(mod_finale))
# Creo una tabella con i risultati
risultati_residui <- data.frame(
Test = c("Breusch-Pagan", "Durbin-Watson", "Shapiro-Wilk"),
Statistica = c(unname(test_bp$statistic),
unname(test_dw$statistic),
unname(test_shapiro$statistic)),
P_value = c(test_bp$p.value,
test_dw$p.value,
test_shapiro$p.value)
)
# Formatto i p-value
risultati_residui$P_value <- format.pval(
risultati_residui$P_value,
digits = 3, eps = 0.001
)
# Mostro la tabella
knitr::kable(risultati_residui,
digits = 3,
row.names = FALSE,
caption = "Risultati dei test sui residui")| Test | Statistica | P_value |
|---|---|---|
| Breusch-Pagan | 11.343 | 0.045 |
| Durbin-Watson | 1.954 | 0.123 |
| Shapiro-Wilk | 0.989 | <0.001 |
# Visualizzo la distribuzione dei residui
plot(density(residuals(mod_finale)),
main = "Distribuzione dei residui",
xlab = "Residui (g)",
ylab = "Densità",
col = "steelblue")
Interpretazione dei risultati
Il test di Breusch-Pagan (p-value = 0,045) indica che la variabilità dei residui non è costante.
Anche il test di Shapiro-Wilk (p-value < 0,001) mostra che i residui non seguono una distribuzione normale.
Il test di Durbin-Watson non è significativo (p-value = 0,123), ma è poco utile perché le nascite non hanno un ordine temporale noto.
Il modello presenta quindi alcuni limiti e le previsioni devono essere interpretate con cautela.
4.3 Individuazione delle osservazioni influenti
Controllo se alcune osservazioni hanno valori particolarmente diversi dalle altre e potrebbero influenzare le stime del modello.
Utilizzo i residui studentizzati, i valori di leverage e la distanza di Cook.
# Calcolo gli indicatori
lev <- hatvalues(mod_finale)
res_student <- rstudent(mod_finale)
cook <- cooks.distance(mod_finale)
# Soglie indicative
n <- nobs(mod_finale)
p <- length(coef(mod_finale))
soglia_lev <- 2 * p / n
soglia_cook <- 4 / n
# Creo i tre grafici
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
plot(lev, main = "Leverage",
xlab = "Osservazione", ylab = "Leverage", pch = 20)
abline(h = soglia_lev, col = "red")
plot(res_student, main = "Residui studentizzati",
xlab = "Osservazione", ylab = "Residuo", pch = 20)
abline(h = c(-3, 3), col = "red")
plot(cook, main = "Distanza di Cook",
xlab = "Osservazione", ylab = "Distanza di Cook", pch = 20)
abline(h = soglia_cook, col = "red")
par(mfrow = c(1, 1))
# Conto le osservazioni oltre le soglie
riepilogo_influenti <- data.frame(
Indicatore = c("Leverage", "Residui studentizzati",
"Distanza di Cook"),
Osservazioni = c(
sum(lev > soglia_lev),
sum(abs(res_student) > 3),
sum(cook > soglia_cook)
)
)
knitr::kable(riepilogo_influenti, row.names = FALSE,
caption = "Osservazioni oltre le soglie")| Indicatore | Osservazioni |
|---|---|
| Leverage | 158 |
| Residui studentizzati | 22 |
| Distanza di Cook | 129 |
# Individuo le 5 osservazioni con Cook maggiore
influenti <- order(cook, decreasing = TRUE)[1:5]
tabella_influenti <- data.frame(
Riga_originale = as.integer(
rownames(model.frame(mod_finale)))[influenti],
Cook = cook[influenti],
Leverage = lev[influenti],
Residuo = res_student[influenti]
)
knitr::kable(tabella_influenti, digits = 3,
row.names = FALSE,
caption = "Cinque osservazioni più influenti")| Riga_originale | Cook | Leverage | Residuo |
|---|---|---|---|
| 310 | 0.065 | 0.029 | -3.620 |
| 1780 | 0.037 | 0.026 | 2.892 |
| 155 | 0.035 | 0.007 | 5.285 |
| 928 | 0.034 | 0.023 | 2.967 |
| 1429 | 0.027 | 0.022 | -2.693 |
Interpretazione dei risultati
Sono state trovate 158 osservazioni con leverage elevato, 22 con residui studentizzati oltre la soglia e 129 con distanza di Cook superiore al limite indicativo.
L’osservazione 310 è quella più influente, con una distanza di Cook di 0,065.
Questi valori non indicano necessariamente errori, ma è importante verificare se le osservazioni più influenti cambiano i risultati del modello.
4.4 Verifica delle osservazioni influenti
L’osservazione 310 è quella che influenza maggiormente il modello.
Provo quindi a ripetere la regressione senza questa osservazione, per controllare se cambiano i coefficienti e i risultati. Non la elimino definitivamente dal dataset.
# Controllo l'osservazione più influente
knitr::kable(dati_puliti[310, ],
row.names = TRUE,
caption = "Osservazione 310")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| 310 | 40 | 3 | No | 28 | 1560 | 42 | 37.9 | Naturale | osp3 | F |
# Creo il modello senza l'osservazione 310
dati_senza_310 <- model.frame(mod_finale)
dati_senza_310 <- dati_senza_310[
rownames(dati_senza_310) != "310", ]
mod_senza_310 <- lm(formula(mod_finale),
data = dati_senza_310)
# Confronto i coefficienti
confronto_coefficienti <- data.frame(
Variabile = names(coef(mod_finale)),
Originale = as.numeric(coef(mod_finale)),
Senza_310 = as.numeric(coef(mod_senza_310))
)
confronto_coefficienti$Differenza <-
confronto_coefficienti$Senza_310 -
confronto_coefficienti$Originale
knitr::kable(confronto_coefficienti,
digits = 2, row.names = FALSE,
caption = "Confronto dei coefficienti")| Variabile | Originale | Senza_310 | Differenza |
|---|---|---|---|
| (Intercept) | -6683.83 | -6662.87 | 20.96 |
| N.gravidanze | 13.14 | 13.26 | 0.11 |
| Gestazione | 29.63 | 28.03 | -1.61 |
| Lunghezza | 108.90 | 108.49 | -0.40 |
| Cranio | 99.19 | 101.02 | 1.83 |
| SessoM | 78.14 | 77.65 | -0.49 |
# Confronto le prestazioni
confronto_prestazioni <- data.frame(
Modello = c("Originale", "Senza osservazione 310"),
R2_aggiustato = c(
summary(mod_finale)$adj.r.squared,
summary(mod_senza_310)$adj.r.squared
),
Errore_residuo = c(
sigma(mod_finale),
sigma(mod_senza_310)
)
)
knitr::kable(confronto_prestazioni,
digits = 3, row.names = FALSE,
caption = "Confronto delle prestazioni")| Modello | R2_aggiustato | Errore_residuo |
|---|---|---|
| Originale | 0.737 | 269.349 |
| Senza osservazione 310 | 0.737 | 268.697 |
Interpretazione dei risultati
L’osservazione 310 è quella che influenza maggiormente il modello aggiornato.
Escludendola temporaneamente, l’R² aggiustato rimane a 0,737 e l’errore residuo diminuisce leggermente, da 269,35 g a 268,70 g. Anche i coefficienti cambiano poco.
Questa osservazione non modifica quindi le conclusioni principali e viene mantenuta nel dataset.
5. Previsione del peso alla nascita per un caso specifico
Utilizzo il modello finale per stimare il peso di una neonata alla 39ª settimana, con madre alla terza gravidanza.
Poiché non conosco la lunghezza e la misura craniale, utilizzo le medie del campione. La previsione è quindi un esempio indicativo e non una stima personalizzata.
# Creo il caso da prevedere
nuovo_caso <- data.frame(
N.gravidanze = 2,
Gestazione = 39,
Lunghezza = mean(dati_puliti$Lunghezza, na.rm = TRUE),
Cranio = mean(dati_puliti$Cranio, na.rm = TRUE),
Sesso = factor("F", levels = levels(dati_puliti$Sesso))
)
# Mostro le caratteristiche utilizzate
knitr::kable(nuovo_caso, digits = 2,
row.names = FALSE,
caption = "Caratteristiche del caso considerato")| N.gravidanze | Gestazione | Lunghezza | Cranio | Sesso |
|---|---|---|---|---|
| 2 | 39 | 49.48 | 34 | F |
# Calcolo la previsione con intervallo al 95%
previsione <- predict(mod_finale,
newdata = nuovo_caso,
interval = "prediction",
level = 0.95)
# Mostro il risultato
risultato_previsione <- data.frame(
Peso_previsto = previsione[1, "fit"],
Limite_inferiore = previsione[1, "lwr"],
Limite_superiore = previsione[1, "upr"]
)
knitr::kable(risultato_previsione,
digits = 1, row.names = FALSE,
caption = "Previsione del peso alla nascita (g)")| Peso_previsto | Limite_inferiore | Limite_superiore |
|---|---|---|
| 3258.9 | 2730.5 | 3787.4 |
Interpretazione dei risultati
Il modello prevede un peso di circa 3.259 g per la neonata considerata. L’intervallo di previsione al 95% va da circa 2.731 g a 3.787 g, quindi il peso effettivo potrebbe essere abbastanza diverso dalla stima.
La previsione è indicativa perché per lunghezza e cranio ho utilizzato i valori medi del campione. Inoltre, i problemi trovati nei residui rendono necessario interpretare l’intervallo con cautela.
6. Visualizzazioni per il personale clinico
6.1 Peso previsto e settimane di gestazione
Mostro come cambia il peso previsto al variare delle settimane di gestazione. Per le altre variabili utilizzo le stesse caratteristiche considerate nella previsione precedente.
# Creo i dati per il grafico
nuovi_dati <- data.frame(
N.gravidanze = 2,
Gestazione = 35:42,
Lunghezza = mean(dati_puliti$Lunghezza, na.rm = TRUE),
Cranio = mean(dati_puliti$Cranio, na.rm = TRUE),
Sesso = factor("F", levels = levels(dati_puliti$Sesso))
)
# Calcolo il peso previsto
nuovi_dati$Peso_previsto <- predict(mod_finale,
newdata = nuovi_dati)
# Creo il grafico
library(ggplot2)
ggplot(nuovi_dati,
aes(x = Gestazione, y = Peso_previsto)) +
geom_line(color = "steelblue", linewidth = 1.2) +
geom_point(color = "steelblue") +
labs(
title = "Peso previsto e settimane di gestazione",
x = "Settimane di gestazione",
y = "Peso previsto (g)"
) +
theme_minimal()
Interpretazione dei risultati
Il grafico mostra che il peso previsto aumenta all’aumentare delle settimane di gestazione. Secondo il modello, ogni settimana in più è associata a circa 30 g in più, mantenendo fisse le altre caratteristiche.
Le previsioni sono indicative perché lunghezza e cranio sono stati mantenuti costanti.
6.2 Peso osservato e fumo materno
Anche se il fumo materno non è stato mantenuto nel modello finale, voglio osservare se il peso alla nascita cambia tra madri fumatrici e non fumatrici.
Utilizzo un boxplot e confronto i pesi medi dei due gruppi.
# Calcolo il peso medio nei due gruppi
media_fumo <- aggregate(Peso ~ Fumatrici,
data = dati_puliti, FUN = mean)
names(media_fumo) <- c("Fumo_materno", "Peso_medio")
# Mostro la tabella
knitr::kable(media_fumo, digits = 2,
row.names = FALSE,
caption = "Peso medio per fumo materno")| Fumo_materno | Peso_medio |
|---|---|
| No | 3286.15 |
| Si | 3236.35 |
# Creo il boxplot
library(ggplot2)
ggplot(dati_puliti, aes(x = Fumatrici, y = Peso)) +
geom_boxplot(fill = "lightblue") +
labs(
title = "Peso alla nascita e fumo materno",
x = "Madre fumatrice",
y = "Peso osservato (g)"
) +
theme_minimal()
Interpretazione dei risultati
I neonati delle madri fumatrici pesano in media circa 50 g in meno rispetto a quelli delle non fumatrici. Il grafico mostra però che i pesi dei due gruppi si sovrappongono molto. Questa differenza è solo descrittiva e non dimostra che sia causata dal fumo.
6.3 Peso previsto e misure antropometriche
Mostro come cambia il peso previsto in base alla lunghezza e al diametro craniale.
Per ciascun grafico modifico una misura alla volta, mantenendo fisse le altre caratteristiche.
# Creo i dati per il grafico della lunghezza
dati_grafici <- data.frame(
N.gravidanze = 2,
Gestazione = 39,
Lunghezza = seq(45, 55, length.out = 100),
Cranio = mean(dati_puliti$Cranio, na.rm = TRUE),
Sesso = factor("F", levels = levels(dati_puliti$Sesso))
)
# Calcolo il peso previsto
dati_grafici$Peso_previsto <- predict(mod_finale,
newdata = dati_grafici)
# Grafico della lunghezza
grafico_lunghezza <- ggplot(dati_grafici,
aes(x = Lunghezza, y = Peso_previsto)) +
geom_line(color = "steelblue", linewidth = 1.2) +
labs(title = "Peso previsto e lunghezza",
x = "Lunghezza (cm)",
y = "Peso previsto (g)") +
theme_minimal()
# Modifico i dati per il grafico del cranio
dati_grafici$Lunghezza <- mean(dati_puliti$Lunghezza,
na.rm = TRUE)
dati_grafici$Cranio <- seq(30, 38, length.out = 100)
# Calcolo le nuove previsioni
dati_grafici$Peso_previsto <- predict(mod_finale,
newdata = dati_grafici)
# Grafico del cranio
grafico_cranio <- ggplot(dati_grafici,
aes(x = Cranio, y = Peso_previsto)) +
geom_line(color = "darkgreen", linewidth = 1.2) +
labs(title = "Peso previsto e diametro craniale",
x = "Diametro craniale (cm)",
y = "Peso previsto (g)") +
theme_minimal()
# Mostro i due grafici
print(grafico_lunghezza)
print(grafico_cranio)
Interpretazione dei risultati
I grafici mostrano che il peso previsto aumenta con la lunghezza e con il diametro craniale.
Secondo il modello, ogni centimetro in più di lunghezza è associato a circa 109 g in più, mentre ogni centimetro in più di diametro craniale è associato a circa 99 g in più.
Le altre caratteristiche rimangono costanti e queste relazioni non indicano necessariamente un rapporto causa-effetto.
7. Conclusioni
L’analisi delle 2.500 nascite ha permesso di studiare le caratteristiche dei neonati e di costruire un modello per stimare il peso alla nascita.
I test non hanno mostrato differenze significative nella frequenza dei cesarei tra i tre ospedali. Il peso medio è risultato compatibile con il riferimento di 3.300 g, mentre la lunghezza media è leggermente inferiore ai 50 cm. I maschi presentano in media peso, lunghezza e misura craniale maggiori rispetto alle femmine.
Il modello finale utilizza cinque variabili: gravidanze precedenti, gestazione, lunghezza, cranio e sesso. Spiega circa il 74% delle differenze di peso osservate e mantiene risultati simili ai modelli più complessi.
Per la neonata considerata, il peso previsto è di circa 3.259 g, con un intervallo di previsione al 95% compreso tra 2.731 e 3.787 g.
Il modello presenta però alcuni limiti perché i residui non rispettano completamente le ipotesi di normalità e varianza costante. Il controllo dell’osservazione 310 non ha modificato sostanzialmente i risultati.
Per una previsione durante la gravidanza servirebbero anche misure ecografiche affidabili di lunghezza e cranio. Nel nostro esempio sono stati utilizzati i valori medi del campione.
Nel complesso, il modello può aiutare a comprendere le relazioni tra le variabili e fornire stime indicative, ma non è sufficientemente verificato per essere utilizzato da solo nelle decisioni cliniche.