knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
set.seed(1234)

# core
library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(janitor)
library(forcats)

# descrittive / test / tabeles
library(moments)
library(knitr)

# model
library(MASS) 
library(car)    # VIF
dati_raw <- read_excel("C:/Users/Trasciati/Desktop/professionAI/250728_neonati.xlsx")
  1. struttura dati
# pulizia nomi + etichette
dati <- dati_raw %>%
janitor::clean_names() %>%
dplyr::rename(
eta_madre = anni_madre,
num_gravidanze = n_gravidanze,
fumo_raw = fumatrici,
settimane_gestazione = gestazione,
peso_neonato = peso,
lunghezza_neonato = lunghezza,
diametro_craniale = cranio,
tipo_parto = tipo_parto,
ospedale = ospedale,
sesso = sesso
) %>%
dplyr::mutate(
# numerici
dplyr::across(c(eta_madre, num_gravidanze, settimane_gestazione,
peso_neonato, lunghezza_neonato, diametro_craniale),
~ suppressWarnings(as.numeric(.))),
# fattoriali
fumo = factor(ifelse(fumo_raw %in% c(1, "1"), "Si", "No"), levels = c("No","Si")),
tipo_parto = forcats::fct_recode(as.factor(tipo_parto), Naturale = "Nat", Cesareo = "Ces"),
ospedale = as.factor(ospedale),
sesso = factor(sesso, levels = c("F","M")),
# specifiche
primigravida = ifelse(num_gravidanze == 0, 1L, 0L), 
grav_tot = num_gravidanze + 1L,
eta_madre_flag_na = ifelse(!is.na(eta_madre) & eta_madre < 13, TRUE, FALSE),
eta_madre = ifelse(eta_madre < 13, NA_real_, eta_madre)
)


list(
n_righe = nrow(dati),
n_variabili = ncol(dati),
livelli_ospedale = levels(dati$ospedale),
livelli_parto = levels(dati$tipo_parto),
livelli_fumo = levels(dati$fumo),
livelli_sesso = levels(dati$sesso)
)
## $n_righe
## [1] 2500
## 
## $n_variabili
## [1] 14
## 
## $livelli_ospedale
## [1] "osp1" "osp2" "osp3"
## 
## $livelli_parto
## [1] "Cesareo"  "Naturale"
## 
## $livelli_fumo
## [1] "No" "Si"
## 
## $livelli_sesso
## [1] "F" "M"
  1. analisi e modellizzazione (preliminare)

1 Statistiche Descrittive

  • per tutte le variabili si calcolano i missing
  • per le variabili numeriche si riportano: min, max, media/mediana, deviazione standard, intervallo interquartile (iqr), asimmetria (skew), curtosi (kurt), valore e p-value dello shapiro–wilk (stimato su campione ≤500 per evitare rifiuti automatici con n elevato) e conteggio degli outlier (regola 1.5×iqr). per le stesse variabili si mostrano istogrammi e box plot.
  • per le variabili categoriali si riportano tabelle di frequenza e grafici a barre in percentuale. I test sono lo χ² su cesarei vs ospedale; t a 1 campione per peso e lunghezza; t di Welch + Wilcoxon per differenze F vs M.

