predizione_neonati

# 1. IMPOSTAZIONI GLOBALI
# Questo blocco viene eseguito per primo. Imposta le regole per TUTTI i blocchi successivi.
# message = FALSE e warning = FALSE nasconderanno i messaggi tecnici in tutto il documento.
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, fig.align = "center")

# Installazione e caricamento pacchetti in un unico blocco ottimizzato
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
library(lmtest) # Diagnostica residui (Breusch-Pagan, Durbin-Watson)
library(MASS)   # Box-Cox

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) 

summary(data)

##    Anni.madre     N.gravidanze     Fumatrici   Gestazione         Peso     
##  Min.   : 0.00   Min.   : 0.0000   0:2396    Min.   :25.00   Min.   : 830  
##  1st Qu.:25.00   1st Qu.: 0.0000   1: 104    1st Qu.:38.00   1st Qu.:2990  
##  Median :28.00   Median : 1.0000             Median :39.00   Median :3300  
##  Mean   :28.16   Mean   : 0.9812             Mean   :38.98   Mean   :3284  
##  3rd Qu.:32.00   3rd Qu.: 1.0000             3rd Qu.:40.00   3rd Qu.:3620  
##  Max.   :46.00   Max.   :12.0000             Max.   :43.00   Max.   :4930  
##    Lunghezza         Cranio    Tipo.parto Ospedale   Sesso   
##  Min.   :310.0   Min.   :235   Ces: 728   osp1:816   F:1256  
##  1st Qu.:480.0   1st Qu.:330   Nat:1772   osp2:849   M:1244  
##  Median :500.0   Median :340              osp3:835           
##  Mean   :494.7   Mean   :340                                 
##  3rd Qu.:510.0   3rd Qu.:350                                 
##  Max.   :565.0   Max.   :390

Nota: grazie al summary possiamo conoscere meglio le nostre variabili, controllando mi rendo conto di alcuni problemi di cleaning nel campo Anni.madre, dato che il minimo rilevato è 0, posso dire con certezza che ci sono delle trasformazioni necessarie da portare a termine, nel prossima cella cercheremo di capire come.

1.1 Analisi approfondita e Pulizia per la variabile Anni.madre

# --- ANALISI SPECIFICA: Controllo Valori Anomali 'Anni.madre' ---
# Obiettivo: Identificare valori errati o poco plausibili (outliers) nella variabile età.

# -----------------------------------------------------------------------------
# 1. TABELLA DI FREQUENZA COMPLETA
# -----------------------------------------------------------------------------
print("--- Tabella Frequenze: Anni Madre ---")

## [1] "--- Tabella Frequenze: Anni Madre ---"

freq_table <- data %>%
  count(Anni.madre) %>%        
  arrange(Anni.madre) %>%             
  mutate(Percentuale = round(n / sum(n) * 100, 2))

print(knitr::kable(freq_table, caption = "Distribuzione Età Madri"))

## 
## 
## Table: Distribuzione Età Madri
## 
## | Anni.madre|   n| Percentuale|
## |----------:|---:|-----------:|
## |          0|   1|        0.04|
## |          1|   1|        0.04|
## |         13|   1|        0.04|
## |         14|   2|        0.08|
## |         15|   6|        0.24|
## |         16|  13|        0.52|
## |         17|  18|        0.72|
## |         18|  24|        0.96|
## |         19|  45|        1.80|
## |         20|  66|        2.64|
## |         21|  74|        2.96|
## |         22| 100|        4.00|
## |         23| 115|        4.60|
## |         24| 131|        5.24|
## |         25| 180|        7.20|
## |         26| 184|        7.36|
## |         27| 197|        7.88|
## |         28| 172|        6.88|
## |         29| 174|        6.96|
## |         30| 200|        8.00|
## |         31| 147|        5.88|
## |         32| 159|        6.36|
## |         33| 110|        4.40|
## |         34|  96|        3.84|
## |         35|  66|        2.64|
## |         36|  64|        2.56|
## |         37|  41|        1.64|
## |         38|  38|        1.52|
## |         39|  27|        1.08|
## |         40|  19|        0.76|
## |         41|  13|        0.52|
## |         42|   8|        0.32|
## |         43|   2|        0.08|
## |         44|   4|        0.16|
## |         45|   1|        0.04|
## |         46|   1|        0.04|

# -----------------------------------------------------------------------------
# 2. IDENTIFICAZIONE VALORI SOSPETTI (Check Biologico)
# -----------------------------------------------------------------------------
soglia_min <- 15
soglia_max <- 55

valori_sospetti <- data %>%
  filter(Anni.madre < soglia_min | Anni.madre > soglia_max) %>%
  dplyr::select(Anni.madre, everything())

if (nrow(valori_sospetti) > 0) {
  print(paste("ATTENZIONE: Trovati", nrow(valori_sospetti), "valori sospetti (Età <", soglia_min, "o >", soglia_max, "):"))
  print(knitr::kable(valori_sospetti))
} else {
  print(paste("Nessun valore sospetto trovato (Tutte le età sono tra", soglia_min, "e", soglia_max, ")."))
}

