knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, fig.align = "center")

# Installazione e caricamento pacchetti in un unico blocco per pulizia
required_packages <- c("dplyr", "moments", "ggplot2", "scales", "car", "lmtest", "MASS")
new_packages <- required_packages[!(required_packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(moments)
library(ggplot2)
library(scales)
library(car)    # VIF
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(lmtest) # Diagnostica residui (Breusch-Pagan, Durbin-Watson)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(MASS)   # Box-Cox
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select

1. Introduzione e Preparazione Dati

In questo documento analizziamo il dataset relativo alle nascite. L’obiettivo è esplorare le variabili che influenzano il peso del neonato e costruire un modello predittivo performante, confrontando approcci lineari classici con trasformazioni logaritmiche e polinomiali.

# Caricamento dati
# Assicurati che il file neonati.csv sia nella working directory
setwd("~/Desktop/Università/Master Data Science/Progetto 3 - Statistica inferenziale")
data <- read.csv("neonati.csv", sep=",")

# Rimozione valori mancanti (CRUCIALE per confronto modelli e ANOVA)
data <- na.omit(data)

# Anteprima rapida
head(data)
##   Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1         26            0         0         42 3380       490    325        Nat
## 2         21            2         0         39 3150       490    345        Nat
## 3         34            3         0         38 3640       500    375        Nat
## 4         28            1         0         41 3690       515    365        Nat
## 5         20            0         0         38 3700       480    335        Nat
## 6         32            0         0         40 3200       495    340        Nat
##   Ospedale Sesso
## 1     osp3     M
## 2     osp1     F
## 3     osp2     M
## 4     osp2     M
## 5     osp3     F
## 6     osp2     F
# Preprocessing: Conversione in Fattori
# Fondamentale per trattare correttamente le variabili categoriche nelle analisi
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)

2. Analisi Descrittiva (EDA)

2.1 Statistiche Variabili Numeriche

Calcoliamo media, deviazione standard, coefficiente di variazione (CV) e indici di forma.

stats_summary <- data.frame()

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

knitr::kable(stats_summary, digits = 2, caption = "Statistiche Descrittive")
Statistiche Descrittive
Variabile Media SD CV_Percent Skewness Kurtosis
Anni.madre 28.16 5.27 18.72 0.04 3.38
N.gravidanze 0.98 1.28 130.51 2.51 13.99
Gestazione 38.98 1.87 4.79 -2.07 11.26
Peso 3284.08 525.04 15.99 -0.65 5.03
Lunghezza 494.69 26.32 5.32 -1.51 9.49
Cranio 340.03 16.43 4.83 -0.79 5.95

2.2 Distribuzione Variabili Continue

Visualizziamo la forma delle distribuzioni (Istogrammi + Densità).

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)
}

2.3 Analisi degli Outliers

Utilizziamo il criterio \(1.5 \times IQR\) per identificare i valori anomali.

df_outliers <- data.frame()

for (var in colnames(data)) {
  colonna_corrente <- data[[var]]
  
  if (is.numeric(colonna_corrente)) {
    # Calcolo soglie
    Q1 <- quantile(colonna_corrente, 0.25, na.rm = TRUE)
    Q3 <- quantile(colonna_corrente, 0.75, na.rm = TRUE)
    IQR_val <- Q3 - Q1
    lower_bound <- Q1 - 1.5 * IQR_val
    upper_bound <- Q3 + 1.5 * IQR_val
    
    # Identificazione indici
    idx_outliers <- which(colonna_corrente < lower_bound | colonna_corrente > upper_bound)
    val_outliers <- colonna_corrente[idx_outliers]
    
    # Salvataggio
    if (length(val_outliers) > 0) {
      temp_df <- data.frame(
        Variabile = var,
        Riga_Indice = idx_outliers,
        Valore = val_outliers
      )
      df_outliers <- rbind(df_outliers, temp_df)
    }
    
    # Plot
    boxplot(colonna_corrente, main = paste("Boxplot Outliers:", var), 
            ylab = var, col = "lightblue", outpch = 19, outcol = "red")
    abline(h = c(lower_bound, upper_bound), col = "red", lty = 2)
  }
}

