# Impostazioni globali: echo=TRUE mostra il codice, warning/message=FALSE puliscono l'output
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, fig.align = "center")
# Caricamento pacchetti
required_packages <- c("dplyr", "moments", "ggplot2", "scales", "car", "lmtest", "MASS", "knitr", "kableExtra")
new_packages <- required_packages[!(required_packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)

library(dplyr)
library(moments)
library(ggplot2)
library(scales)
library(car)    # VIF e diagnostica
library(lmtest) # Test sui residui
library(MASS)   # Box-Cox
library(knitr)  # Tabelle

1. Introduzione e preparazione del dataset

L’obiettivo di questa analisi è identificare i determinanti statisticamente significativi del peso alla nascita e costruire un modello di regressione lineare parsimonioso e robusto.

# Caricamento e pulizia preliminare
setwd("~/Desktop/Università/Master Data Science/Progetto 3 - Statistica inferenziale") # Impostare directory se necessario
data <- read.csv("neonati.csv", sep=",")

# Rimozione casi incompleti (Listwise Deletion) per garantire comparabilità dei modelli
data <- na.omit(data)

# Conversione delle variabili categoriche in fattori
data$Tipo.parto <- as.factor(data$Tipo.parto)
data$Ospedale <- as.factor(data$Ospedale)
data$Sesso <- as.factor(data$Sesso)
data$Fumatrici <- as.factor(data$Fumatrici)

Eseguiamo un controllo preliminare sulla variabile Anni.madre per identificare valori non plausibili (es. età < 15 anni) che potrebbero rappresentare errori di registrazione.

# Filtro biologico: Esclusione madri con età < 15 anni
n_pre <- nrow(data)
data <- data %>% filter(Anni.madre >= 15)
n_post <- nrow(data)

if(n_pre > n_post) {
  cat("Rimossi", n_pre - n_post, "record con età materna non plausibile (< 15 anni).\n")
}
## Rimossi 5 record con età materna non plausibile (< 15 anni).

2. Analisi Esplorativa (EDA)

2.1 Statistiche Descrittive Univariate

Analizziamo la tendenza centrale e la dispersione delle variabili quantitative. Utilizziamo il coefficiente di variazione (CV) per confrontare la variabilità relativa tra grandezze diverse.

stats_summary <- data.frame()

for (var in colnames(data)) {
  vec <- data[[var]]
  if (is.numeric(vec)) {
    temp <- data.frame(
      Variabile = var,
      Media = mean(vec),
      SD = sd(vec),
      CV_Perc = (sd(vec) / mean(vec)) * 100,
      Skewness = skewness(vec),
      Kurtosis = kurtosis(vec)
    )
    stats_summary <- rbind(stats_summary, temp)
  }
}

kable(stats_summary, digits = 2, caption = "Statistiche Descrittive Variabili Numeriche")
Statistiche Descrittive Variabili Numeriche
Variabile Media SD CV_Perc Skewness Kurtosis
Anni.madre 28.20 5.20 18.42 0.17 2.87
N.gravidanze 0.98 1.28 130.44 2.51 13.97
Gestazione 38.98 1.87 4.80 -2.06 11.25
Peso 3284.20 525.39 16.00 -0.65 5.03
Lunghezza 494.71 26.34 5.32 -1.52 9.48
Cranio 340.02 16.43 4.83 -0.79 5.95 
var_continue <- c("Anni.madre", "Peso", "Lunghezza", "Cranio")

for (var in var_continue) {
  mean_val <- mean(data[[var]], na.rm = TRUE)
  p <- ggplot(data, aes_string(x = var)) +
    geom_histogram(aes(y = ..density..), bins = 30, fill = "lightblue", color = "white") +
    geom_density(alpha = 0.2, fill = "blue") +
    geom_vline(aes(xintercept = mean_val), color = "red", linetype = "dashed", size = 1) +
    labs(title = paste("Distribuzione di:", var),
         subtitle = paste("Linea rossa = Media (", round(mean_val, 2), ")"),
         x = var, y = "Densità") +
    theme_minimal()
  print(p)
}

Distribuzione più asimmetrica: Numero gravidanze e Gestazione, alta asimettria ma a senso opposto per queste due variabili, infatti vediamo come in N.gravidanze la maggioranza delle osservazioni si trovano nella parte sinistra (quindi a valori bassi 0-1) mentre per Gestazione si hanno frequenze assolute maggiori a valori più alti (da 30 sett in poi).

Distribuzione più normale: Anni.madre risulta avere un valore di Skewness molto prossimo allo 0 che indica normalità in distribuzione, tuttavia il test di Shapiro effettuato su questa variabile ci dice che la distribuzione è significativamente non normale, lo stesso per la variabile Peso, che a vista sembrerebbe abbastanza normale con un valore di skewness anch’esso prossimo allo 0 ma comunque il test di normalità ci suggerisce che è significativamente distribuita in modo non normale.

2.2 Identificazione degli outliers

Utilizziamo il criterio dell’intervallo interquartile (\(1.5 \times IQR\)) per identificare i valori anomali nelle variabili quantitative.

df_outliers <- data.frame()

for (var in colnames(data)) {
  colonna <- data[[var]]
  if (is.numeric(colonna)) {
    Q1 <- quantile(colonna, 0.25)
    Q3 <- quantile(colonna, 0.75)
    IQR_val <- Q3 - Q1
    
    # Identificazione indici outlier
    idx <- which(colonna < (Q1 - 1.5 * IQR_val) | colonna > (Q3 + 1.5 * IQR_val))
    
    if (length(idx) > 0) {
      df_outliers <- rbind(df_outliers, data.frame(Variabile = var, Riga_Indice = idx, Valore = colonna[idx]))
    }
    
    # Visualizzazione Boxplot
    boxplot(colonna, main = paste("Boxplot:", var), horizontal = TRUE, col = "lightblue", pch = 19)
  }
}

# Analisi specifica della composizione degli outliers
if (nrow(df_outliers) > 0) {
  indici_out <- unique(df_outliers$Riga_Indice)
  subset_out <- data[indici_out, ]
  subset_clean <- data[-indici_out, ]
  
  # Confronto medie chiave
  mean_grav_out <- mean(subset_out$N.gravidanze)
  mean_grav_clean <- mean(subset_clean$N.gravidanze)
  
  cat("\nAnalisi Comparativa Outliers vs Popolazione:\n")
  cat("Media N.Gravidanze (outliers):", round(mean_grav_out, 2), "\n")
  cat("Media N.Gravidanze (totale):  ", round(mean(data$N.gravidanze), 2), "\n")
}
## 
## Analisi Comparativa Outliers vs Popolazione:
## Media N.Gravidanze (outliers): 2.89 
## Media N.Gravidanze (totale):   0.98

Commento sull’analisi degli outliers: l’analisi grafica e descrittiva evidenzia che la presenza di outliers è particolarmente rilevante per la variabile N.gravidanze, mentre nel dataset completo la media è di circa 1 gravidanza pregressa, nel sottogruppo identificato come “anomalo” la media sale drasticamente a circa 3. Questo indica che le madri pluripare con elevato numero di parti costituiscono un cluster statistico distinto. Per le variabili biometriche (“Lunghezza”, “Cranio”), gli outlier rappresentano casi estremi fisiologici che è opportuno mantenere nel modello per garantirne la generalizzabilità.

2.3 Relazioni bivariate (variabili qualitative)

Esaminiamo visivamente la relazione tra le variabili categoriche e la variabile target (Peso) tramite Boxplot.

vars_cat <- c("Fumatrici", "Sesso", "Tipo.parto", "Ospedale")

par(mfrow=c(2,2))
for (var in vars_cat) {
  boxplot(as.formula(paste("Peso ~", var)), data = data, 
          main = paste("Peso vs", var), col = "lightgreen")
}

par(mfrow=c(1,1))

t_test_sex <- t.test(Peso ~ Sesso, data = data, var.equal = FALSE)
print(t_test_sex)
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.077, df = 2486.1, 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.0099 -206.8273
## sample estimates:
## mean in group F mean in group M 
##        3161.381        3408.300

Osservazione: Si nota visivamente una differenza di peso mediana a favore del sesso maschile. Le differenze legate all’ospedale o al tipo di parto appaiono trascurabili, suggerendo una loro scarsa rilevanza predittiva. Il t-test sulle medie quantifica a livello significativo la differenza tra le medie delle due modalità di “Sesso” in relazione alla variabile “Peso”, per questo decido di lasciare la variabile “Sesso” come campo di controllo durante la fase di modellazione.

3. Analisi di correlazione e multicollinearità

Esaminiamo le relazioni lineari tra le variabili quantitative per prevenire problemi di multicollinearità nei modelli di regressione.

# Matrice di correlazione
numeric_cols <- sapply(data, is.numeric)
cor_matrix <- cor(data[, numeric_cols])

kable(cor_matrix, digits = 2, caption = "Matrice di Correlazione di Pearson")
Matrice di Correlazione di Pearson
Anni.madre N.gravidanze Gestazione Peso Lunghezza Cranio
Anni.madre 1.00 0.38 -0.14 -0.02 -0.07 0.02
N.gravidanze 0.38 1.00 -0.10 0.00 -0.06 0.04
Gestazione -0.14 -0.10 1.00 0.59 0.62 0.46
Peso -0.02 0.00 0.59 1.00 0.80 0.70
Lunghezza -0.07 -0.06 0.62 0.80 1.00 0.60
Cranio 0.02 0.04 0.46 0.70 0.60 1.00 
# Focus: Età Madre vs Peso (Test di Spearman per non linearità/non normalità)
cor_test <- cor.test(data$Anni.madre, data$Peso, method = "spearman")
print(cor_test)
## 
##  Spearman's rank correlation rho
## 
## data:  data$Anni.madre and data$Peso
## S = 2605601257, p-value = 0.7426
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##          rho 
## -0.006578443
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, use = "complete.obs")
    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 * abs(r) + 0.5) 
}