## [1] "ATTENZIONE: Trovati 5 valori sospetti (Età < 15 o > 55 ):"
## 
## 
## | Anni.madre| N.gravidanze|Fumatrici | Gestazione| Peso| Lunghezza| Cranio|Tipo.parto |Ospedale |Sesso |
## |----------:|------------:|:---------|----------:|----:|---------:|------:|:----------|:--------|:-----|
## |         13|            0|0         |         38| 2760|       470|    325|Nat        |osp2     |F     |
## |         14|            1|0         |         39| 3510|       490|    365|Nat        |osp2     |M     |
## |          1|            1|0         |         41| 3250|       490|    350|Nat        |osp2     |F     |
## |          0|            0|0         |         39| 3060|       490|    330|Nat        |osp3     |M     |
## |         14|            0|0         |         39| 3550|       500|    355|Ces        |osp1     |M     |

# -----------------------------------------------------------------------------
# 3. VISUALIZZAZIONE GRAFICA (Dot Plot)
# -----------------------------------------------------------------------------
ggplot(data, aes(x = Anni.madre)) +
  geom_dotplot(binwidth = 1, fill = "orange", color = "black", dotsize = 0.8) +
  geom_vline(xintercept = c(soglia_min, soglia_max), color = "red", linetype = "dashed", size = 1) +
  scale_x_continuous(breaks = seq(min(data$Anni.madre, na.rm=TRUE), 
                                  max(data$Anni.madre, na.rm=TRUE), by = 2)) +
  labs(title = "Distribuzione Anni Madre (Controllo Errori)",
       subtitle = "I punti oltre le linee rosse tratteggiate potrebbero essere errori.",
       x = "Età Madre", y = "Conteggio (Proporzione)") +
  theme_minimal() +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

Pulizia e aggiornamento dataset originale

Nota: qui identifichiamo i valori estremi e quanti sono, ipotizzando che il periodo in cui è possibile rimanere fecondi è intorno ai 12 anni, ma secondo alcune ricerche condotte esternamente il periodo di maggiore fecondazione è tra 20 e 25 anni, decido di troncare i valori sotto ai 15 anni. Dal momento che nel dataset si riscontra una triplicazione di casi con divario di 1 anno (da 14 a 15 i casi di parto triplicano suggerendoci che sia meno dovuto al caso e che non siano presenti errori di cleaning).

# 2. Analisi Pre-Filtraggio
n_totale <- nrow(data)
n_da_rimuovere <- sum(data$Anni.madre < 15, na.rm = TRUE)

print(paste("Totale osservazioni:", n_totale))

## [1] "Totale osservazioni: 2500"

print(paste("Osservazioni con Anni.madre < 15:", n_da_rimuovere))

## [1] "Osservazioni con Anni.madre < 15: 5"

if (n_da_rimuovere > 0) {
  # 3. Filtraggio e Riassegnazione
  # Sovrascriviamo la variabile 'data' mantenendo solo chi ha >= 15 anni
  data <- data %>% filter(Anni.madre >= 15)
  
  # 4. Verifica Post-Filtraggio
  n_rimasti <- nrow(data)
  print("--- Situazione Finale ---")
  print(paste("Osservazioni rimanenti:", n_rimasti))
  print(paste("Età minima attuale nel dataset:", min(data$Anni.madre, na.rm=TRUE)))
  
  if ((n_totale - n_rimasti) == n_da_rimuovere) {
    print("SUCCESSO: Il dataset 'data' è stato aggiornato ed è pronto per le analisi successive.")
  } 
} else {
  print("Nessuna osservazione da rimuovere. Il dataset 'data' è rimasto invariato.")
}

## [1] "--- Situazione Finale ---"
## [1] "Osservazioni rimanenti: 2495"
## [1] "Età minima attuale nel dataset: 15"
## [1] "SUCCESSO: Il dataset 'data' è stato aggiornato ed è pronto per le analisi successive."

# Ristampiamo il summary per verificare la modifica.
summary(data)

##    Anni.madre    N.gravidanze     Fumatrici   Gestazione         Peso     
##  Min.   :15.0   Min.   : 0.0000   0:2391    Min.   :25.00   Min.   : 830  
##  1st Qu.:25.0   1st Qu.: 0.0000   1: 104    1st Qu.:38.00   1st Qu.:2990  
##  Median :28.0   Median : 1.0000             Median :39.00   Median :3300  
##  Mean   :28.2   Mean   : 0.9824             Mean   :38.98   Mean   :3284  
##  3rd Qu.:32.0   3rd Qu.: 1.0000             3rd Qu.:40.00   3rd Qu.:3620  
##  Max.   :46.0   Max.   :12.0000             Max.   :43.00   Max.   :4930  
##    Lunghezza         Cranio    Tipo.parto Ospedale   Sesso   
##  Min.   :310.0   Min.   :235   Ces: 727   osp1:815   F:1254  
##  1st Qu.:480.0   1st Qu.:330   Nat:1768   osp2:846   M:1241  
##  Median :500.0   Median :340              osp3:834           
##  Mean   :494.7   Mean   :340                                 
##  3rd Qu.:510.0   3rd Qu.:350                                 
##  Max.   :565.0   Max.   :390

Come si può vedere ora abbiamo un range di età materna che varia da 15 anni ad un massimo di 46, un range di 31 anni.

2. Analisi Descrittiva (EDA)

2.1 Statistiche Variabili Numeriche