if (nrow(df_outliers) > 0) {
  print(paste("Totale osservazioni classificate come outlier:", nrow(df_outliers)))
  print(table(df_outliers$Variabile))
}
## [1] "Totale osservazioni classificate come outlier: 502"
## 
##   Anni.madre       Cranio   Gestazione    Lunghezza N.gravidanze         Peso 
##           13           48           67           59          246           69

Nota: Quasi il 20% del campione presenta outlier. Sarà necessario monitorare i residui dei modelli per capire se questi valori influenzano negativamente la regressione. Gli outliers tendono ad essere asimettrici negativi, quindi con frequenze maggiori a valori alti.

2.4 Relazione Variabili Discrete vs Peso

Visualizziamo come il peso varia in base alle categorie qualitative.

vars_discrete <- c("N.gravidanze", "Fumatrici", "Tipo.parto", "Ospedale", "Sesso")

for (var in vars_discrete) {
  if (var %in% colnames(data)) {
    p <- ggplot(data, aes_string(x = paste0("as.factor(", var, ")"), 
                                 y = "Peso", 
                                 fill = paste0("as.factor(", var, ")"))) +
      geom_boxplot(alpha = 0.6, outlier.colour = "red") +
      labs(title = paste("Distribuzione Peso per:", var),
           x = var, y = "Peso (g)", fill = var) +
      theme_minimal() +
      theme(legend.position = "none")
    print(p)
  }
}

3. Analisi di Correlazione

3.1 Matrice di Scatterplot e Correlazioni (Pairs Plot)

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))

3.2 Heatmap delle Correlazioni

cor_matrix <- cor(data_numeric, method = "pearson")
melted_cormat <- as.data.frame(as.table(cor_matrix))
names(melted_cormat) <- c("Var1", "Var2", "Correlazione")

ggplot(melted_cormat, aes(x = Var1, y = Var2, fill = Correlazione)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(-1,1), name="Pearson\nCorr") +
  geom_text(aes(label = round(Correlazione, 2)), color = "black", size = 3) +
  theme_minimal() + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  coord_fixed() +
  labs(title = "Matrice di Correlazione Completa")

3.3 Focus: Età Madre vs Peso Neonato

var_x <- "Anni.madre"
var_y <- "Peso"
cor_test <- cor.test(data[[var_x]], data[[var_y]], method = "spearman")

print(cor_test)
## 
##  Spearman's rank correlation rho
## 
## data:  data[[var_x]] and data[[var_y]]
## S = 2619756994, p-value = 0.7648
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##          rho 
## -0.005986847
ggplot(data, aes_string(x = var_x, y = var_y)) +
  geom_point(alpha = 0.6, color = "darkblue", size = 2) +
  geom_smooth(method = "loess", color = "red", se = TRUE) +
  labs(title = "Relazione Età Madre vs Peso",
       subtitle = paste("Corr (Spearman):", round(cor_test$estimate, 3), 
                        "| P-value:", format.pval(cor_test$p.value, digits=3))) +
  theme_minimal()

Interpretazione: La correlazione non è statisticamente significativa (p > 0.05), suggerendo che l’età della madre non ha una relazione lineare diretta evidente con il peso del nascituro. Si usa in questo caso spearman dato che i dati sono particolarmente asimettrici.

4. Statistica Inferenziale (Test di Ipotesi)

4.1 T-Test e Chi-Quadro

Eseguiamo i test per verificare le differenze tra i gruppi prima della modellazione.

# 1. T-Test: Peso ~ Fumatrici
print("--- T-Test: Peso ~ Fumatrici ---")
## [1] "--- T-Test: Peso ~ Fumatrici ---"
print(t.test(Peso ~ Fumatrici, data = data, var.equal = FALSE))
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Fumatrici
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -45.61354 145.22674
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.153        3236.346
# Conclusione: Non rifiutiamo H0. Fumo non significativo univariatamente.