Cosa si nota: - missing trascurabili (non c’è bisogno di fare imputazione) - x numeriche: distribuzioni asimmetriche (Shapiro spesso < 0.05. si useranno per questo test come Welch/Wilcoxon - x categoriali: sesso è abbastanza bilanciato; ospedali con numerosità comparabili; fumo circa il 4% - misure antropometriche in mm - outlier (1.5×IQR, solo continue): presenti su settimane, peso, lunghezza, cranio e gravidanze; sono clinicamente plausibili ((prematurità, SGA/LGA, multiparità - dalle mie nozioni, lavorando in ospedale), per cui non h effettuato rimozioni. solo valori implausibili sono stati esclusi a priori (es. età materna <13).

# missing, descrittive, normalità, outlier, grafici

# 1) percentuali di missing
miss_tbl <- sapply(dati, function(x) mean(is.na(x))) |>
  tibble::enframe(name = "variabile", value = "perc_missing") |>
  dplyr::arrange(dplyr::desc(perc_missing))

knitr::kable(
  miss_tbl,
  digits = 3,
  caption = "percentuali di missing per variabile"
)
percentuali di missing per variabile
variabile perc_missing
eta_madre 0.001
num_gravidanze 0.000
fumo_raw 0.000
settimane_gestazione 0.000
peso_neonato 0.000
lunghezza_neonato 0.000
diametro_craniale 0.000
tipo_parto 0.000
ospedale 0.000
sesso 0.000
fumo 0.000
primigravida 0.000
grav_tot 0.000
eta_madre_flag_na 0.000 
#2) variabili numeriche e continue (magg uguale 10 valori distinti)
num_vars <- dati |> dplyr::select(where(is.numeric)) |> names()

num_vars_cont <- dati |>
  dplyr::select(dplyr::all_of(num_vars)) |>
  dplyr::select(where(~ dplyr::n_distinct(., na.rm = TRUE) >= 10)) |>
  names()

# 3) statistiche descrittive per numeriche (tabella)
desc_num <- dati |>
  dplyr::select(all_of(num_vars)) |>
  dplyr::summarise(dplyr::across(
    dplyr::everything(),
    list(
      n      = ~sum(!is.na(.)),
      mean   = ~mean(., na.rm = TRUE),
      sd     = ~sd(., na.rm = TRUE),
      q1     = ~quantile(., .25, na.rm = TRUE),
      median = ~median(., na.rm = TRUE),
      q3     = ~quantile(., .75, na.rm = TRUE),
      min    = ~min(., na.rm = TRUE),
      max    = ~max(., na.rm = TRUE),
      skew   = ~moments::skewness(., na.rm = TRUE),
      kurt   = ~moments::kurtosis(., na.rm = TRUE)
    ),
    .names = "{.col}_{.fn}"
  ))

desc_num_tidy <- desc_num |>
  tidyr::pivot_longer(
    cols = dplyr::everything(),
    names_to = c("variabile", "stat"),
    names_sep = "_(?=[^_]+$)",
    values_to = "val"
  ) |>
  tidyr::pivot_wider(
    names_from = stat,
    values_from = val
  ) |>
  dplyr::select(
    variabile, n, mean, sd, q1, median, q3, min, max, skew, kurt
  ) |>
  dplyr::arrange(variabile)

knitr::kable(
  desc_num_tidy,
  digits = 2,
  caption = "statistiche descrittive variabili numeriche (riassunto)"
)
statistiche descrittive variabili numeriche (riassunto)
variabile n mean sd q1 median q3 min max skew kurt
diametro_craniale 2500 340.03 16.43 330 340 350 235 390 -0.79 5.95
eta_madre 2498 28.19 5.22 25 28 32 13 46 0.15 2.89
fumo_raw 2500 0.04 0.20 0 0 0 0 1 4.59 22.08
grav_tot 2500 1.98 1.28 1 2 2 1 13 2.51 13.99
lunghezza_neonato 2500 494.69 26.32 480 500 510 310 565 -1.51 9.49
num_gravidanze 2500 0.98 1.28 0 1 1 0 12 2.51 13.99
peso_neonato 2500 3284.08 525.04 2990 3300 3620 830 4930 -0.65 5.03
primigravida 2500 0.44 0.50 0 0 1 0 1 0.25 1.06
settimane_gestazione 2500 38.98 1.87 38 39 40 25 43 -2.07 11.26 
# 4) shapiro–wilk (campione minore uguale 500 per variabile)
set.seed(1)
shapiro_tbl <- lapply(num_vars_cont, function(v){
  x <- stats::na.omit(dati[[v]])
  if (length(x) < 3) {
    return(data.frame(variabile = v, shapiro_W = NA, shapiro_p = NA))
  }
  xs <- if (length(x) > 500) sample(x, 500) else x
  sw <- shapiro.test(xs)
  data.frame(
    variabile  = v,
    shapiro_W  = unname(sw$statistic),
    shapiro_p  = sw$p.value
  )
}) |>
  dplyr::bind_rows()

knitr::kable(
  shapiro_tbl,
  digits = 4,
  caption = "shapiro–wilk (campione ≤500) per variabili continue"
)
shapiro–wilk (campione ≤500) per variabili continue
variabile shapiro_W shapiro_p
eta_madre 0.9936 0.0329
num_gravidanze 0.7431 0.0000
settimane_gestazione 0.7987 0.0000
peso_neonato 0.9900 0.0017
lunghezza_neonato 0.9118 0.0000
diametro_craniale 0.9557 0.0000
grav_tot 0.7125 0.0000 
# 5) istogrammi (solo continue)
col_fill <- "#6EC5FF"
col_line <- "#1F6FEB"

dati |>
  tidyr::pivot_longer(all_of(num_vars_cont),
                      names_to = "variabile",
                      values_to = "val") |>
  ggplot2::ggplot(ggplot2::aes(x = val)) +
  ggplot2::geom_histogram(bins = 30,
                          fill = col_fill,
                          color = col_line) +
  ggplot2::facet_wrap(~ variabile, scales = "free") +
  ggplot2::labs(title = "distribuzioni variabili continue",
                x = "",
                y = "conteggio") +
  ggplot2::theme_minimal() +
  ggplot2::theme(
    panel.grid.minor = ggplot2::element_blank(),
    strip.background = ggplot2::element_rect(fill = "white", color = NA),
    panel.background = ggplot2::element_rect(fill = "white", color = NA)
  )

# 6) outlier (regola 1.5 × IQR) — solo continue
iqr_flag <- function(x){
  x <- x[!is.na(x)]
  if (length(x) < 4 || diff(range(x))==0) return(0L)
  q1  <- stats::quantile(x, 0.25)
  q3  <- stats::quantile(x, 0.75)
  iqr <- q3 - q1
  low <- q1 - 1.5 * iqr
  up  <- q3 + 1.5 * iqr
  sum(x < low | x >up)
}

out_tbl <- tibble::tibble(
  variabile = num_vars_cont,
  n_outlier_IQR = sapply(dati[num_vars_cont], iqr_flag)
)

knitr::kable(
  out_tbl,
  caption = "outlier conteggiati (1.5×IQR) — variabili continue"
)
outlier conteggiati (1.5×IQR) — variabili continue
variabile n_outlier_IQR
eta_madre 11
num_gravidanze 246
settimane_gestazione 67
peso_neonato 69
lunghezza_neonato 59
diametro_craniale 48
grav_tot 246 
# 7) boxplot (solo continue)
dati |>
  dplyr::select(all_of(num_vars_cont)) |>
  tidyr::pivot_longer(dplyr::everything(),
                      names_to = "variabile",
                      values_to = "val") |>
  ggplot2::ggplot(ggplot2::aes(x = variabile, y = val)) +
  ggplot2::geom_boxplot(
    fill = col_fill,
    color = col_line,
    outlier.colour = col_line,
    outlier.alpha = 0.5
  ) +
  ggplot2::coord_flip() +
  ggplot2::labs(
    title = "boxplot variabili continue",
    x = "",
    y = ""
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(
    panel.grid.minor = ggplot2::element_blank(),
    panel.background = ggplot2::element_rect(fill = "white", color = NA)
  )

# 8) variabili categoriali: tabelle di frequenza + barre
cat_vars <- dati |>
  dplyr::select(where(is.factor)) |>
  names()

for (v in cat_vars) {
  tab <- dati |>
    dplyr::count(.data[[v]], name = "n") |>
    dplyr::mutate(pct = n / sum(n))

  cat("\n\n###", v, "\n")
  print(
    knitr::kable(
      tab,
      digits = 3,
      caption = paste("distribuzione di", v)
    )
  )

  p <- ggplot2::ggplot(tab, ggplot2::aes(x = .data[[v]], y = pct)) +
    ggplot2::geom_col(fill = col_fill, color = col_line) +
    ggplot2::scale_y_continuous(labels = scales::percent_format()) +
    ggplot2::labs(
      title = paste("distribuzione di", v),
      x = v,
      y = "percentuale"
    ) +
    ggplot2::coord_flip() +
    ggplot2::theme_minimal() +
    ggplot2::theme(
      panel.grid.minor = ggplot2::element_blank(),
      panel.background = ggplot2::element_rect(fill = "white", color = NA)
    )

  print(p)
}
## 
## 
## ### tipo_parto 
## 
## 
## Table: distribuzione di tipo_parto
## 
## |tipo_parto |    n|   pct|
## |:----------|----:|-----:|
## |Cesareo    |  728| 0.291|
## |Naturale   | 1772| 0.709|

## 
## 
## ### ospedale 
## 
## 
## Table: distribuzione di ospedale
## 
## |ospedale |   n|   pct|
## |:--------|---:|-----:|
## |osp1     | 816| 0.326|
## |osp2     | 849| 0.340|
## |osp3     | 835| 0.334|

## 
## 
## ### sesso 
## 
## 
## Table: distribuzione di sesso
## 
## |sesso |    n|   pct|
## |:-----|----:|-----:|
## |F     | 1256| 0.502|
## |M     | 1244| 0.498|

## 
## 
## ### fumo 
## 
## 
## Table: distribuzione di fumo
## 
## |fumo |    n|   pct|
## |:----|----:|-----:|
## |No   | 2396| 0.958|
## |Si   |  104| 0.042|

TEST IPOTESI

# 9.1) in alcuni ospedali si fanno più parti cesarei (chi yest + cramér's v + post hoc)
tab_parti <- table(dati$ospedale, dati$tipo_parto)
chi_res   <- chisq.test(tab_parti, correct = FALSE)
N <- sum(tab_parti); r <- nrow(tab_parti); c <- ncol(tab_parti)
cramers_v <- as.numeric(sqrt((chi_res$statistic / N) / (min(r - 1, c - 1))))

prop_ces <- prop.table(tab_parti, 1)[, "Cesareo"]
prop_tbl <- data.frame(
  ospedale = names(prop_ces),
  proporzione_cesarei = as.numeric(prop_ces)
)
knitr::kable(prop_tbl, digits = 3, caption = "proporzione di cesarei per ospedale")
proporzione di cesarei per ospedale
ospedale proporzione_cesarei
osp1 0.297
osp2 0.299
osp3 0.278 
hosp <- rownames(tab_parti); pairs <- t(combn(hosp, 2))
p_raw <- apply(pairs, 1, function(idx){
  a <- idx[1]; b <- idx[2]
  prop.test(x = c(tab_parti[a, "Cesareo"], tab_parti[b, "Cesareo"]),
            n = c(sum(tab_parti[a, ]),       sum(tab_parti[b, ])))$p.value
})
p_holm <- p.adjust(p_raw, method = "holm")
posthoc_tbl <- data.frame(
  ospedale_1 = pairs[,1],
  ospedale_2 = pairs[,2],
  p_value = round(p_raw, 5),
  p_adj_holm = round(p_holm, 5)
)
knitr::kable(posthoc_tbl, caption = "confronti a coppie sulle proporzioni di cesarei (holm)")
confronti a coppie sulle proporzioni di cesarei (holm)
ospedale_1 ospedale_2 p_value p_adj_holm
osp1 osp2 0.95003 1
osp1 osp3 0.43164 1
osp2 osp3 0.36170 1 
data.frame(
  X2 = as.numeric(chi_res$statistic),
  df = unname(chi_res$parameter),
  p_value = as.numeric(chi_res$p.value),
  cramers_v = cramers_v
) |>
  knitr::kable(digits = 4, caption = "chi-quadro: cesarei ~ ospedale (con cramér's v)")
chi-quadro: cesarei ~ ospedale (con cramér’s v)
X2 df p_value cramers_v
1.0972 2 0.5778 0.0209

2.1) TEST 9.1) hp: in alcuni ospedali si fanno più parti cesarei nel campione analizzato la proporzione di parti cesarei è simile nei tre ospedali (osp1 29.7%, osp2 29.9%, osp3 27.8%); il test χ² non evidenzia differenze (χ²(2)=1.10, p=0.578) e l’effetto è trascurabile (Cramér’s V=0.02). I confronti a coppie, corretti per Holm, confermano l’assenza di differenze significative.

9.2) hp: La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione come popolazione sono stati utilizzati gli standard internazionali INTERGROWTH-21st per peso e lunghezza neonatale per età gestazionale e sesso (villar et al., lancet 2014; intergrowth-21st newborn size standards) assumendo per i nati un valore medio di riferimento di circa 3300 g e 50 cm. nel mio campione il peso medio alla nascita (3284 g) è molto vicino al valore di riferimento considerato e le differenze non risultano significative, per cui c’è coerenza con la popolazione attesa. per la lunghezza invece osservo una media leggermente più bassa (49.5 cm vs. 50 cm), con differenza statisticamente significativa (di 5 mm). i dati descrivono una popolazione di neonati in linea con gli standard internazionali per il peso, con una liev riduzione della lunghezza.

