Modello Statistico per la Previsione del Peso Neonatale

Obiettivo del progetto:

L’azienda Neonatal Health Solutions desidera creare un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali. Il progetto mira a migliorare la gestione delle gravidanze ad alto rischio, ottimizzare le risorse ospedaliere e garantire migliori risultati per la salute neonatale.

Il progetto si inserisce all’interno di un contesto di crescente attenzione verso la prevenzione delle complicazioni neonatali. La possibilità di prevedere il peso alla nascita dei neonati rappresenta un’opportunità fondamentale per migliorare la pianificazione clinica e ridurre i rischi associati a nascite problematiche, come parti prematuri o neonati con basso peso. Di seguito, i principali benefici che questo progetto porterà all’azienda e al settore sanitario:

Miglioramento delle previsioni cliniche.
Ottimizzazione delle risorse ospedaliere.
Prevenzione e identificazione dei fattori di rischio.
Valutazione delle pratiche ospedaliere.
Supporto alla pianificazione strategica.

1. Raccolta dei Dati e Struttura del Dataset:

Per costruire il modello predittivo, abbiamo raccolto dati su 2500 neonati provenienti da tre ospedali. Le variabili raccolte includono:

Età della madre: Misura dell’età in anni.
Numero di gravidanze: Quante gravidanze ha avuto la madre.
Fumo materno: Un indicatore binario (0=non fumatrice, 1=fumatrice).
Durata della gravidanza: Numero di settimane di gestazione.
Peso del neonato: Peso alla nascita in grammi.
Lunghezza e diametro del cranio: Lunghezza del neonato e diametro craniale, misurabili anche durante la gravidanza tramite ecografie.
Tipo di parto: Naturale o cesareo.
Ospedale di nascita: Ospedale 1, 2 o 3.
Sesso del neonato: Maschio (M) o femmina (F).

L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature.

Caricamento del dataset “neonati.csv” in un dataframe R denominato df e visualizzazione dell’intestazione in formato tabella per verificarne il caricamento.

# Caricamento del dataset
df <- read.csv("neonati.csv", sep = ",", stringsAsFactors = FALSE)

# Visualizzazione dell'intestazione
knitr::kable(
  head(df),
  caption = "Prime 6 osservazioni del dataset neonati"
)

Prime 6 osservazioni del dataset neonati
Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
26	0	42	3380	490	325	Nat	osp3	M
21	2	39	3150	490	345	Nat	osp1	F
34	3	38	3640	500	375	Nat	osp2	M
28	1	41	3690	515	365	Nat	osp2	M
20	0	38	3700	480	335	Nat	osp3	F
32	0	40	3200	495	340	Nat	osp2	F

Esaminazione della struttura del dataframe per identificare: il tipo di ciascuna variabile, la presenza di eventuali anomalie o codifiche non appropriate e di valori mancanti o incoerenti al fine di assicurare la qualità dei dati utilizzati.

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   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ 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  : chr  "Nat" "Nat" "Nat" "Nat" ...
##  $ Ospedale    : chr  "osp3" "osp1" "osp2" "osp2" ...
##  $ Sesso       : chr  "M" "F" "M" "M" ...

Il dataset contiene 2500 osservazioni.

colSums(is.na(df))

##   Anni.madre N.gravidanze    Fumatrici   Gestazione         Peso    Lunghezza 
##            0            0            0            0            0            0 
##       Cranio   Tipo.parto     Ospedale        Sesso 
##            0            0            0            0

Il dataset non contiene valori nulli.

df[duplicated(df), ]

##  [1] Anni.madre   N.gravidanze Fumatrici    Gestazione   Peso        
##  [6] Lunghezza    Cranio       Tipo.parto   Ospedale     Sesso       
## <0 rows> (or 0-length row.names)

Il dataset non contiene valori duplicati.

Caricamento delle librerie necessarie per condurre l’analisi.

library(dplyr)
library(kableExtra)
library(tidyr)
library(corrplot)
library(car)
library(broom)
library(ggplot2)
library(plotly)
library(MASS)
library(lmtest)
library(gridExtra)

Imputazione dei valori anomali di Anni.madre:

# Summary prima dell'imputazione
summary(df$Anni.madre)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   25.00   28.00   28.16   32.00   46.00

Il summary() della variabile Anni.madre rivela la presenza di valori pari a 0 e 1, chiaramente non plausibili per l’età di una madre. Questi valori vengono imputati prima dell’analisi inferenziale con la media della variabile calcolata escludendoli, in modo che tutte le statistiche e i modelli successivi operino su dati coerenti.

# Identificazione valori anomali
cat("Valori anomali di Anni.madre (<=1):",
    df$Anni.madre[df$Anni.madre <= 1], "\n")

## Valori anomali di Anni.madre (<=1): 1 0

# Calcolo della media (escludendo i valori anomali)
media_anni <- mean(df$Anni.madre[df$Anni.madre > 1])

# Imputazione
df$Anni.madre[df$Anni.madre <= 1] <- round(media_anni, 0)

cat("Media utilizzata per l'imputazione:", round(media_anni, 1), "anni\n")

## Media utilizzata per l'imputazione: 28.2 anni

# Summary dopo l'imputazione
summary(df$Anni.madre)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13.00   25.00   28.00   28.19   32.00   46.00

2. Analisi e Modellizzazione:

2.1 Analisi Preliminare: nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.

Inoltre si saggeranno le seguenti ipotesi con i test adatti:

In alcuni ospedali si fanno più parti cesarei.
La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione.
Le misure antropometriche sono significativamente diverse tra i due sessi.

# Statistiche descrittive per le variabili numeriche
knitr::kable(
  summary(df[, c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")]),
  caption = "Statistiche descrittive per le variabili numeriche"
)

Statistiche descrittive per le variabili numeriche
Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio
Min. :13.00	Min. : 0.0000	Min. :25.00	Min. : 830	Min. :310.0	Min. :235
1st Qu.:25.00	1st Qu.: 0.0000	1st Qu.:38.00	1st Qu.:2990	1st Qu.:480.0	1st Qu.:330
Median :28.00	Median : 1.0000	Median :39.00	Median :3300	Median :500.0	Median :340
Mean :28.19	Mean : 0.9812	Mean :38.98	Mean :3284	Mean :494.7	Mean :340
3rd Qu.:32.00	3rd Qu.: 1.0000	3rd Qu.:40.00	3rd Qu.:3620	3rd Qu.:510.0	3rd Qu.:350
Max. :46.00	Max. :12.0000	Max. :43.00	Max. :4930	Max. :565.0	Max. :390

# Distribuzioni grafiche
df_hist <- df[, c("Peso", "Gestazione", "Lunghezza", "Cranio", "Anni.madre", "N.gravidanze")] %>%
  pivot_longer(cols = everything(), names_to = "Variabile", values_to = "Valore") %>%
  mutate(Variabile = factor(Variabile,
         levels = c("Peso", "Gestazione", "Lunghezza", "Cranio", "Anni.madre", "N.gravidanze"),
         labels = c("Distribuzione del peso", "Settimane di gestazione", "Lunghezza",
                     "Diametro craniale", "Età della madre", "Numero gravidanze")))

ggplot(df_hist, aes(x = Valore)) +
  geom_histogram(fill = "#1f78b4", color = "white", bins = 20) +
  facet_wrap(~ Variabile, scales = "free", ncol = 3) +
  labs(title = "Distribuzioni grafiche delle variabili numeriche",
       x = "", y = "Frequenza") +
  theme_minimal(base_size = 10) +
  theme(
    plot.title = element_text(face = "bold"),
    strip.text = element_text(face = "plain"),
    legend.position = "none"
  )

La tabella seguente riassume frequenze assolute e proporzioni per le quattro variabili categoriche del dataset.

# Frequenze e proporzioni per le variabili categoriche
cat_vars <- list(
  data.frame(Variabile = "Fumatrici",  Categoria = names(table(df$Fumatrici)),  Frequenza = as.integer(table(df$Fumatrici))),
  data.frame(Variabile = "Sesso",      Categoria = names(table(df$Sesso)),      Frequenza = as.integer(table(df$Sesso))),
  data.frame(Variabile = "Tipo.parto", Categoria = names(table(df$Tipo.parto)), Frequenza = as.integer(table(df$Tipo.parto))),
  data.frame(Variabile = "Ospedale",   Categoria = names(table(df$Ospedale)),   Frequenza = as.integer(table(df$Ospedale)))
)

cat_df <- bind_rows(cat_vars)
cat_df$Proporzione <- round(cat_df$Frequenza / nrow(df), 4)

kable(cat_df,
      caption = "Frequenze e proporzioni per le variabili categoriche",
      col.names = c("Variabile", "Categoria", "Frequenza", "Proporzione")) |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE) |>
  collapse_rows(columns = 1, valign = "top")

Frequenze e proporzioni per le variabili categoriche
Variabile	Categoria	Frequenza	Proporzione
Fumatrici	0	2396	0.9584
Fumatrici	1	104	0.0416
Sesso	F	1256	0.5024
Sesso	M	1244	0.4976
Tipo.parto	Ces	728	0.2912
Tipo.parto	Nat	1772	0.7088
Ospedale	osp1	816	0.3264
	osp2	849	0.3396
	osp3	835	0.3340

Si osserva come il campione sia bilanciato per sesso (circa 50/50) e distribuito in modo uniforme tra i tre ospedali (circa 33% ciascuno). I parti naturali rappresentano la netta maggioranza (70.9%) rispetto ai cesarei (29.1%), mentre le madri fumatrici costituiscono solo il 4.2% del campione.

l grafico seguente mostra la distribuzione delle frequenze assolute dei parti cesarei e naturali nei tre ospedali.

# Barplot tipo parto per ospedale (ggplot2)
tab_bar <- as.data.frame(table(df$Ospedale, df$Tipo.parto))
colnames(tab_bar) <- c("Ospedale", "Tipo_parto", "Frequenza")

ggplot(tab_bar, aes(x = Ospedale, y = Frequenza, fill = Tipo_parto)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.8)) +
  geom_text(aes(label = Frequenza),
            position = position_dodge(width = 0.8),
            vjust = -0.3,
            size = 3) +
  scale_fill_manual(values = c(
    "Ces" = "#33a02c",
    "Nat" = "#1f78b4"
  )) +
  labs(
    title = "Distribuzione del tipo di parto per ospedale",
    x = "Ospedale",
    y = "Frequenza",
    fill = "Tipo di parto"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "right"
  )