# 2. T-Test: Peso ~ Sesso
print("--- T-Test: Peso ~ Sesso ---")
## [1] "--- T-Test: Peso ~ Sesso ---"
print(t.test(Peso ~ Sesso, data = data, var.equal = FALSE))
## 
##  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
# Conclusione: Rifiutiamo H0. Sesso significativo.

# 3. Chi-Quadro: Ospedale vs Tipo di Parto
print("--- Chi-Quadro: Ospedale vs Tipo Parto ---")
## [1] "--- Chi-Quadro: Ospedale vs Tipo Parto ---"
print(chisq.test(table(data$Ospedale, data$Tipo.parto)))
## 
##  Pearson's Chi-squared test
## 
## data:  table(data$Ospedale, data$Tipo.parto)
## X-squared = 1.0972, df = 2, p-value = 0.5778
# Conclusione: Non rifiutiamo H0. Indipendenza tra ospedale e tipo parto.

4.2 Focus Specifico: Numero di Gravidanze

Approfondiamo se essere al primo parto (Primipara, 0) o ai successivi (Pluripara, 1+) influenza il peso.

data <- data %>%
  mutate(Gruppo_Gravidanze = as.factor(ifelse(N.gravidanze == 0, "Primipara", "Pluripara")))

t_test_grav <- t.test(Peso ~ Gruppo_Gravidanze, data = data)
print(t_test_grav)
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Gruppo_Gravidanze
## t = 0.43065, df = 2375.8, p-value = 0.6668
## alternative hypothesis: true difference in means between group Pluripara and group Primipara is not equal to 0
## 95 percent confidence interval:
##  -32.30465  50.48640
## sample estimates:
## mean in group Pluripara mean in group Primipara 
##                3288.066                3278.975
ggplot(data, aes(x = Gruppo_Gravidanze, y = Peso, fill = Gruppo_Gravidanze)) +
  geom_boxplot(alpha = 0.6) +
  labs(title = "Peso Neonato: Primipare vs Pluripare",
       subtitle = paste("P-value T-test:", format.pval(t_test_grav$p.value, digits=3))) +
  theme_minimal()

Conclusione: Il P-value > 0.05 indica che, a livello univariato, non c’è differenza significativa di peso medio tra chi è al primo parto e chi ha già partorito.

5. Modellazione (Regressione Lineare)

Procediamo con la costruzione dei modelli, partendo da quello completo e rimuovendo via via le variabili non significative (Stepwise manuale).

5.1 Selezione Variabili (Modelli Lineari Standard)

# Modello 1: Tutte le variabili (Baseline)
mod1 <- lm(Peso ~ ., data = data)

# Modello 2: Rimozione Fumatrici (non significativa)
mod2 <- update(mod1, . ~ . - Fumatrici)

# Modello 3: Reset su data pulita
mod3 <- lm(Peso ~ ., data = data)

# Modello 4: Rimozione Fumatrici e Ospedale
mod4 <- update(mod3, . ~ . - Fumatrici - Ospedale)

# Modello 5: Rimozione Lunghezza (Test feature importance)
mod5 <- update(mod4, . ~ . - Lunghezza)

# Modello 6: Reintroduzione Lunghezza + Interazione (Cranio*Lunghezza)
mod6 <- update(mod5, . ~ . + Lunghezza + Cranio:Lunghezza)

# Modelli successivi: Rimozione variabili meno significative
mod7 <- update(mod6, . ~ . - Lunghezza - Cranio) 
mod8 <- update(mod7, . ~ . - Tipo.parto)
mod9 <- update(mod8, . ~ . - Anni.madre) # Candidato "Simple & Robust"

# Modello 10: Reintroduzione Lunghezza per verifica stabilità
mod10 <- update(mod9, . ~ . + Lunghezza)

# Modello 11: Rimozione N.gravidanze
mod11 <- update(mod10, . ~ . - N.gravidanze)

# Altri Modelli Test
mod5_rev = update(mod5, ~ . - N.gravidanze - Tipo.parto - Anni.madre)

