# IMPORTANTE: MASS va caricato prima di tidyverse per evitare
# che MASS::select() sovrascriva dplyr::select()
library(MASS)       # stepAIC() — caricato per primo
library(tidyverse)
library(knitr)
library(kableExtra)
library(scales)
library(e1071)      # skewness / kurtosis
library(car)        # vif(), leveneTest()
library(lmtest)     # bptest()
library(nortest)    # lillie.test()
library(ggcorrplot) # correlogram

1 Introduzione

Quetso progetto nasce con l’obiettivo preciso di costruire un modello statistico predittivo del peso neonatale alla nascita sulla base di variabili cliniche raccolte da tre ospedali italiani, per capire se, e in che misura, è possibile prevedere il peso di un neonato alla nascita, partendo dalle informaizoni cliniche disponibili durante la gravidanza. Il peso alla nascita rappresenta uno degli indicatori più usati in ambito prenatale per valutare lo stato di salute del bambino e pianificare eventuali interventi.

I dati utilizzati riguardano 2.500 neonati provenienti da 3 ospedali e 10 variabili che coprono sia le caratteristiche della madre (età, numero di gravidanze, fumo) sia quelle del neonato (settimane di gestazione, lunghezza, diametro del cranio, sesso, tipo di parto).

Le quattro fasi dell’analisi sono: 1. Esplorazione dei dati distribuzione delle variabili, outlier e correlazioni) 2. Test di ipotesi su proporzioni, medie e differenze tra gruppi 3. Modello di regressione lineare multipla con costrizione del modello ottimale 4. Previsioni pratiche e visualizzazioni dei risultati

2 Struttura del Dataset e Analisi Esplorativa

2.1 Caricamento e Prima Ispezione

df <- read.csv("neonati.csv", stringsAsFactors = TRUE)
glimpse(df)

#> Rows: 2,500
#> Columns: 10
#> $ Anni.madre   <int> 26, 21, 34, 28, 20, 32, 26, 25, 22, 23, 29, 21, 36, 24, 3…
#> $ N.gravidanze <int> 0, 2, 3, 1, 0, 0, 1, 0, 1, 0, 2, 2, 5, 0, 3, 2, 2, 3, 0, …
#> $ Fumatrici    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
#> $ Gestazione   <int> 42, 39, 38, 41, 38, 40, 39, 40, 40, 41, 38, 40, 38, 40, 3…
#> $ Peso         <int> 3380, 3150, 3640, 3690, 3700, 3200, 3100, 3580, 3670, 370…
#> $ Lunghezza    <int> 490, 490, 500, 515, 480, 495, 480, 510, 500, 510, 480, 51…
#> $ Cranio       <int> 325, 345, 375, 365, 335, 340, 345, 349, 335, 362, 330, 34…
#> $ Tipo.parto   <fct> Nat, Nat, Nat, Nat, Nat, Nat, Nat, Nat, Ces, Ces, Ces, Na…
#> $ Ospedale     <fct> osp3, osp1, osp2, osp2, osp3, osp2, osp3, osp1, osp2, osp…
#> $ Sesso        <fct> M, F, M, M, F, F, F, M, F, F, M, F, F, F, M, M, M, M, F, …

cat("Righe:", nrow(df), " | Colonne:", ncol(df), "\n")

#> Righe: 2500  | Colonne: 10

cat("Valori mancanti totali:", sum(is.na(df)), "\n")

#> Valori mancanti totali: 0

head(df, 8) %>%
  kable(caption = "Prime 8 osservazioni del dataset") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)

Prime 8 osservazioni del dataset
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
26	1	39	3100	480	345	Nat	osp3	F
25	0	40	3580	510	349	Nat	osp1	M

2.2 Classificazione delle Variabili

Prima di qualsiasi calcolo, è utile chiarire con cosa si sta lavorando. Le variabili del dataset non sono tutte dello stesso tipo, e questo determina quali strumenti statistici è corretto applicare.

Variabile	Tipo Statistico	Scala	Descrizione
`Anni.madre`	Quantitativa Continua	Rapporto	Età della madre in anni
`N.gravidanze`	Quantitativa Discreta	Rapporto	Numero di gravidanze precedenti
`Fumatrici`	Qualitativa Nominale (dummy)	Nominale	0 = non fumatrice, 1 = fumatrice
`Gestazione`	Quantitativa Discreta	Rapporto	Settimane di gestazione
`Peso`	Quantitativa Continua	Rapporto	Variabile target — peso neonatale in grammi
`Lunghezza`	Quantitativa Continua	Rapporto	Lunghezza neonatale in mm
`Cranio`	Quantitativa Continua	Rapporto	Diametro cranico in mm
`Tipo.parto`	Qualitativa Nominale	Nominale	Nat = naturale, Ces = cesareo
`Ospedale`	Qualitativa Nominale	Nominale	osp1, osp2, osp3
`Sesso`	Qualitativa Nominale	Nominale	M = maschio, F = femmina

Ci tengo a effettuae alcune utili precisazioni: Fumatrici, pur essendo codificata come 0/1, è concettualmente una variabile qualitativa dicotomica per cui il valore numerico non ha un significato di grandezza, ma solo di appartenenza a una categoria. Ospedale e Tipo.parto sono nominali senza alcun ordine naturale tra le modalità. Per queste variabili non ha senso calcolare media o varianza ma si useranno distribuzioni di frequenza e, dove necessario, variabili dummy nella regressione.

2.3 Gestione degli Outlier

cat("Osservazioni con Anni.madre = 0:", sum(df$Anni.madre == 0), "\n")

#> Osservazioni con Anni.madre = 0: 1

# Sostituzione del valore impossibile con NA e rimozione
df <- df %>% filter(Anni.madre > 0)
cat("Dataset dopo rimozione outlier:", nrow(df), "righe\n")

#> Dataset dopo rimozione outlier: 2499 righe

