neonati_analysis

1. Caricamento dei pacchetti e importazione dei dati

Per iniziare l’analisi, sono stati caricati i principali pacchetti R necessari per la manipolazione, la visualizzazione e la modellazione statistica dei dati. In particolare:

• tidyverse è usato per l’importazione e la trasformazione dei dati;

• ggplot2 e cowplot per la visualizzazione grafica;

• car per i test di omogeneità delle varianze e diagnostica del modello;

• broom per organizzare i risultati della regressione in formato tabellare;

• ggcorrplot per la costruzione di matrici di correlazione.

Successivamente, è stato importato il dataset contenente le informazioni relative a 2500 neonati e alle caratteristiche cliniche e anagrafiche delle loro madri. I nomi delle colonne sono stati convertiti nel formato snake_case per migliorarne la leggibilità e l’utilizzo nel codice.

Questa fase è fondamentale per impostare correttamente il flusso dell’analisi, garantendo che i dati siano facilmente gestibili e che le variabili siano denominate in modo coerente e intuitivo.

library(tidyverse)
library(reshape2)
library(cowplot)
library(car)
library(broom)
library(ggcorrplot)

# 1. Importazione del dataset
data <- read_csv("neonati.csv")

# 2. Rinomino colonne per usare nomi in snake_case (più leggibili)
data <- data %>%
  rename_with(~ str_replace_all(., "\\.", "_")) %>%
  rename(
    eta_madre = Anni_madre,
    num_gravidanze = N_gravidanze,
    fumo = Fumatrici,
    settimane = Gestazione,
    peso = Peso,
    lunghezza = Lunghezza,
    cranio = Cranio,
    tipo_parto = Tipo_parto,
    ospedale = Ospedale,
    sesso = Sesso
  )

2. Esplorazione del dataset e pulizia

Dopo l’importazione dei dati, è stato effettuato un primo controllo sulla struttura del dataset e sulla tipologia delle variabili. Le variabili categoriche (fumo, sesso, tipo_parto, ospedale) sono state convertite in fattori, in modo da poter essere correttamente interpretate durante l’analisi statistica e nella successiva modellazione.

La funzione str() ha confermato che il dataset è composto da 2500 osservazioni e 10 variabili. Le variabili numeriche come peso, lunghezza, cranio e settimane sono state rilevate come variabili di tipo numerico continuo, mentre le variabili categoriche presentano un numero di livelli coerente con quanto atteso (es. 2 livelli per sesso, 3 per ospedale, ecc.).

Successivamente, è stata eseguita una descrizione statistica delle variabili numeriche, da cui si osserva:

• una media del peso alla nascita di circa 3284 grammi,

• una durata media della gestazione di 38.98 settimane,

• una lunghezza media dei neonati di circa 494.7 mm,

• una circonferenza cranica media di circa 340 mm.

Infine, è stata verificata la presenza di valori mancanti. La funzione colSums(is.na()) ha confermato che non ci sono valori mancanti nel dataset, pertanto non è necessario applicare metodi di imputazione o rimozione.

Questa fase preliminare garantisce la qualità dei dati su cui si baseranno le analisi successive.

Successivamente sono stati utilizzati i boxplot per esplorare le sei variabili numeriche principali: eta_madre, num_gravidanze, settimane, peso, lunghezza e cranio.

L’obiettivo è stato quello di individuare rapidamente la presenza di valori anomali (outlier) potenzialmente errati o estremi. Come visibile nei grafici, sono emersi valori sospetti in tutte le variabili, evidenziati come puntini rossi esterni alle “whiskers” del boxplot. In particolare:

• Per eta_madre, si notano alcuni valori estremamente bassi (0 e 1 anni), clinicamente non plausibili;

• num_gravidanze mostra valori elevati fino a 12, rari ma teoricamente possibili;

• settimane presenta valori inferiori a 30 settimane, compatibili con gravidanze molto premature;

• peso, lunghezza e cranio presentano valori molto bassi compatibili con nati pretermine o patologici, ma che meritano comunque un’osservazione.

In seguito a questa analisi grafica, si è deciso di intervenire solo sui valori chiaramente impossibili, come l’età materna pari a 0 o 1 anno.

Rimozione di valori non plausibili

È stato applicato un filtro per rimuovere tutte le osservazioni con eta_madre ≤ 1, considerate errori di inserimento o registrazione. Dopo questa operazione, il dataset contiene solo età materne uguali o superiori a 13 anni, che rappresentano una soglia più realistica dal punto di vista biologico e sociale.

Questa fase di pulizia preliminare dei dati consente di migliorare l’affidabilità delle successive analisi statistiche e di evitare che i modelli predittivi siano influenzati da casi palesemente errati.