Si può osservare come le proporzioni siano sostanzialmente omogenee tra le strutture: i parti cesarei oscillano tra 232 e 254, mentre quelli naturali tra 574 e 603, confermando visivamente l’assenza di associazione tra ospedale e tipo di parto che verrà poi verificata formalmente con il test chi-quadrato.

l grafico seguente confronta la distribuzione del peso neonatale tra madri fumatrici e non fumatrici tramite curve di densità sovrapposte.

# Density plot fumatrici vs non fumatrici
ggplot(df, aes(x = Peso, fill = factor(Fumatrici, labels = c("Non fumatrice", "Fumatrice")))) +
  geom_density(alpha = 0.5) +
  scale_fill_manual(values = c("Non fumatrice" = "#66c2a5", "Fumatrice" = "#fc8d62")) +
  labs(title = "Distribuzione del peso: fumatrici vs non fumatrici",
       x = "Peso (g)", y = "Densità", fill = "") +
  theme_minimal(base_size = 10) +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "right")

Si può osservare come la curva delle fumatrici risulti leggermente spostata verso sinistra rispetto a quella delle non fumatrici, suggerendo un peso tendenzialmente inferiore nei neonati di madri fumatrici. Questo effetto verrà poi quantificato nel modello di regressione.

Vengono ora identificati gli outlier tramite il metodo IQR. La tabella seguente riporta, per ogni variabile numerica, le soglie IQR, il numero di outlier e la percentuale sul totale.

# Funzione per identificare outlier con metodo IQR
find_outliers_iqr <- function(x) {
  Q1 <- quantile(x, 0.25, na.rm = TRUE)
  Q3 <- quantile(x, 0.75, na.rm = TRUE)
  IQR_val <- Q3 - Q1
  lower <- Q1 - 1.5 * IQR_val
  upper <- Q3 + 1.5 * IQR_val
  list(idx = which(x < lower | x > upper),
       lower = lower, upper = upper, Q1 = Q1, Q3 = Q3, IQR = IQR_val)
}

# Variabili numeriche da analizzare
vars_num <- c("Peso", "Lunghezza", "Cranio", "Gestazione", "Anni.madre")

# Creazione tabella riassuntiva
out_table <- lapply(vars_num, function(v) {
  res <- find_outliers_iqr(df[[v]])
  data.frame(
    Variabile      = v,
    Q1             = round(res$Q1, 2),
    Q3             = round(res$Q3, 2),
    IQR            = round(res$IQR, 2),
    Soglia_inf     = round(res$lower, 2),
    Soglia_sup     = round(res$upper, 2),
    N_outlier      = length(res$idx),
    Perc_outlier   = paste0(round(length(res$idx) / nrow(df) * 100, 2), "%")
  )
}) |> bind_rows()

# Tabella formattata
out_table |>
  kable("html", caption = "Outlier individuati con metodo IQR — dettaglio per variabile") |>
  kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE) |>
  column_spec(9, width = "25em")

Outlier individuati con metodo IQR — dettaglio per variabile
	Variabile	Q1	Q3	IQR	Soglia_inf	Soglia_sup	N_outlier	Perc_outlier
25%…1	Peso	2990	3620	630	2045.0	4565.0	69	2.76%
25%…2	Lunghezza	480	510	30	435.0	555.0	59	2.36%
25%…3	Cranio	330	350	20	300.0	380.0	48	1.92%
25%…4	Gestazione	38	40	2	35.0	43.0	67	2.68%
25%…5	Anni.madre	25	32	7	14.5	42.5	11	0.44%

# Boxplot outliers (ggplot2)

df_long <- df[, vars_num] %>%
  pivot_longer(cols = everything(), names_to = "Variabile", values_to = "Valore")

ggplot(df_long, aes(x = Variabile, y = Valore)) +
  geom_boxplot(fill = "lightblue", width = 0.4, 
               outlier.colour = "red", outlier.shape = 16, outlier.size = 2) +
  facet_wrap(~ Variabile, scales = "free", ncol = 3) +
  labs(title = "Boxplot delle variabili numeriche") +
  theme_minimal(base_size = 10) +
  theme(
    plot.title = element_text(face = "bold"),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    strip.text = element_text(face = "bold")
  )

Gli outlier identificati vengono mantenuti nel dataset in questa fase. Eventuali osservazioni influenti verranno studiate nella diagnostica del modello (Sezione 2.3).

Si può notare quindi come le variabili numeriche del dataset neonatale mostrano distribuzioni complessivamente regolari e coerenti con dati clinici reali, pur presentando un numero non trascurabile di outlier (soprattutto per peso, lunghezza, cranio e settimane di gestazione) che riflettono la naturale variabilità biologica e la presenza di casi estremi come prematurità o neonati macrosomici. L’assenza di valori mancanti o duplicati conferma inoltre una buona qualità del dataset, rendendolo idoneo per le successive analisi inferenziali e per la costruzione del modello predittivo del peso neonatale.

In seguito a quanto emerso dall’analisi preliminare possiamo saggiare l’ipotesi che in alcuni ospedali si fanno più parti cesarei utlizzando il test di chi-quadrato. Considerate le due variabili categoriche: Tipo.parto (Nat / Ces) e Ospedale (Osp1 / osp2 / osp3), l’obiettivo è verificare se la frequenza dei parti cesarei dipende dall’ospedale.
L’ipotesi nulla \(H_0\) è che la proporzione di parti cesarei è la stessa nei tre ospedali e quindi non c’è associazione tra ospedale e tipo di parto.
L’ipotesi alternativa \(H_1\) è che la proporzione di parti cesarei è diversa in almeno un ospedale e quindi esiste un’associazione tra ospedale e tipo di parto.

# Tabella di contingenza
tab <- table(df$Ospedale, df$Tipo.parto)
tab_df <- as.data.frame.matrix(tab)

# Test chi-quadrato
chisq.test(tab)

## 
##  Pearson's Chi-squared test
## 
## data:  tab
## X-squared = 1.0972, df = 2, p-value = 0.5778

# Proporzioni per ospedale
prop_tab <- round(prop.table(tab, margin = 1), 3)

# Creazione di un data frame
combined_df <- data.frame(
  Cesareo_freq = tab[, "Ces"],
  Naturale_freq = tab[, "Nat"],
  Cesareo_prop = prop_tab[, "Ces"],
  Naturale_prop = prop_tab[, "Nat"]
)

# Tabella
kable(combined_df,
      caption = "Frequenze e proporzioni per tipo di parto e ospedale",
      col.names = c("Ospedale",
                    "Cesareo (n)", "Naturale (n)",
                    "Cesareo (%)", "Naturale (%)")) |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE)

Frequenze e proporzioni per tipo di parto e ospedale
Ospedale	Cesareo (n)	Naturale (n)	Cesareo (%)	Naturale (%)
osp1	242	574	0.297	0.703
osp2	254	595	0.299	0.701
osp3	232	603	0.278	0.722

In base ai risultati del test del chi-quadrato, in particolare al p-value = 0.5778 > 0.05, non si rifiuta l’ipotesi nulla H_0 e quindi non ci sono differenze significative tra gli ospedali.

Saggiamo l’ipotesi successiva, ovvero che la media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione, utilizzando il Test-t per singolo campione e assumendo (da linee guida cliniche) che il peso medio della popolazione sia pari a 3300 g e che la lunghezza media della popolazione sia pari a 500 mm. L’ipotesi nulla H_0 prevede che il peso campione sia uguale al peso della popolazione e che la lunghezza campione sia uguale alla lunghezza della popolazione.

# Valori medi della popolazione
mu_peso <- 3300
mu_lunghezza <- 500

# Test t per un campione
t_peso <- t.test(df$Peso, mu = mu_peso)
t_lung <- t.test(df$Lunghezza, mu = mu_lunghezza)

t_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.081

t_lung

## 
##  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.692

# Creazione tabella
summary_df <- data.frame(
  Variabile = c("Peso", "Lunghezza"),
  Media_campione = c(mean(df$Peso), mean(df$Lunghezza)),
  Media_popolazione = c(mu_peso, mu_lunghezza),
  Differenza = c(mean(df$Peso) - mu_peso,
                 mean(df$Lunghezza) - mu_lunghezza),
  T_value = c(t_peso$statistic, t_lung$statistic),
  P_value = c(t_peso$p.value, t_lung$p.value),
  Decisione = c(
    ifelse(t_peso$p.value < 0.05, "Rifiuto H0", "Non rifiuto H0"),
    ifelse(t_lung$p.value < 0.05, "Rifiuto H0", "Non rifiuto H0")
  )
)