9.3) hp: le misure antropometriche sono diverse tra i sessi nel campione analizzato, i maschi risultano mediamente più pesanti e più lunghi, con diametro craniale maggiore rispetto alle femmine. in dettaglio: peso 3408 g vs 3161 g (Δ=+247 g); lunghezza 499.7 mm vs 489.8 mm (delta=+9.9 mm, quindi circa 1 cm); diametro craniale 342.4 mm vs 337.6 mm (delta=+4.8 mm, circa 0.48 cm). le differenze sono statisticamente molto solide (welch e wilcoxon: p < 0.001 per tutte) e hanno un’entità piccola–media (cohen’s d: 0.48 per il peso, 0.38 per la lunghezza, 0.30 per il diametro craniale).

# tabella ospedale × tipo di parto
tab <- table(dati$ospedale, dati$tipo_parto)
tab
##       
##        Cesareo Naturale
##   osp1     242      574
##   osp2     254      595
##   osp3     232      603
chi <- chisq.test(tab, correct = FALSE)
chi
## 
##  Pearson's Chi-squared test
## 
## data:  tab
## X-squared = 1.0972, df = 2, p-value = 0.5778
N <- sum(tab); r <- nrow(tab); c <- ncol(tab)
cramers_v <- as.numeric(sqrt((chi$statistic/N) / (min(r-1, c-1))))
cramers_v
## [1] 0.02094948
prop_ces <- prop.table(tab, 1)[, "Cesareo"]
prop_ces
##      osp1      osp2      osp3 
## 0.2965686 0.2991755 0.2778443
hosp <- rownames(tab); pairs <- t(combn(hosp, 2))
p_raw <- apply(pairs, 1, function(idx){
  a <- idx[1]; b <- idx[2]
  prop.test(x = c(tab[a,"Cesareo"], tab[b,"Cesareo"]),
            n = c(sum(tab[a,]),       sum(tab[b,])) )$p.value
})
p_holm <- p.adjust(p_raw, method = "holm")
cbind(pairs, p_value = round(p_raw,5), p_adj_holm = round(p_holm,5))
##                    p_value   p_adj_holm
## [1,] "osp1" "osp2" "0.95003" "1"       
## [2,] "osp1" "osp3" "0.43164" "1"       
## [3,] "osp2" "osp3" "0.3617"  "1"
# 9.2) media del peso e della lunghezza vs. valori di riferimento di popolazione (letteratura citata sorpa)