5.2 Modellazione Avanzata (Trasformazioni e GLM)

Esploriamo trasformazioni della variabile target e termini non lineari.

# --- MODELLO 12: Log-Level Regression ---
# Trasformazione log(Peso) per stabilizzare varianza
mod12 <- lm(log(Peso) ~ . - Fumatrici - N.gravidanze - Ospedale + Lunghezza*Cranio, data = data)

# --- MODELLO 13: Log + Gestazione Quadrata ---
# Aggiungiamo I(Gestazione^2) per catturare l'accelerazione della crescita fetale
mod13 <- update(mod12, . ~ . + I(Gestazione^2))

# --- ANALISI BOX-COX ---
# Verifica lambda ottimale
mod_temp <- lm(Peso ~ . - Fumatrici - N.gravidanze - Ospedale + Lunghezza*Cranio + I(Gestazione^2), data = data)
bc <- boxcox(mod_temp, plotit = FALSE, lambda = seq(-2, 2, by = 0.1))
lambda_opt <- bc$x[which.max(bc$y)]
cat("Lambda Ottimale Box-Cox:", round(lambda_opt, 4), "\n")
## Lambda Ottimale Box-Cox: 0.3
# --- MODELLO 14: Box-Cox ---
if(abs(lambda_opt) > 0.1) {
  data$Peso_BoxCox <- (data$Peso^lambda_opt - 1) / lambda_opt
  mod14 <- lm(Peso_BoxCox ~ . - Fumatrici - N.gravidanze - Ospedale + Lunghezza*Cranio + I(Gestazione^2) - Peso, data = data)
}

# --- MODELLO 15: Lineare + Gestazione Quadrata ---
mod15 <- lm(Peso ~ . - Fumatrici - N.gravidanze - Ospedale + Lunghezza*Cranio + I(Gestazione^2), data = data)

# --- NUOVI MODELLI (Log Predittori) ---
# Modello 999: Mod9 ma con log(Gestazione)
mod999 <- update(mod9, . ~ . - Gestazione + log(Gestazione))

# Modello 1000: Mod999 esteso (Log target e Log gestazione)
mod1000 <- lm(log(Peso) ~ . - N.gravidanze - Tipo.parto - Anni.madre - Lunghezza - Fumatrici - Ospedale - Gestazione + log(Gestazione), data=data)

6. Confronto Globale e Diagnostica

6.1 Classifica Prestazioni (Tutti i Modelli)

Confrontiamo tutti i modelli generati in termini di \(R^2\) Adjusted, AIC, BIC e diagnostica dei residui.

# Recupero lista modelli automatici
auto_models <- ls(pattern = "^mod[0-9]+$")
manual_models <- c("mod5_rev", "mod999", "mod1000")
all_models_names <- unique(c(auto_models, manual_models))

global_results <- data.frame()

for (m_name in all_models_names) {
  if(exists(m_name)) {
    curr_mod <- get(m_name)
    if(inherits(curr_mod, "lm")) {
      summ <- summary(curr_mod)
      
      # Test Residui (su campione se N > 5000)
      res <- resid(curr_mod)
      if(length(res) > 5000) { set.seed(123); res <- sample(res, 5000) }
      
      # P-Values Test Diagnostici (con gestione errori)
      shapiro_p <- tryCatch(shapiro.test(res)$p.value, error=function(e) NA)
      bp_val <- tryCatch(bptest(curr_mod)$p.value, error=function(e) NA)
      dw_val <- tryCatch(dwtest(curr_mod)$p.value, error=function(e) NA)
      
      # Identifica Target
      target_var <- as.character(formula(curr_mod)[[2]])
      if(length(target_var) > 1) target_var <- paste(target_var, collapse=" ") 
      
      temp <- data.frame(
        Modello = m_name, Target = target_var,
        Adj_R2 = summ$adj.r.squared,
        AIC = AIC(curr_mod), BIC = BIC(curr_mod),
        F_Stat = summ$fstatistic[1],
        Shapiro_P = shapiro_p, BP_P = bp_val, DW_P = dw_val
      )
      global_results <- rbind(global_results, temp)
    }
  }
}