Calcoliamo media, deviazione standard, coefficiente di variazione (CV) e indici di forma (asimmetria e curtosi) per avere un quadro generale della distribuzione dei dati.

# Calcolo specifico degli Indici di Variazione
variation_summary <- data.frame()

for (var in colnames(data)) {
  vec <- data[[var]]
  
  if (is.numeric(vec)) {
    # Calcolo metriche di dispersione
    varianza <- var(vec, na.rm = TRUE)
    dev_std <- sd(vec, na.rm = TRUE)
    cv_pct <- (dev_std / mean(vec, na.rm = TRUE)) * 100
    range_val <- diff(range(vec, na.rm = TRUE))
    iqr_val <- IQR(vec, na.rm = TRUE)
    
    temp <- data.frame(
      Variabile = var,
      Varianza = varianza,
      Dev_Std = dev_std,
      CV_Percent = cv_pct,
      Range = range_val,
      IQR = iqr_val
      
    )
    variation_summary <- rbind(variation_summary, temp)
  }
}

knitr::kable(variation_summary, digits = 2, caption = "Indici di Variazione (Dispersione)")

Indici di Variazione (Dispersione)

Variabile	Varianza	Dev_Std	CV_Percent	Range	IQR
Anni.madre	27.00	5.20	18.42	31	7
N.gravidanze	1.64	1.28	130.44	12	1
Gestazione	3.50	1.87	4.80	18	2
Peso	276038.78	525.39	16.00	4100	630
Lunghezza	693.78	26.34	5.32	255	30
Cranio	269.82	16.43	4.83	155	20

shapiro.test(data$Anni.madre)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$Anni.madre
## W = 0.99414, p-value = 2.063e-08

shapiro.test(data$Peso)

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

Di seguito riportiamo il quadro completo che include anche gli 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 Complete")

Statistiche Descrittive Complete

Variabile	Media	SD	CV_Percent	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

Analisi univariate per ogni variabile:

Caratteristica con tasso di variazione più alto: Numero gravidanze <- dai plot vediamo meglio il grande divario tra le diverse modalità

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 Distribuzione Variabili Continue

Visualizziamo la forma delle distribuzioni. Notiamo come alcune variabili (es. Peso) possano approssimare una normale, mentre altre potrebbero essere asimmetriche.

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

Nota: vediamo come la variabile Gestazione va da un minimo di 25 settimane ad un max di 43, questo arco di tempo è idoneo per permettere al feto di crescere e di essere partorito. Quindi non andiamo a cambiare il dataset. Anche dopo la trasformazione del dataset eseguiamo un test di Shapiro Wilk per capire se la distribuzione degli Anni.madre segue un’andamento normale.

P-value <<< 0,05 –> quindi si rigetta l’ipotesi nulla nonostante il valore del test è molto prossimo all’1 e che suggerirebbe quindi un andamento normale della distribuzione.

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)) {
    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
    
    idx_outliers <- which(colonna_corrente < lower_bound | colonna_corrente > upper_bound)
    val_outliers <- colonna_corrente[idx_outliers]
    
    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)
    }
    
    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: 497"
## 
##   Anni.madre       Cranio   Gestazione    Lunghezza N.gravidanze         Peso 
##            8           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.

Analisi approfondita degli Outliers

Vediamo se ci sono differenze significative tra outliers e non:

if (exists("df_outliers") && nrow(df_outliers) > 0) {
  indici_outliers <- unique(df_outliers$Riga_Indice)
  data_outliers_subset <- data[indici_outliers, ]
  data_clean_subset <- data[-indici_outliers, ]
  
  print(paste("Numero di osservazioni (righe) con almeno un outlier:", nrow(data_outliers_subset)))
  print(paste("Percentuale sul totale:", round(nrow(data_outliers_subset)/nrow(data)*100, 2), "%"))
  
  # Statistiche Descrittive Outliers
  stats_outliers <- data.frame()
  for (var in colnames(data_outliers_subset)) {
    vec <- data_outliers_subset[[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_outliers <- rbind(stats_outliers, temp)
    }
  }
  print("--- Statistiche Descrittive del Sottoinsieme Outliers ---")
  print(knitr::kable(stats_outliers, digits = 2, caption = "Statistiche: Subset Outliers"))
  
  # Generazione Grafici Distribuzione Outliers
  print("--- Generazione Grafici Distribuzione Outliers ---")
  for (var in colnames(data_outliers_subset)) {
    if (is.numeric(data_outliers_subset[[var]])) {
      mean_val <- mean(data_outliers_subset[[var]], na.rm = TRUE)
      p <- ggplot(data_outliers_subset, aes_string(x = var)) +
        geom_histogram(aes(y = ..density..), bins = 20, fill = "salmon", color = "white", alpha = 0.7) +
        geom_density(alpha = 0.2, fill = "red") +
        geom_vline(aes(xintercept = mean_val), color = "black", linetype = "dashed", size = 1) +
        labs(title = paste("Distribuzione (Solo Outliers):", var), x = var, y = "Densità") +
        theme_minimal()
      print(p)
    } else {
      p <- ggplot(data_outliers_subset, aes_string(x = var)) +
        geom_bar(fill = "darkred", alpha = 0.7) +
        labs(title = paste("Frequenza (Solo Outliers):", var), x = var, y = "Conteggio") +
        theme_minimal()
      print(p)
    }
  }
  
  # TEST DI CONFRONTO MEDIE: N.GRAVIDANZE
  print("--- Test T per N.gravidanze: Outliers vs Resto del Dataset ---")
  grav_outliers <- data_outliers_subset$N.gravidanze
  grav_clean <- data_clean_subset$N.gravidanze
  t_test_res <- t.test(grav_outliers, grav_clean)
  print(t_test_res)
  
  if(t_test_res$p.value < 0.05) {
    print("RISULTATO SIGNIFICATIVO: Il gruppo degli Outliers ha un numero medio di gravidanze DIVERSO dal resto della popolazione.")
  } else {
    print("NESSUNA DIFFERENZA SIGNIFICATIVA: Essere un outlier (in generale) non sembra correlato al numero di gravidanze.")
  }
  
  # Visualizzazione Comparativa
  df_compare <- data.frame(
    N_Gravidanze = c(grav_outliers, grav_clean),
    Gruppo = factor(c(rep("Outliers", length(grav_outliers)), rep("Non-Outliers", length(grav_clean))))
  )
  p_comp <- ggplot(df_compare, aes(x = Gruppo, y = N_Gravidanze, fill = Gruppo)) +
    geom_boxplot(alpha = 0.6) +
    stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "black") + 
    labs(title = "Confronto N. Gravidanze: Outliers vs Non-Outliers", y = "Numero di Gravidanze") +
    theme_minimal() + scale_fill_manual(values = c("lightgreen", "salmon"))
  print(p_comp)
}

