Modello statistico per la previsione del peso neonatale
Neonatal Health Solutions
Introduzione
Questo report mostra l’analisi esplorativa, i test statistici, la modellizzazione e una previsione di esempio del peso neonatale.
1) Dati e descrizione
Dimensioni del dataset
tibble(
n_osservazioni = nrow(df),
n_variabili = ncol(df)
)## # A tibble: 1 × 2
## n_osservazioni n_variabili
## <int> <int>
## 1 2500 10Statistiche riassuntive (variabili chiave)
summary(df[, c("peso_neonato","settimane_gravidanza","lunghezza","diametro_cranio")])## peso_neonato settimane_gravidanza lunghezza diametro_cranio
## Min. : 830 Min. :25.00 Min. :310.0 Min. :235
## 1st Qu.:2990 1st Qu.:38.00 1st Qu.:480.0 1st Qu.:330
## Median :3300 Median :39.00 Median :500.0 Median :340
## Mean :3284 Mean :38.98 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:40.00 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :43.00 Max. :565.0 Max. :390Cardinalità delle categoriche
sapply(df[, intersect(cats, names(df))], nlevels)## fumo_materno tipo_parto ospedale sesso
## 2 1 3 22) Grafici esplorativi
ggplot(df, aes(peso_neonato)) +
geom_histogram(bins=30) +
labs(x="Peso (g)", y="Frequenza", title="Distribuzione del peso alla nascita")
ggplot(df, aes(settimane_gravidanza, peso_neonato, color=fumo_materno)) +
geom_point(alpha=.5) + geom_smooth(method="lm") +
labs(x="Settimane di gestazione", y="Peso (g)", color="Fumo")
ggplot(df, aes(sesso, peso_neonato, fill=sesso)) +
geom_boxplot() + labs(x="Sesso", y="Peso (g)") + guides(fill="none")
# Composizione percentuale tipo di parto per ospedale
df %>%
mutate(cesareo = if_else(tipo_parto == "cesareo", "Cesareo", "Naturale")) %>%
ggplot(aes(ospedale, fill=cesareo)) +
geom_bar(position="fill") +
scale_y_continuous(labels=scales::percent) +
labs(x="Ospedale", y="Composizione", fill="Tipo parto",
title="Composizione percentuale del tipo di parto per ospedale")
3) Test statistici (ipotesi del progetto)
a) Parti cesarei diversi tra ospedali (Chi-quadrato)
if (all(c("tipo_parto","ospedale") %in% keep_cats)) {
chisq.test(table(df$ospedale, df$tipo_parto))
} else {
"Test non eseguibile: 'tipo_parto' o 'ospedale' hanno un solo livello nel dataset."
}## [1] "Test non eseguibile: 'tipo_parto' o 'ospedale' hanno un solo livello nel dataset."b) Medie campione vs popolazione (aggiorna i riferimenti clinici)
t.test(df$peso_neonato, mu = 3500) # <-- aggiorna 3500 con media di riferimento reale (g)##
## One Sample t-test
##
## data: df$peso_neonato
## t = -20.562, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 3500
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081t.test(df$lunghezza, mu = 50) # <-- aggiorna 50 con media di riferimento reale (cm)##
## One Sample t-test
##
## data: df$lunghezza
## t = 844.82, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 50
## 95 percent confidence interval:
## 493.6598 495.7242
## sample estimates:
## mean of x
## 494.692c) Differenze antropometriche tra i sessi
if ("sesso" %in% keep_cats) {
list(
t_peso_sesso = t.test(peso_neonato ~ sesso, data=df),
t_lung_sesso = t.test(lunghezza ~ sesso, data=df),
t_cranio_sesso = t.test(diametro_cranio ~ sesso, data=df)
)
} else {
"Test non eseguibile: 'sesso' ha un solo livello."
}## [1] "Test non eseguibile: 'sesso' ha un solo livello."4) Modellizzazione (train/test, selezione AIC/BIC, diagnostica)
Split train/test e modello pieno
set.seed(42)
idx <- sample(seq_len(nrow(df)), floor(0.8*nrow(df)))
train <- df[idx, ]; test <- df[-idx, ]
mod_full <- lm(form, data = train)
summary(mod_full)##
## Call:
## lm(formula = form, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1131.02 -179.34 -13.44 172.91 1456.98
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6887.2986 156.1477 -44.108 < 2e-16 ***
## eta_madre 0.2204 1.2647 0.174 0.86164
## numero_gravidanze 14.8134 5.2415 2.826 0.00476 **
## settimane_gravidanza 33.3256 4.3035 7.744 1.52e-14 ***
## lunghezza 10.8107 0.3410 31.706 < 2e-16 ***
## diametro_cranio 10.3079 0.4825 21.363 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 276.1 on 1994 degrees of freedom
## Multiple R-squared: 0.731, Adjusted R-squared: 0.7303
## F-statistic: 1084 on 5 and 1994 DF, p-value: < 2.2e-16Selezione modello (AIC/BIC)
mod_aic <- stepAIC(mod_full, direction="both", trace=FALSE)
mod_bic <- stepAIC(mod_full, direction="both", k=log(nrow(train)), trace=FALSE)
# Scegli il migliore secondo AIC (puoi cambiare in mod_bic)
mod <- mod_aic
summary(mod)##
## Call:
## lm(formula = peso_neonato ~ numero_gravidanze + settimane_gravidanza +
## lunghezza + diametro_cranio, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1128.48 -179.29 -13.54 172.80 1457.63
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6879.9132 150.2531 -45.789 < 2e-16 ***
## numero_gravidanze 15.1342 4.9066 3.084 0.00207 **
## settimane_gravidanza 33.2328 4.2695 7.784 1.12e-14 ***
## lunghezza 10.8127 0.3407 31.739 < 2e-16 ***
## diametro_cranio 10.3112 0.4820 21.391 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 276 on 1995 degrees of freedom
## Multiple R-squared: 0.731, Adjusted R-squared: 0.7304
## F-statistic: 1355 on 4 and 1995 DF, p-value: < 2.2e-16Prestazioni su test
pred <- predict(mod, newdata = test)
RMSE <- sqrt(mean((test$peso_neonato - pred)^2))
R2 <- 1 - sum((test$peso_neonato - pred)^2) / sum((test$peso_neonato - mean(test$peso_neonato))^2)
tibble(RMSE = round(RMSE,1), R2 = round(R2,3))## # A tibble: 1 × 2
## RMSE R2
## <dbl> <dbl>
## 1 282. 0.676Diagnostica residui e controlli
par(mfrow=c(2,2)); plot(mod); par(mfrow=c(1,1))
list(
breusch_pagan = lmtest::bptest(mod),
VIF = tryCatch(car::vif(mod), error = function(e) "VIF non disponibile per questa specifica.")
)## $breusch_pagan
##
## studentized Breusch-Pagan test
##
## data: mod
## BP = 16.732, df = 4, p-value = 0.002179
##
##
## $VIF
## numero_gravidanze settimane_gravidanza lunghezza
## 1.024340 1.689451 2.127292
## diametro_cranio
## 1.6643955) Previsione di esempio
nuovo <- tibble(
eta_madre = 30,
numero_gravidanze = 3,
fumo_materno = if ("fumo_materno" %in% keep_cats) factor("no", levels=c("no","si")) else NULL,
settimane_gravidanza = 39,
lunghezza = 50,
diametro_cranio = 34,
tipo_parto = if ("tipo_parto" %in% keep_cats) factor("naturale", levels=c("naturale","cesareo")) else NULL,
ospedale = if ("ospedale" %in% keep_cats) factor(levels(df$ospedale)[1], levels=levels(df$ospedale)) else NULL,
sesso = if ("sesso" %in% keep_cats) factor("F", levels=c("F","M")) else NULL
) %>% dplyr::select(any_of(all_terms))
predict(mod, nuovo, interval="prediction", level=.95)## fit lwr upr
## 1 -4647.215 -5267.195 -4027.2366) Visualizzazione effetto settimane × fumo
if ("fumo_materno" %in% keep_cats) {
newgrid <- expand.grid(
eta_madre = median(df$eta_madre, na.rm=TRUE),
numero_gravidanze = median(df$numero_gravidanze, na.rm=TRUE),
fumo_materno = factor(c("no","si"), levels=c("no","si")),
settimane_gravidanza = seq(min(df$settimane_gravidanza), max(df$settimane_gravidanza), length.out=50),
lunghezza = median(df$lunghezza, na.rm=TRUE),
diametro_cranio = median(df$diametro_cranio, na.rm=TRUE),
tipo_parto = if ("tipo_parto" %in% keep_cats) factor("naturale", levels=c("naturale","cesareo")) else NULL,
ospedale = if ("ospedale" %in% keep_cats) factor(levels(df$ospedale)[1], levels=levels(df$ospedale)) else NULL,
sesso = if ("sesso" %in% keep_cats) factor("F", levels=c("F","M")) else NULL
) %>% dplyr::select(any_of(all_terms))
newgrid$yhat <- predict(mod, newdata=newgrid)
ggplot() +
geom_point(data=df, aes(settimane_gravidanza, peso_neonato, color=fumo_materno), alpha=.3) +
geom_line(data=newgrid, aes(settimane_gravidanza, yhat, color=fumo_materno), linewidth=1.1) +
labs(x="Settimane di gestazione", y="Peso previsto (g)", color="Fumo")
} else {
plot.new(); title("Grafico non disponibile: 'fumo_materno' ha un solo livello nel dataset.")
}
7) Conclusioni
##- Le settimane di gestazione e (quando presente) lo stato di fumo
materno risultano driver rilevanti del peso previsto.
##- Il modello selezionato (AIC/BIC) raggiunge le prestazioni riportate
(RMSE, R²) sul set di test.
##- Le previsioni vanno interpretate come supporto al
giudizio clinico e monitorate nel tempo per eventuale drift.