mu_peso_pop_g    <- 3300   # 3.3 kg
mu_lung_pop_mm   <- 500    # 50 cm

x_peso <- dati$peso_neonato
x_lung <- dati$lunghezza_neonato

# t-test a 1 camp.
t_peso <- t.test(x_peso, mu = mu_peso_pop_g)
t_lung <- t.test(x_lung, mu = mu_lung_pop_mm)

# wilcoxon 
w_peso <- suppressWarnings(wilcox.test(x_peso, mu = mu_peso_pop_g, exact = FALSE))
w_lung <- suppressWarnings(wilcox.test(x_lung, mu = mu_lung_pop_mm, exact = FALSE))

one_sample_tbl <- data.frame(
  variabile       = c("peso (g)", "lunghezza (mm)"),
  media_campione  = c(mean(x_peso, na.rm = TRUE),
                      mean(x_lung, na.rm = TRUE)),
  target_pop      = c(mu_peso_pop_g, mu_lung_pop_mm),
  diff            = c(mean(x_peso, na.rm = TRUE) - mu_peso_pop_g,
                      mean(x_lung, na.rm = TRUE) - mu_lung_pop_mm),
  t_stat          = c(unname(t_peso$statistic),
                      unname(t_lung$statistic)),
  t_p_value       = c(t_peso$p.value, t_lung$p.value),
  wilcox_p_value  = c(w_peso$p.value, w_lung$p.value)
)