## [1] "Numero di osservazioni (righe) con almeno un outlier: 351"
## [1] "Percentuale sul totale: 14.07 %"
## [1] "--- Statistiche Descrittive del Sottoinsieme Outliers ---"
## 
## 
## Table: Statistiche: Subset Outliers
## 
## |Variabile    |   Media|     SD| CV_Percent| Skewness| Kurtosis|
## |:------------|-------:|------:|----------:|--------:|--------:|
## |Anni.madre   |   31.86|   5.47|      17.15|    -0.14|     2.78|
## |N.gravidanze |    2.89|   2.06|      71.15|     0.96|     5.36|
## |Gestazione   |   37.53|   3.34|       8.91|    -1.22|     4.08|
## |Peso         | 3051.91| 878.53|      28.79|    -0.40|     2.74|
## |Lunghezza    |  478.03|  45.17|       9.45|    -1.15|     4.49|
## |Cranio       |  334.28|  26.85|       8.03|    -0.75|     3.80|
## [1] "--- Generazione Grafici Distribuzione Outliers ---"

## [1] "--- Test T per N.gravidanze: Outliers vs Resto del Dataset ---"
## 
##  Welch Two Sample t-test
## 
## data:  grav_outliers and grav_clean
## t = 20.019, df = 364.45, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  2.000705 2.436588
## sample estimates:
## mean of x mean of y 
## 2.8888889 0.6702425 
## 
## [1] "RISULTATO SIGNIFICATIVO: Il gruppo degli Outliers ha un numero medio di gravidanze DIVERSO dal resto della popolazione."

Commento: Il gruppo degli Outliers ha un numero medio di gravidanze DIVERSO dal resto della popolazione. Quindi deduciamo con significatività statistica che le medie di n.gravidanze tra outliers e non è diversa.

Distribuzioni degli ouliers:

Anni.madre –> si distribuisce come una gaussiana, con valori molto alti di Shapiro-Wilk e p-value significativo per non rifiutare l’ipotesi nulla di normalità, a differenza della distribuzione della stessa variabile ma per i non outliers.

N.gravidanze –> tra tutto il dataset questa è la variabile che cambia molto rispetto all’intero dataset, mentre nell’intero dataset abbiamo una media pari ad una gravidanza circa, qui notiamo come aumenta fino ad arrivare a 3 gravidanze arrotondando.

Il resto delle variabili rimangono molto vicini ai valori dei non outliers.

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

Dall’analisi visiva dei boxplot, possiamo vedere come le medie delle coppie di variabili sono pressappoco le stesse, le uniche variabili a suggerirci una lieve differenza sono i grafici della variabile Sesso (e le sue 2 modalità) in relazione al Peso, inoltre anche Fumatrici e non in relazione al Peso. Ma prima di saltare a conclusioni, termineremo l’analisi nelle celle successive in cui vedrai la significatività statistica dei risultati di cui sopra.

3. Analisi di Correlazione

3.1 Matrice di Scatterplot e Correlazioni (Pairs Plot)

Utilizziamo una visualizzazione avanzata per vedere simultaneamente la distribuzione grafica e il coefficiente di Pearson tra tutte le variabili numeriche.

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

Interpretazione: notiamo visivamente le correlazioni tra le varie coppie di variabili. Le variabili più correlate fra loro sono:

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.

3.2 Heatmap delle Correlazioni

Una visione alternativa per individuare rapidamente la multicollinearità (rosso scuro o blu scuro)

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 = 2605601257, p-value = 0.7426
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##          rho 
## -0.006578443

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 in questo campione.

Nota: si usa spearman in questo caso perché le variabili non si distribuiscono come una normale.

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