# Ordinamento e Formattazione
if(nrow(global_results) > 0) {
  global_results <- global_results[order(global_results$Adj_R2, decreasing = TRUE), ]
  
  global_results_display <- global_results
  fmt_p <- function(x) format.pval(x, digits = 3, eps = 0.001, na.form = "-")
  
  global_results_display$Shapiro_P <- fmt_p(global_results$Shapiro_P)
  global_results_display$BP_P <- fmt_p(global_results$BP_P)
  global_results_display$DW_P <- fmt_p(global_results$DW_P)
  
  print(knitr::kable(global_results_display, digits = 4, row.names = FALSE, 
                     caption = "Performance Globali Modelli"))
}
## 
## 
## Table: Performance Globali Modelli
## 
## |Modello  |Target      | Adj_R2|        AIC|        BIC|      F_Stat|Shapiro_P |BP_P   |DW_P   |
## |:--------|:-----------|------:|----------:|----------:|-----------:|:---------|:------|:------|
## |mod1000  |log Peso    | 0.9959| -15152.141| -15111.373| 122793.1819|<0.001    |<0.001 |0.6224 |
## |mod15    |Peso        | 0.9957|  24782.409|  24852.298|  58328.7523|<0.001    |<0.001 |<0.001 |
## |mod13    |log Peso    | 0.7915|  -5297.276|  -5233.212|   1054.9819|<0.001    |<0.001 |0.0674 |
## |mod12    |log Peso    | 0.7891|  -5269.775|  -5211.535|   1169.7498|<0.001    |<0.001 |0.0600 |
## |mod14    |Peso_BoxCox | 0.7711|   6809.709|   6873.773|    936.1481|<0.001    |<0.001 |0.0830 |
## |mod6     |Peso        | 0.7293|  35156.925|  35220.990|    749.1705|<0.001    |<0.001 |0.1128 |
## |mod10    |Peso        | 0.7287|  35159.623|  35206.216|   1119.7590|<0.001    |<0.001 |0.1224 |
## |mod11    |Peso        | 0.7285|  35160.443|  35201.212|   1342.1720|<0.001    |<0.001 |0.1274 |
## |mod1     |Peso        | 0.7277|  35173.678|  35249.390|    608.1984|<0.001    |<0.001 |0.1170 |
## |mod3     |Peso        | 0.7277|  35173.678|  35249.390|    608.1984|<0.001    |<0.001 |0.1170 |
## |mod2     |Peso        | 0.7277|  35172.931|  35242.819|    668.8271|<0.001    |<0.001 |0.1126 |
## |mod4     |Peso        | 0.7269|  35178.137|  35236.377|    832.4856|<0.001    |<0.001 |0.1142 |
## |mod7     |Peso        | 0.7255|  35190.295|  35242.712|    944.4085|<0.001    |0.295  |0.1080 |
## |mod9     |Peso        | 0.7252|  35191.065|  35231.834|   1319.7599|<0.001    |0.499  |0.1164 |
## |mod8     |Peso        | 0.7251|  35192.670|  35239.262|   1099.5984|<0.001    |0.218  |0.1119 |
## |mod999   |Peso        | 0.7251|  35191.793|  35232.561|   1319.2310|<0.001    |0.467  |0.1200 |
## |mod5_rev |Peso        | 0.5992|  36133.564|  36168.508|    934.8531|<0.001    |0.128  |0.0234 |
## |mod5     |Peso        | 0.5988|  36138.572|  36190.988|    533.9125|<0.001    |0.144  |0.0215 |

6.2 Analisi VIF (Multicollinearità)

vif_data_list <- list()
for (m_name in sort(all_models_names)) {
  if (exists(m_name)) {
    model_obj <- get(m_name)
    if (inherits(model_obj, "lm")) {
      tryCatch({
        vif_res <- vif(model_obj)
        vals <- if (is.matrix(vif_res)) vif_res[, 1] else vif_res
        temp_list <- as.list(vals)
        temp_list$Modello <- m_name
        vif_data_list[[m_name]] <- temp_list
      }, error = function(e) {})
    }
  }
}

