neonati
Chiara Tombolini
2025-10-01
knitr::opts_chunk$set(echo = TRUE)
df <- read.csv("/Users/chiaratombolini/Desktop/RSTUDIO/neonati.csv", stringsAsFactors = FALSE)
library(janitor)## Warning: package 'janitor' was built under R version 4.3.3##
## Attaching package: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.testdf <- janitor::clean_names(df)
# Variabili categoriche come fattori
df$tipo_parto <- as.factor(df$tipo_parto)
df$ospedale <- as.factor(df$ospedale)
df$sesso <- as.factor(df$sesso)
df$fumatrici <- as.factor(df$fumatrici)
# 1) Analisi descrittiva ----------------------------------
# Struttura + statistiche sintetiche
str(df)## 'data.frame': 2500 obs. of 10 variables:
## $ anni_madre : int 26 21 34 28 20 32 26 25 22 23 ...
## $ n_gravidanze: int 0 2 3 1 0 0 1 0 1 0 ...
## $ fumatrici : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ gestazione : int 42 39 38 41 38 40 39 40 40 41 ...
## $ peso : int 3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
## $ lunghezza : int 490 490 500 515 480 495 480 510 500 510 ...
## $ cranio : int 325 345 375 365 335 340 345 349 335 362 ...
## $ tipo_parto : Factor w/ 2 levels "Ces","Nat": 2 2 2 2 2 2 2 2 1 1 ...
## $ ospedale : Factor w/ 3 levels "osp1","osp2",..: 3 1 2 2 3 2 3 1 2 2 ...
## $ sesso : Factor w/ 2 levels "F","M": 2 1 2 2 1 1 1 2 1 1 ...summary(df)## anni_madre n_gravidanze fumatrici gestazione peso
## Min. : 0.00 Min. : 0.0000 0:2396 Min. :25.00 Min. : 830
## 1st Qu.:25.00 1st Qu.: 0.0000 1: 104 1st Qu.:38.00 1st Qu.:2990
## Median :28.00 Median : 1.0000 Median :39.00 Median :3300
## Mean :28.16 Mean : 0.9812 Mean :38.98 Mean :3284
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:40.00 3rd Qu.:3620
## Max. :46.00 Max. :12.0000 Max. :43.00 Max. :4930
## lunghezza cranio tipo_parto ospedale sesso
## Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1256
## 1st Qu.:480.0 1st Qu.:330 Nat:1772 osp2:849 M:1244
## Median :500.0 Median :340 osp3:835
## Mean :494.7 Mean :340
## 3rd Qu.:510.0 3rd Qu.:350
## Max. :565.0 Max. :390# Distribuzioni
hist(df$peso, main="Distribuzione Peso", xlab="Peso (g)")
hist(df$lunghezza, main="Distribuzione Lunghezza", xlab="Lunghezza (cm)")
hist(df$anni_madre, main="Distribuzione Età madre", xlab="Età (anni)")
hist(df$gestazione, main="Distribuzione Settimane", xlab="Settimane di gestazione")
# Outlier per peso e lunghezza
outlier_iqr <- function(x) {
x <- x[!is.na(x)]
q1 <- quantile(x, 0.25); q3 <- quantile(x, 0.75)
iqr <- q3 - q1
low <- q1 - 1.5 * iqr; high <- q3 + 1.5 * iqr
list(bounds=c(low=low, high=high))
}
oi_peso <- outlier_iqr(df$peso)
oi_lunghezza <- outlier_iqr(df$lunghezza)
boxplot(df$peso, main="Boxplot Peso", ylab="Peso (g)")
boxplot(df$lunghezza, main="Boxplot Lunghezza", ylab="Lunghezza (cm)")
oi_peso## $bounds
## low.25% high.75%
## 2045 4565oi_lunghezza## $bounds
## low.25% high.75%
## 435 555# 2) Ipotesi ------------------------------------------------
# (a) In alcuni ospedali si fanno più parti cesarei
tab_parti <- table(df$ospedale, df$tipo_parto)
tab_parti##
## Ces Nat
## osp1 242 574
## osp2 254 595
## osp3 232 603chisq.test(tab_parti)##
## Pearson's Chi-squared test
##
## data: tab_parti
## X-squared = 1.0972, df = 2, p-value = 0.5778Conclusione: non si rifiuta H₀ → nessuna evidenza che la frequenza di cesarei differisca tra ospedali nel campione.
# (b) Media del peso e della lunghezza uguali alla popolazione(presi valori medi su internet)
mu_peso <- 3300
mu_lunghezza <- 500
t.test(df$peso, mu = mu_peso)##
## One Sample t-test
##
## data: df$peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081t.test(df$lunghezza, mu = mu_lunghezza)##
## One Sample t-test
##
## data: df$lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 493.6598 495.7242
## sample estimates:
## mean of x
## 494.692Peso medio: coerente con il riferimento (≈ 3300 g), differenza non significativa e di entità pratica minima. Lunghezza media: significativamente più bassa del riferimento (≈ −5 mm)
# (c) Misure antropometriche diverse tra i sessi
t.test(peso ~ sesso, data = df) # peso M vs F##
## Welch Two Sample t-test
##
## data: peso by sesso
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -287.1051 -207.0615
## sample estimates:
## mean in group F mean in group M
## 3161.132 3408.215t.test(lunghezza ~ sesso, data = df) # lunghezza M vs F##
## Welch Two Sample t-test
##
## data: lunghezza by sesso
## t = -9.582, df = 2459.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -11.929470 -7.876273
## sample estimates:
## mean in group F mean in group M
## 489.7643 499.6672t.test(cranio ~ sesso, data = df) # cranio M vs F##
## Welch Two Sample t-test
##
## data: cranio by sesso
## t = -7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -6.089912 -3.541270
## sample estimates:
## mean in group F mean in group M
## 337.6330 342.4486differenze significative tra i sessi per peso, lunghezza e cranio (M > F).
library(moments)
# Indici di forma (asimmetria e curtosi) per il peso
skewness(df$peso, na.rm = TRUE) # asimmetria## [1] -0.6470308kurtosis(df$peso, na.rm = TRUE) - 3 # curtosi centrata## [1] 2.031532L’asimmetria è negativa (skewness ≈ −0,65) e la curtosi in eccesso è ≈ +2,03. Interpretazione: c’è una coda a sinistra (alcuni pesi più bassi) e code più pesanti rispetto a una normale (presenza di casi non comuni).
# Scatterplot peso ~ settimane di gestazione
plot(df$gestazione, df$peso,
xlab = "Settimane di gestazione",
ylab = "Peso (g)",
main = "Relazione Peso ~ Settimane Gestazione",
pch = 19, col = "darkblue")
# Retta di regressione
abline(lm(peso ~ gestazione, data = df), col = "red")
#regressione multipla
mod_full <- lm(
peso ~ anni_madre + n_gravidanze + fumatrici +
gestazione + lunghezza + cranio +
tipo_parto + ospedale + sesso,
data = df
)
summary(mod_full)##
## Call:
## lm(formula = peso ~ anni_madre + n_gravidanze + fumatrici + gestazione +
## lunghezza + cranio + tipo_parto + ospedale + sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1124.40 -181.66 -14.42 160.91 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6738.4762 141.3087 -47.686 < 2e-16 ***
## anni_madre 0.8921 1.1323 0.788 0.4308
## n_gravidanze 11.2665 4.6608 2.417 0.0157 *
## fumatrici1 -30.1631 27.5386 -1.095 0.2735
## gestazione 32.5696 3.8187 8.529 < 2e-16 ***
## lunghezza 10.2945 0.3007 34.236 < 2e-16 ***
## cranio 10.4707 0.4260 24.578 < 2e-16 ***
## tipo_partoNat 29.5254 12.0844 2.443 0.0146 *
## ospedaleosp2 -11.2095 13.4379 -0.834 0.4043
## ospedaleosp3 28.0958 13.4957 2.082 0.0375 *
## sessoM 77.5409 11.1776 6.937 5.08e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 669.2 on 10 and 2489 DF, p-value: < 2.2e-16mod_step <- step(mod_full, direction = "both", trace = FALSE)
summary(mod_step)##
## Call:
## lm(formula = peso ~ n_gravidanze + gestazione + lunghezza + cranio +
## tipo_parto + ospedale + sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1113.18 -181.16 -16.58 161.01 2620.19
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.4293 135.9438 -49.340 < 2e-16 ***
## n_gravidanze 12.3619 4.3325 2.853 0.00436 **
## gestazione 31.9909 3.7896 8.442 < 2e-16 ***
## lunghezza 10.3086 0.3004 34.316 < 2e-16 ***
## cranio 10.4922 0.4254 24.661 < 2e-16 ***
## tipo_partoNat 29.2803 12.0817 2.424 0.01544 *
## ospedaleosp2 -11.0227 13.4363 -0.820 0.41209
## ospedaleosp3 28.6408 13.4886 2.123 0.03382 *
## sessoM 77.4412 11.1756 6.930 5.36e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared: 0.7287, Adjusted R-squared: 0.7278
## F-statistic: 836.3 on 8 and 2491 DF, p-value: < 2.2e-16BIC(mod_full, mod_step)## df BIC
## mod_full 12 35241.84
## mod_step 10 35228.03R2 <- summary(mod_step)$r.squared
res <- residuals(mod_step)
RMSE <- sqrt(mean(res^2, na.rm = TRUE))
R2; RMSE## [1] 0.7286934## [1] 273.4227par(mfrow = c(2, 2))
plot(mod_step)
# Diagnostica test
lmtest::bptest(mod_step)
##
## studentized Breusch-Pagan test
##
## data: mod_step
## BP = 91.768, df = 8, p-value < 2.2e-16shapiro.test(residuals(mod_step))##
## Shapiro-Wilk normality test
##
## data: residuals(mod_step)
## W = 0.97408, p-value < 2.2e-16car::outlierTest(mod_step)## rstudent unadjusted p-value Bonferroni p
## 1551 10.004278 3.9642e-23 9.9104e-20
## 155 5.046640 4.8210e-07 1.2053e-03
## 1306 4.872419 1.1716e-06 2.9291e-03# Punti influenti
par(mfrow = c(1, 1), mar=c(5,4,2,1))
plot(cooks.distance(mod_step), main = "Distanza di Cook")
lev <- hatvalues(mod_step)
par(mar=c(5,4,2,1))
plot(lev, main = "Leverage"); abline(h = 2 * mean(lev), col = "red")
Bontà di adattamento: R² ≈ 0,729, RMSE ≈ 273 g → ~73% della variabilità
del peso spiegata; errore medio delle previsioni ~±273 g. Effetti (a
parità delle altre variabili): Gestazione: circa +32 g per settimana.
Sesso (M vs riferimento F): circa +77 g. Lunghezza: circa +10,3 g per
mm. Cranio: circa +10,5 g per mm. n_gravidanze: effetto positivo, ≈ +12
g per gravidanza.
# PREVISIONI
moda <- function(x) { ux <- na.omit(x); if (!length(ux)) return(NA); names(sort(table(ux), decreasing = TRUE))[1] }
typ_num <- function(x) median(x, na.rm = TRUE)
newx <- data.frame(
anni_madre = if ("anni_madre" %in% names(df)) typ_num(df$anni_madre) else NA_real_,
n_gravidanze = if ("n_gravidanze" %in% names(df)) 3 else NA_real_,
fumatrici = if ("fumatrici" %in% names(df)) factor(moda(df$fumatrici), levels = levels(df$fumatrici)) else NA,
gestazione = if ("gestazione" %in% names(df)) 39 else NA_real_,
lunghezza = if ("lunghezza" %in% names(df)) typ_num(df$lunghezza) else NA_real_,
cranio = if ("cranio" %in% names(df)) typ_num(df$cranio) else NA_real_,
tipo_parto = if ("tipo_parto" %in% names(df)) factor(moda(df$tipo_parto), levels = levels(df$tipo_parto)) else NA,
ospedale = if ("ospedale" %in% names(df)) factor(moda(df$ospedale), levels = levels(df$ospedale)) else NA,
sesso = if ("sesso" %in% names(df)) factor("F", levels = levels(df$sesso)) else NA
)
nm_needed <- all.vars(formula(mod_step))[-1]
missing_cols <- setdiff(nm_needed, colnames(newx))
for (mc in missing_cols) {
if (mc %in% names(df)) {
newx[[mc]] <- if (is.numeric(df[[mc]])) typ_num(df[[mc]]) else factor(moda(df[[mc]]), levels = levels(df[[mc]]))
} else {
newx[[mc]] <- NA
}
}
newx <- newx[, nm_needed, drop = FALSE]
pred_point <- predict(mod_step, newdata = newx)
pred_pred <- predict(mod_step, newdata = newx, interval = "prediction") # singola osservazione
pred_conf <- predict(mod_step, newdata = newx, interval = "confidence") # media attesa
pred_point; pred_pred; pred_conf## 1
## 3317.222## fit lwr upr
## 1 3317.222 2779.3 3855.145## fit lwr upr
## 1 3317.222 3287.977 3346.467# === VISUALIZZAZIONI--------------
# Modello: peso ~ gestazione * fumatrici
m_vis <- lm(peso ~ gestazione * fumatrici, data = df)
summary(m_vis) ##
## Call:
## lm(formula = peso ~ gestazione * fumatrici, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1609.03 -289.03 -11.54 280.97 1898.44
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3240.771 178.861 -18.119 <2e-16 ***
## gestazione 167.495 4.585 36.534 <2e-16 ***
## fumatrici1 1343.426 1164.746 1.153 0.249
## gestazione:fumatrici1 -36.764 29.646 -1.240 0.215
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 422.9 on 2496 degrees of freedom
## Multiple R-squared: 0.352, Adjusted R-squared: 0.3513
## F-statistic: 452 on 3 and 2496 DF, p-value: < 2.2e-16# Predizioni su una griglia di gestazione × fumatrici
gest_vals <- seq(min(df$gestazione, na.rm = TRUE),
max(df$gestazione, na.rm = TRUE), by = 1)
fum_lvls <- levels(df$fumatrici)
grid <- expand.grid(gestazione = gest_vals, fumatrici = fum_lvls)
grid$yhat <- predict(m_vis, newdata = grid)
# Grafico con due linee (una per livello di fumo)
par(mfrow = c(1,1), mar = c(4,4,2,1))
plot(NA,
xlim = range(gest_vals),
ylim = range(grid$yhat, na.rm = TRUE),
xlab = "Settimane di gestazione",
ylab = "Peso previsto (g)",
main = "Peso previsto vs Gestazione per livello di fumo (modello semplice)")
for (fv in fum_lvls) {
sel <- grid$fumatrici == fv
lines(grid$gestazione[sel], grid$yhat[sel], lwd = 2)
}
legend("topleft", legend = fum_lvls, lwd = 2, bty = "n", title = "Fumatrici")
commento al grafico:Le due linee (non fumatrice/fumatrice) crescono con
le settimane (più gestazione ⇒ più peso). Le linee non parallele e
l’incrocio indicano interazione: l’effetto del fumo dipende dalla durata
della gestazione.
CONCLUSIONI: nessuna differenza di cesarei tra ospedali; peso campione allineato a 3300 g; differenze tra sessi nette e significative per tutte le misure. Modello: buona capacità esplicativa/predittiva (R²~0,73; RMSE~273 g). Variabili chiave: gestazione, misure antropometriche (lunghezza, cranio) e sesso. Diagnostica: presenza di eteroschedasticità → preferibili SE robusti; residui non normali (attesi con n grande)..