t_test_smoke <- t.test(Peso ~ Fumatrici, data = data, var.equal = FALSE)
print(t_test_smoke)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Fumatrici
## t = 1.0365, df = 114.14, p-value = 0.3021
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -45.49575 145.36053
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.279        3236.346

# 2. T-Test: Peso ~ Sesso
print("--- T-Test: Peso ~ Sesso ---")

## [1] "--- T-Test: Peso ~ Sesso ---"

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

# 3. Chi-Quadro: Ospedale vs Tipo di Parto
print("--- Chi-Quadro: Ospedale vs Tipo Parto ---")

## [1] "--- Chi-Quadro: Ospedale vs Tipo Parto ---"

chisq_res <- chisq.test(table(data$Ospedale, data$Tipo.parto))
print(chisq_res)

## 
##  Pearson's Chi-squared test
## 
## data:  table(data$Ospedale, data$Tipo.parto)
## X-squared = 1.0993, df = 2, p-value = 0.5772

ggplot(data, aes(x = Ospedale, fill = Tipo.parto)) +
  geom_bar(position = "fill") +
  labs(title = "Tipo di Parto per Ospedale", subtitle = paste("P-value:", format.pval(chisq_res$p.value, digits=3)), y="Proporzione") +
  theme_minimal() +
  scale_y_continuous(labels = percent)

Senza conoscere l’incidenza delle altre variabili e analizzandole singolarmente, possiamo commentare i risultati dei test indipendenti appena lanciati.

Commento T-Test: Ospedale ~ Tipo Parto: Non rifiutiamo l’ipotesi nulla di indipendenza. Non c’è evidenza statistica che la frequenza dei tipi di parto (Cesareo vs Naturale) vari in base all’ospedale. Questo suggerisce che le pratiche o le casistiche sono omogenee tra le diverse strutture sanitarie considerate

Commento T-Test: Peso ~ Sesso: Il test conferma una differenza altamente significativa tra i sessi. I maschi pesano mediamente di più delle femmine (3408g contro 3161g). L’intervallo di confidenza al 95% ci dice che la differenza media stimata oscilla tra i 206g e i 287g a favore dei maschi. Questa variabile dovrà essere inclusa necessariamente nel modello predittivo.

Commento T-Test: Peso ~ Fumatrici: Il test non mostra una differenza statisticamente significativa tra le medie dei due gruppi (p > 0.05). Non abbiamo prove sufficienti per affermare che, in questo campione, lo status di fumatrice influenzi il peso alla nascita. È interessante notare i gradi di libertà bassi (df ≈ 114), che suggeriscono un numero esiguo di fumatrici nel dataset, il che potrebbe limitare la potenza del test nel rilevare differenze minori.

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.40349, df = 2369.6, p-value = 0.6866
## alternative hypothesis: true difference in means between group Pluripara and group Primipara is not equal to 0
## 95 percent confidence interval:
##  -32.93583  50.00089
## sample estimates:
## mean in group Pluripara mean in group Primipara 
##                3287.935                3279.403

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, ma comunque manteniamo come variabile di controllo come facciamo per la variabile sesso.

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)
print(summary(mod1))

## 
## Call:
## lm(formula = Peso ~ ., data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1123.73  -180.56   -14.42   163.17  2609.45 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                -6729.3344   142.2196 -47.317  < 2e-16 ***
## Anni.madre                     0.7550     1.1581   0.652   0.5145    
## N.gravidanze                   9.3864     5.9938   1.566   0.1175    
## Fumatrici1                   -30.9120    27.5973  -1.120   0.2628    
## Gestazione                    32.7399     3.8349   8.537  < 2e-16 ***
## Lunghezza                     10.2920     0.3012  34.167  < 2e-16 ***
## Cranio                        10.4545     0.4282  24.417  < 2e-16 ***
## Tipo.partoNat                 29.8576    12.1096   2.466   0.0137 *  
## Ospedaleosp2                 -11.0253    13.4697  -0.819   0.4131    
## Ospedaleosp3                  28.1308    13.5237   2.080   0.0376 *  
## SessoM                        77.5648    11.2010   6.925 5.54e-12 ***
## Gruppo_GravidanzePrimipara    -8.0889    15.2002  -0.532   0.5947    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.2 on 2483 degrees of freedom
## Multiple R-squared:  0.7288, Adjusted R-squared:  0.7276 
## F-statistic: 606.5 on 11 and 2483 DF,  p-value: < 2.2e-16

par(mfrow=c(2,2)); plot(mod1); par(mfrow=c(1,1))

# Modello 2: Rimozione Fumatrici (non significativa al T-Test e nel mod1)
mod2 <- update(mod1, . ~ . - Fumatrici)

# Modello 3: Reset su data pulita (equivalente a mod1 ma esplicito)
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)
# L'interazione cattura l'effetto congiunto delle dimensioni corporee
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 - Gruppo_Gravidanze) # Questo sarà il nostro candidato "Simple & Robust"
print(summary(mod9))

## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Sesso + Cranio:Lunghezza, 
##     data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1119.04  -185.48    -9.86   166.68  2489.98 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.784e+03  1.169e+02 -23.822  < 2e-16 ***
## N.gravidanze      1.051e+01  4.336e+00   2.424   0.0154 *  
## Gestazione        4.133e+01  3.686e+00  11.213  < 2e-16 ***
## SessoM            7.671e+01  1.125e+01   6.820 1.14e-11 ***
## Cranio:Lunghezza  2.617e-02  4.663e-04  56.122  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275.5 on 2490 degrees of freedom
## Multiple R-squared:  0.7256, Adjusted R-squared:  0.7251 
## F-statistic:  1646 on 4 and 2490 DF,  p-value: < 2.2e-16

par(mfrow=c(2,2)); plot(mod9); par(mfrow=c(1,1))

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

# Modello 11: Rimozione N.gravidanze (basandosi sul T-test precedente)
mod11 <- update(mod10, . ~ . - N.gravidanze)

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

Dai risultati del summary notiamo come la riduzioni del numero di variabili per il mod9 abbia migliorato le performance, anche se Rquadro adjusted è peggiorato lievemente, vediamo come la varianza spiegata è molto più alta e le variabli sono tutte significative.

5.2 Modellazione Avanzata (Trasformazioni GLM/Box-Cox)

Esploriamo se trasformare la variabile target o aggiungere termini quadratici migliora il fit.

Nel grafico Pairs si potevano notare le diverse correlazioni tra coppie di variabili, nel relazioni non lineari discusse mi suggerisce di provare a trasformare le variabili e testare modelli non lineari e capire se possono spiegare meglio e sono più affidabili di altri.

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

## 
## Call:
## lm(formula = log(Peso) ~ . - Fumatrici - N.gravidanze - Ospedale + 
##     Lunghezza * Cranio - Gruppo_Gravidanze, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.52396 -0.05330  0.00042  0.05153  0.66295 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.460e-01  3.143e-01   2.056   0.0399 *  
## Anni.madre        5.238e-04  3.296e-04   1.589   0.1121    
## Gestazione        1.408e-02  1.228e-03  11.469  < 2e-16 ***
## Lunghezza         1.161e-02  6.795e-04  17.090  < 2e-16 ***
## Cranio            1.525e-02  9.841e-04  15.496  < 2e-16 ***
## Tipo.partoNat     7.929e-03  3.721e-03   2.131   0.0332 *  
## SessoM            2.188e-02  3.455e-03   6.332 2.86e-10 ***
## Lunghezza:Cranio -2.419e-05  2.014e-06 -12.009  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08432 on 2487 degrees of freedom
## Multiple R-squared:  0.7892, Adjusted R-squared:  0.7886 
## F-statistic:  1330 on 7 and 2487 DF,  p-value: < 2.2e-16

par(mfrow=c(2,2)); plot(mod12); par(mfrow=c(1,1))

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

## 
## Call:
## lm(formula = log(Peso) ~ Anni.madre + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Sesso + I(Gestazione^2) + Lunghezza:Cranio, 
##     data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.38728 -0.05401 -0.00030  0.05111  0.70989 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       1.466e-01  3.252e-01   0.451   0.6523    
## Anni.madre        4.548e-04  3.279e-04   1.387   0.1655    
## Gestazione        1.270e-01  2.044e-02   6.210 6.19e-10 ***
## Lunghezza         8.290e-03  9.040e-04   9.170  < 2e-16 ***
## Cranio            1.049e-02  1.303e-03   8.053 1.24e-15 ***
## Tipo.partoNat     7.926e-03  3.699e-03   2.143   0.0322 *  
## SessoM            2.263e-02  3.437e-03   6.583 5.61e-11 ***
## I(Gestazione^2)  -1.484e-03  2.683e-04  -5.531 3.52e-08 ***
## Lunghezza:Cranio -1.455e-05  2.654e-06  -5.484 4.58e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08382 on 2486 degrees of freedom
## Multiple R-squared:  0.7917, Adjusted R-squared:  0.7911 
## F-statistic:  1181 on 8 and 2486 DF,  p-value: < 2.2e-16

par(mfrow=c(2,2)); plot(mod13); par(mfrow=c(1,1))

# --- ANALISI BOX-COX ---
# Verifica se il Logaritmo è la trasformazione ideale
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: Trasformazione Box-Cox ---
# Se lambda non è vicino a 0 o 1, applichiamo la trasformazione specifica
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 - Gruppo_Gravidanze, data = data)
} 

# --- MODELLO 15: Lineare + Gestazione Quadrata ---
# Manteniamo la struttura complessa ma con target Peso originale (non trasformato)
mod15 <- lm(Peso ~ . - Fumatrici - N.gravidanze - Ospedale + Lunghezza*Cranio + I(Gestazione^2) - Gruppo_Gravidanze, data = data)
print(summary(mod15))

## 
## Call:
## lm(formula = Peso ~ . - Fumatrici - N.gravidanze - Ospedale + 
##     Lunghezza * Cranio + I(Gestazione^2) - Gruppo_Gravidanze, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -200.05  -15.33   -6.63    8.16  469.62 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.146e+03  1.349e+02  45.578   <2e-16 ***
## Anni.madre        3.240e-01  1.343e-01   2.412   0.0159 *  
## Gestazione       -1.888e+02  8.404e+00 -22.467   <2e-16 ***
## Lunghezza        -1.847e+01  3.724e-01 -49.599   <2e-16 ***
## Cranio           -2.632e+01  5.356e-01 -49.142   <2e-16 ***
## Tipo.partoNat     1.682e+00  1.516e+00   1.109   0.2674    
## SessoM            6.659e-01  1.420e+00   0.469   0.6391    
## Peso_BoxCox       2.876e+02  7.292e-01 394.485   <2e-16 ***
## I(Gestazione^2)   2.412e+00  1.102e-01  21.894   <2e-16 ***
## Lunghezza:Cranio  5.350e-02  1.087e-03  49.199   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 34.32 on 2485 degrees of freedom
## Multiple R-squared:  0.9957, Adjusted R-squared:  0.9957 
## F-statistic: 6.467e+04 on 9 and 2485 DF,  p-value: < 2.2e-16