if(length(vif_data_list) > 0) {
  vif_table <- bind_rows(vif_data_list) %>% dplyr::select(Modello, everything())
  vif_table <- vif_table %>% mutate(across(everything(), as.character))
  vif_table[is.na(vif_table)] <- "-"
  print(knitr::kable(vif_table, caption = "Livelli di VIF per Variabile"))
}
## 
## 
## Table: Livelli di VIF per Variabile
## 
## |Modello  |Anni.madre       |N.gravidanze     |Fumatrici        |Gestazione       |Lunghezza        |Cranio           |Tipo.parto       |Ospedale         |Sesso            |Gruppo_Gravidanze |Lunghezza:Cranio |Peso_BoxCox      |log(Gestazione)  |I(Gestazione^2)  |Cranio:Lunghezza |
## |:--------|:----------------|:----------------|:----------------|:----------------|:----------------|:----------------|:----------------|:----------------|:----------------|:-----------------|:----------------|:----------------|:----------------|:----------------|:----------------|
## |mod1     |1.1968201538097  |1.95582453346558 |1.00933686150721 |1.70522390768659 |2.0857584335596  |1.63906792154197 |1.0045793702441  |1.00485881581363 |1.04070540206104 |1.88667406791593  |-                |-                |-                |-                |-                |
## |mod10    |-                |1.85340935779114 |-                |1.66084433899482 |5.82461147434374 |-                |-                |-                |1.04107405815535 |1.86406185487713  |5.57242325807231 |-                |-                |-                |-                |
## |mod1000  |-                |-                |-                |-                |-                |2.06891280773069 |-                |-                |1.0539563477396  |1.02848646951495  |-                |2.82599083806204 |1.75839682559755 |-                |-                |
## |mod11    |-                |-                |-                |1.6603678017212  |5.81472521896731 |-                |-                |-                |1.04018284153286 |1.02850740974693  |5.56680403535622 |-                |-                |-                |-                |
## |mod12    |1.13415730702627 |-                |-                |1.86792257450727 |112.732196800078 |92.29337395777   |1.00426397570473 |-                |1.0474708556533  |1.14234398497377  |314.885580676347 |-                |-                |-                |-                |
## |mod13    |1.13464296664542 |-                |-                |519.447570346048 |201.328091704502 |162.919910630903 |1.00426588599791 |-                |1.04912321478107 |1.14489876958434  |551.449429174834 |-                |-                |493.321640208192 |-                |
## |mod14    |1.13464296664547 |-                |-                |519.447570345897 |201.328091704475 |162.919910630886 |1.00426588599791 |-                |1.04912321478106 |1.14489876958428  |551.449429174748 |-                |-                |493.321640208045 |-                |
## |mod15    |1.13495123979753 |-                |-                |523.51444886488  |203.826819243249 |164.211817583939 |1.00629739472524 |-                |1.06750480909462 |1.14714466424463  |552.248192952987 |4.38366767832181 |-                |496.116636034339 |-                |
## |mod2     |1.19668596797496 |1.95555320019302 |-                |1.69823984366408 |2.08192212730182 |1.6388135237447  |1.00417970761843 |1.00443359439492 |1.04047231364652 |1.88303793549059  |-                |-                |-                |-                |-                |
## |mod3     |1.1968201538097  |1.95582453346558 |1.00933686150721 |1.70522390768659 |2.0857584335596  |1.63906792154197 |1.0045793702441  |1.00485881581363 |1.04070540206104 |1.88667406791593  |-                |-                |-                |-                |-                |
## |mod4     |1.19604481566846 |1.95455392409214 |-                |1.69693305264703 |2.08004067893252 |1.63852253342956 |1.00372757631702 |-                |1.04019433488858 |1.88165140158832  |-                |-                |-                |-                |-                |
## |mod5     |1.19604481555786 |1.95195939804381 |-                |1.32924622676014 |-                |1.30547395790098 |1.00114490896784 |-                |1.02870213071338 |1.88164980401082  |-                |-                |-                |-                |-                |
## |mod5_rev |-                |-                |-                |1.30982216532279 |-                |1.30188225663671 |-                |-                |1.02789997739555 |1.02961563638688  |-                |-                |-                |-                |-                |
## |mod6     |1.19607778768903 |1.95544419836685 |-                |1.86817469757152 |112.808960046059 |92.3256269268869 |1.00434770517865 |-                |1.04811608407097 |1.88827610763836  |-                |-                |-                |-                |315.029007024637 |
## |mod7     |1.19529418055755 |1.95241393285817 |-                |1.59141588235532 |-                |-                |1.00120061848155 |-                |1.04005775452098 |1.87732965209732  |-                |-                |-                |-                |1.58535417797033 |
## |mod8     |1.19527476361522 |1.9523193302883  |-                |1.59136863677715 |-                |-                |-                |-                |1.04005775248542 |1.87690167739332  |-                |-                |-                |-                |1.58503406904003 |
## |mod9     |-                |1.85026352080116 |-                |1.57326334943691 |-                |-                |-                |-                |1.03998615519062 |1.86167865038122  |-                |-                |-                |-                |1.58220975631089 |
## |mod999   |-                |1.85078114373214 |-                |-                |-                |-                |-                |-                |1.03967241387812 |1.85903333522014  |-                |-                |1.59379288422143 |-                |1.6062966801985  |