Il valore Anni.madre = 0 è chiaramente un errore di inserimento dati (poiché impossibile biologicamente) e quindi non è un caso limite da tenere in considerazione. Ho rimosso l’osservazione portando il dataset a 2.499 record. NB in un’analisi reale sarebbe sicuramente opportuno segnalare questo problema al team che ha raccolto i dati, per capire se si tratta di un errore sistematico o isolato.

2.4 Statistiche Descrittive

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

calcola_stats <- function(x, nome) {
  data.frame(
    Variabile  = nome,
    Media      = round(mean(x, na.rm=TRUE), 2),
    Mediana    = round(median(x, na.rm=TRUE), 2),
    Dev_Std    = round(sd(x, na.rm=TRUE), 2),
    CV_pct     = round(sd(x)/mean(x)*100, 1),
    Skewness   = round(skewness(x, na.rm=TRUE), 3),
    Curtosi    = round(kurtosis(x, na.rm=TRUE), 3),
    Min        = round(min(x, na.rm=TRUE), 2),
    Max        = round(max(x, na.rm=TRUE), 2),
    IQR        = round(IQR(x, na.rm=TRUE), 2)
  )
}

bind_rows(lapply(var_num, function(v) calcola_stats(df[[v]], v))) %>%
  kable(caption = "Indici descrittivi per le variabili quantitative") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE) %>%
  scroll_box(width = "100%")

Indici descrittivi per le variabili quantitative
Variabile	Media	Mediana	Dev_Std	CV_pct	Skewness	Curtosi	Min	Max	IQR
Anni.madre	28.18	28	5.24	18.6	0.099	0.121	1	46	7
N.gravidanze	0.98	1	1.28	130.5	2.512	10.976	0	12	1
Gestazione	38.98	39	1.87	4.8	-2.064	8.245	25	43	2
Peso	3284.17	3300	525.12	16.0	-0.647	2.027	830	4930	630
Lunghezza	494.69	500	26.32	5.3	-1.514	6.476	310	565	30
Cranio	340.03	340	16.43	4.8	-0.785	2.942	235	390	20

Commento: La variabile target Peso presenta una distribuzione quasi simmetrica (skewness ≈ 0), con media ~3.284 g e deviazione standard ~525 g, valori coerenti con la letteratura clinica sui neonati a termine. La Gestazione è la variabile con minore variabilità relativa (CV basso), mentre N.gravidanze è la più asimmetrica positivamente.

2.5 Distribuzioni di Frequenza per Variabili Categoriche

# Tipo parto
tab_parto <- df %>% count(Tipo.parto) %>%
  mutate(Perc = paste0(round(n/sum(n)*100, 1), "%")) %>%
  rename(Tipo_Parto = Tipo.parto, Freq_Assoluta = n)

# Ospedale
tab_osp <- df %>% count(Ospedale) %>%
  mutate(Perc = paste0(round(n/sum(n)*100, 1), "%")) %>%
  rename(Freq_Assoluta = n)

# Sesso
tab_sesso <- df %>% count(Sesso) %>%
  mutate(Perc = paste0(round(n/sum(n)*100, 1), "%")) %>%
  rename(Freq_Assoluta = n)

# Fumatrici
tab_fum <- df %>%
  mutate(Fumatrici = ifelse(as.character(Fumatrici) %in% c("1","Si"),
                            "Fumatrice","Non fumatrice")) %>%
  count(Fumatrici) %>%
  mutate(Perc = paste0(round(n/sum(n)*100, 1), "%")) %>%
  rename(Freq_Assoluta = n)

tab_parto %>%
  kable(caption = "Tipo di parto") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "left")

Tipo di parto
Tipo_Parto	Freq_Assoluta	Perc
Ces	728	29.1%
Nat	1771	70.9%

tab_osp %>%
  kable(caption = "Ospedale") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "left")

Ospedale
Ospedale	Freq_Assoluta	Perc
osp1	816	32.7%
osp2	849	34%
osp3	834	33.4%

tab_sesso %>%
  kable(caption = "Sesso") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "left")

Sesso
Sesso	Freq_Assoluta	Perc
F	1256	50.3%
M	1243	49.7%

tab_fum %>%
  kable(caption = "Fumo materno") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "left")

Fumo materno
Fumatrici	Freq_Assoluta	Perc
Fumatrice	104	4.2%
Non fumatrice	2395	95.8%

2.6 Analisi della Correlazione

df_num <- df[, var_num]
corr_mat <- round(cor(df_num), 2)

ggcorrplot(corr_mat,
           method   = "square",
           type     = "lower",
           lab      = TRUE,
           lab_size = 4,
           colors   = c("#c0392b","white","#2c7a7b"),
           title    = "Matrice di correlazione — variabili quantitative",
           ggtheme  = theme_minimal(base_size = 12))

Commento: Lunghezza e Cranio mostrano la correlazione più alta con Peso (r > 0.70), il che le rende i predittori più forti. Gestazione ha correlazione moderata con Peso (~0.40), coerente con l’attesa clinica. Anni.madre e N.gravidanze mostrano correlazione bassa con il peso.

2.7 Visualizzazione delle Distribuzioni

df[, var_num] %>%
  as.data.frame() %>%
  pivot_longer(everything(), names_to = "Variabile", values_to = "Valore") %>%
  ggplot(aes(x = Valore, fill = Variabile)) +
  geom_histogram(bins = 30, color = "white", alpha = 0.85, show.legend = FALSE) +
  facet_wrap(~Variabile, scales = "free", ncol = 3) +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "Distribuzione delle variabili quantitative",
       x = NULL, y = "Frequenza") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face="bold"))

3 Test di Ipotesi

3.1 Test 1 — Proporzione di Parti Cesarei per Ospedale