# Tabella finale
kable(summary_df,
      digits = 4,
      caption = "Risultati dei test t per il confronto con le medie della popolazione") |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE)

Risultati dei test t per il confronto con le medie della popolazione
Variabile	Media_campione	Media_popolazione	Differenza	T_value	P_value	Decisione
Peso	3284.081	3300	-15.9192	-1.5160	0.1296	Non rifiuto H0
Lunghezza	494.692	500	-5.3080	-10.0841	0.0000	Rifiuto H0

In base ai risultati ottenuti dal t-test è possibile concludere che per la media del peso dei neonati non si rifiuta l’ipotesi nulla, con un p-value = 0.1296 > 0.05 e quindi la media del campione non differisce significativamente da quella della popolazione. Invece per la lunghezza media del neonati l’ipotesi nulla viene rifutata, con un p-value = 0 < 0.05 e quindi la media del campione è significativamente diversa da quella della popolazione.

La terza e ultima ipotesi da saggiare dice che le misure antropometriche (Peso, Lunghezza, Cranio) sono significativamente diverse tra i due sessi. Si utilizza anche in questo caso un Test-t per campioni indipendenti. La formulazione delle ipotesi ne prevede una nulla H_0: La media della variabile è uguale nei due sessi e una alternativa H_1: La media della variabile è diversa nei due sessi.

# Test t per campioni indipendenti
t_peso  <- t.test(Peso ~ Sesso, data = df)
t_lung  <- t.test(Lunghezza ~ Sesso, data = df)
t_cranio <- t.test(Cranio ~ Sesso, data = df)

t_peso

## 
##  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.215

t_lung

## 
##  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.6672

t_cranio

## 
##  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.4486

summary_sex <- data.frame(
  Variabile = c("Peso", "Lunghezza", "Cranio"),
  Media_M = c(mean(df$Peso[df$Sesso=="M"]),
              mean(df$Lunghezza[df$Sesso=="M"]),
              mean(df$Cranio[df$Sesso=="M"])),
  Media_F = c(mean(df$Peso[df$Sesso=="F"]),
              mean(df$Lunghezza[df$Sesso=="F"]),
              mean(df$Cranio[df$Sesso=="F"])),
  T_value = c(t_peso$statistic,
              t_lung$statistic,
              t_cranio$statistic),
  P_value = c(t_peso$p.value,
              t_lung$p.value,
              t_cranio$p.value),
  Decisione = c(
    ifelse(t_peso$p.value < 0.05, "Rifiuto H0", "Non rifiuto H0"),
    ifelse(t_lung$p.value < 0.05, "Rifiuto H0", "Non rifiuto H0"),
    ifelse(t_cranio$p.value < 0.05, "Rifiuto H0", "Non rifiuto H0")
  )
)

kable(summary_sex,
      digits = 4,
      caption = "Confronto delle misure antropometriche tra i due sessi (t-test)") |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE)

Confronto delle misure antropometriche tra i due sessi (t-test)
Variabile	Media_M	Media_F	T_value	Decisione
Peso	3408.2154	3161.1322	-12.1061	Rifiuto H0
Lunghezza	499.6672	489.7643	-9.5820	Rifiuto H0
Cranio	342.4486	337.6330	-7.4102	Rifiuto H0

A supporto dei t-test effettuati tra i sessi vengono rappresentati i risultati tramite violin plot con boxplot sovrapposti, che permettono di visualizzare sia la distribuzione complessiva dei dati sia i principali indicatori di posizione.

# Boxplot confronto M vs F con violin + jitter
p1 <- ggplot(df, aes(x = Sesso, y = Peso, fill = Sesso)) +
  geom_violin(alpha = 0.3, width = 0.9, trim = FALSE) +
  geom_boxplot(width = 0.15, alpha = 0.7, outlier.colour = "red", outlier.size = 1.5) +
  scale_color_viridis_d(option = "H")+
  labs(title = "Peso per sesso", y = "Peso (g)", x = "") +
  theme_minimal(base_size = 10) +
  theme(plot.title = element_text(face = "bold"), legend.position = "none")

p2 <- ggplot(df, aes(x = Sesso, y = Lunghezza, fill = Sesso)) +
  geom_violin(alpha = 0.3, width = 0.9, trim = FALSE) +
  geom_boxplot(width = 0.15, alpha = 0.7, outlier.colour = "red", outlier.size = 1.5) +
  scale_color_viridis_d(option = "H")+
  labs(title = "Lunghezza per sesso", y = "Lunghezza (mm)", x = "") +
  theme_minimal(base_size = 10) +
  theme(plot.title = element_text(face = "bold"), legend.position = "none")

p3 <- ggplot(df, aes(x = Sesso, y = Cranio, fill = Sesso)) +
  geom_violin(alpha = 0.3, width = 0.9, trim = FALSE) +
  geom_boxplot(width = 0.15, alpha = 0.7, outlier.colour = "red", outlier.size = 1.5) +
  scale_color_viridis_d(option = "H") +
  labs(title = "Cranio per sesso", y = "Cranio (mm)", x = "") +
  theme_minimal(base_size = 10) +
  theme(plot.title = element_text(face = "bold"), legend.position = "none")

grid.arrange(p1, p2, p3, ncol = 3)

Si può osservare come per tutte e tre le misure antropometriche — peso, lunghezza e diametro cranico — la distribuzione dei maschi risulti sistematicamente spostata verso valori più alti rispetto a quella delle femmine, con mediane visibilmente differenti. La forma dei violin plot evidenzia inoltre una maggiore concentrazione attorno ai valori centrali per entrambi i sessi, con code più pronunciate verso i valori estremi inferiori (neonati prematuri o sottopeso).

In seguito all’analisi condotta è emerso che l’ipotesi nulla viene rifiutata in tutti e tre i casi essendo il p-value < 0.05 per tutte e tre le misure. Quindi le misure antropometriche sono significativamente diverse tra maschi e femmine. In particolare i maschi hanno un peso medio (differenza media di circa 247 g), una lunghezza media e una circonferenza cranica media significativamente maggiori rispetto alle femmine.

2.2 Creazione del modello di regressione: verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo, potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato ed eventuali interazioni. Ad esempio, ci aspettiamo che una maggiore durata della gestazione aumenterebbe in media il peso del neonato.

Il processo di modellazione, come prima cosa, richiede la predisposizione di un modello completo ma solo con variabili che hanno senso logico in ottica predittiva. Per questo veongono escluse a priori le vairiabili Tipo.parto e Ospedale:

Tipo.parto viene deciso al momento del parto (spesso in emergenza); non è una variabile nota in anticipo e quindi non ha valore predittivo ex-ante.
Ospedale è una variabile logistica che non ha alcuna influenza biologica sul peso del neonato. Inoltre il test chi-quadrato (p = 0.578) ha confermato che la distribuzione dei parti non differisce tra ospedali.

Calcolo della correlazione tra le variabili e visualizzazione della matrice al fine di identificare la presenza di eventuali forti correlazioni. Si considerano a tal fine le variabili numeriche (Anni.madre, N.gravidanze, Gestazione, Peso, Lunghezza, Cranio).

# Selezione variabili numeriche
vars_num <- df[, c("Anni.madre", "N.gravidanze", "Gestazione",
                   "Peso", "Lunghezza", "Cranio")]

# Matrice di correlazione
cor_mat <- round(cor(vars_num), 3)

# Tabella kable
kable(cor_mat,
      caption = "Matrice di correlazione tra le variabili numeriche",
      col.names = colnames(cor_mat)) |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE)

Matrice di correlazione tra le variabili numeriche
	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio
Anni.madre	1.000	0.383	-0.135	-0.024	-0.065	0.016
N.gravidanze	0.383	1.000	-0.101	0.002	-0.060	0.039
Gestazione	-0.135	-0.101	1.000	0.592	0.619	0.461
Peso	-0.024	0.002	0.592	1.000	0.796	0.705
Lunghezza	-0.065	-0.060	0.619	0.796	1.000	0.603
Cranio	0.016	0.039	0.461	0.705	0.603	1.000

Per analizzare le relazioni lineari tra le variabili numeriche si calcola la matrice di correlazione di Pearson, rappresentata graficamente tramite una heatmap con i coefficienti sovrapposti.

corrplot(cor_mat, method = "color", type = "upper",
         addCoef.col = "black", tl.col = "black", number.cex = 0.8)

Dalla matrice emergono correlazioni forti e positive tra le misure antropometriche: Peso–Lunghezza (0.80), Peso–Cranio (0.70) e Lunghezza–Cranio (0.60). La Gestazione mostra correlazioni moderate con Peso (0.59), Lunghezza (0.62) e Cranio (0.46), confermando che una maggiore durata della gravidanza si associa a una crescita fetale più avanzata. Le variabili materne (Anni.madre e N.gravidanze) presentano invece correlazioni molto deboli o nulle con le misure neonatali, suggerendo un impatto trascurabile. L’unica correlazione apprezzabile tra variabili materne è Anni.madre–N.gravidanze (0.38), fisiologicamente attesa. Nel complesso, non emergono valori superiori a 0.80 tali da indicare problemi di multicollinearità per il modello di regressione.