7. Approfondimenti e Conclusioni Finali

Dopo aver analizzato vari modelli, il Modello 9 e il Modello 13 emergono come i più interessanti.

7.1 Confronto Modello 9 vs Modello 1000

Verifichiamo se la trasformazione logaritmica della gestazione (mod1000) è preferibile alla lineare (mod9).

if(exists("mod9") && exists("mod1000")) {
  # Confronto Metriche
  comp_df <- data.frame(
    Metric = c("RSS (Devianza)", "AIC", "Adj. R2"),
    Mod9_Lineare = c(deviance(mod9), AIC(mod9), summary(mod9)$adj.r.squared),
    Mod1000_Log = c(deviance(mod1000), AIC(mod1000), summary(mod1000)$adj.r.squared)
  )
  print(knitr::kable(comp_df, digits = 4, caption = "Confronto Mod9 vs Mod1000"))
  
  # Logica di Decisione
  diff_aic <- AIC(mod1000) - AIC(mod9)
  if (diff_aic > 2) {
    print("Vince MOD9 (AIC più basso). La relazione lineare fitta meglio i dati e ha più senso biologico.")
  } else {
    print("I modelli sono simili, si preferisce Mod9 per semplicità interpretativa.")
  }
  
  # Confronto Grafico Residui
  plot(density(resid(mod9)), col="blue", lwd=2, main="Confronto Residui: Mod9 vs Mod1000")
  lines(density(resid(mod1000)), col="green", lwd=2, lty=2)
  legend("topright", legend=c("Mod9 (Lineare)", "Mod1000 (Log)"), col=c("blue", "green"), lty=c(1, 2))
}
## 
## 
## Table: Confronto Mod9 vs Mod1000
## 
## |Metric         | Mod9_Lineare| Mod1000_Log|
## |:--------------|------------:|-----------:|
## |RSS (Devianza) | 1.889504e+08|      0.3395|
## |AIC            | 3.519107e+04| -15152.1411|
## |Adj. R2        | 7.252000e-01|      0.9959|
## [1] "I modelli sono simili, si preferisce Mod9 per semplicità interpretativa."

7.2 Inferenza Finale (Modello 9)

Il modello finale scelto è il Modello 9, che bilancia accuratezza e interpretabilità. Ecco i risultati inferenziali definitivi.