knitr::kable(one_sample_tbl, digits = 4,
             caption = "confronto con valori di riferimento di popolazione (t a 1 campione + wilcoxon)")
confronto con valori di riferimento di popolazione (t a 1 campione + wilcoxon)
variabile media_campione target_pop diff t_stat t_p_value wilcox_p_value
peso (g) 3284.081 3300 -15.9192 -1.5160 0.1296 0.9612
lunghezza (mm) 494.692 500 -5.3080 -10.0841 0.0000 0.0000 
# 9.3) misure antropometriche diverse tra sessi (welch + wilcoxon + cohen)

cohen_d <- function(x, y){
  nx <- sum(!is.na(x)); ny <- sum(!is.na(y))
  sp <- sqrt(((nx - 1) * stats::var(x, na.rm = TRUE) + (ny - 1) * stats::var(y, na.rm = TRUE)) / (nx + ny - 2))
  (mean(y, na.rm = TRUE) - mean(x, na.rm = TRUE)) / sp
}

vars_test <- c("peso_neonato", "lunghezza_neonato", "diametro_craniale")

sex_tbl <- lapply(vars_test, function(v){
  x <- dati[dati$sesso == "F", v, drop = TRUE]
  y <- dati[dati$sesso == "M", v, drop = TRUE]
  t_res <- t.test(y, x, var.equal = FALSE)         # welch
  w_res <- suppressWarnings(wilcox.test(y, x, exact = FALSE))
  data.frame(
    variabile = v,
    media_F = mean(x, na.rm = TRUE),
    media_M = mean(y, na.rm = TRUE),
    diff_M_meno_F = mean(y, na.rm = TRUE) - mean(x, na.rm = TRUE),
    welch_p = as.numeric(t_res$p.value),
    wilcoxon_p = as.numeric(w_res$p.value),
    cohen_d = cohen_d(x, y)
  )
}) |>
  dplyr::bind_rows()

knitr::kable(sex_tbl, digits = 3,
             caption = "differenze f vs m: medie,p-value (welch/wilcoxon) ed effect size (cohen)")
differenze f vs m: medie,p-value (welch/wilcoxon) ed effect size (cohen)
variabile media_F media_M diff_M_meno_F welch_p wilcoxon_p cohen_d
peso_neonato 3161.132 3408.215 247.083 0 0 0.484
lunghezza_neonato 489.764 499.667 9.903 0 0 0.383
diametro_craniale 337.633 342.449 4.816 0 0 0.296