Per visualizzare nel dettaglio la forma delle relazioni bivariate tra le variabili numeriche si produce una scatterplot matrix, con scatter plot e smoother loess nel triangolo superiore e coefficienti di correlazione nel triangolo inferiore.

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) {
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- cor(x, y)
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8 / strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor, col = "red")
}

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

pairs(vars_num,
      upper.panel = panel.smooth,
      lower.panel = panel.cor,
      main = "Scatterplot matrix")

La scatterplot matrix conferma visivamente quanto evidenziato dalla matrice di correlazione. Le coppie Peso–Lunghezza, Peso–Cranio e Lunghezza–Cranio mostrano nuvole di punti con trend ascendente ben definito e smoother lineari, coerenti con le correlazioni elevate. Per Gestazione–Peso e Gestazione–Lunghezza si osserva una dispersione crescente con possibile lieve curvatura, che sarà indagata nella fase di modellazione tramite effetti quadratici. Le variabili materne (Anni.madre, N.gravidanze) non mostrano trend apprezzabili con le misure neonatali, confermando il loro scarso potere predittivo sul peso alla nascita.

Queste osservazioni confermano la coerenza fisiologica del dataset. Nel complesso, non emergono problemi di multicollinearità tali da compromettere la costruzione del modello di regressione multipla.

Al fine di verificare la normalità della variabile Peso si effettua uno Shapiro-Wilk test.

shapiro.test(df$Peso)

## 
##  Shapiro-Wilk normality test
## 
## data:  df$Peso
## W = 0.97066, p-value < 2.2e-16

Finora abbiamo riscontrato una distribuzione del peso asimmetrica, la presenza di outlier e una coda sinistra più lunga (ovvero con valori molto bassi). E’ quindi corretto aspettarsi come esito del test un p-value < 0.05 e il risultato lo conferma. Si rifiuta quindi l’ipotesi nulla H_0 per cui il peso sarebbe normalmente distribuito: il peso non è normalmente distribuito. Questo, tuttavia, non è un problema per la regressione lineare, perché la regressione richiede la normalità dei residui, non della variabile dipendente e con n = 2500, il teorema del limite centrale rende il modello robusto.

Creazione del modello di regressione lineare multipla con le sole variabili che hanno senso predittivo: Gestazione, Lunghezza, Cranio, Anni.madre, N.gravidanze, Fumatrici, Sesso.

mod_full <- lm(Peso ~ Gestazione + Lunghezza + Cranio +
                 Anni.madre + N.gravidanze + Fumatrici + Sesso,
               data = df)

summary(mod_full)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Anni.madre + 
##     N.gravidanze + Fumatrici + Sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1160.62  -181.17   -15.91   163.47  2631.35 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6711.5440   141.2543 -47.514  < 2e-16 ***
## Gestazione      32.8936     3.8259   8.598  < 2e-16 ***
## Lunghezza       10.2348     0.3009  34.010  < 2e-16 ***
## Cranio          10.5192     0.4268  24.644  < 2e-16 ***
## Anni.madre       0.8772     1.1487   0.764   0.4452    
## N.gravidanze    11.4029     4.6745   2.439   0.0148 *  
## Fumatrici      -30.2865    27.5981  -1.097   0.2726    
## SessoM          78.0898    11.2042   6.970 4.05e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2492 degrees of freedom
## Multiple R-squared:  0.7272, Adjusted R-squared:  0.7264 
## F-statistic:   949 on 7 and 2492 DF,  p-value: < 2.2e-16

tab_mod <- tidy(mod_full)

kable(tab_mod,
      digits = 4,
      caption = "Coefficienti del modello completo (senza Tipo.parto e Ospedale)") |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE)

Coefficienti del modello completo (senza Tipo.parto e Ospedale)
term	estimate	std.error	statistic	p.value
(Intercept)	-6711.5440	141.2543	-47.5139	0.0000
Gestazione	32.8936	3.8259	8.5976	0.0000
Lunghezza	10.2348	0.3009	34.0103	0.0000
Cranio	10.5192	0.4268	24.6440	0.0000
Anni.madre	0.8772	1.1487	0.7636	0.4452
N.gravidanze	11.4029	4.6745	2.4394	0.0148
Fumatrici	-30.2865	27.5981	-1.0974	0.2726
SessoM	78.0898	11.2042	6.9697	0.0000

Una volta creato il modello completo calcoliamo il VIF per la verifica della multicollinearità:

# Calcolo del VIF
vif_values <- vif(mod_full)

# Conversione in data frame
vif_df <- as.data.frame(vif_values)

# Tabella kable
kable(vif_df,
      digits = 4,
      caption = "Valori del VIF per la verifica della multicollinearità") |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE)

Valori del VIF per la verifica della multicollinearità
	vif_values
Gestazione	1.6937
Lunghezza	2.0787
Cranio	1.6289
Anni.madre	1.1892
N.gravidanze	1.1874
Fumatrici	1.0067
Sesso	1.0404

Da questi risultati emerge che tutti i VIF sono ben al di sotto della soglia critica di 5-10: non ci sono problemi di multicollinearità.

2.2 Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC), selezioneremo il modello più parsimonioso, eliminando le variabili non significative. Verranno considerati anche modelli con interazioni tra le variabili e possibili effetti non lineari.

Selezione del modello tramite AIC (stepwise), usando il metodo stepAIC, che prova aggiunte e rimozioni di variabili per minimizzare l’AIC.

mod_step_AIC <- stepAIC(mod_full,
                        direction = "both",
                        trace = FALSE)
summary(mod_step_AIC)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze + 
##     Sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.44  -180.81   -15.58   163.64  2639.72 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.1445   135.7229 -49.226  < 2e-16 ***
## Gestazione      32.3321     3.7980   8.513  < 2e-16 ***
## Lunghezza       10.2486     0.3006  34.090  < 2e-16 ***
## Cranio          10.5402     0.4262  24.728  < 2e-16 ***
## N.gravidanze    12.4750     4.3396   2.875  0.00408 ** 
## SessoM          77.9927    11.2021   6.962 4.26e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1328 on 5 and 2494 DF,  p-value: < 2.2e-16

Secondo AIC, queste variabili aggiungono un piccolo miglioramento predittivo, sufficiente a giustificare la loro presenza.

Selezione del modello tramite BIC. Questo modello penalizza di più i modelli complessi, quindi tende a selezionare modelli più parsimoniosi.

mod_step_BIC <- stepAIC(mod_full,
                        direction = "both",
                        k = log(nrow(df)),
                        trace = FALSE)
summary(mod_step_BIC)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze + 
##     Sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.44  -180.81   -15.58   163.64  2639.72 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.1445   135.7229 -49.226  < 2e-16 ***
## Gestazione      32.3321     3.7980   8.513  < 2e-16 ***
## Lunghezza       10.2486     0.3006  34.090  < 2e-16 ***
## Cranio          10.5402     0.4262  24.728  < 2e-16 ***
## N.gravidanze    12.4750     4.3396   2.875  0.00408 ** 
## SessoM          77.9927    11.2021   6.962 4.26e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1328 on 5 and 2494 DF,  p-value: < 2.2e-16

Secondo BIC, queste variabili non migliorano abbastanza il modello da giustificare la loro inclusione.

#Variabili selezionate da AIC e BIC
vars_AIC <- names(coef(mod_step_AIC))[-1]  # escludi intercetta
vars_BIC <- names(coef(mod_step_BIC))[-1]

sel_df <- data.frame(
  Variabile = unique(c(vars_AIC, vars_BIC))
) %>%
  mutate(
    AIC = ifelse(Variabile %in% vars_AIC, "✔", ""),
    BIC = ifelse(Variabile %in% vars_BIC, "✔", "")
  )

kable(sel_df,
      caption = "Variabili selezionate dai criteri AIC e BIC",
      col.names = c("Variabile", "Selezionata da AIC", "Selezionata da BIC")) |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE)

Variabili selezionate dai criteri AIC e BIC
Variabile	Selezionata da AIC	Selezionata da BIC
Gestazione	✔	✔
Lunghezza	✔	✔
Cranio	✔	✔
N.gravidanze	✔	✔
SessoM	✔	✔

Entrambi i modelli confermano che le misure antropometriche (Lunghezza, Cranio) e la durata della gestazione sono i principali predittori del peso neonatale, con un contributo significativo anche del sesso del neonato.

2.3 Una volta ottenuto il modello finale, valuteremo la sua capacità predittiva utilizzando metriche come R² e il Root Mean Squared Error (RMSE). Un’attenzione particolare sarà rivolta all’analisi dei residui e alla presenza di valori influenti, che potrebbero distorcere le previsioni, indagando su di essi.

Prima di saggiare effetti quadratici o di interazione, si esamina graficamente la relazione tra ciascun predittore e il peso, con smoother loess, per individuare eventuali non-linearità. Si saggiano solo gli effetti che si ipotizzano dai grafici, non tutti quanti.

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

par(mfrow = c(2, 3))
for (v in vars_plot) {
  plot(df[[v]], df$Peso,
       xlab = v, ylab = "Peso (g)",
       main = paste("Peso vs", v),
       pch = 16, col = rgb(0, 0, 1, 0.2), cex = 0.6)
  lines(lowess(df[[v]], df$Peso), col = "red", lwd = 2)
}
par(mfrow = c(1, 1))