if(exists("mod9")) {
  summary_mod9 <- summary(mod9)
  print(summary_mod9$coefficients)
  
  # Intervalli di Confidenza
  conf_intervals <- confint(mod9, level = 0.95)
  
  # Forest Plot
  coef_data <- as.data.frame(summary_mod9$coefficients)
  coef_data$Term <- rownames(coef_data)
  coef_data$CI_Low <- conf_intervals[,1]
  coef_data$CI_High <- conf_intervals[,2]
  colnames(coef_data)[1] <- "Estimate"
  coef_data_plot <- coef_data[coef_data$Term != "(Intercept)", ]
  
  ggplot(coef_data_plot, aes(x = Term, y = Estimate)) +
    geom_point(color = "blue", size = 3) +
    geom_errorbar(aes(ymin = CI_Low, ymax = CI_High), width = 0.2, color = "darkblue") +
    geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
    coord_flip() +
    labs(title = "Intervalli di Confidenza (Modello 9)",
         subtitle = "Variabili significative: intervallo non tocca lo zero",
         y = "Stima Coefficiente", x = "Variabile") +
    theme_minimal()
}
##                                 Estimate   Std. Error     t value      Pr(>|t|)
## (Intercept)                -2.781147e+03 1.167771e+02 -23.8158625 3.913020e-113
## N.gravidanze                8.367012e+00 5.848578e+00   1.4306061  1.526684e-01
## Gestazione                  4.146580e+01 3.695879e+00  11.2194682  1.559309e-28
## SessoM                      7.677906e+01 1.122805e+01   6.8381461  1.005353e-11
## Gruppo_GravidanzePrimipara -8.314177e+00 1.513768e+01  -0.5492373  5.828918e-01
## Cranio:Lunghezza            2.615358e-02 4.667644e-04  56.0316500  0.000000e+00

8. Conclusioni Finali

L’analisi condotta ha permesso di identificare i principali determinanti del peso alla nascita attraverso un processo iterativo di selezione delle variabili e confronto tra diversi modelli di regressione.

1. Selezione del Modello Migliore: Dopo aver testato modelli lineari semplici, trasformazioni logaritmiche e polinomiali, il Modello 9 è stato selezionato come il migliore compromesso.

Sebbene modelli più complessi (come quelli con trasformazione logaritmica o termini quadratici) offrissero un leggero incremento dell’R-quadro o un AIC inferiore, la loro complessità aggiuntiva e la difficoltà di interpretazione non giustificavano il guadagno marginale.

Il Modello 9 spiega circa il 72.5% della varianza del peso (Adj R2 ~ 0.725) utilizzando solo variabili lineari e un’interazione significativa.

2. Variabili Significative: Le analisi inferenziali confermano che il peso del neonato è influenzato positivamente e significativamente da:

Dimensioni Corporee: Lunghezza e Cranio (e la loro interazione).

Età Gestazionale: Ogni settimana in più aumenta il peso atteso, sembra seguire una curva non lineare e un aumento esponenziale tendendo verso destra, quindi avvicinandoci al terzo quartile.

Sesso: I maschi tendono a pesare di più delle femmine a parità di altre condizioni.

Numero di Gravidanze: Un fattore emerso solo nell’analisi multivariata.

3. Il Paradosso delle Gravidanze (Effetto di Mascheramento): Un risultato chiave di questo studio è la discrepanza tra l’analisi univariata (T-Test) e quella multivariata riguardo al numero di gravidanze.

Mentre il T-Test non mostrava differenze significative tra primipare e pluripare, il modello di regressione ha rivelato che, a parità di gestazione e dimensioni fisiche, aver avuto gravidanze precedenti contribuisce positivamente al peso. Questo suggerisce che l’effetto biologico dell’esperienza uterina era “mascherato” da altre variabili di confusione nell’analisi semplice.

4. Validità del Modello: L’analisi dei residui ha mostrato una distribuzione approssimativamente normale, con una leggera deviazione sulle code dovuta alla presenza di outlier (circa il 20% del dataset). Tuttavia, i test diagnostici confermano l’assenza di eteroschedasticità grave e di autocorrelazione, rendendo il modello robusto per scopi inferenziali e predittivi.