numeric_cols <- sapply(data, is.numeric)
data_numeric <- na.omit(data[, numeric_cols])

pairs(data_numeric, 
      lower.panel = panel.cor, 
      upper.panel = panel.smooth, 
      main = "Matrice di Scatterplot e Correlazioni",
      pch = 20, col = alpha("black", 0.4))

Commento: esiste una forte correlazione positiva tra le misure biometriche (“Lunghezza”, “Cranio”) e il “Peso”, come atteso biologicamente.

La correlazione tra “Anni.madre” e “Peso” non è statisticamente significativa (\(p > 0.05\)), suggerendo che l’età materna non ha un effetto lineare diretto rilevante sul peso del nascituro in questo campione.

Peso e Gestazione: giustificato dal fatto che più aumenta il numero di settimane di gestazione e maggiore sarà il peso, ed ha senso dal momento che l’embrione cresce sempre di più al passare delle settimane. Ma questa non sembra avere un andamente proprio lineare, ma anzi logaritmico, in cui si arriva quasi a saturazione dopo molto tempo di gestazione (intorno a 40 sett).

Lunghezza e Crano: notiamo come la correlazione è abbastanza alta tra queste due misure, e che a loro volta sono correlate anche alla gestazione, anche loro non esattamente in modo lineare ma leggeremente curvo. Questo ha senso dal mmomento che il feto cresce all’aumentare delle settimane di gestazione e quindi crescerà il diametro del cranio e la sua lunghezza fino a saturarsi e quindi fermarsi di crescere.