Osservando i grafici si nota una curvatura nel rapporto tra Gestazione e Peso e tra Lunghezza e Peso. Per Cranio, N.gravidanze e Anni.madre la relazione appare sostanzialmente lineare (o assente). Si saggia pertanto l’effetto quadratico solo per Gestazione e Lunghezza.

Secondo il principio di marginalità, quando si include un termine quadratico I(x^2), si mantiene sempre anche il corrispondente effetto principale x. Lo stesso vale per le interazioni: se si include A:B, si mantengono sia A sia B.

# Modello con effetti quadratici (solo dove suggerito dai grafici)
# Si parte dalle variabili del modello AIC/BIC e si aggiungono i termini quadratici

mod_quad <- lm(Peso ~ Gestazione + I(Gestazione^2) +
                 Lunghezza + I(Lunghezza^2) +
                 Cranio + N.gravidanze + Fumatrici + Sesso,
               data = df)

summary(mod_quad)

## 
## Call:
## lm(formula = Peso ~ Gestazione + I(Gestazione^2) + Lunghezza + 
##     I(Lunghezza^2) + Cranio + N.gravidanze + Fumatrici + Sesso, 
##     data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1191.75  -179.83   -13.33   165.58  1402.42 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -2.379e+03  9.056e+02  -2.627 0.008663 ** 
## Gestazione       3.370e+02  6.270e+01   5.374 8.40e-08 ***
## I(Gestazione^2) -3.877e+00  8.245e-01  -4.703 2.71e-06 ***
## Lunghezza       -3.210e+01  4.038e+00  -7.950 2.80e-15 ***
## I(Lunghezza^2)   4.364e-02  4.141e-03  10.539  < 2e-16 ***
## Cranio           1.044e+01  4.192e-01  24.910  < 2e-16 ***
## N.gravidanze     1.468e+01  4.252e+00   3.451 0.000568 ***
## Fumatrici       -2.456e+01  2.699e+01  -0.910 0.362981    
## SessoM           7.277e+01  1.100e+01   6.616 4.49e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 268.5 on 2491 degrees of freedom
## Multiple R-squared:  0.7393, Adjusted R-squared:  0.7385 
## F-statistic: 883.1 on 8 and 2491 DF,  p-value: < 2.2e-16

Si verificano anche possibili interazioni suggerite dal contesto clinico:

Gestazione e Fumatrici: il fumo potrebbe avere un effetto diverso a diverse età gestazionali.
Gestazione e Sesso: la crescita potrebbe differire tra i sessi nelle ultime settimane.

mod_quad_int <- lm(Peso ~ Gestazione + I(Gestazione^2) +
                     Lunghezza + I(Lunghezza^2) +
                     Cranio + N.gravidanze + Fumatrici + Sesso +
                     Gestazione:Fumatrici + Gestazione:Sesso,
                   data = df)

summary(mod_quad_int)

## 
## Call:
## lm(formula = Peso ~ Gestazione + I(Gestazione^2) + Lunghezza + 
##     I(Lunghezza^2) + Cranio + N.gravidanze + Fumatrici + Sesso + 
##     Gestazione:Fumatrici + Gestazione:Sesso, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1188.10  -181.69   -13.02   164.20  1403.40 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -2.398e+03  9.130e+02  -2.627 0.008677 ** 
## Gestazione            3.387e+02  6.362e+01   5.324 1.10e-07 ***
## I(Gestazione^2)      -3.922e+00  8.431e-01  -4.653 3.45e-06 ***
## Lunghezza            -3.200e+01  4.041e+00  -7.919 3.57e-15 ***
## I(Lunghezza^2)        4.354e-02  4.145e-03  10.505  < 2e-16 ***
## Cranio                1.044e+01  4.193e-01  24.887  < 2e-16 ***
## N.gravidanze          1.468e+01  4.253e+00   3.451 0.000568 ***
## Fumatrici             6.784e+02  7.449e+02   0.911 0.362481    
## SessoM               -1.176e+02  2.353e+02  -0.500 0.617286    
## Gestazione:Fumatrici -1.792e+01  1.896e+01  -0.945 0.344549    
## Gestazione:SessoM     4.894e+00  6.029e+00   0.812 0.417031    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 268.5 on 2489 degrees of freedom
## Multiple R-squared:  0.7395, Adjusted R-squared:  0.7384 
## F-statistic: 706.5 on 10 and 2489 DF,  p-value: < 2.2e-16

Di seguito si conforntano i modelli ottenuti:

comp_df <- data.frame(
  Modello = c("Completo lineare", "AIC", "BIC", "Quadratico", "Quadratico + Interazioni"),
  AIC     = c(AIC(mod_full), AIC(mod_step_AIC), AIC(mod_step_BIC), AIC(mod_quad), AIC(mod_quad_int)),
  BIC     = c(BIC(mod_full), BIC(mod_step_AIC), BIC(mod_step_BIC), BIC(mod_quad), BIC(mod_quad_int)),
  R2      = c(glance(mod_full)$r.squared, glance(mod_step_AIC)$r.squared,
              glance(mod_step_BIC)$r.squared, glance(mod_quad)$r.squared,
              glance(mod_quad_int)$r.squared),
  R2_adj  = c(glance(mod_full)$adj.r.squared, glance(mod_step_AIC)$adj.r.squared,
              glance(mod_step_BIC)$adj.r.squared, glance(mod_quad)$adj.r.squared,
              glance(mod_quad_int)$adj.r.squared),
  RMSE    = c(glance(mod_full)$sigma, glance(mod_step_AIC)$sigma,
              glance(mod_step_BIC)$sigma, glance(mod_quad)$sigma,
              glance(mod_quad_int)$sigma)
)

kable(comp_df,
      digits = 3,
      caption = "Confronto tra modelli (AIC, BIC, R², R² adj, RMSE)") |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE)

Confronto tra modelli (AIC, BIC, R², R² adj, RMSE)
Modello	AIC	BIC	R2	R2_adj	RMSE
Completo lineare	35181.52	35233.94	0.727	0.726	274.615
AIC	35179.33	35220.10	0.727	0.726	274.604
BIC	35179.33	35220.10	0.727	0.726	274.604
Quadratico	35069.90	35128.14	0.739	0.738	268.498
Quadratico + Interazioni	35072.32	35142.21	0.739	0.738	268.521

Da questa tabella emerge il modello migliore sulla base del trade-off tra AIC, BIC, R² aggiustato e RMSE. Se le interazioni non risultano significative, si seleziona il modello quadratico senza interazioni.

# Selezione del modello finale
# (aggiornare se mod_quad_int risulta significativamente migliore)
mod_finale <- mod_quad

Lo step successivo del processo di modellazione prevede, una volta individuato il modello finale ovvero quello più parsimonioso, un commento rigoroso di tutte le sue stime.

tab_finale <- tidy(mod_finale)

kable(tab_finale,
      digits = 4,
      caption = "Coefficienti del modello finale") |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE)

Coefficienti del modello finale
term	estimate	std.error	statistic	p.value
(Intercept)	-2379.0609	905.5721	-2.6271	0.0087
Gestazione	336.9969	62.7042	5.3744	0.0000
I(Gestazione^2)	-3.8773	0.8245	-4.7026	0.0000
Lunghezza	-32.0986	4.0375	-7.9501	0.0000
I(Lunghezza^2)	0.0436	0.0041	10.5394	0.0000
Cranio	10.4425	0.4192	24.9105	0.0000
N.gravidanze	14.6751	4.2525	3.4509	0.0006
Fumatrici	-24.5584	26.9912	-0.9099	0.3630
SessoM	72.7731	10.9991	6.6163	0.0000

cat("R² =", round(glance(mod_finale)$r.squared, 4), "\n")

## R² = 0.7393

cat("R² adj =", round(glance(mod_finale)$adj.r.squared, 4), "\n")

## R² adj = 0.7385

cat("RMSE =", round(glance(mod_finale)$sigma, 2), "grammi\n")

## RMSE = 268.5 grammi

Vengono quindi interpretati i coefficienti con particolare attenzione ai punti seguenti:

Gestazione + I(Gestazione²): Effetto non lineare — il peso aumenta con le settimane di gestazione, ma con rendimenti decrescenti (o crescenti) a seconda del segno del termine quadratico. L’effetto principale è mantenuto in accordo con il principio di marginalità.
Lunghezza + I(Lunghezza²): Effetto non lineare — la lunghezza ha un impatto positivo sul peso, la curvatura cattura il fatto che a lunghezze molto elevate l’incremento marginale di peso si modifica. Anche qui l’effetto principale è mantenuto.
Cranio: Effetto lineare positivo — un cranio più grande è associato a un peso maggiore.
N.gravidanze: Effetto della parità — madri con più gravidanze tendono ad avere neonati leggermente più pesanti.
Fumatrici: Effetto negativo — il fumo materno riduce il peso neonatale.
Sesso (M): I maschi tendono a pesare di più rispetto alle femmine.

Il modello spiega circa il 74% della variabilità del peso neonatale e l’errore medio di previsione è di circa 268 grammi. I termini quadratici per gestazione e lunghezza sono risultati altamente significativi, indicando relazioni non strettamente lineari tra queste variabili e il peso alla nascita.