2.2) Creazione del Modello di Regressione

  1. preparazione dati per il modeling

  2. modelli candidati: ols_full: lineare multiplo “tutte le variabili”, ols_int: aggiungo interazione settimane × fumo (hp clinica), ols_ns: 1consento non-linearità di settimane (spline naturale df=3) - modello alternativo in più per non-linearità

  3. selezione modello (aic/bic + stepwise)

  • in-sample (aic/bic): migliore ols_ns (aic 28091.2; bic 28169.6). stepAIC vicino (aic 28096.1).
  • scelta predittiva (hold-out): migliore ols_int (rmse 279.96 g, r² 0.672). differenze tra modelli < 1 g di rmse => prestazioni quasi identiche.
  • scelta operativa: ols_int (priorità alla predizione); ols_ns preferibile solo se si privilegia parsimonia/aic.
  1. qualità predizione (rmse/mae/r2) sul test:diagnostica dei residui e cook’s distance per individuare osservazioni influenti; collinearità (vif) -test set: ols_int rmse 279.96, mae 202.14, r² 0.672 ols_full rmse 280.14, mae 202.23, r² 0.671 ols_stepAIC rmse 280.37, mae 202.35, r² 0.671 ols_ns rmse 280.69, mae 201.44, r² 0.670
  • diagnostica residui: lieve eteroschedasticità, qq-plot con code un po’ pesanti: ok per n grande.

  1. interpretazione (coef + vif) settimane: +31.6 g/sett nei non fumatori (p < 1e-12). fumo × settimane: −13.1 g/sett: nei fumatori la pendenza è circa 18.5 g/sett. lunghezza: +11.0 g per mm (p << 0.001). cranio: +9.7 g per mm (p << 0.001). num gravidanze: +11.5 g per gravidanza (p = 0.025). parto cesareo: −26.6 g (p = 0.048, borderline). sesso (m): +76 g. ospedale: differenze piccole, al limite (osp2 −26 g p=0.081; osp3 +30 g p=0.050). età madre: non significativa. vif: molto alto per fumo e fumo:settimane (gvif^(1/(2*df)) curca 27)

  2. diagnostica residui e valori influenti

  • residui vs fitted / scale-location: piccola eteroschedasticità.
  • leverage/cook: 100 osservazioni > soglia 4/n; max cook circa 0.043, nessun caso “critico”, per cui si può non rimuovere
# 10) preparazione: target, predittori

vars_model <- c(
  "peso_neonato","eta_madre","num_gravidanze","fumo",
  "settimane_gestazione","lunghezza_neonato","diametro_craniale",
  "tipo_parto","ospedale","sesso"
)

mod_df <- dati |>
  dplyr::select(all_of(vars_model)) |>
  dplyr::mutate(
    fumo       = droplevels(stats::relevel(fumo, ref = "No")),
    tipo_parto = droplevels(stats::relevel(tipo_parto, ref = "Naturale")),
    ospedale   = droplevels(stats::relevel(ospedale, ref = "osp1")),
    sesso      = droplevels(stats::relevel(sesso, ref = "F"))
  ) |>
  stats::na.omit()  
set.seed(1234)
n   <- nrow(mod_df)
idx <- sample.int(n, size = floor(0.8*n))
train <- mod_df[idx, ]
test  <- mod_df[-idx, ]
# utility metriche
rmse <- function(obs, pred) sqrt(mean((obs - pred)^2))
mae  <- function(obs, pred) mean(abs(obs - pred))
r2   <- function(obs, pred) 1 - sum((obs - pred)^2)/sum((obs - mean(obs))^2)





# 11) fit modelli 
form_full <- peso_neonato ~ eta_madre + num_gravidanze + fumo +
  settimane_gestazione + lunghezza_neonato + diametro_craniale +
  tipo_parto + ospedale + sesso

form_int  <- update(form_full, . ~ . + settimane_gestazione:fumo)

form_ns   <- peso_neonato ~ eta_madre + num_gravidanze + fumo +
  splines::ns(settimane_gestazione, df = 3) +
  lunghezza_neonato + diametro_craniale + tipo_parto + ospedale + sesso

fit_full <- lm(form_full, data = train)
fit_int  <- lm(form_int,  data = train)
fit_ns   <- lm(form_ns,   data = train)






# 12) confronto aic/bic
comp_ic <- data.frame(
  modello = c("ols_full","ols_int","ols_ns"),
  AIC = c(AIC(fit_full), AIC(fit_int), AIC(fit_ns)),
  BIC = c(BIC(fit_full), BIC(fit_int), BIC(fit_ns))
)
knitr::kable(comp_ic, digits = 1, caption = "confronto AIC/BIC tra modelli candidati")
confronto AIC/BIC tra modelli candidati
modello AIC BIC
ols_full 28099.6 28166.8
ols_int 28101.2 28174.0
ols_ns 28091.2 28169.6 
# 12) stepAIC con MASS 
if (requireNamespace("MASS", quietly = TRUE)) {
  fit_step <- MASS::stepAIC(fit_full, direction = "both", trace = FALSE)
  comp_ic <- rbind(
    comp_ic,
    data.frame(modello = "ols_stepAIC", AIC = AIC(fit_step), BIC = BIC(fit_step))
  )
  knitr::kable(comp_ic[order(comp_ic$AIC),], digits = 1,
               caption = "confronto AIC/BIC (incluso stepAIC se disponibile)")
} else {
  fit_step <- NULL
}
confronto AIC/BIC (incluso stepAIC se disponibile)
modello AIC BIC
3 ols_ns 28091.2 28169.6
4 ols_stepAIC 28096.1 28152.1
1 ols_full 28099.6 28166.8
2 ols_int 28101.2 28174.0 
# 13) hold-out su test
pred_full <- predict(fit_full, newdata = test)
pred_int  <- predict(fit_int,  newdata = test)
pred_ns   <- predict(fit_ns,   newdata = test)