Peso e lunghezza: andamento più lineare per questo, ma comunque altissima correlazione (0,8), quindi all’aumentare della lunghezza aumenta anche il peso del nascituro.

4. Modellazione Statistica

Procediamo alla costruzione dei modelli adottando un approccio logico (“foward selection” ragionata) piuttosto che puramente algoritmico.

4.1 Selezione delle Variabili

Escludiamo a priori le variabili “Ospedale” e “Tipo.parto”.

Motivazione logica: in un contesto di predizione prenatale o neonatale, l’ospedale di nascita è una variabile logistica che non dovrebbe influenzare biologicamente il peso (salvo bias di selezione non oggetto di studio). Il tipo di parto è spesso una conseguenza delle dimensioni del feto (es. cesareo per macrosomia), non una causa, e includerlo potrebbe introdurre endogeneità.

4.2 Costruzione dei Modelli Lineari

Costruiamo tre modelli principali per confrontare la complessità:

  1. Modello Base: Solo variabili biometriche essenziali.
  2. Modello Intermedio: Aggiunta di variabili materne e sesso.
  3. Modello con Interazioni (Mod9): Inclusione di termini di interazione per catturare la tridimensionalità del feto.
# Modello Base: solo "Gestazione" e "Sesso" (variabili note prima del parto)
mod_base <- lm(Peso ~ Gestazione + Sesso, data = data)