A questo punto dell’analisi si effettua una diagnostica dei residui e, poichè si è ottenuto il modello finale, si valuta la sua capacità predittiva e il rispetto delle assunzioni.

Si analizzano i residui per verificare che: la varianza sia costante (omoschedasticità), i residui siano distribuiti normalmente, non ci siano pattern sistematici, non ci siano osservazioni influenti che distorcono il modello.

par(mfrow = c(2, 2))
plot(mod_finale)

par(mfrow = c(1, 1))

I grafici ottenuti si interpretano come segue:

Residuals vs Fitted: Verifica linearità e omoschedasticità. La dispersione dovrebbe essere casuale e uniforme attorno allo 0.
Normal Q-Q: Verifica normalità dei residui. I punti dovrebbero seguire la retta diagonale.
Scale-Location: Verifica omoschedasticità. La linea rossa dovrebbe essere approssimativamente orizzontale.
Residuals vs Leverage: Identifica osservazioni influenti (Cook’s distance).

residui <- residuals(mod_finale)

par(mfrow = c(1, 2))
hist(residui, breaks = 40, col = "#1f78b4",
     main = "Distribuzione dei residui", xlab = "Residui")
abline(v = 0, col = "red", lwd = 2)

qqnorm(residui, main = "QQ-Plot dei residui")
qqline(residui, col = "red", lwd = 2)

par(mfrow = c(1, 1))

Si eseguono i test formali di Breusch-Pagan, Durbin-Watson e Shapiro-Wilk sui residui e si riassumono i risultati ottenuti in una tabella.

# Test di Breusch-Pagan per l'omoschedasticità
bp_test <- bptest(mod_finale)
bp_test

## 
##  studentized Breusch-Pagan test
## 
## data:  mod_finale
## BP = 98.19, df = 8, p-value < 2.2e-16

# Test di Durbin-Watson per l'autocorrelazione
dw_test <- dwtest(mod_finale)
dw_test

## 
##  Durbin-Watson test
## 
## data:  mod_finale
## DW = 1.9499, p-value = 0.1052
## alternative hypothesis: true autocorrelation is greater than 0

# Test di Shapiro-Wilk sui residui
sw_test <- shapiro.test(residui)
sw_test

## 
##  Shapiro-Wilk normality test
## 
## data:  residui
## W = 0.98911, p-value = 7.151e-13

diag_summary <- data.frame(
  Test = c("Breusch-Pagan (omoschedasticità)",
           "Durbin-Watson (autocorrelazione)",
           "Shapiro-Wilk (normalità residui)"),
  Statistica = c(bp_test$statistic, dw_test$statistic, sw_test$statistic),
  P_value = c(bp_test$p.value, dw_test$p.value, sw_test$p.value),
  Decisione = c(
    ifelse(bp_test$p.value < 0.05, "Eteroschedasticità", "Omoschedasticità OK"),
    ifelse(dw_test$p.value < 0.05, "Autocorrelazione presente", "No autocorrelazione"),
    ifelse(sw_test$p.value < 0.05, "Non normalità", "Normalità OK")
  )
)

kable(diag_summary,
      digits = 4,
      caption = "Riepilogo test diagnostici sui residui") |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE)

Riepilogo test diagnostici sui residui
	Test	Statistica	P_value	Decisione
BP	Breusch-Pagan (omoschedasticità)	98.1900	0.0000	Eteroschedasticità
DW	Durbin-Watson (autocorrelazione)	1.9499	0.1052	No autocorrelazione
W	Shapiro-Wilk (normalità residui)	0.9891	0.0000	Non normalità

Essendo i test formali molto sensibili, possono rifiutare H_0 anche per deviazioni minime dalla distribuzione teorica. È importante quindi valutare i risultati congiuntamente ai grafici diagnostici. Con n = 2500 il teorema del limite centrale garantisce comunque la robustezza dell’inferenza sui coefficienti.

Si passa allo studio degli outlier e alle osservazioni influenti:

Cook’s Distance:

cooks_d <- cooks.distance(mod_finale)
soglia_cook <- 4 / nrow(df)

plot(cooks_d, type = "h", main = "Cook's Distance",
     ylab = "Cook's Distance", xlab = "Osservazione",
     col = ifelse(cooks_d > soglia_cook, "red", "grey"))
abline(h = soglia_cook, col = "red", lty = 2)
legend("topright", legend = paste("Soglia =", round(soglia_cook, 5)),
       col = "red", lty = 2)

influenti_cook <- which(cooks_d > soglia_cook)
cat("Numero di osservazioni con Cook's D > soglia:", length(influenti_cook), "\n")

## Numero di osservazioni con Cook's D > soglia: 131

DFFITS:

dffits_val <- dffits(mod_finale)
p <- length(coef(mod_finale))
n <- nrow(df)
soglia_dffits <- 2 * sqrt(p / n)

cat("Soglia DFFITS:", round(soglia_dffits, 4), "\n")

## Soglia DFFITS: 0.12

influenti_dffits <- which(abs(dffits_val) > soglia_dffits)
cat("Osservazioni con |DFFITS| > soglia:", length(influenti_dffits), "\n")

## Osservazioni con |DFFITS| > soglia: 131

Leverage (hat values):

hat_val <- hatvalues(mod_finale)
soglia_hat <- 2 * p / n

plot(hat_val, type = "h", main = "Leverage (hat values)",
     ylab = "Leverage", xlab = "Osservazione",
     col = ifelse(hat_val > soglia_hat, "red", "grey"))
abline(h = soglia_hat, col = "red", lty = 2)

influenti_hat <- which(hat_val > soglia_hat)
cat("Osservazioni con leverage > soglia:", length(influenti_hat), "\n")

## Osservazioni con leverage > soglia: 213

Tabella riepilogativa osservazioni influenti:

# Osservazioni flaggate da almeno 2 criteri su 3
flag_cook   <- cooks_d > soglia_cook
flag_dffits <- abs(dffits_val) > soglia_dffits
flag_hat    <- hat_val > soglia_hat

flag_sum <- flag_cook + flag_dffits + flag_hat
idx_multi <- which(flag_sum >= 2)

if (length(idx_multi) > 0) {
  infl_df <- data.frame(
    Osservazione = idx_multi,
    Cook_D       = round(cooks_d[idx_multi], 5),
    DFFITS       = round(dffits_val[idx_multi], 4),
    Leverage     = round(hat_val[idx_multi], 4),
    Peso         = df$Peso[idx_multi],
    Gestazione   = df$Gestazione[idx_multi],
    Lunghezza    = df$Lunghezza[idx_multi]
  )

  kable(infl_df,
        caption = "Osservazioni influenti (flaggate da almeno 2 criteri su 3)") |>
    kable_styling(full_width = FALSE,
                  bootstrap_options = c("striped", "hover")) |>
    column_spec(1, bold = TRUE)
} else {
  cat("Nessuna osservazione flaggata da almeno 2 criteri su 3.\n")
}

Osservazioni influenti (flaggate da almeno 2 criteri su 3)
	Osservazione	Cook_D	DFFITS	Leverage	Peso	Gestazione	Lunghezza