sc <- data.frame(
  modello = c("ols_full","ols_int","ols_ns"),
  RMSE = c(rmse(test$peso_neonato, pred_full),
           rmse(test$peso_neonato, pred_int),
           rmse(test$peso_neonato, pred_ns)),
  MAE  = c(mae(test$peso_neonato, pred_full),
           mae(test$peso_neonato, pred_int),
           mae(test$peso_neonato, pred_ns)),
  R2   = c(r2(test$peso_neonato, pred_full),
           r2(test$peso_neonato, pred_int),
           r2(test$peso_neonato, pred_ns))
)

if (!is.null(fit_step)) {
  pred_step <- predict(fit_step, newdata = test)
  sc <- rbind(
    sc,
    data.frame(
      modello = "ols_stepAIC",
      RMSE = rmse(test$peso_neonato, pred_step),
      MAE  = mae(test$peso_neonato, pred_step),
      R2   = r2(test$peso_neonato, pred_step)
    )
  )
}

knitr::kable(
  sc[order(sc$RMSE),],
  digits = 3,
  caption = "performance sul test (più basso RMSE/MAE, più alto R² = modello migliore)"
)
performance sul test (più basso RMSE/MAE, più alto R² = modello migliore)
modello RMSE MAE R2
2 ols_int 279.957 202.143 0.672
1 ols_full 280.144 202.227 0.671
4 ols_stepAIC 280.373 202.351 0.671
3 ols_ns 280.693 201.442 0.670 
# scelta automatica del modello migliore in base al RMSE
best <- sc$modello[which.min(sc$RMSE)]
fit_final <- switch(
  best,
  "ols_full"    = fit_full,
  "ols_int"     = fit_int,
  "ols_ns"      = fit_ns,
  "ols_stepAIC" = if (!is.null(fit_step)) fit_step else fit_full
)

best
## [1] "ols_int"
# 14) coefficienti stimati + vif
coef_tbl <- data.frame(
  term = names(coef(fit_final)),
  estimate = unname(coef(fit_final))
)
knitr::kable(coef_tbl, digits = 3, caption = "coefficienti del modello finale")
coefficienti del modello finale
term estimate
(Intercept) -6729.759
eta_madre 0.727
num_gravidanze 11.483
fumoSi 501.021
settimane_gestazione 31.595
lunghezza_neonato 10.994
diametro_craniale 9.657
tipo_partoCesareo -26.633
ospedaleosp2 -26.017
ospedaleosp3 29.574
sessoM 76.040
fumoSi:settimane_gestazione -13.101 
# multicollinearità
if (requireNamespace("car", quietly = TRUE)) {
  car::vif(fit_final)
}
##                                 GVIF Df GVIF^(1/(2*Df))
## eta_madre                   1.198218  1        1.094631
## num_gravidanze              1.200279  1        1.095573
## fumo                      741.274498  1       27.226357
## settimane_gestazione        1.759611  1        1.326503
## lunghezza_neonato           2.195305  1        1.481656
## diametro_craniale           1.691468  1        1.300564
## tipo_parto                  1.007786  1        1.003886
## ospedale                    1.008175  2        1.002037
## sesso                       1.048562  1        1.023993
## fumo:settimane_gestazione 741.766928  1       27.235398
# 15) grafica
par(mfrow = c(2,2)); plot(fit_final); par(mfrow = c(1,1))

# cook's distance (sservazioni molto influenti)
cook  <- cooks.distance(fit_final)
thr   <- 4 / nrow(train)
infl  <- which(cook > thr)
length(infl)  # quante sopra soglia 4/n
## [1] 100
head(infl, 15)
##  21  23  66  70  73  80  95 161 222 258 289 313 326 336 386 
##  21  23  66  70  73  80  95 161 222 258 289 313 326 336 386
infl_tbl <- data.frame(id_train = infl, cooks = cook[infl]) |>
  dplyr::arrange(dplyr::desc(cooks))