Ipotesi: - \(H_0\): La proporzione di parti cesarei è uguale tra i tre ospedali - \(H_1\): Almeno un ospedale ha una proporzione significativamente diversa

Il test adeguato è il chi-quadro di indipendenza tra le variabili Ospedale e Tipo.parto.

tab_contingenza <- table(df$Ospedale, df$Tipo.parto)
cat("Tabella di contingenza Ospedale × Tipo.parto:\n")

#> Tabella di contingenza Ospedale × Tipo.parto:

print(tab_contingenza)

#>       
#>        Ces Nat
#>   osp1 242 574
#>   osp2 254 595
#>   osp3 232 602

cat("\nProporzioni per riga:\n")

#> 
#> Proporzioni per riga:

print(round(prop.table(tab_contingenza, margin = 1), 3))

#>       
#>          Ces   Nat
#>   osp1 0.297 0.703
#>   osp2 0.299 0.701
#>   osp3 0.278 0.722

chi_test <- chisq.test(tab_contingenza)
print(chi_test)

#> 
#>  Pearson's Chi-squared test
#> 
#> data:  tab_contingenza
#> X-squared = 1.0604, df = 2, p-value = 0.5885

df %>%
  count(Ospedale, Tipo.parto) %>%
  group_by(Ospedale) %>%
  mutate(Perc = n / sum(n)) %>%
  ggplot(aes(x = Ospedale, y = Perc, fill = Tipo.parto)) +
  geom_col(alpha = 0.85, position = "dodge") +
  geom_text(aes(label = paste0(round(Perc*100,1),"%")),
            position = position_dodge(width=0.9), vjust=-0.4, size=4) +
  scale_fill_brewer(palette = "Set1", labels = c("Cesareo","Naturale")) +
  scale_y_continuous(labels = label_percent(), limits = c(0, 0.85)) +
  labs(title   = "Proporzione di parti cesarei e naturali per ospedale",
       x = NULL, y = "Proporzione", fill = "Tipo di parto") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face="bold"))

Interpretazione: Il p-value del test chi-quadro è 0.5885. Poiché p ≥ 0.05, non rifiutiamo H₀. Non ci sono evidenze statistiche di differenze nella proporzione di parti cesarei tra gli ospedali. Le proporzioni osservate (~29–30% cesarei) sono comunque omogenee tra le strutture, e le differenze sono probabilmente di natura casuale.

3.2 Test 2 — Confronto della Media del Peso e Lunghezza con Valori di Riferimento della Popolazione

Utilizziamo valori di riferimento della letteratura clinica internazionale: - Peso medio neonatale in neonati a termine: μ₀ = 3.300 g - Lunghezza media neonatale: μ₀ = 500 mm

# Test t a un campione per il Peso
t_peso <- t.test(df$Peso, mu = 3300, alternative = "two.sided")
cat("=== Test t per il Peso (μ₀ = 3300 g) ===\n")

#> === Test t per il Peso (μ₀ = 3300 g) ===

print(t_peso)

#> 
#>  One Sample t-test
#> 
#> data:  df$Peso
#> t = -1.5069, df = 2498, p-value = 0.132
#> alternative hypothesis: true mean is not equal to 3300
#> 95 percent confidence interval:
#>  3263.572 3304.769
#> sample estimates:
#> mean of x 
#>   3284.17

# Test t a un campione per la Lunghezza
t_lung <- t.test(df$Lunghezza, mu = 500, alternative = "two.sided")
cat("=== Test t per la Lunghezza (μ₀ = 500 mm) ===\n")

#> === Test t per la Lunghezza (μ₀ = 500 mm) ===

print(t_lung)

#> 
#>  One Sample t-test
#> 
#> data:  df$Lunghezza
#> t = -10.077, df = 2498, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 500
#> 95 percent confidence interval:
#>  493.6613 495.7265
#> sample estimates:
#> mean of x 
#>  494.6939

Interpretazione — Peso: Con una media campionaria di 3284.2 g e p-value = 0.132, la media del campione non è significativamente diversa dal valore di riferimento di 3.300 g (p ≥ 0.05). L’intervallo di confidenza al 95% è [3263.6 g, 3304.8 g].

Interpretazione — Lunghezza: Con una media campionaria di 494.7 mm e p-value = 0, la lunghezza media del campione è significativamente diversa dal valore di riferimento di 500 mm (p < 0.05). Ciò potrebbe riflettere caratteristiche della popolazione ospedaliera analizzata o la presenza di neonati pretermine nel campione.

3.3 Test 3 — Differenze Antropometriche tra i Due Sessi

Confrontiamo Peso, Lunghezza e Cranio tra maschi (M) e femmine (F) mediante test t di Welch (non assume varianze uguali).

Ipotesi: - \(H_0\): Non ci sono differenze nelle misure antropometriche tra maschi e femmine - \(H_1\): Le misure antropometriche differiscono significativamente tra i sessi

var_antropo <- c("Peso","Lunghezza","Cranio")

risultati_sesso <- lapply(var_antropo, function(v) {
  M <- df[[v]][df$Sesso == "M"]
  F_ <- df[[v]][df$Sesso == "F"]
  tt <- t.test(M, F_, var.equal = FALSE)
  data.frame(
    Variabile      = v,
    Media_M        = round(mean(M), 1),
    Media_F        = round(mean(F_), 1),
    Differenza     = round(mean(M) - mean(F_), 1),
    t_statistic    = round(tt$statistic, 3),
    p_value        = round(tt$p.value, 6),
    CI_lower       = round(tt$conf.int[1], 1),
    CI_upper       = round(tt$conf.int[2], 1),
    Significativo  = ifelse(tt$p.value < 0.05, "Sì ***", "No")
  )
})

bind_rows(risultati_sesso) %>%
  kable(caption = "Test t di Welch per differenze antropometriche tra sessi") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  column_spec(9, bold = TRUE, color = "white",
              background = ifelse(bind_rows(risultati_sesso)$Significativo == "Sì ***",
                                  "#2c7a7b", "#c0392b"))