par(mfrow=c(2,2)); plot(mod15); par(mfrow=c(1,1))

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

Per questo compito ho usato il metodo Box-Cox per trovare il parametro gamma preferibile.

Possiamo vedere subito che graficamente abbiamo un peggioramento dell’efficacia del modello, anche se l’Rquadro adjusted è molto alto, vediamo che i residui non soddisfano le assunzioni di base per questo tipo di modello, come l’omoschedasticità dei residui e la normalità in distribuzione e i grafici suggeriscono anche la presenza di autocorrelazione dei residui, il che ci fa capire che il modello non è affidabile per l’inferenza.

Inoltre usando modelli trasformati (GLM) si notano dei valori leva maggiori e più vicini al limite della distanza di Cook (guarda grafico: Residuals vs Leverage).

Ma andiamo a confrontare tutti insieme.

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| -15117.936| -15077.182| 122538.3625|<0.001    |<0.001 |0.6302 |
## |mod15    |Peso        | 0.9957|  24735.768|  24799.810|  64666.4550|<0.001    |<0.001 |<0.001 |
## |mod13    |log Peso    | 0.7911|  -5279.003|  -5220.783|   1181.3931|<0.001    |<0.001 |0.0688 |
## |mod12    |log Peso    | 0.7886|  -5250.487|  -5198.089|   1329.9681|<0.001    |<0.001 |0.0616 |
## |mod14    |Peso_BoxCox | 0.7706|   6803.874|   6862.095|   1048.1125|<0.001    |<0.001 |0.0848 |
## |mod6     |Peso        | 0.7292|  35091.346|  35155.388|    747.1257|<0.001    |<0.001 |0.1186 |
## |mod10    |Peso        | 0.7286|  35092.576|  35133.330|   1340.1788|<0.001    |<0.001 |0.1275 |
## |mod11    |Peso        | 0.7278|  35099.198|  35134.130|   1667.9673|<0.001    |<0.001 |0.1392 |
## |mod1     |Peso        | 0.7276|  35108.049|  35183.736|    606.5443|<0.001    |<0.001 |0.1207 |
## |mod3     |Peso        | 0.7276|  35108.049|  35183.736|    606.5443|<0.001    |<0.001 |0.1207 |
## |mod2     |Peso        | 0.7276|  35107.310|  35177.174|    667.0049|<0.001    |<0.001 |0.1163 |
## |mod4     |Peso        | 0.7268|  35112.507|  35170.727|    830.2138|<0.001    |<0.001 |0.1192 |
## |mod7     |Peso        | 0.7253|  35124.465|  35176.864|    941.9196|<0.001    |0.337  |0.1147 |
## |mod9     |Peso        | 0.7251|  35123.495|  35158.427|   1645.7709|<0.001    |0.370  |0.1223 |
## |mod999   |Peso        | 0.7251|  35124.132|  35159.064|   1645.1916|<0.001    |0.343  |0.1258 |
## |mod8     |Peso        | 0.7250|  35126.906|  35173.482|   1096.6558|<0.001    |0.255  |0.1184 |
## |mod5_rev |Peso        | 0.5992|  36064.388|  36099.320|    933.1699|<0.001    |0.128  |0.0238 |
## |mod5     |Peso        | 0.5988|  36069.624|  36122.022|    532.8694|<0.001    |0.154  |0.0222 |

6.2 Analisi VIF (Multicollinearità)