99	99	0.00304	-0.1655	0.0106	3350	40	520
101	101	0.00324	-0.1709	0.0173	1370	34	390
119	119	0.00952	-0.2935	0.0062	3410	40	550
130	130	0.00224	0.1422	0.0018	4240	39	485
134	134	0.00294	0.1627	0.0078	3950	37	500
140	140	0.00460	0.2035	0.0125	4420	41	530
146	146	0.00174	0.1252	0.0023	3820	40	500
155	155	0.02847	0.5081	0.0128	3610	36	410
161	161	0.00421	-0.1946	0.0233	3760	42	540
220	220	0.00255	0.1516	0.0077	3520	40	445
295	295	0.00323	-0.1708	0.0042	1850	40	460
304	304	0.00163	-0.1212	0.0088	2060	37	420
310	310	0.07921	-0.8453	0.0926	1560	28	420
312	312	0.00325	-0.1709	0.0386	1280	32	360
322	322	0.00227	-0.1431	0.0059	1750	37	430
329	329	0.00246	0.1489	0.0040	4560	40	540
335	335	0.00201	0.1344	0.0112	3150	38	465
375	375	0.00237	0.1461	0.0037	4270	38	510
378	378	0.00706	0.2522	0.0586	1285	28	400
383	383	0.00234	-0.1452	0.0077	3600	39	550
390	390	0.00309	0.1671	0.0029	3700	40	470
391	391	0.00252	-0.1507	0.0112	3400	40	525
413	413	0.00228	-0.1433	0.0106	2590	40	480
424	424	0.00200	0.1342	0.0113	3820	41	500
442	442	0.00260	-0.1530	0.0165	2430	38	460
471	471	0.00599	-0.2323	0.0101	3680	40	560
472	472	0.00169	0.1235	0.0025	3990	41	495
473	473	0.00271	-0.1563	0.0118	3090	41	495
478	478	0.00251	0.1504	0.0058	4310	42	505
492	492	0.00346	-0.1764	0.0205	1410	33	380
516	516	0.00320	0.1698	0.0135	3520	38	470
567	567	0.00382	0.1854	0.0106	3280	38	475
582	582	0.00440	-0.1991	0.0123	2220	35	470
615	615	0.00189	-0.1306	0.0047	2100	35	440
616	616	0.00366	-0.1817	0.0047	3540	42	540
656	656	0.00225	-0.1423	0.0062	2320	41	450
657	657	0.00217	-0.1399	0.0056	2970	41	520
684	684	0.00427	-0.1962	0.0089	3000	39	475
729	729	0.00264	-0.1543	0.0096	3900	39	555
750	750	0.00182	-0.1279	0.0129	1450	35	405
791	791	0.00272	0.1569	0.0019	4440	41	510
828	828	0.00292	0.1621	0.0073	4200	40	510
890	890	0.00173	0.1248	0.0023	3660	39	490
906	906	0.00161	0.1203	0.0046	3790	42	500
991	991	0.00308	-0.1666	0.0074	4050	40	550
1008	1008	0.00165	-0.1220	0.0057	2400	36	470
1036	1036	0.00196	0.1329	0.0017	4330	40	500
1109	1109	0.00180	-0.1275	0.0044	3740	41	545
1181	1181	0.00235	0.1456	0.0061	4070	36	500
1194	1194	0.00192	-0.1316	0.0037	2160	39	450
1215	1215	0.00330	-0.1726	0.0057	3630	40	550
1230	1230	0.00378	0.1846	0.0042	4010	41	470
1253	1253	0.00171	-0.1239	0.0042	3450	39	523
1268	1268	0.00363	0.1811	0.0026	3790	40	460
1273	1273	0.00230	-0.1440	0.0121	2040	33	480
1293	1293	0.00958	0.2944	0.0063	4600	38	485
1306	1306	0.00593	0.2322	0.0022	4900	41	510
1323	1323	0.00505	-0.2134	0.0085	3660	39	555
1341	1341	0.00249	0.1498	0.0039	3780	36	480
1357	1357	0.00286	0.1606	0.0130	2340	32	445
1385	1385	0.00465	0.2046	0.0405	1620	29	410
1395	1395	0.00440	0.1994	0.0054	3790	39	505
1398	1398	0.00192	0.1314	0.0109	3180	40	470
1399	1399	0.00608	-0.2348	0.0028	2560	38	525
1402	1402	0.00186	-0.1293	0.0049	2660	39	460
1426	1426	0.00428	-0.1964	0.0133	2250	39	460
1428	1428	0.01875	-0.4112	0.0295	1280	36	385
1429	1429	0.04106	-0.6088	0.0446	1280	29	390
1433	1433	0.00373	0.1834	0.0047	4810	41	530
1450	1450	0.00553	0.2232	0.0155	3730	41	480
1472	1472	0.00190	0.1307	0.0044	3720	37	480
1505	1505	0.00330	-0.1724	0.0135	2860	39	490
1525	1525	0.00319	-0.1696	0.0105	2920	39	490
1541	1541	0.00360	0.1802	0.0035	4540	38	530
1551	1551	1.32448	3.4773	0.2474	4370	38	315
1553	1553	0.01181	0.3267	0.0088	4520	35	520
1556	1556	0.00200	-0.1344	0.0062	2420	41	490
1583	1583	0.00287	0.1607	0.0127	4220	40	510
1588	1588	0.00260	-0.1532	0.0049	3100	42	510
1593	1593	0.00537	-0.2201	0.0075	1500	35	420
1619	1619	0.01674	-0.3882	0.0638	990	31	340
1626	1626	0.00198	-0.1335	0.0112	3030	38	500
1635	1635	0.00424	0.1956	0.0041	3430	39	445
1639	1639	0.00258	0.1525	0.0078	4760	40	550
1644	1644	0.00175	-0.1257	0.0052	3850	41	550
1691	1691	0.00241	-0.1472	0.0113	3220	41	500
1694	1694	0.00450	0.2018	0.0022	3850	36	460
1701	1701	0.00317	-0.1690	0.0188	1430	32	380
1712	1712	0.00545	0.2217	0.0073	3800	39	520
1718	1718	0.00391	-0.1877	0.0086	2660	42	500
1743	1743	0.00202	-0.1349	0.0050	1950	38	445
1780	1780	0.04858	0.6615	0.1300	900	25	325
1838	1838	0.00206	0.1363	0.0023	4580	40	515
1856	1856	0.00325	0.1712	0.0047	3110	36	465
1868	1868	0.00480	-0.2080	0.0055	3470	40	525
1893	1893	0.00347	0.1768	0.0060	3030	34	470
1915	1915	0.00163	-0.1212	0.0031	2480	38	500
1920	1920	0.02230	0.4489	0.0178	4930	39	550
1937	1937	0.00161	-0.1206	0.0023	2940	41	500
1962	1962	0.00459	0.2036	0.0063	4370	38	530
1963	1963	0.00293	0.1625	0.0054	4700	42	540
1971	1971	0.00262	0.1535	0.0124	3530	38	495
2003	2003	0.00265	-0.1544	0.0112	2740	40	500
2016	2016	0.00183	-0.1282	0.0162	3870	42	530
2023	2023	0.00189	0.1307	0.0012	4650	39	510
2040	2040	0.00676	0.2467	0.0189	3240	38	410
2046	2046	0.00180	0.1273	0.0150	4300	40	530
2049	2049	0.00271	0.1563	0.0106	4180	40	520
2076	2076	0.00290	0.1618	0.0042	4720	40	540
2086	2086	0.00347	-0.1768	0.0133	3250	40	500
2087	2087	0.00188	0.1301	0.0083	4140	43	520
2110	2110	0.00243	-0.1478	0.0110	2910	41	490
2115	2115	0.01506	-0.3686	0.0223	1890	32	500
2120	2120	0.01611	0.3809	0.0697	1140	27	370
2132	2132	0.00178	0.1266	0.0034	3710	36	480
2135	2135	0.00189	0.1306	0.0024	3950	39	480
2175	2175	0.00613	0.2349	0.0625	930	28	355
2195	2195	0.00193	0.1321	0.0015	3980	40	480
2219	2219	0.00169	-0.1237	0.0014	2500	39	490
2221	2221	0.00471	-0.2059	0.0218	2950	39	495
2225	2225	0.00842	0.2760	0.0062	3140	35	465
2245	2245	0.00205	0.1358	0.0140	4050	37	530
2287	2287	0.00186	0.1297	0.0031	3640	40	470
2315	2315	0.00437	-0.1988	0.0037	2800	42	520
2317	2317	0.00209	-0.1372	0.0081	2530	38	460
2337	2337	0.00230	0.1437	0.0145	3440	37	460
2392	2392	0.00206	0.1364	0.0033	4720	40	540
2412	2412	0.00433	0.1975	0.0105	4030	39	510
2422	2422	0.00590	-0.2305	0.0217	3090	40	485
2452	2452	0.04212	0.6160	0.0938	930	26	345
2471	2471	0.00362	-0.1805	0.0214	2880	38	470

Considerazioni: Con 2500 osservazioni è normale avere alcuni punti con valori superiori alle soglie. Le osservazioni influenti non vengono rimosse automaticamente: si valuta caso per caso se sono errori di registrazione o valori clinicamente plausibili (es. neonati molto prematuri o macrosomici, già identificati come outlier nelle analisi preliminari). In assenza di errori evidenti, vengono mantenute nel modello. Nessuno di questi punti compromette la stabilità del modello, come confermato dai valori di R² e dalla diagnostica dei residui.

3. Previsioni e Risultati:

STEP 7 — Forecasting:

Una volta validato il modello, lo useremo per fare previsioni pratiche. Invece di una singola previsione puntuale, si costruisce una tabella con diversi scenari clinicamente rilevanti, variando settimane di gestazione e condizione di fumatrice/non fumatrice.

Previsione multi-scenario:

media_lunghezza <- mean(df$Lunghezza)
media_cranio    <- mean(df$Cranio)

scenari <- expand.grid(
  Gestazione   = c(36, 37, 38, 39, 40, 41),
  Fumatrici    = c(0, 1),
  Lunghezza    = media_lunghezza,
  Cranio       = media_cranio,
  N.gravidanze = 2,
  Sesso        = "F"
)

pred_scenari <- predict(mod_finale, newdata = scenari, interval = "prediction")
scenari_result <- cbind(scenari[, c("Gestazione", "Fumatrici")],
                        round(pred_scenari, 1))
scenari_result$Fumatrici <- factor(scenari_result$Fumatrici,
                                    labels = c("No", "Sì"))

kable(scenari_result,
      col.names = c("Gestazione (sett.)", "Fumatrice", "Peso stimato (g)",
                     "IC inf (g)", "IC sup (g)"),
      caption = "Previsione del peso neonatale per diversi scenari (femmina, 2ª gravidanza, misure medie)") |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover")) |>
  column_spec(1, bold = TRUE)

Previsione del peso neonatale per diversi scenari (femmina, 2ª gravidanza, misure medie)
Gestazione (sett.)	Fumatrice	Peso stimato (g)	IC inf (g)	IC sup (g)
36	No	3109.2	2581.9	3636.5
37	No	3163.2	2636.2	3690.2
38	No	3209.4	2682.5	3736.3
39	No	3247.8	2721.0	3774.6
40	No	3278.5	2751.7	3805.4
41	No	3301.5	2774.3	3828.6
36	Sì	3084.7	2554.9	3614.5
37	Sì	3138.6	2609.1	3668.1
38	Sì	3184.8	2655.5	3714.1
39	Sì	3223.3	2694.1	3752.5
40	Sì	3254.0	2724.7	3783.2
41	Sì	3276.9	2747.4	3806.4