knitr::kable(head(infl_tbl, 20), digits = 4, caption = "osservazioni più influenti (cook > 4/n)")
osservazioni più influenti (cook > 4/n)
id_train cooks
1041 1041 0.0429
1506 1506 0.0251
1719 1719 0.0235
1520 1520 0.0230
439 439 0.0220
95 95 0.0177
222 222 0.0165
1135 1135 0.0164
1578 1578 0.0144
386 386 0.0127
1171 1171 0.0118
557 557 0.0117
886 886 0.0099
1193 1193 0.0084
66 66 0.0072
1952 1952 0.0068
1343 1343 0.0068
648 648 0.0068
1240 1240 0.0066
956 956 0.0061
  1. previsioni e risultati validato il modello lineare multiplo, l’ho utilizzato per produrre previsioni in scenari realistici. (es. neonata da madre alla 3a gravidanza), non fumatrice, parto naturale in ospedale 1, alla 39ª settimana di gestazione, con parametri antropometrici attesi (lunghezza 500 mm, diametro craniale 340 mm). applicando il modello finale, si ottiene: -peso previsto (fit): 3327.4 g (il modello si aspetta un peso dicirca 3.3 kg. -lwr / upr (2791 g – 3863.8 g): intervallo di previsione al 95% (cioè per una neonata con quelle caratteristiche, ci si può attendere un peso in questo range
# 3) previsione esempio con modello finale

# neonata, madre alla 3a gravidanza, 39 settimane,
# non fumatrice, parto naturale, ospedale 1
new_case <- data.frame(
  eta_madre            = 30,
  num_gravidanze       = 2,                                  
  fumo                 = factor("No", levels = levels(mod_df$fumo)),
  settimane_gestazione = 39,
  lunghezza_neonato    = 500,
  diametro_craniale    = 340,
  tipo_parto           = factor("Naturale", levels = levels(mod_df$tipo_parto)),
  ospedale             = factor("osp1",     levels = levels(mod_df$ospedale)),
  sesso                = factor("F",        levels = levels(mod_df$sesso))
)

# previsione  peso con intervallo di previsione aL 95%
pred_example <- predict(fit_final, newdata = new_case,
                        interval = "prediction", level = 0.95)

knitr::kable(pred_example, digits = 1,
             caption = "stima del peso alla nascita")
stima del peso alla nascita
fit lwr upr
3327.4 2791 3863.8
  1. Visualizzazioni dal grafico settimane × fumo si vede che più la gravidanza prosegue, più cresce il peso atteso. la crescita è più ripida per le non fumatrici, mentre nei fumatori la curva è più piatta. oltre circa 38 settimane le differenze diventano visibili, anche se vanno lette con prudenza perché la popolazione fumatrici è poco nuemrosa.

il confronto tra peso osservato e peso predetto mostra i punti concentrati intorno alla diagonale: il modello “ci prende” abbastanza bene.

il grafico dei residui rispetto ai valori predetti mostra residui distribuiti intorno allo 0, senza pattern particolari, c’è solo un po’ più di dispersione ai pesi più alti, ma non sembra mettere in crisi il modello.

# 4) visualizzazioni: effetto fumo×settimane, performance e residui

# 4.1) effetto settimane × fumo sul peso previsto

newdf <- expand.grid(
eta_madre            = median(train$eta_madre, na.rm = TRUE),
num_gravidanze       = median(train$num_gravidanze, na.rm = TRUE),
fumo                 = factor(c("No","Si"), levels = levels(train$fumo)),
settimane_gestazione = seq(35, 41, by = 0.5),
lunghezza_neonato    = median(train$lunghezza_neonato, na.rm = TRUE),
diametro_craniale    = median(train$diametro_craniale, na.rm = TRUE),
tipo_parto           = factor("Naturale", levels = levels(train$tipo_parto)),
ospedale             = factor("osp1",      levels = levels(train$ospedale)),
sesso                = factor("F",         levels = levels(train$sesso))
)

newdf$pred <- predict(fit_final, newdata = newdf)

p1 <- ggplot(newdf, aes(x = settimane_gestazione, y = pred, linetype = fumo)) +
geom_line() +
labs(
title = "peso previsto - settimane di gestazione per stato di fumo",
x = "settimane di gestazione",
y = "peso previsto (g)",
linetype = "fumo materno"
) +
theme_minimal()



#4.2) peso osservato vs peso predetto (test set)

test$pred <- predict(fit_final, newdata = test)

p2 <- ggplot(test, aes(x = pred, y = peso_neonato)) +
geom_point(alpha = 0.4) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
labs(
title = "peso osservato vs peso predetto (test set)",
x = "peso predetto (g)",
y = "peso osservato (g)"
) +
theme_minimal()



# 4.3) residui vs. peso predetto (test set)

test$resid <- test$peso_neonato - test$pred

p3 <- ggplot(test, aes(x = pred, y = resid)) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_point(alpha = 0.4) +
labs(
title = "residui vs peso predetto",
x = "peso predetto (g)",
y = "residuo (g)"
) +
theme_minimal()



print(p1)

print(p2)

print(p3)