Controlliamo se le variabili indipendenti sono troppo correlate tra loro.

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.20097975995677 |1.95647387861763 |1.0094089304997 |1.70552805704238 |2.08784319721726 |1.64049031837929 |1.00460794133163 |1.00505972001199 |1.04064385264289 |1.88707451936217  |-                |-                |-                |-                |-                |
## |mod10    |-                |1.02272335019553 |-               |1.65154876690726 |5.82777785933241 |-                |-                |-                |1.04095454510511 |-                 |5.55691708530652 |-                |-                |-                |-                |
## |mod1000  |-                |-                |-               |-                |-                |2.06871094520832 |-                |-                |1.05380182365294 |1.02852290795519  |-                |2.82647687685971 |1.75926046311684 |-                |-                |
## |mod11    |-                |-                |-               |1.63740962869979 |5.78270419084527 |-                |-                |-                |1.0397478802436  |-                 |5.49867898818683 |-                |-                |-                |-                |
## |mod12    |1.02868106690689 |-                |-               |1.84935482027669 |112.379453972833 |91.665323470413  |1.00346633891906 |-                |1.04717416355096 |-                 |313.587147992162 |-                |-                |-                |-                |
## |mod13    |1.03016915021577 |-                |-               |518.722541611669 |201.260349962752 |162.555712838952 |1.00346635442542 |-                |1.04880654988698 |-                 |550.938249728183 |-                |-                |492.376318754492 |-                |
## |mod14    |1.03016915021575 |-                |-               |518.722541611641 |201.260349962761 |162.555712838966 |1.00346635442553 |-                |1.0488065498871  |-                 |550.938249728217 |-                |-                |492.376318754492 |-                |
## |mod15    |1.03111171286199 |-                |-               |522.943229811226 |203.773794191495 |163.874528648497 |1.00540638412497 |-                |1.06733459883679 |-                 |551.762598209403 |4.37284793368075 |-                |495.317064564015 |-                |
## |mod2     |1.20078827619692 |1.95618985298894 |-               |1.69854345627303 |2.08393047801562 |1.64025449035963 |1.00420062366655 |1.00463190210113 |1.04040911808096 |1.88338442289539  |-                |-                |-                |-                |-                |
## |mod3     |1.20097975995677 |1.95647387861763 |1.0094089304997 |1.70552805704238 |2.08784319721726 |1.64049031837929 |1.00460794133163 |1.00505972001199 |1.04064385264289 |1.88707451936217  |-                |-                |-                |-                |-                |
## |mod4     |1.20002687839735 |1.95515448639891 |-               |1.69720865766338 |2.08185611080799 |1.63990831694497 |1.00375229507407 |-                |1.04016461144497 |1.88195958594633  |-                |-                |-                |-                |-                |
## |mod5     |1.19999475926698 |1.95262405687324 |-               |1.33012961766111 |-                |1.30552271621648 |1.00121622834664 |-                |1.02867994571879 |1.88195940687469  |-                |-                |-                |-                |-                |
## |mod5_rev |-                |-                |-               |1.31032476006067 |-                |1.30161462784384 |-                |-                |1.02786106968386 |1.02966548136286  |-                |-                |-                |-                |-                |
## |mod6     |1.20004504434146 |1.95607082039983 |-               |1.86875748668551 |112.792927267316 |92.1791139892427 |1.00437837051613 |-                |1.04814209080905 |1.88864152988099  |-                |-                |-                |-                |314.828793244112 |
## |mod7     |1.19889039329512 |1.95306376501972 |-               |1.59190154081061 |-                |-                |1.00124505446471 |-                |1.03997812849454 |1.87789147606498  |-                |-                |-                |-                |1.58508091405061 |
## |mod8     |1.19883461791227 |1.95296210753078 |-               |1.59184742729161 |-                |-                |-                |-                |1.03997789262826 |1.87739964901294  |-                |-                |-                |-                |1.58478750444519 |
## |mod9     |-                |1.01481332096081 |-               |1.56163215293438 |-                |-                |-                |-                |1.03982772845012 |-                 |-                |-                |-                |-                |1.57509719092921 |
## |mod999   |-                |1.01448682337985 |-               |-                |-                |-                |-                |-                |1.03948858273923 |-                 |-                |-                |1.58424058543783 |-                |1.60000778714964 |

Conclusione: dall’analisi VIF (multicollinearità) vediamo quali sono i valori di multicollinearità dei vari regressori per ogni modello creato. Anche se di poco il mod9 è migliore in termini di semplicità e non collinearità rispetto a mod7.

Quindi si preferisce continuare con mod9.

PS. : qui abbiamo la certezza che i modelli trasformati eseguiti prima non sono affidabile dato il loro alto valore VIF tra i regressori, questo indica multicollinearità e quindi inaffidabilità.

7. Approfondimenti e Conclusioni Finali

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.889344e+08|      0.3393|
## |AIC            | 3.512349e+04| -15117.9360|
## |Adj. R2        | 7.251000e-01|      0.9959|
## [1] "I modelli sono simili, si preferisce Mod9 per semplicità interpretativa."

Commento: in questa cella abbiamo comparato i risultati di AIC e BIC test dei due modelli per capire il migliore modello con trade-off semplicità e spiegabilità.

Il risultato è che i due modelli differiscono davvero di poco, ma si preferisce comunque usare mod9 per semplicità dal momento che una relazione lineare è meglio interpretabile rispetto a una non lineare.

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)
  print("Intervalli di Confidenza (95%):")
  print(conf_intervals)
  
  # 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)", ]
  
  
}

##                       Estimate   Std. Error    t value      Pr(>|t|)
## (Intercept)      -2.784045e+03 1.168670e+02 -23.822336 3.589374e-113
## N.gravidanze      1.050948e+01 4.336208e+00   2.423658  1.543589e-02
## Gestazione        4.133342e+01 3.686057e+00  11.213449  1.667693e-28
## SessoM            7.670764e+01 1.124701e+01   6.820269  1.136562e-11
## Cranio:Lunghezza  2.616832e-02 4.662765e-04  56.121896  0.000000e+00
## [1] "Intervalli di Confidenza (95%):"
##                          2.5 %        97.5 %
## (Intercept)      -3.013211e+03 -2.554878e+03
## N.gravidanze      2.006541e+00  1.901243e+01
## Gestazione        3.410536e+01  4.856147e+01
## SessoM            5.465318e+01  9.876210e+01
## Cranio:Lunghezza  2.525399e-02  2.708265e-02

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.