# 1. Conversione delle variabili categoriche in fattori
data <- data %>%
  mutate(
    tipo_parto = as.factor(tipo_parto),
    ospedale = as.factor(ospedale),
    sesso = as.factor(sesso),
    fumo = as.factor(fumo)
  )

# 2. Esplorazione del dataset

# Struttura del dataset: tipo di ogni colonna
str(data)

## tibble [2,500 × 10] (S3: tbl_df/tbl/data.frame)
##  $ eta_madre     : num [1:2500] 26 21 34 28 20 32 26 25 22 23 ...
##  $ num_gravidanze: num [1:2500] 0 2 3 1 0 0 1 0 1 0 ...
##  $ fumo          : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ settimane     : num [1:2500] 42 39 38 41 38 40 39 40 40 41 ...
##  $ peso          : num [1:2500] 3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
##  $ lunghezza     : num [1:2500] 490 490 500 515 480 495 480 510 500 510 ...
##  $ cranio        : num [1:2500] 325 345 375 365 335 340 345 349 335 362 ...
##  $ tipo_parto    : Factor w/ 2 levels "Ces","Nat": 2 2 2 2 2 2 2 2 1 1 ...
##  $ ospedale      : Factor w/ 3 levels "osp1","osp2",..: 3 1 2 2 3 2 3 1 2 2 ...
##  $ sesso         : Factor w/ 2 levels "F","M": 2 1 2 2 1 1 1 2 1 1 ...

# Statistiche descrittive per le variabili numeriche
summary(data)

##    eta_madre     num_gravidanze    fumo       settimane          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

# 3. Verifica dei valori mancanti in ogni variabile
# Mostra il numero di NA presenti per ciascuna colonna
colSums(is.na(data))

##      eta_madre num_gravidanze           fumo      settimane           peso 
##              0              0              0              0              0 
##      lunghezza         cranio     tipo_parto       ospedale          sesso 
##              0              0              0              0              0

colSums(is.na(data)) %>% .[. > 0]

## named numeric(0)

# 5.Boxplot per identificare outlier nelle variabili continue principali
# Lista delle variabili continue
vars <- c("eta_madre", "num_gravidanze", "settimane", "peso", "lunghezza", "cranio")

# Creazione di un boxplot per ciascuna variabile
boxplot_list <- lapply(vars, function(var) {
  ggplot(data, aes(y = .data[[var]])) +
    geom_boxplot(outlier.colour = "red", outlier.shape = 8) +
    theme_minimal() +
    labs(title = paste("Boxplot di", var), y = NULL, x = NULL)
})

# Visualizzazione affiancata (modifica ncol per il layout)
plot_grid(plotlist = boxplot_list, ncol = 3)

# 6.Verifica manuale di valori sospetti per anni madre
sort(unique(data$eta_madre))

##  [1]  0  1 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
## [26] 36 37 38 39 40 41 42 43 44 45 46

# 7.Rimozione dei valori non plausibili dell'età madre
data <- data %>%
  filter(eta_madre > 1)

3.Visualizzazione iniziale

In questa fase vengono esplorate graficamente le principali caratteristiche delle variabili numeriche e categoriali presenti nel dataset. L’obiettivo è comprendere la distribuzione dei dati, verificare la presenza di pattern rilevanti e individuare eventuali anomalie che possano influenzare la modellazione.

3.1 Distribuzione delle variabili numeriche.

Il primo grafico mostra gli istogrammi delle variabili continue, ovvero: età della madre, numero di gravidanze, settimane di gestazione, peso, lunghezza e circonferenza cranica. La maggior parte di queste variabili si distribuisce in modo approssimativamente simmetrico, ad eccezione del numero di gravidanze, che è fortemente asimmetrica sinistra(la maggior parte delle madri ha meno di 3 gravidanze, ma ci sono alcuni casi isolati con più di 10), il peso che mostra una distribuzione simile a quella normale ma con lieve asimmetria destra, e infine il numero di settimane che presenta una forte asimmetria sinistra.