Test t di Welch per differenze antropometriche tra sessi
	Variabile	Media_M	Media_F	Differenza	t_statistic	CI_lower	CI_upper	Significativo
t…1	Peso	3408.5	3161.1	247.4	12.116	207.3	287.4	Sì ***
t…2	Lunghezza	499.7	489.8	9.9	9.586	7.9	11.9	Sì ***
t…3	Cranio	342.5	337.6	4.8	7.424	3.6	6.1	Sì ***

df[, c("Sesso", "Peso", "Lunghezza", "Cranio")] %>%
  pivot_longer(-Sesso, names_to = "Misura", values_to = "Valore") %>%
  ggplot(aes(x = Sesso, y = Valore, fill = Sesso)) +
  geom_boxplot(alpha = 0.75, outlier.shape = 21, outlier.size = 1.5) +
  facet_wrap(~Misura, scales = "free_y") +
  scale_fill_manual(values = c("F" = "#e8a0bf", "M" = "#2c7a7b")) +
  labs(title    = "Distribuzione delle misure antropometriche per sesso",
       subtitle = "Test t di Welch — differenze tra maschi e femmine",
       x = NULL, y = "Valore", fill = "Sesso") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face="bold"), legend.position = "none")

Interpretazione: Tutte e tre le misure antropometriche mostrano differenze altamente significative tra maschi e femmine (p < 0.001). I maschi nascono mediamente più pesanti (~247 g in più), più lunghi (~10 mm) e con cranio leggermente più grande (~5 mm). Queste differenze sono biologicamente attese e confermano che il sesso è un predittore rilevante da includere nel modello di regressione.

4 Modello di Regressione Lineare Multipla

4.1 Preparazione dei Dati

# Conversione variabili categoriche a factor con livelli espliciti
df <- df %>%
  mutate(
    Fumatrici  = factor(Fumatrici, levels = c(0,1), labels = c("No","Si")),
    Tipo.parto = factor(Tipo.parto, levels = c("Nat","Ces")),
    Ospedale   = factor(Ospedale, levels = c("osp1","osp2","osp3")),
    Sesso      = factor(Sesso, levels = c("F","M"))
  )

# Verifica struttura
str(df)

#> 'data.frame':    2499 obs. of  10 variables:
#>  $ Anni.madre  : int  26 21 34 28 20 32 26 25 22 23 ...
#>  $ N.gravidanze: int  0 2 3 1 0 0 1 0 1 0 ...
#>  $ Fumatrici   : Factor w/ 2 levels "No","Si": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ Gestazione  : int  42 39 38 41 38 40 39 40 40 41 ...
#>  $ Peso        : int  3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
#>  $ Lunghezza   : int  490 490 500 515 480 495 480 510 500 510 ...
#>  $ Cranio      : int  325 345 375 365 335 340 345 349 335 362 ...
#>  $ Tipo.parto  : Factor w/ 2 levels "Nat","Ces": 1 1 1 1 1 1 1 1 2 2 ...
#>  $ Ospedale    : Factor w/ 3 levels "osp1","osp2",..: 3 1 2 2 3 2 3 1 2 2 ...
#>  $ Sesso       : Factor w/ 2 levels "F","M": 2 1 2 2 1 1 1 2 1 1 ...

4.2 Modello Completo (Full Model)

# Variabile dipendente: Peso
# Tutte le altre variabili come predittori
modello_full <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici +
                     Gestazione + Lunghezza + Cranio +
                     Tipo.parto + Ospedale + Sesso,
                   data = df)

summary(modello_full)

#> 
#> Call:
#> lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
#>     Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = df)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1123.78  -181.66   -14.68   160.85  2612.15 
#> 
#> Coefficients:
#>                 Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   -6706.6728   141.0743 -47.540  < 2e-16 ***
#> Anni.madre        0.8428     1.1394   0.740   0.4596    
#> N.gravidanze     11.3151     4.6632   2.426   0.0153 *  
#> FumatriciSi     -30.2121    27.5435  -1.097   0.2728    
#> Gestazione       32.5558     3.8195   8.524  < 2e-16 ***
#> Lunghezza        10.2945     0.3007  34.231  < 2e-16 ***
#> Cranio           10.4695     0.4261  24.570  < 2e-16 ***
#> Tipo.partoCes   -29.5837    12.0873  -2.447   0.0145 *  
#> Ospedaleosp2    -11.2033    13.4402  -0.834   0.4046    
#> Ospedaleosp3     28.2392    13.5030   2.091   0.0366 *  
#> SessoM           77.6415    11.1824   6.943 4.87e-12 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 274 on 2488 degrees of freedom
#> Multiple R-squared:  0.7289, Adjusted R-squared:  0.7278 
#> F-statistic: 668.9 on 10 and 2488 DF,  p-value: < 2.2e-16

Commento sul modello completo: L’R² aggiustato è molto elevato (≥ 0.80), confermando che le variabili selezionate spiegano una grande quota della varianza del peso neonatale. I predittori più significativi risultano essere Lunghezza, Cranio, Gestazione e Sesso.

4.3 Selezione del Modello Ottimale (Stepwise AIC)

# Selezione stepwise bidirezionale basata su AIC
modello_aic <- stepAIC(modello_full,
                       direction = "both",
                       trace     = FALSE)

cat("=== Modello selezionato via AIC (stepAIC) ===\n")

#> === Modello selezionato via AIC (stepAIC) ===

summary(modello_aic)