Grafico: Peso previsto vs Gestazione con banda di confidenza:

# Griglia fine per il grafico
grid_pred <- expand.grid(
  Gestazione   = seq(min(df$Gestazione), max(df$Gestazione), by = 0.5),
  Fumatrici    = c(0, 1),
  Lunghezza    = media_lunghezza,
  Cranio       = media_cranio,
  N.gravidanze = 2,
  Sesso        = "F"
)

pred_grid <- predict(mod_finale, newdata = grid_pred, interval = "confidence")
grid_pred <- cbind(grid_pred, pred_grid)
grid_pred$Fumo <- factor(grid_pred$Fumatrici, labels = c("Non fumatrice", "Fumatrice"))

ggplot(grid_pred, aes(x = Gestazione, y = fit, color = Fumo, fill = Fumo)) +
  geom_line(linewidth = 1) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.2, color = NA) +
  labs(title = "Peso previsto vs Settimane di gestazione",
       subtitle = "Banda di confidenza al 95% — Femmina, 2ª gravidanza, misure medie",
       x = "Settimane di gestazione",
       y = "Peso previsto (g)",
       color = "", fill = "") +
  scale_color_manual(values = c("Non fumatrice" = "#66c2a5", "Fumatrice" = "#fc8d62")) +
  scale_fill_manual(values = c("Non fumatrice" = "#66c2a5", "Fumatrice" = "#fc8d62")) +
  theme_minimal(base_size = 13)

Esempio di previsione puntuale:

Si calcolano le medie di lunghezza e cranio e si crea il nuovo caso da prevedere effettuando una previsione con intervallo di previsione. Stimiamo il peso di una neonata considerando una madre non fumatrice alla terza gravidanza che partorirà alla 39esima settimana:

nuovo_neonato <- data.frame(
  Gestazione   = 39,
  Lunghezza    = media_lunghezza,
  Cranio       = media_cranio,
  N.gravidanze = 3,
  Fumatrici    = 0,
  Sesso        = "F"
)

pred_punto <- predict(mod_finale, newdata = nuovo_neonato, interval = "prediction")

pred_df <- data.frame(
  Stima    = round(pred_punto[1], 2),
  IC_lower = round(pred_punto[2], 2),
  IC_upper = round(pred_punto[3], 2)
)

kable(pred_df,
      col.names = c("Peso stimato (g)", "IC inf (g)", "IC sup (g)"),
      caption = "Previsione puntuale: neonata, 39 sett., 3ª gravidanza, non fumatrice, misure medie") |>
  kable_styling(full_width = FALSE,
                bootstrap_options = c("striped", "hover"))

Previsione puntuale: neonata, 39 sett., 3ª gravidanza, non fumatrice, misure medie
Peso stimato (g)	IC inf (g)	IC sup (g)
3262.51	2735.5	3789.53

Questa previsione rappresenta un’applicazione pratica del modello e dimostra la sua capacità di fornire stime clinicamente plausibili e coerenti con la fisiologia neonatale.

4. Visualizzazioni:

Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili. Ad esempio, potremmo visualizzare l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.

Si usa plotly per generare un grafico 3D, interattivo e ruotabile, che mostra l’effetto combinato di Gestazione × Fumo sul Peso previsto, utilizzando il modello quadratico.

Questo tipo di grafico è molto utile perché visualizza l’effetto principale della gestazione, mostra la differenza tra fumatrici e non fumatrici, evidenzia eventuali divergenze nelle pendenze (interazioni) e comunica in modo intuitivo l’impatto del fumo sul peso neonatale.

# Sequenza di settimane di gestazione
gest_seq <- seq(min(df$Gestazione), max(df$Gestazione), length.out = 50)

# Data frame per fumatrici e non fumatrici
df_3d <- expand.grid(
  Gestazione = gest_seq,
  Fumatrici = c(0, 1)
)

# Aggiunta delle altre variabili ai valori medi
df_3d$Lunghezza    <- mean(df$Lunghezza)
df_3d$Cranio       <- mean(df$Cranio)
df_3d$N.gravidanze <- mean(df$N.gravidanze)
df_3d$Sesso        <- "F"

# Predizione del modello
df_3d$Peso_pred <- predict(mod_finale, newdata = df_3d)

# Conversione del fumo in etichette
df_3d$Fumo_label <- factor(df_3d$Fumatrici, labels = c("Non fumatrici", "Fumatrici"))

# Grafico 3D
plot_ly(
  df_3d,
  x = ~Gestazione,
  y = ~Fumo_label,
  z = ~Peso_pred,
  type = "scatter3d",
  mode = "markers",
  color = ~Fumo_label,
  colors = c("blue", "red"),
  marker = list(size = 4)
) %>%
  layout(
    title = "Effetto combinato di Gestazione e Fumo sul Peso previsto",
    scene = list(
      xaxis = list(title = "Settimane di gestazione"),
      yaxis = list(title = "Fumo materno"),
      zaxis = list(title = "Peso previsto (g)")
    )
  )

Si può interpretare il grafico dicendo che la curva ascendente lungo l’asse della gestazione mostra l’aumento del peso previsto con l’avanzare delle settimane.
Le due “nuvole” di punti (fumatrici vs non fumatrici) sono separate verticalmente, quindi le fumatrici hanno un peso previsto inferiore a parità di gestazione.
La distanza tra le due superfici indica l’effetto del fumo: se il modello avesse un’interazione forte, le due superfici avrebbero pendenze diverse; nel tuo caso l’effetto è moderato, quindi le superfici sono quasi parallele.

In conclusione, per visualizzare l’effetto combinato delle settimane di gestazione e del fumo materno sul peso previsto, è stato prodotto un grafico 3D interattivo. Il grafico mostra due superfici distinte: una per le madri fumatrici e una per le non fumatrici. In entrambi i gruppi il peso previsto aumenta con la gestazione, ma la superficie relativa alle fumatrici risulta sistematicamente più bassa, indicando un effetto negativo del fumo sul peso neonatale. Questa rappresentazione consente di cogliere in modo immediato sia l’effetto principale della gestazione sia la differenza tra i due gruppi.

Si genera un secondo grafico 3D, sempre usando plotly, per visualizzare la relazione tra: Gestazione, Lunghezza e Peso previsto. Tale grafico permette di visualizzare i punti simulati dal modello, la superficie di regressione stimata e l’effetto combinato di due variabili chiave sul peso previsto.

# Griglia di valori per le due variabili
gest_seq <- seq(min(df$Gestazione), max(df$Gestazione), length.out = 40)
lung_seq <- seq(min(df$Lunghezza), max(df$Lunghezza), length.out = 40)

grid <- expand.grid(
  Gestazione = gest_seq,
  Lunghezza = lung_seq,
  Cranio = mean(df$Cranio),
  N.gravidanze = mean(df$N.gravidanze),
  Fumatrici = 0,
  Sesso = "F"
)

# Predizioni sulla griglia
grid$Peso_pred <- predict(mod_finale, newdata = grid)

# Grafico 3D
plot_ly(
  x = grid$Gestazione,
  y = grid$Lunghezza,
  z = grid$Peso_pred,
  type = "scatter3d",
  mode = "markers",
  marker = list(size = 3, color = grid$Peso_pred, colorscale = "Viridis")
) %>%
  layout(
    title = "Scatter plot 3D del peso previsto",
    scene = list(
      xaxis = list(title = "Settimane di gestazione"),
      yaxis = list(title = "Lunghezza (mm)"),
      zaxis = list(title = "Peso previsto (g)")
    )
  )

La sua interpretazione consente di notare, osservando la superficie, come gestazione e lunghezza influenzano il peso previsto. La forma leggermente curva riflette i termini quadratici del modello, mentre le sue zone più alte corrispondono a neonati più maturi e più lunghi.

In conclusione, si può dire che per visualizzare l’effetto combinato delle variabili più rilevanti, è stato prodotto uno scatter plot 3D che mostra la relazione tra settimane di gestazione, lunghezza del neonato e peso previsto dal modello quadratico. La superficie risultante evidenzia un incremento non lineare del peso con l’aumentare della gestazione e della lunghezza, confermando la presenza di effetti quadratici e l’adeguatezza del modello nel catturare la fisiologia della crescita fetale.

Conclusioni:

Il modello quadratico finale spiega circa il 74% della variabilità del peso neonatale. Le variabili più rilevanti sono la durata della gestazione e le misure antropometriche (lunghezza e cranio), con effetti non lineari per gestazione e lunghezza. Il fumo materno ha un effetto negativo significativo. Le variabili Tipo.parto e Ospedale sono state escluse a priori in quanto prive di senso predittivo (il tipo di parto si decide al momento, l’ospedale non influenza il peso — confermato dal chi-quadro p = 0.578). Gli effetti quadratici sono stati saggiati solo dove i grafici loess suggerivano curvatura, e i relativi effetti principali sono stati mantenuti in accordo con il principio di marginalità. La diagnostica completa sui residui (Breusch-Pagan, Durbin-Watson, Shapiro-Wilk, QQ-plot) e lo studio delle osservazioni influenti (Cook’s distance, DFFITS, leverage) confermano la robustezza del modello. Le previsioni multi-scenario dimostrano l’applicabilità clinica del modello per supportare la pianificazione perinatale.