# 1. Istogrammi delle variabili continue
data %>%
  select(eta_madre, num_gravidanze, settimane, peso, lunghezza, cranio) %>%
  pivot_longer(cols = everything(), names_to = "variabile", values_to = "valore") %>%
  ggplot(aes(x = valore)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  facet_wrap(~ variabile, scales = "free") +
  theme_minimal() +
  labs(title = "Distribuzione delle variabili numeriche", x = NULL, y = "Frequenza")

3.2 Peso in base a variabili categoriche.

I tre boxplot successivi mostrano la distribuzione del peso neonatale in funzione di:

• Sesso del neonato: si osserva che i maschi tendono ad avere un peso medio superiore rispetto alle femmine. La mediana del peso è infatti più alta per i maschi.

• Fumo materno: il peso dei neonati le cui madri fumano appare tendenzialmente inferiore rispetto a quello dei neonati da madri non fumatrici. Anche in questo caso si nota una differenza a livello di mediana e range interquartile.

• Tipo di parto: il confronto tra parti cesarei e naturali non mostra una differenza evidente, anche se i neonati da parto naturale sembrano avere una leggera tendenza ad avere pesi superiori.

# 2. Boxplot del peso in base a sesso, fumo, tipo di parto
ggplot(data, aes(x = sesso, y = peso, fill = sesso)) +
  geom_boxplot() +
  scale_fill_manual(values = c("F" = "#FF9999", "M" = "#9999FF")) +
  theme_minimal() +
  labs(title = "Distribuzione del peso per sesso", x = "Sesso", y = "Peso (g)")

ggplot(data, aes(x = fumo, y = peso, fill = fumo)) +
  geom_boxplot() +
  scale_fill_manual(values = c("0" = "#66C2A5", "1" = "#FC8D62")) +
  theme_minimal() +
  labs(title = "Peso in base al fumo materno", x = "Fumatrice", y = "Peso (g)")

ggplot(data, aes(x = tipo_parto, y = peso, fill = tipo_parto)) +
  geom_boxplot() +
  scale_fill_manual(values = c("Ces" = "#8DA0CB", "Nat" = "#E78AC3")) +
  theme_minimal() +
  labs(title = "Peso in base al tipo di parto", x = "Tipo di parto", y = "Peso (g)")

3.3 Relazione tra settimane di gestazione e peso.

Il grafico a dispersione mostra una chiara relazione positiva tra le settimane di gestazione e il peso alla nascita: all’aumentare delle settimane, il peso aumenta in modo pressoché lineare. Questo conferma l’importanza della durata della gestazione come predittore chiave del peso.

# 3. Relazione tra settimane di gestazione e peso (scatter)
ggplot(data, aes(x = settimane, y = peso)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE, color = "darkred") +
  theme_minimal() +
  labs(title = "Relazione tra settimane di gestazione e peso", x = "Settimane", y = "Peso (g)")

3.4 Correlazione tra variabili numeriche.

La heatmap delle correlazioni tra variabili continue evidenzia forti correlazioni positive tra peso e le seguenti variabili:

• Lunghezza (r = 0.80)

• Cranio (r = 0.70)

• Settimane di gestazione (r = 0.59)

Questo conferma che peso, lunghezza e cranio sono tra loro fortemente collegati, suggerendo che queste variabili misurano, almeno in parte, la crescita fisica complessiva del neonato. La bassa correlazione con l’età materna e il numero di gravidanze suggerisce che queste ultime variabili potrebbero avere un impatto più limitato o indiretto.

# 4. Heatmap di correlazione tra variabili numeriche
numeric_vars <- data %>%
  select(eta_madre, num_gravidanze, settimane, peso, lunghezza, cranio)

cor_matrix <- round(cor(numeric_vars), 2)

ggplot(melt(cor_matrix), aes(Var1, Var2, fill = value)) +
  geom_tile(color = "white") +
  geom_text(aes(label = value), color = "black", size = 4) +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0) +
  theme_minimal() +
  labs(title = "Heatmap di correlazione tra variabili numeriche", x = NULL, y = NULL)

3.5 Analisi delle interazioni tramite matrice di scatterplot.

Il grafico pairs() permette di osservare la relazione tra coppie di variabili continue distinguendo i punti per sesso. Le correlazioni osservate visivamente confermano i risultati della heatmap. In particolare:

• Si conferma la forte relazione lineare tra lunghezza, cranio e peso.

• Si osserva che le interazioni tra variabili come settimane e cranio, oppure peso e lunghezza, sono abbastanza lineari.

• Non sembrano emergere interazioni complesse (ad esempio non-lineari o a “cuneo”) tra le variabili, ma la correlazione forte tra alcune coppie (es. lunghezza-cranio) potrebbe giustificare l’inclusione di qualche effetto di interazione nella fase successiva.

pairs(data[, c("peso", "settimane", "cranio", "lunghezza")], col = data$sesso)

4.Analisi inferenziale

In questa fase sono stati condotti test inferenziali per verificare l’esistenza di differenze statisticamente significative tra gruppi per alcune variabili biometriche, e per confrontare il campione con valori medi attesi in letteratura.

4.1 Confronto del peso in base al sesso.

Il test di Levene non ha rilevato differenze significative nelle varianze del peso tra maschi e femmine (p = 0.3646), permettendo l’utilizzo di un t-test classico. Il test t ha evidenziato una differenza altamente significativa (p < 2.2e-16), con i neonati maschi che presentano un peso medio superiore di circa 247 grammi rispetto alle femmine (3408.5g vs 3161.1g). Questo risultato è coerente con la letteratura clinica.

# 1. Verifica dell'omogeneità delle varianze (test di Levene)
leveneTest(peso ~ sesso, data = data)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value Pr(>F)
## group    1  0.8222 0.3646
##       2496