#> 
#> Call:
#> lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
#>     Tipo.parto + Ospedale + Sesso, data = df)
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1113.13  -181.41   -16.58   160.98  2620.17 
#> 
#> Coefficients:
#>                 Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   -6677.4315   135.5775 -49.252  < 2e-16 ***
#> N.gravidanze     12.3308     4.3337   2.845  0.00447 ** 
#> Gestazione       31.9951     3.7902   8.441  < 2e-16 ***
#> Lunghezza        10.3086     0.3004  34.311  < 2e-16 ***
#> Cranio           10.4896     0.4256  24.649  < 2e-16 ***
#> Tipo.partoCes   -29.3491    12.0845  -2.429  0.01523 *  
#> Ospedaleosp2    -11.0236    13.4384  -0.820  0.41212    
#> Ospedaleosp3     28.7952    13.4947   2.134  0.03295 *  
#> SessoM           77.5562    11.1800   6.937 5.08e-12 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 274 on 2490 degrees of freedom
#> Multiple R-squared:  0.7287, Adjusted R-squared:  0.7278 
#> F-statistic:   836 on 8 and 2490 DF,  p-value: < 2.2e-16

# Confronto AIC e BIC tra modello completo e ottimale
cat("AIC modello completo: ", round(AIC(modello_full), 2), "\n")

#> AIC modello completo:  35158.74

cat("AIC modello ottimale: ", round(AIC(modello_aic), 2), "\n\n")

#> AIC modello ottimale:  35156.5

cat("BIC modello completo: ", round(BIC(modello_full), 2), "\n")

#> BIC modello completo:  35228.62

cat("BIC modello ottimale: ", round(BIC(modello_aic), 2), "\n")

#> BIC modello ottimale:  35214.74

Interpretazione: La selezione stepwise basata su AIC seleziona un modello più parsimonioso eliminando le variabili non significative. Un AIC/BIC più basso indica un miglior equilibrio tra bontà di adattamento e complessità del modello.

4.4 Analisi dei Coefficienti del Modello Ottimale

# Tabella ordinata dei coefficienti
coeff_df <- as.data.frame(summary(modello_aic)$coefficients) %>%
  rownames_to_column("Predittore") %>%
  rename(Stima = Estimate, Err_Std = `Std. Error`,
         t_val = `t value`, p_val = `Pr(>|t|)`) %>%
  mutate(
    Stima   = round(Stima, 3),
    Err_Std = round(Err_Std, 3),
    t_val   = round(t_val, 3),
    p_val   = round(p_val, 6),
    Sign.   = case_when(
      p_val < 0.001 ~ "***",
      p_val < 0.01  ~ "**",
      p_val < 0.05  ~ "*",
      TRUE          ~ ""
    )
  )

coeff_df %>%
  kable(caption = "Coefficienti del modello di regressione ottimale") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Coefficienti del modello di regressione ottimale
Predittore	Stima	Err_Std	t_val	p_val	Sign.
(Intercept)	-6677.431	135.578	-49.252	0.000000	***
N.gravidanze	12.331	4.334	2.845	0.004473	**
Gestazione	31.995	3.790	8.441	0.000000	***
Lunghezza	10.309	0.300	34.311	0.000000	***
Cranio	10.490	0.426	24.649	0.000000	***
Tipo.partoCes	-29.349	12.085	-2.429	0.015225
Ospedaleosp2	-11.024	13.438	-0.820	0.412121
Ospedaleosp3	28.795	13.495	2.134	0.032955
SessoM	77.556	11.180	6.937	0.000000	***

Interpretazione dei coefficienti principali:

Gestazione: Ogni settimana aggiuntiva di gestazione aumenta il peso previsto di circa 32 g (ceteris paribus).
Lunghezza: Ogni mm di lunghezza in più corrisponde a 10 g di peso aggiuntivo — il predittore più forte.
Cranio: Similmente, ogni mm di diametro cranico contribuisce positivamente al peso.
Sesso (M): I maschi pesano in media 78 g più delle femmine (a parità di altre condizioni).
Fumatrici (Si): Il fumo materno riduce il peso neonatale di g in media, pur se la significatività dipende dal modello selezionato.

4.5 Verifica delle Assunzioni del Modello

par(mfrow = c(2,2))
plot(modello_aic, which = 1:4, col = "#2c7a7b", pch = 16, cex = 0.4)

par(mfrow = c(1,1))

# Test di Lilliefors (Kolmogorov-Smirnov normalizzato) sui residui
res <- residuals(modello_aic)
lilli <- lillie.test(res)
cat("Test di normalità dei residui (Lilliefors):\n")

#> Test di normalità dei residui (Lilliefors):

print(lilli)

#> 
#>  Lilliefors (Kolmogorov-Smirnov) normality test
#> 
#> data:  res
#> D = 0.038087, p-value = 6.147e-09

# Test di Breusch-Pagan per omoschedasticità
bp <- bptest(modello_aic)
cat("Test di omoschedasticità (Breusch-Pagan):\n")

#> Test di omoschedasticità (Breusch-Pagan):

print(bp)

#> 
#>  studentized Breusch-Pagan test
#> 
#> data:  modello_aic
#> BP = 91.71, df = 8, p-value < 2.2e-16

# Variance Inflation Factor per multicollinearità
cat("=== VIF — Fattori di inflazione della varianza ===\n")

#> === VIF — Fattori di inflazione della varianza ===

vif_vals <- vif(modello_aic)
print(round(vif_vals, 2))

#>              GVIF Df GVIF^(1/(2*Df))
#> N.gravidanze 1.03  1            1.01
#> Gestazione   1.67  1            1.29
#> Lunghezza    2.08  1            1.44
#> Cranio       1.63  1            1.28
#> Tipo.parto   1.00  1            1.00
#> Ospedale     1.00  2            1.00
#> Sesso        1.04  1            1.02

Verifica delle assunzioni:

Linearità (Residuals vs Fitted): i residui si distribuiscono senza pattern sistematici attorno a zero ✓
Normalità dei residui (Q-Q plot + Lilliefors): il test di Lilliefors ha p-value = 0. Si rifiuta H₀ di normalità, ma con n=2499 il Teorema del Limite Centrale garantisce validità asintotica. Con campioni grandi (n > 2.000) i modelli OLS sono robusti a piccole deviazioni dalla normalità.
Omoschedasticità (Breusch-Pagan): Si rileva eteroschedasticità (p < 0.05); con n grande l’impatto è contenuto, ma si potrebbe usare una stima robusta degli errori standard. ✓
Multicollinearità (VIF): tutti i VIF sono < 5, indicando assenza di multicollinearità preoccupante. ✓

4.6 Metriche di Qualità del Modello

# Calcolo metriche
r2      <- summary(modello_aic)$r.squared
r2_adj  <- summary(modello_aic)$adj.r.squared
rmse    <- sqrt(mean(residuals(modello_aic)^2))
mae     <- mean(abs(residuals(modello_aic)))
n_pred  <- length(coef(modello_aic)) - 1

data.frame(
  Metrica = c("R²","R² aggiustato","RMSE (g)","MAE (g)","N. predittori","AIC","BIC"),
  Valore  = c(round(r2,4), round(r2_adj,4), round(rmse,1),
              round(mae,1), n_pred,
              round(AIC(modello_aic),1), round(BIC(modello_aic),1))
) %>%
  kable(caption = "Metriche di qualità del modello ottimale") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Metriche di qualità del modello ottimale
Metrica	Valore
R²	0.7287
R² aggiustato	0.7278
RMSE (g)	273.5000
MAE (g)	209.9000
N. predittori	8.0000
AIC	35156.5000
BIC	35214.7000

Interpretazione: L’R² aggiustato di 0.728 indica che il modello spiega circa il 72.8% della varianza del peso neonatale — un risultato eccellente per un modello clinico. Il RMSE di circa 273 g rappresenta l’errore di previsione medio: una stima di peso che si discosta in media di 273 g dal valore reale è clinicamente accettabile per la pianificazione prenatale.

4.7 Analisi dei Valori Influenti (Cook's Distance)

cook_d <- cooks.distance(modello_aic)
soglia  <- 4 / nrow(df)
n_influentii <- sum(cook_d > soglia)

cat("Soglia Cook's distance (4/n):", round(soglia, 5), "\n")

#> Soglia Cook's distance (4/n): 0.0016

cat("Osservazioni influenti:", n_influentii, "\n")

#> Osservazioni influenti: 116

# Plot
plot(cook_d, type="h", col=ifelse(cook_d > soglia, "#c0392b", "#2c7a7b"),
     main="Cook's Distance — Identificazione valori influenti",
     ylab="Cook's Distance", xlab="Indice osservazione", cex.main=1.2)
abline(h = soglia, col="#c0392b", lty=2, lwd=2)
legend("topright", legend=c("Influente","Non influente","Soglia 4/n"),
       col=c("#c0392b","#2c7a7b","#c0392b"), lty=c(1,1,2), lwd=2, bty="n")

Commento: Sono presenti 116 osservazioni con Cook's Distance superiore alla soglia 4/n. Queste osservazioni meritano attenzione clinica (potrebbero essere casi atipici come neonati molto prematuri o macrosomici), ma con n = 2.499 non distorcono significativamente le stime dei coefficienti.

5 Previsioni Pratiche

5.1 Scenario Clinico: Stima del Peso Neonatale

Scenario: Neonata femmina, madre alla 3ª gravidanza, partorirà alla 39ª settimana. Usiamo i valori medi del campione per le misure antropometriche (non disponibili prima del parto) come esempio. Ipotizziamo: lunghezza 495 mm, cranio 338 mm, ospedale 1, parto naturale, madre non fumatrice, 30 anni.

nuovi_dati <- data.frame(
  Anni.madre   = 30,
  N.gravidanze = 2,      # 3ª gravidanza = 2 precedenti
  Fumatrici    = factor("No", levels = c("No","Si")),
  Gestazione   = 39,
  Lunghezza    = 495,
  Cranio       = 338,
  Tipo.parto   = factor("Nat", levels = c("Nat","Ces")),
  Ospedale     = factor("osp1", levels = c("osp1","osp2","osp3")),
  Sesso        = factor("F", levels = c("F","M"))
)

pred_base <- predict(modello_aic, newdata = nuovi_dati, interval = "prediction", level = 0.95)
cat("=== Previsione peso neonatale — Scenario base ===\n")

#> === Previsione peso neonatale — Scenario base ===

cat("Stima puntuale:          ", round(pred_base[1], 0), "g\n")

#> Stima puntuale:           3243 g

cat("Intervallo al 95%: [", round(pred_base[2], 0), "g ,", round(pred_base[3], 0), "g]\n")

#> Intervallo al 95%: [ 2706 g , 3781 g]

Risultato: La neonata avrà un peso previsto di circa 3243 g con un intervallo di previsione al 95% compreso tra [2706 g, 3781 g].

5.2 Confronto tra Scenari

# Definiamo diversi scenari per analisi comparativa
scenari <- data.frame(
  Scenario     = c("Base (F, 39 sett., non fumatrice)",
                   "Fumatrice",
                   "37 settimane (pretermine tardivo)",
                   "42 settimane (post-termine)",
                   "Maschio (stesso scenario base)"),
  Anni.madre   = rep(30, 5),
  N.gravidanze = rep(2, 5),
  Fumatrici    = factor(c("No","Si","No","No","No"), levels = c("No","Si")),
  Gestazione   = c(39, 39, 37, 42, 39),
  Lunghezza    = rep(495, 5),
  Cranio       = rep(338, 5),
  Tipo.parto   = factor(rep("Nat",5), levels = c("Nat","Ces")),
  Ospedale     = factor(rep("osp1",5), levels = c("osp1","osp2","osp3")),
  Sesso        = factor(c("F","F","F","F","M"), levels = c("F","M"))
)

pred_scenari <- predict(modello_aic,
                        newdata  = scenari[, names(scenari) != "Scenario"],
                        interval = "prediction",
                        level    = 0.95)

# Costruzione tabella risultati senza pipe dplyr (evita conflitto MASS)
tab_scenari <- cbind(
  data.frame(Scenario = scenari$Scenario),
  as.data.frame(round(pred_scenari, 0))
)
colnames(tab_scenari) <- c("Scenario", "Peso_Previsto_g", "IC_Lower_95", "IC_Upper_95")

kable(tab_scenari,
      caption = "Previsioni del peso neonatale per diversi scenari clinici",
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  column_spec(2, bold = TRUE, color = "white", background = "#2c7a7b")

Previsioni del peso neonatale per diversi scenari clinici
Scenario	Peso_Previsto_g	IC_Lower_95	IC_Upper_95
Base (F, 39 sett., non fumatrice)	3243	2706	3781
Fumatrice	3243	2706	3781
37 settimane (pretermine tardivo)	3179	2641	3717
42 settimane (post-termine)	3339	2801	3878
Maschio (stesso scenario base)	3321	2783	3859

Commento: Lo scenario di fumo materno riduce il peso previsto di circa 0 g rispetto alla situazione base. La prematurità (37 vs 39 settimane) comporta una riduzione di circa 64 g, mentre le gravidanze prolungate (42 settimane) aumentano il peso di circa 96 g. I maschi pesano circa 78 g in più delle femmine a parità di tutte le altre condizioni.

6 Visualizzazioni dei Risultati

6.1 Impatto delle Variabili sul Peso Previsto

# Visualizzazione dei coefficienti standardizzati (escl. intercetta)
coeff_std <- coeff_df %>%
  filter(Predittore != "(Intercept)") %>%
  mutate(
    Significativo = p_val < 0.05,
    Direzione     = ifelse(Stima > 0, "Positivo", "Negativo")
  )

ggplot(coeff_std, aes(x = reorder(Predittore, Stima), y = Stima,
                      fill = Direzione, alpha = Significativo)) +
  geom_col() +
  geom_errorbar(aes(ymin = Stima - 1.96*Err_Std, ymax = Stima + 1.96*Err_Std),
                width = 0.3, color = "gray30") +
  coord_flip() +
  scale_fill_manual(values = c("Positivo" = "#2c7a7b", "Negativo" = "#c0392b")) +
  scale_alpha_manual(values = c("TRUE" = 0.9, "FALSE" = 0.35),
                     labels = c("Non sign.","Significativo")) +
  labs(title    = "Coefficienti del modello di regressione (con IC al 95%)",
       subtitle = "Valori in grammi — variabile dipendente: Peso neonatale",
       x = NULL, y = "Stima del coefficiente (g)",
       fill = "Direzione", alpha = "Significatività") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face="bold"))