# Modello Intermedio: aggiunta misure biometriche ("Cranio", "Lunghezza") e "N.gravidanze"
# rimuoviamo "Fumatrici" e "Anni.madre" in quanto risultati non significativi in EDA
mod_inter <- lm(Peso ~  Gestazione + Sesso + Lunghezza + Cranio + N.gravidanze, data = data)
plot(mod_inter, caption="Plot mod_inter")

# Modello Avanzato (Mod9): aggiunta interazione Cranio:Lunghezza
# L'interazione è giustificata fisicamente: il peso è funzione del volume (Lunghezza * Circonferenza)
mod9 <- lm(Peso ~  Gestazione + Sesso + N.gravidanze + Lunghezza * Cranio, data = data)
plot(mod9, caption="Plot mod9")

# Confronto rapido tramite ANOVA (test F per modelli annidati)
anova_res <- anova(mod_inter, mod9)
kable(anova_res, digits = 3, caption = "Test ANOVA: Modello Intermedio vs Modello con Interazione")
Test ANOVA: Modello Intermedio vs Modello con Interazione
Res.Df RSS Df Sum of Sq F Pr(>F)
2489 188039612 NA NA NA NA
2488 186292253 1 1747359 23.337 0

Esito selezione: Il test ANOVA mostra un p-value estremamente significativo (\(< 1.442e-06\)) per l’aggiunta dell’interazione. Questo conferma che mod9 cattura meglio la variabilità dei dati rispetto al modello additivo semplice.

4.3 Modellazione non lineare (trasformazioni)

Verifichiamo se una trasformazione della variabile target o dei regressori migliora le performance, dato che la crescita biologica è spesso non lineare.

# Modello 12: trasformazione logaritmica della target (log-Level)
# Utile per stabilizzare la varianza dei residui
mod12 <- lm(log(Peso) ~ Gestazione + Sesso + N.gravidanze + Lunghezza * Cranio, data = data)
plot(mod12, caption = "Plot mod12")

# Modello 13: log-Level + gestazione quadratica
# Ipotizziamo una crescita esponenziale del feto nelle ultime settimane
mod13 <- update(mod12, . ~ . + I(Gestazione^2))
plot(mod13, caption = "Plot mod13")

# Modello 15: target Lineare + gestazione quadratica
# Manteniamo la scala originale del peso ma permettiamo curvatura sulla gestazione
mod15 <- lm(Peso ~ Sesso + N.gravidanze + Lunghezza * Cranio + I(Gestazione^2), data = data)
plot(mod15, caption = "Plot mod15")