# T-test a due campioni (maschi vs femmine)
t.test(peso ~ sesso, data = data, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  peso by sesso
## t = -12.111, df = 2496, 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.4963 -207.3721
## sample estimates:
## mean in group F mean in group M 
##        3161.061        3408.496

4.2 Confronto della lunghezza in base al sesso.

Il test di Levene ha indicato varianze disomogenee tra i due gruppi (p = 0.0012), motivo per cui è stato impiegato un t-test di Welch. Il test ha mostrato una differenza significativa (p < 2.2e-16), con una lunghezza media superiore nei maschi (499.7 mm) rispetto alle femmine (489.8 mm), con un intervallo di confidenza che esclude lo zero.

# Test varianze
leveneTest(lunghezza ~ sesso, data = data)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value   Pr(>F)   
## group    1  10.571 0.001164 **
##       2496                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# T-test lunghezza
t.test(lunghezza ~ sesso, data = data, var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  lunghezza by sesso
## t = -9.5823, df = 2457.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -11.939001  -7.882672
## sample estimates:
## mean in group F mean in group M 
##        489.7641        499.6750

4.3 Confronto della circonferenza cranica in base al sesso.

In questo caso, le varianze sono risultate omogenee (p = 0.2715) e il t-test ha confermato una differenza significativa tra i due gruppi (p = 1.4e-13), con i maschi che presentano un cranio mediamente più grande di circa 4.8 mm rispetto alle femmine.

# Test varianze
leveneTest(cranio ~ sesso, data = data)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value Pr(>F)
## group    1  1.2098 0.2715
##       2496

# T-test cranio
t.test(cranio ~ sesso, data = data, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  cranio by sesso
## t = -7.4344, df = 2496, p-value = 1.436e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.110877 -3.560044
## sample estimates:
## mean in group F mean in group M 
##        337.6231        342.4586

4.4 Confronto del peso in base al fumo materno.

Non sono state rilevate differenze significative nella varianza del peso tra figli di madri fumatrici e non fumatrici (p = 0.2324), e il t-test ha confermato l’assenza di differenze statisticamente significative nelle medie (p = 0.3428). Sebbene il peso medio dei neonati da madri fumatrici sia leggermente inferiore (3236g contro 3286g), la differenza non è sufficiente per essere considerata significativa.

# Verifica omogeneità varianze
leveneTest(peso ~ fumo, data = data)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value Pr(>F)
## group    1  1.4267 0.2324
##       2496

# T-test (peso in base al fumo)
t.test(peso ~ fumo, data = data, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  peso by fumo
## t = 0.94878, df = 2496, p-value = 0.3428
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -53.24919 153.08153
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.262        3236.346

4.5 Associazione tra tipo di parto e ospedale.

È stata costruita una tabella di contingenza tra tipo di parto (cesareo o naturale) e ospedale. Il test chi-quadro ha restituito un p-value pari a 0.5819, indicando che non esistono differenze statisticamente significative tra gli ospedali nella distribuzione dei tipi di parto. In altre parole, la frequenza di parti cesarei e naturali è simile nei tre ospedali analizzati.

# Tabella di contingenza tra tipo di parto e ospedale
tabella_parto_ospedale <- table(data$tipo_parto, data$ospedale)
print(tabella_parto_ospedale)

##      
##       osp1 osp2 osp3
##   Ces  242  254  232
##   Nat  574  594  602

# Test chi-quadro
chi_test <- chisq.test(tabella_parto_ospedale)
print(chi_test)

## 
##  Pearson's Chi-squared test
## 
## data:  tabella_parto_ospedale
## X-squared = 1.083, df = 2, p-value = 0.5819

4.6 Confronto con i valori medi della popolazione (OMS).

Infine, sono stati effettuati due test t per campione singolo per confrontare il campione con i valori medi di riferimento della WHO Child Growth Standards:

• Peso alla nascita: la media osservata (3284g) non differisce significativamente dal valore di riferimento di 3300g (p = 0.1324).

• Lunghezza alla nascita: la media osservata (494.7 mm) è significativamente inferiore al valore di riferimento di 500 mm (p < 2.2e-16).

📌 Fonte dei valori di riferimento: https://www.who.int/tools/child-growth-standards

#fonte valori medi Tabella di crescita OMS https://www.who.int/tools/child-growth-standards
t.test(data$peso, mu = 3300)

## 
##  One Sample t-test
## 
## data:  data$peso
## t = -1.505, df = 2497, p-value = 0.1324
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.577 3304.791
## sample estimates:
## mean of x 
##  3284.184

t.test(data$lunghezza, mu = 500)

## 
##  One Sample t-test
## 
## data:  data$lunghezza
## t = -10.069, df = 2497, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6628 495.7287
## sample estimates:
## mean of x 
##  494.6958

5.Costruzione del modello di regressione

In questa fase si è proceduto alla costruzione di un modello di regressione lineare multipla per prevedere il peso del neonato alla nascita a partire da un insieme di variabili cliniche e anagrafiche.

5.1 Modello completo.

Nel modello completo iniziale sono state incluse tutte le variabili disponibili: età della madre, numero di gravidanze, settimane di gestazione, lunghezza, cranio, sesso del neonato e abitudine al fumo. Ma sono state escluse a priori il tipo di parto e l’ospedale, in quanto non predittive per logica. I risultati del modello mostrano un buon livello di adattamento con un R² pari a 0.7272, indicando che circa il 72.7% della variabilità del peso neonatale è spiegata dal modello.

Tra i coefficienti stimati, le variabili che risultano statisticamente significative (p < 0.05) sono:

• Numero di gravidanze: ogni gravidanza in più è associata a un aumento medio di circa 11.4 grammi nel peso del neonato.

• Settimane di gestazione: ogni settimana aggiuntiva è associata a un incremento medio di circa 32.9 grammi.

• Lunghezza e cranio: ogni millimetro in più di lunghezza e cranio è associato rispettivamente a un aumento medio di 10.2 g e 10.5 g nel peso.

• Sesso maschile: i neonati maschi pesano mediamente 78 g in più rispetto alle femmine, a parità di altre condizioni.

L’età della madre e il fumo non risultano significativi nel modello completo, pur essendo variabili teoricamente rilevanti dal punto di vista clinico.

# 1. Modello di regressione lineare multipla completo
modello_completo <- lm(peso ~ eta_madre + num_gravidanze + settimane +
                         lunghezza + cranio + sesso + fumo, data = data)

# 2. Sommario del modello e valutazione
summary(modello_completo)

## 
## Call:
## lm(formula = peso ~ eta_madre + num_gravidanze + settimane + 
##     lunghezza + cranio + sesso + fumo, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1160.6  -181.3   -15.7   163.6  2630.7 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -6712.2405   141.3339 -47.492  < 2e-16 ***
## eta_madre          0.8803     1.1491   0.766    0.444    
## num_gravidanze    11.3789     4.6767   2.433    0.015 *  
## settimane         32.9472     3.8288   8.605  < 2e-16 ***
## lunghezza         10.2316     0.3011  33.979  < 2e-16 ***
## cranio            10.5198     0.4271  24.633  < 2e-16 ***
## sessoM            78.0787    11.2132   6.963 4.24e-12 ***
## fumo1            -30.3958    27.6080  -1.101    0.271    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2490 degrees of freedom
## Multiple R-squared:  0.7272, Adjusted R-squared:  0.7264 
## F-statistic: 948.3 on 7 and 2490 DF,  p-value: < 2.2e-16

5.2 Selezione automatica del modello.

Per ottenere un modello più parsimonioso, si è proceduto a una selezione automatica delle variabili secondo i criteri AIC (Akaike Information Criterion) e BIC (Bayesian Information Criterion), utilizzando la funzione step().

Entrambe le strategie di selezione (AIC e BIC) hanno condotto allo stesso modello ridotto, in cui sono state eliminate le variabili eta_madre e fumo, non risultate significative.

Il modello selezionato ha mantenuto:

• num_gravidanze

• settimane

• lunghezza

• cranio

• sesso

Con un R² identico (0.727), ma meno gradi di libertà occupati, quindi più efficiente e interpretabile. L’AIC e il BIC sono entrambi più bassi rispetto al modello completo:

Modello	AIC	BIC
Modello completo	35155.07	35207.48
Modello step_AIC	35152.89	35193.65
Modello step_BIC	35152.89	35193.65

Questo conferma che il modello selezionato è più parsimonioso e ugualmente performante, risultando quindi la scelta migliore per la previsione del peso neonatale.

modello_step_AIC <- step(modello_completo, direction = "both", trace = FALSE)
summary(modello_step_AIC)

## 
## Call:
## lm(formula = peso ~ num_gravidanze + settimane + lunghezza + 
##     cranio + sesso, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.37  -180.98   -15.57   163.69  2639.09 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -6681.7251   135.8036 -49.201  < 2e-16 ***
## num_gravidanze    12.4554     4.3416   2.869  0.00415 ** 
## settimane         32.3827     3.8008   8.520  < 2e-16 ***
## lunghezza         10.2455     0.3008  34.059  < 2e-16 ***
## cranio            10.5410     0.4265  24.717  < 2e-16 ***
## sessoM            77.9807    11.2111   6.956 4.47e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1327 on 5 and 2492 DF,  p-value: < 2.2e-16

n <- nrow(data)  # numero osservazioni
modello_step_BIC <- step(modello_completo, direction = "both", k = log(n), trace = FALSE)
summary(modello_step_BIC)

## 
## Call:
## lm(formula = peso ~ num_gravidanze + settimane + lunghezza + 
##     cranio + sesso, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.37  -180.98   -15.57   163.69  2639.09 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -6681.7251   135.8036 -49.201  < 2e-16 ***
## num_gravidanze    12.4554     4.3416   2.869  0.00415 ** 
## settimane         32.3827     3.8008   8.520  < 2e-16 ***
## lunghezza         10.2455     0.3008  34.059  < 2e-16 ***
## cranio            10.5410     0.4265  24.717  < 2e-16 ***
## sessoM            77.9807    11.2111   6.956 4.47e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1327 on 5 and 2492 DF,  p-value: < 2.2e-16

AIC(modello_completo, modello_step_AIC, modello_step_BIC)

##                  df      AIC
## modello_completo  9 35155.07
## modello_step_AIC  7 35152.89
## modello_step_BIC  7 35152.89

BIC(modello_completo, modello_step_AIC, modello_step_BIC)

##                  df      BIC
## modello_completo  9 35207.48
## modello_step_AIC  7 35193.65
## modello_step_BIC  7 35193.65

Analisi dei residui

A partire dal modello selezionato con BIC, sono stati poi verificati i presupposti classici della regressione lineare:

• Omogeneità della varianza: il grafico dei residui rispetto ai valori predetti non mostra pattern sistematici, suggerendo che l’ipotesi di omoscedasticità può essere considerata plausibile.

• Normalità dei residui: il QQ plot mostra un buon allineamento dei quantili osservati con quelli teorici, con qualche deviazione alle code, ma complessivamente accettabile.

• Individuazione di outlier influenti: la Cook’s distance ha evidenziato l’osservazione 1549 come particolarmente influente, da considerare con attenzione nell’interpretazione finale.

• Multicollinearità: i valori del VIF sono tutti inferiori a 2.1, indicando un’assenza di multicollinearità problematica tra i predittori.

Infine, si è valutata la possibilità di migliorare il modello includendo un termine di interazione tra lunghezza e cranio, suggerito da un’osservazione del grafico pairs e della heatmap. L’aggiunta di questo termine comporta un leggero aumento di R² (da 0.727 a 0.7296), ma la maggiore complessità del modello non appare giustificata. Pertanto, si conferma la preferenza per il modello semplificato selezionato tramite BIC.

# 3. Residui vs valori predetti
plot(modello_step_BIC$fitted.values, resid(modello_step_BIC),
     xlab = "Valori predetti", ylab = "Residui",
     main = "Residui vs Valori predetti")
abline(h = 0, col = "red", lty = 2)

# 6. QQ plot dei residui
qqnorm(resid(modello_step_BIC))
qqline(resid(modello_step_BIC), col = "blue", lwd = 2)

# 7. Scale-Location plot
plot(modello_step_BIC, which = 3)

# 8. Cook’s distance
plot(modello_step_BIC, which = 4)

# 9. Multicollinearità (VIF)
vif(modello_step_BIC)

## num_gravidanze      settimane      lunghezza         cranio          sesso 
##       1.023462       1.669779       2.075747       1.624568       1.040184

# 10. Calcolo dell'RMSE
rmse <- sqrt(mean(residuals(modello_step_BIC)^2))
print(rmse)

## [1] 274.3666

# 11. Modello con tutte le interazioni a 2 vie tra i predittori del modello BIC esteso
modello_interazioni <- lm(peso ~ num_gravidanze + settimane + lunghezza +
                                         cranio + sesso + lunghezza:cranio, data = data)
summary(modello_interazioni) #Vado a complicare il modello rispetto a quello che ho trovato prima per migliorare di pochissimo R2, scelgo il modello precedente perche piu semplice

## 
## Call:
## lm(formula = peso ~ num_gravidanze + settimane + lunghezza + 
##     cranio + sesso + lunghezza:cranio, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1150.65  -180.93   -13.48   165.99  2865.46 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.803e+03  1.018e+03  -1.771   0.0767 .  
## num_gravidanze    1.293e+01  4.323e+00   2.991   0.0028 ** 
## settimane         3.815e+01  3.967e+00   9.616  < 2e-16 ***
## lunghezza        -3.060e-01  2.203e+00  -0.139   0.8895    
## cranio           -4.755e+00  3.192e+00  -1.490   0.1365    
## sessoM            7.324e+01  1.120e+01   6.537 7.59e-11 ***
## lunghezza:cranio  3.157e-02  6.531e-03   4.835 1.41e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.5 on 2491 degrees of freedom
## Multiple R-squared:  0.7296, Adjusted R-squared:  0.7289 
## F-statistic:  1120 on 6 and 2491 DF,  p-value: < 2.2e-16

Identificazione delle osservazioni influenti

L’analisi della Cook’s distance segnala l’osservazione 1549 come fortemente influente, con una distanza superiore a tutte le altre. Per investigare il motivo di tale anomalia, è stata estratta la riga corrispondente. L’osservazione in questione rappresenta un neonato con un peso molto alto (4370g), ma con lunghezza decisamente inferiore alla media (315mm) e circonferenza cranica molto elevata (374mm). Il boxplot con evidenziazione dell’osservazione mostra che, pur essendo compreso nel range dei dati, questo valore risulta potenzialmente anomalo e meritevole di esclusione a scopo esplorativo.

data[1549, ]

## # A tibble: 1 × 10
##   eta_madre num_gravidanze fumo  settimane  peso lunghezza cranio tipo_parto
##       <dbl>          <dbl> <fct>     <dbl> <dbl>     <dbl>  <dbl> <fct>     
## 1        35              1 0            38  4370       315    374 Nat       
## # ℹ 2 more variables: ospedale <fct>, sesso <fct>

ggplot(data, aes(x = "", y = peso)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(alpha = 0.3, width = 0.1) +
  geom_point(data = data[1549, ], aes(x = "", y = peso), color = "red", size = 3) +
  theme_minimal() +
  labs(title = "Posizione dell'osservazione 1549 rispetto alla distribuzione del peso")

Modellazione senza osservazione influente

Per verificare l’impatto dell’osservazione 1549, è stato stimato un nuovo modello escludendo questo dato. Il modello aggiornato presenta un R² leggermente superiore (da 0.727 a 0.737), una riduzione dell’errore standard dei residui (da 274.7 a 269.3) e una stabilità dei coefficienti stimati, suggerendo che l’osservazione 1549 stesse introducendo rumore nella stima. Anche i grafici diagnostici relativi al nuovo modello (ultime tre immagini) mostrano un miglioramento generale nella distribuzione dei residui e nella riduzione dell’influenza di singoli punti.

Conclusione: Sebbene l’osservazione 1549 sia clinicamente plausibile, il suo peso e la combinazione di misure associate la rendono un valore potenzialmente anomalo per il contesto del modello. Si è deciso di escluderla nella versione finale del modello per ottenere stime più stabili e una diagnostica migliore. Tale scelta viene giustificata in un’ottica di robustezza del modello predittivo.

Inoltre, dalla nuova analisi dei residui il modello risulta valido, mentre i punti identificati dal grafico della Cook’s Distance hanno dei valori molto bassi, per cui non è necessario indagare ulteriormente.

modello_senza_outlier <- lm(peso ~ num_gravidanze + settimane + lunghezza +
                              cranio + sesso, data = data[-1549, ])
summary(modello_senza_outlier)

## 
## Call:
## lm(formula = peso ~ num_gravidanze + settimane + lunghezza + 
##     cranio + sesso, data = data[-1549, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1165.68  -179.74   -12.42   162.92  1410.68 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -6683.8326   133.1602 -50.194  < 2e-16 ***
## num_gravidanze    13.1448     4.2576   3.087  0.00204 ** 
## settimane         29.6341     3.7369   7.930 3.27e-15 ***
## lunghezza         10.8899     0.3019  36.077  < 2e-16 ***
## cranio             9.9192     0.4227  23.465  < 2e-16 ***
## sessoM            78.1376    10.9929   7.108 1.53e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.3 on 2491 degrees of freedom
## Multiple R-squared:  0.7372, Adjusted R-squared:  0.7367 
## F-statistic:  1398 on 5 and 2491 DF,  p-value: < 2.2e-16

# Nuovo modello
modello_finale <- lm(peso ~ num_gravidanze + settimane + lunghezza + cranio + sesso, data = data[-1549, ])

# Residui vs predetti
plot(modello_finale$fitted.values, resid(modello_finale),
     xlab = "Valori predetti", ylab = "Residui",
     main = "Residui vs Valori predetti")
abline(h = 0, col = "red", lty = 2)

# QQ plot
qqnorm(resid(modello_finale))
qqline(resid(modello_finale), col = "blue", lwd = 2)

# Cook’s distance
plot(modello_finale, which = 4)

6. Previsioni

Per illustrare la capacità predittiva del modello finale, è stata effettuata una stima del peso alla nascita in un caso ipotetico di interesse pratico. In particolare, si è voluto prevedere il peso di una neonata (sesso = F), alla terza gravidanza e partorita alla 39ª settimana di gestazione.

Poiché il modello finale include altre variabili (lunghezza, cranio, tipo di parto), queste sono state impostate ai valori medi o modali osservati nel dataset:

• Lunghezza: media del campione;

• Circonferenza cranica: media del campione;

• Tipo di parto: modalità più frequente tra cesareo e naturale.

Il modello ha fornito una stima del peso pari a circa 3271 grammi, con un intervallo di predizione compreso tra 2743g e 3800g (a livello di confidenza del 95%).

Questo intervallo riflette la variabilità attesa nel peso neonatale, anche a parità di caratteristiche osservabili. Il risultato è coerente con quanto atteso per un parto a termine e conferma la buona capacità del modello di fornire stime utili in contesti clinici e demografici realistici.

# 1. Calcola valori medi/modali dalle osservazioni reali
lunghezza_media <- mean(data$lunghezza)
cranio_medio <- mean(data$cranio)
tipo_piu_frequente <- names(sort(table(data$tipo_parto), decreasing = TRUE))[1]

# 2. Nuovo caso con solo sesso, settimane, n_gravidanze (gli altri completati automaticamente)
nuovo_caso_corretto <- data.frame(
  num_gravidanze = 3,
  settimane = 39,
  lunghezza = lunghezza_media,
  cranio = cranio_medio,
  sesso = factor("F", levels = levels(data$sesso)),
  tipo_parto = factor(tipo_piu_frequente, levels = levels(data$tipo_parto))
)

# 3. Previsione
predict(modello_finale, newdata = nuovo_caso_corretto, interval = "prediction")

##        fit      lwr      upr
## 1 3271.329 2742.661 3799.997

7. Visualizzazioni finali

Infine, sono stati costruiti diversi grafici per meglio interpretare e comunicare i risultati:

• Un scatter plot mostra la relazione tra settimane di gestazione e peso, distinguendo per sesso: la crescita ponderale risulta più marcata nei maschi.

• Un boxplot per il fumo materno suggerisce visivamente un leggero calo nel peso dei neonati di madri fumatrici, sebbene non significativo nel modello finale.

• Il grafico degli effetti stimati (con IC95%) evidenzia visivamente il contributo relativo di ciascuna variabile nel modello finale.

• Le due matrici di correlazione (visualizzate con ggcorrplot e con geom_tile) confermano che le correlazioni più forti si concentrano tra peso, lunghezza e cranio, a supporto delle scelte modellistiche.

ggplot(data, aes(x = settimane, y = peso, color = sesso)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal() +
  labs(title = "Peso vs Settimane di gestazione per sesso",
       x = "Settimane di gestazione", y = "Peso (g)")

ggplot(data, aes(x = fumo, y = peso, fill = fumo)) +
  geom_boxplot() +
  scale_fill_manual(values = c("0" = "#66C2A5", "1" = "#FC8D62")) +
  theme_minimal() +
  labs(title = "Peso alla nascita in base al fumo materno",
       x = "Fumo (0 = no, 1 = sì)", y = "Peso (g)")

coeff_df <- tidy(modello_finale, conf.int = TRUE) %>%
  filter(term != "(Intercept)")

ggplot(coeff_df, aes(x = reorder(term, estimate), y = estimate)) +
  geom_col(fill = "steelblue") +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.25) +
  coord_flip() +
  theme_minimal() +
  labs(title = "Effetti stimati delle variabili sul peso",
       x = NULL, y = "Stima (g)")

ggcorrplot(cor_matrix,
           method = "square",
           type = "lower", # o "full"
           lab = TRUE,
           lab_size = 3,
           colors = c("blue", "white", "red"),
           title = "Matrice di correlazione",
           ggtheme = theme_minimal())

cor_matrix <- round(cor(data %>%
                          select(eta_madre, num_gravidanze, settimane, peso, lunghezza, cranio)), 2)

melted_cor <- melt(cor_matrix)
melted_cor$Var1 <- factor(melted_cor$Var1, levels = colnames(cor_matrix))
melted_cor$Var2 <- factor(melted_cor$Var2, levels = colnames(cor_matrix))

ggplot(melted_cor, aes(Var1, Var2, fill = value)) +
  geom_tile(color = "white") +
  geom_text(aes(label = value), color = "black", size = 4) +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0) +
  theme_minimal() +
  labs(title = "Correlazione tra variabili numeriche", x = NULL, y = NULL)

summary(modello_finale)

## 
## Call:
## lm(formula = peso ~ num_gravidanze + settimane + lunghezza + 
##     cranio + sesso, data = data[-1549, ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1165.68  -179.74   -12.42   162.92  1410.68 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -6683.8326   133.1602 -50.194  < 2e-16 ***
## num_gravidanze    13.1448     4.2576   3.087  0.00204 ** 
## settimane         29.6341     3.7369   7.930 3.27e-15 ***
## lunghezza         10.8899     0.3019  36.077  < 2e-16 ***
## cranio             9.9192     0.4227  23.465  < 2e-16 ***
## sessoM            78.1376    10.9929   7.108 1.53e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.3 on 2491 degrees of freedom
## Multiple R-squared:  0.7372, Adjusted R-squared:  0.7367 
## F-statistic:  1398 on 5 and 2491 DF,  p-value: < 2.2e-16