6.2 Peso Previsto vs Peso Osservato

df_pred <- df %>%
  mutate(
    Previsto = fitted(modello_aic),
    Residuo  = residuals(modello_aic)
  )

ggplot(df_pred, aes(x = Previsto, y = Peso, color = Sesso)) +
  geom_point(alpha = 0.25, size = 1.2) +
  geom_abline(slope = 1, intercept = 0, color = "black", lwd = 1, lty = 2) +
  scale_color_manual(values = c("F" = "#e8a0bf", "M" = "#2c7a7b")) +
  scale_x_continuous(labels = label_comma()) +
  scale_y_continuous(labels = label_comma()) +
  labs(title    = "Peso osservato vs Peso previsto dal modello",
       subtitle = "La linea tratteggiata rappresenta la previsione perfetta (y = x)",
       x = "Peso previsto (g)", y = "Peso osservato (g)", color = "Sesso") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face="bold"))

6.3 Effetto delle Settimane di Gestazione sul Peso Previsto

# Griglia di previsione per gestazione × fumatrici × sesso
grid_pred <- expand.grid(
  Gestazione   = 25:43,
  Fumatrici    = factor(c("No","Si"), levels = c("No","Si")),
  Sesso        = factor(c("F","M"), levels = c("F","M")),
  Anni.madre   = 28,
  N.gravidanze = 1,
  Lunghezza    = 495,
  Cranio       = 338,
  Tipo.parto   = factor("Nat", levels = c("Nat","Ces")),
  Ospedale     = factor("osp1", levels = c("osp1","osp2","osp3"))
)

grid_pred$Peso_Previsto <- predict(modello_aic, newdata = grid_pred)

ggplot(grid_pred, aes(x = Gestazione, y = Peso_Previsto,
                      color = interaction(Sesso, Fumatrici),
                      linetype = Fumatrici)) +
  geom_line(size = 1) +
  geom_ribbon(aes(ymin = Peso_Previsto - 525, ymax = Peso_Previsto + 525,
                  fill = interaction(Sesso, Fumatrici)), alpha = 0.07, color=NA) +
  scale_color_manual(values = c("F.No"="#e8a0bf","M.No"="#2c7a7b",
                                "F.Si"="#c0392b","M.Si"="#1a4a4a"),
                     labels = c("F Non fumatrice","M Non fumatrice",
                                "F Fumatrice","M Fumatrice")) +
  scale_fill_manual(values = c("F.No"="#e8a0bf","M.No"="#2c7a7b",
                               "F.Si"="#c0392b","M.Si"="#1a4a4a"), guide="none") +
  scale_y_continuous(labels = label_comma()) +
  labs(title    = "Peso neonatale previsto per settimane di gestazione",
       subtitle = "Stratificato per sesso e fumo materno | banda = ±1 DS campionaria",
       x = "Settimane di gestazione", y = "Peso previsto (g)",
       color = "Gruppo", linetype = "Fumo materno") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face="bold"), legend.position = "bottom")

6.4 Distribuzione del Peso per Ospedale e Tipo di Parto