shapiro.test(mod13$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod13$residuals
## W = 0.98042, p-value < 2.2e-16
bptest(mod13)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod13
## BP = 249.73, df = 7, p-value < 2.2e-16
dwtest(mod13)
## 
##  Durbin-Watson test
## 
## data:  mod13
## DW = 1.9424, p-value = 0.07532
## alternative hypothesis: true autocorrelation is greater than 0
shapiro.test(mod9$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod9$residuals
## W = 0.96978, p-value < 2.2e-16
bptest(mod9)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod9
## BP = 140.72, df = 6, p-value < 2.2e-16
dwtest(mod9)
## 
##  Durbin-Watson test
## 
## data:  mod9
## DW = 1.9544, p-value = 0.1275
## alternative hypothesis: true autocorrelation is greater than 0

I test di normalità, omoschedasticità e indipendenza tra il modello lineare mod9 e quello trasformato logaritmico mod13 sono praticamente li stessi. Quindi procedo con un confrondo RMSE per controllare l’errore sulla predizione dei due modelli.

5. Confronto e diagnostica finale

Confrontiamo i modelli più promettenti.

Nota metodologica: non è possibile confrontare direttamente l’AIC/BIC tra modelli con variabile dipendente trasformata (mod12, mod13 usano \(log(Y)\)) e modelli lineari (mod9, mod15 usano \(Y\)), poiché la scala della verosimiglianza cambia. Ci baseremo quindi sull’\(R^2\) corretto e, soprattutto, sull’analisi dei residui.

# Estrazione metriche per modelli comparabili (stessa Y)
get_metrics <- function(model) {
  s <- summary(model)
  c(R2_Adj = s$adj.r.squared, 
    Sigma = s$sigma, 
    AIC = AIC(model))
}

res_linear <- rbind(
  Mod9 = get_metrics(mod9),
  Mod15 = get_metrics(mod15)
)

kable(res_linear, digits = 4, caption = "Confronto modelli lineari (stessa scala Y)")
Confronto modelli lineari (stessa scala Y)
R2_Adj Sigma AIC
Mod9 0.7287 273.6354 35092.36
Mod15 0.7286 273.7259 35094.01 
# --- Modello 9 (Lineare) ---
# Le predizioni sono già in grammi
pred_mod9 <- predict(mod9, newdata = data)
residuals_mod9 <- data$Peso - pred_mod9
rmse_mod9 <- sqrt(mean(residuals_mod9^2))

# --- Modello 13 (Logaritmico) ---
# Le predizioni sono in log(grammi). Dobbiamo retro-trasformarle con exp()
pred_mod13_log <- predict(mod13, newdata = data)
pred_mod13 <- exp(pred_mod13_log) # Back-transformation
residuals_mod13 <- data$Peso - pred_mod13
rmse_mod13 <- sqrt(mean(residuals_mod13^2))

# 2. Creazione Tabella di Confronto
rmse_results <- data.frame(
  Modello = c("Modello 9 (Lineare)", "Modello 13 (Log-Quadratico)"),
  RMSE = c(rmse_mod9, rmse_mod13),
  Tipo_Target = c("Peso (g)", "log(Peso)")
)

kable(rmse_results, digits = 2, caption = "Root mean squared error a confronto")
Root mean squared error a confronto
Modello RMSE Tipo_Target
Modello 9 (Lineare) 273.25 Peso (g)
Modello 13 (Log-Quadratico) 269.26 log(Peso) 
# Calcolo differenza percentuale
diff_pct <- (rmse_mod13 - rmse_mod9) / rmse_mod9 * 100
if (diff_pct < 0) {
  print(paste("Il Modello 13 riduce l'errore del", round(abs(diff_pct), 2), "% rispetto al Modello 9."))
} else {
  print(paste("Il Modello 9 ha un errore inferiore del", round(diff_pct, 2), "% rispetto al Modello 13."))
}
## [1] "Il Modello 13 riduce l'errore del 1.46 % rispetto al Modello 9."
# 3. Visualizzazione Grafica (barplot RMSE)
p1 <- ggplot(rmse_results, aes(x = Modello, y = RMSE, fill = Modello)) +
  geom_col(width = 0.5, alpha = 0.8) +
  geom_text(aes(label = round(RMSE, 2)), vjust = -0.5, fontface = "bold") +
  labs(title = "Confronto errore di predizione (RMSE)",
       subtitle = "Valori espressi in grammi (scala originale)",
       y = "RMSE (g)",
       x = "") +
  theme_minimal() +
  scale_fill_manual(values = c("steelblue", "orange")) +
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(min(rmse_results$RMSE) * 0.9, max(rmse_results$RMSE) * 1.1)) # Zoom sull'asse Y

print(p1)

# 4. Visualizzazione distribuzione errori assoluti
# Vediamo dove sbagliano i modelli
errors_df <- data.frame(
  Errore_Assoluto = c(abs(residuals_mod9), abs(residuals_mod13)),
  Modello = factor(c(rep("Mod9", length(residuals_mod9)), rep("Mod13", length(residuals_mod13))))
)

p2 <- ggplot(errors_df, aes(x = Errore_Assoluto, fill = Modello)) +
  geom_density(alpha = 0.4) +
  labs(title = "Distribuzione degli errori assoluti di predizione",
       subtitle = "Confronto densità degli scostamenti (valore reale - predetto)",
       x = "Errore assoluto (g)",
       y = "Densità") +
  theme_minimal() +
  scale_fill_manual(values = c("steelblue", "orange")) +
  xlim(0, 1000) # Limitiamo l'asse X per vedere meglio la parte centrale

print(p2)

Per il confronto tra “Modello 9” e “Modello 13” (trasformato log(Y)) usiamo la metrica RMSE per valutare il potere predittivo. Notiamo come il Modello 13 riduce l’errore del 1.46% rispetto al Modello 9. Perciò possiamo

Il “Modello 9” e il “Modello 15” mostrano performance pressoché identiche. Dai grafici si nota un miglioramento in termini di residui più sparpagliati che suggeriscono una varianza non costante ma comunque dal test di Breusch-Pagan si rileva un \(p-value < 2.2e-16\) ciò indica la presenza di eteroschedasticità dei residui. Mentre dal QQ-Plot notiamo un miglioramento nelle code ma non sufficiente per preferire un modello trasformato su scala logaritmica rispetto ad un modello lineare. Per questo proseguimamo con il modello creato mod9.

5.1 Diagnostica dei residui (Modello 9)

Verifichiamo le assunzioni di validità del modello scelto.

# Test di normalità (Shapiro-Wilk)
# Eseguiamo su un campione casuale se N > 5000 per limiti computazionali di shapiro.test
res <- resid(mod15)
if(length(res) > 5000) res <- sample(res, 5000)
shapiro_res <- shapiro.test(res)

# Test di Omoschedasticità (Breusch-Pagan)
bp_res <- bptest(mod9)

# Test di Autocorrelazione (Durbin-Watson)
dw_res <- dwtest(mod9)

# Tabella risultati diagnostici
diag_table <- data.frame(
  Test = c("Shapiro-Wilk (normalità)", "Breusch-Pagan (omoschedasticità)", "Durbin-Watson (indipendenza)"),
  P_Value = c(format.pval(shapiro_res$p.value, digits=3), 
              format.pval(bp_res$p.value, digits=3), 
              format.pval(dw_res$p.value, digits=3))
)
kable(diag_table, caption = "Diagnostica dei residui (Mod9)")
Diagnostica dei residui (Mod9)
Test P_Value
Shapiro-Wilk (normalità) <2e-16
Breusch-Pagan (omoschedasticità) <2e-16
Durbin-Watson (indipendenza) 0.128 
# Visualizzazione grafica
par(mfrow=c(2,2))
plot(mod9)

par(mfrow=c(1,1))

Analisi dei residui: Sebbene il test di Shapiro-Wilk indichi una deviazione dalla normalità (comune in dataset numerosi con outliers), il Q-Q plot mostra che i residui seguono la linea teorica per la maggior parte della distribuzione, deviando solo sulle code. Il test di Breusch-Pagan e Durbin-Watson confermano la presenza di eteroschedasticità e assenza di autocorrelazione grave. Possiamo considerare il modello non proprio valido per l’inferenza dato l’andamento constante della varianza dei residui, ma nella cella successiva incroceremo queste considerazioni con i risultati dei test d’inferenza sui regressori.

6. Inferenza e conclusioni (Modello 9)

Presentiamo le stime finali dei coefficienti per il modello selezionato.

summary_9 <- summary(mod9)
coefs <- summary_9$coefficients
conf_int <- confint(mod9, level=0.95)

# Uniamo stime e intervalli di confidenza
final_results <- cbind(coefs[,1:2], conf_int, P_Value = coefs[,4])
kable(final_results, digits = 4, caption = "Stime dei Parametri (Modello 9)")
Stime dei Parametri (Modello 9)
Estimate Std. Error 2.5 % 97.5 % P_Value
(Intercept) -1803.6673 1018.7773 -3801.4060 194.0714 0.0768
Gestazione 38.1508 3.9699 30.3661 45.9355 0.0000
SessoM 73.2227 11.2165 51.2280 95.2173 0.0000
N.gravidanze 12.9344 4.3272 4.4491 21.4197 0.0028
Lunghezza -0.3049 2.2043 -4.6274 4.0175 0.8900
Cranio -4.7520 3.1943 -11.0158 1.5117 0.1370
Lunghezza:Cranio 0.0316 0.0065 0.0188 0.0444 0.0000

6. Inferenza e conclusioni (Modello 13)

Presentiamo le stime finali dei coefficienti per il modello selezionato.

summary_13 <- summary(mod13)
coefs <- summary_13$coefficients
conf_int <- confint(mod13, level=0.95)

# Uniamo stime e intervalli di confidenza
final_results <- cbind(coefs[,1:2], conf_int, P_Value = coefs[,4])
kable(final_results, digits = 4, caption = "Stime dei Parametri (Modello 9)")
Stime dei Parametri (Modello 9)
Estimate Std. Error 2.5 % 97.5 % P_Value
(Intercept) 0.1882 0.3248 -0.4487 0.8251 0.5623
Gestazione 0.1286 0.0204 0.0885 0.1686 0.0000
SessoM 0.0224 0.0034 0.0157 0.0291 0.0000
N.gravidanze 0.0037 0.0013 0.0011 0.0063 0.0052
Lunghezza 0.0082 0.0009 0.0064 0.0099 0.0000
Cranio 0.0103 0.0013 0.0078 0.0129 0.0000
I(Gestazione^2) -0.0015 0.0003 -0.0020 -0.0010 0.0000
Lunghezza:Cranio 0.0000 0.0000 0.0000 0.0000 0.0000

Interpretazione dei risultati mod9 Gestazione: la variabile è altamente significativa (\(p < 0.001\)). Il coefficiente positivo (~41.3) indica che, mantenendo costanti le altre variabili, ogni settimana aggiuntiva di gestazione è associata a un aumento medio di peso di circa 41 grammi.

Sesso: essere maschio è associato a un aumento di peso significativo (~76.7 grammi) rispetto alle femmine, a parità di dimensioni e gestazione.

Il “paradosso” delle gravidanze: risultato cruciale di questa analisi è l’effetto della variabile N.gravidanze. Nelle analisi preliminari (T-Test univariato), non emergeva differenza significativa di peso tra primipare e pluripare. Tuttavia, nel modello multivariato, il coefficiente è positivo (~10.5) e significativo (\(p < 0.05\)). Questo fenomeno si spiega con l’effetto di “mascheramento”: l’effetto positivo dell’esperienza uterina (precedenti gravidanze) era nascosto dalla variabilità di altri fattori (come la lunghezza o l’età gestazionale).

Controllando per questi fattori nel modello di regressione, l’effetto biologico netto emerge: a parità di struttura fisica e durata della gestazione, i figli di madri pluripare tendono a pesare leggermente di più.

Interazione biometrica: il termine “Lunghezza:Cranio” è significativo, confermando che il peso non è una semplice somma di lunghezze, ma è meglio rappresentato da un volume corporeo approssimato dall’interazione di queste misure, anche se la stima è molto bassa (pari a 0).

Conclusione operativa: il modello mod9 rappresenta uno strumento affidabile per stimare il peso alla nascita basandosi su parametri ecografici (“Cranio”, “Lunghezza”) e clinici (“Gestazione”,“Sesso”,“Parità”), offrendo un’interpretazione chiara dei meccanismi biologici sottostanti.

Interpretazione dei risultati mod13 (scala log) Per quanto riguarda il Modello 13, osserviamo un andamento dei residui simile al modello lineare mod9 ciò significa che il modello log non peggiora le assunzioni.

Ma confrontando l’RMSE osserviamo come si riduce l’errore sulle previsioni in maniera significativa. Inoltre i coefficienti di regressione, dopo il test d’inferenza, sono risultati tutti molto significativi. Giustificando l’utilizzo del modello mod13 anzichè il mod9 lineare, spiega meglio il modello trasformato con log(Y) e questo ha senso dal momento in cui la crescita del feto e il peso hanno andamento curvo da quanto riscontrato a livello visuale dai plot.

Conclusione finale: tra i due modelli, il modello con variabile dipendente trasformata logaritmicamente (mod13) mostra una migliore capacità predittiva, come indicato da un RMSE inferiore calcolato sulla scala originale di Y. Inoltre, il comportamento dei residui è comparabile a quello del modello lineare (mod9). Pertanto, mod13 viene preferito.