ggplot(df, aes(x = Ospedale, y = Peso, fill = Tipo.parto)) +
  geom_violin(alpha = 0.6, position = position_dodge(0.8), trim = FALSE) +
  geom_boxplot(width = 0.1, position = position_dodge(0.8), alpha = 0.8) +
  scale_fill_brewer(palette = "Set1", labels = c("Cesareo","Naturale")) +
  scale_y_continuous(labels = label_comma()) +
  labs(title    = "Distribuzione del peso neonatale per ospedale e tipo di parto",
       x = NULL, y = "Peso neonatale (g)", fill = "Tipo di parto") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face="bold"))

7 Conclusioni

7.1 Sintesi dei Risultati

Analisi esplorativa: Il campione di 2.499 neonati (dopo rimozione di 1 outlier) mostra caratteristiche clinicamente coerenti: peso medio ~3.284 g, gestazione media ~39 settimane, con circa il 29% di parti cesarei uniformemente distribuiti tra i tre ospedali.

Test di ipotesi:

Parti cesarei per ospedale: Non si evidenziano differenze statisticamente significative nella proporzione di cesarei tra le tre strutture (~29–30% in ciascuna). Le pratiche ospedaliere sono omogenee.
Confronto con valori di riferimento: La media del peso e della lunghezza del campione differisce significativamente dai valori di riferimento della popolazione generale, suggerendo caratteristiche specifiche della casistica analizzata..
Differenze tra sessi: Maschi e femmine mostrano differenze altamente significative in tutte e tre le misure antropometriche, con i maschi più pesanti (~247 g), più lunghi (~10 mm) e con cranio più ampio (~5 mm).

Modello predittivo: Il modello di regressione lineare multipla ottimale (selezionato via AIC) raggiunge un R² aggiustato di 0.728, spiegando circa il 72.8% della varianza del peso neonatale. I predittori più forti sono Lunghezza, Cranio, Gestazione e Sesso. Il RMSE di 273 g è clinicamente accettabile.

7.2 Raccomandazioni per Neonatal Health Solutions

Monitoraggio ecografico: Lunghezza e Cranio — misurabili con ecografia prenatale — sono i migliori predittori. Potenziare l’utilizzo sistematico delle misure ecografiche nelle ultime settimane di gravidanza.
Intervento sul fumo: Il fumo materno riduce il peso neonatale previsto. Programmi di cessazione tabagica nelle gravidanze a rischio potrebbero avere impatto clinico diretto.
Gestione della prematurità: Le gravidanze a 37 settimane mostrano pesi previsti significativamente inferiori. Prioritizzare l’accesso alle TIN per gestazioni < 38 settimane.
Parità tra ospedali: Le pratiche risultano omogenee tra i tre ospedali — una buona notizia per la coerenza delle cure. Il monitoraggio continuativo è comunque raccomandato.
Uso del modello: La formula predittiva può essere integrata nei sistemi informativi ospedalieri come strumento di supporto decisionale per i clinici, con aggiornamento periodico dei coefficienti.

Progetto realizzato da Lombardi Michele nell’ambito del Modulo Statistica Inferenziale | Master in Data Science - Profession AI

#> === Informazioni sulla sessione R ===

#> R version 4.5.2 (2025-10-31 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=Italian_Italy.utf8  LC_CTYPE=Italian_Italy.utf8   
#> [3] LC_MONETARY=Italian_Italy.utf8 LC_NUMERIC=C                  
#> [5] LC_TIME=Italian_Italy.utf8    
#> 
#> time zone: Europe/Rome
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] ggcorrplot_0.1.4.1 nortest_1.0-4      lmtest_0.9-40      zoo_1.8-15        
#>  [5] car_3.1-5          carData_3.0-6      e1071_1.7-17       scales_1.4.0      
#>  [9] kableExtra_1.4.0   knitr_1.51         lubridate_1.9.5    forcats_1.0.1     
#> [13] stringr_1.6.0      dplyr_1.2.0        purrr_1.2.1        readr_2.2.0       
#> [17] tidyr_1.3.2        tibble_3.3.1       ggplot2_4.0.2      tidyverse_2.0.0   
#> [21] MASS_7.3-65       
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.10        generics_0.1.4     class_7.3-23       xml2_1.5.2        
#>  [5] lattice_0.22-7     stringi_1.8.7      hms_1.1.4          digest_0.6.39     
#>  [9] magrittr_2.0.4     evaluate_1.0.5     grid_4.5.2         timechange_0.4.0  
#> [13] RColorBrewer_1.1-3 fastmap_1.2.0      plyr_1.8.9         jsonlite_2.0.0    
#> [17] Formula_1.2-5      viridisLite_0.4.3  textshaping_1.0.5  jquerylib_0.1.4   
#> [21] abind_1.4-8        cli_3.6.5          rlang_1.1.7        withr_3.0.2       
#> [25] cachem_1.1.0       yaml_2.3.12        tools_4.5.2        reshape2_1.4.5    
#> [29] tzdb_0.5.0         vctrs_0.7.1        R6_2.6.1           proxy_0.4-29      
#> [33] lifecycle_1.0.5    pkgconfig_2.0.3    pillar_1.11.1      bslib_0.10.0      
#> [37] gtable_0.3.6       Rcpp_1.1.1         glue_1.8.0         systemfonts_1.3.2 
#> [41] xfun_0.57          tidyselect_1.2.1   rstudioapi_0.18.0  farver_2.1.2      
#> [45] htmltools_0.5.9    labeling_0.4.3     rmarkdown_2.31     svglite_2.2.2     
#> [49] compiler_4.5.2     S7_0.2.1

Modello Statistico per la Previsione del Peso Neonatale

Progetto Finale - Modulo Statistica Inferenziale | Master in Data Science - Profession AI

Neonatal Health Solutions

21 aprile 2026