Modello Statistico per la Previsione del Peso Neonatale

Obiettivo del progetto:

Creare un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita.

Benefici:

Miglioramento delle previsioni cliniche
Ottimizzazione delle risorse ospedaliere
Prevenzione e identificazione dei fattori di rischio
Valutazione delle pratiche ospedaliere
Supporto alla pianificazione strategica

knitr::opts_chunk$set(echo = TRUE)
options(warn = -1)

#installa se necessario i packages che verranno usati nel codice

if (!require("ggplot2")) {
  install.packages("ggplot2")
}
if (!require("car")) {
  install.packages("car")
}
if (!require("plotly")) {
  install.packages("plotly")
}
if (!require("TeachingDemos")) {
  install.packages("TeachingDemos")
}
if (!require("ggrepel")) {
  install.packages("ggrepel")
}
if (!require("ggrepel")) {
  install.packages("ggpubr", type = "binary")
}
if (!require("ggrepel")) {
  install.packages("lmtest")
}
if (!require("dplyr")) {
  install.packages("dplyr")
}
if (!require("kableExtra")) {
  install.packages("kableExtra")
}
if (!require("e1071")) {
  install.packages("e1071")
}
if (!require("gridExtra")) {
  install.packages("gridExtra")
}


#carica i packages
library(ggplot2)
library(car)
library(plotly)
library(TeachingDemos)
library(ggrepel)
library(ggpubr)
library(lmtest)
library(dplyr)
library(kableExtra)
library(purrr)
library(e1071)
library(gridExtra)

#Utility -> funzioni per il calcolo della Moda, Indici, Frequenza, GINI Index

# Funzione per il calcolo della moda
get_mode <- function(v) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}


# Funzione per calcolare tutti gli indici su un singolo vettore
# per variabili Quantitative, la funzione ha due paramtri:
# x -> variabile su cui calcolare gli indici
# variabile_continua -> indica se la variabile è continua o discreta
# se la variabile è continua non viene calcolato l'indice Moda
calcola_indici <- function(x,variabile_continua = FALSE) {
  if (variabile_continua)
    {
       data.frame(
                  media = mean(x, na.rm = TRUE), # Media aritmetica
                  mediana = median(x, na.rm = TRUE), # Mediana
                  minimo = min(x, na.rm = TRUE), # Minimo
                  massimo = max(x, na.rm = TRUE), # Massimo
                  Q1 = quantile(x, 0.25, na.rm = TRUE), # Primo quartile 
                  Q3 = quantile(x, 0.75, na.rm = TRUE), # Terzo quartile
                  varianza = var(x, na.rm = TRUE), # Varianza
                  deviazione_std = sd(x, na.rm = TRUE), # Deviazione standard
                  coeff_var = sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE), # Coefficiente di variazione
                  skewness = skewness(x, na.rm = TRUE), # Asimmetria (skewness)
                  kurtosis = kurtosis(x, na.rm = TRUE) # Curtosi
                )
  }
  else
  {
    
    data.frame(
                  media = mean(x, na.rm = TRUE),
                  mediana = median(x, na.rm = TRUE),
                  moda = get_mode(x), # Moda da funziona 
                  minimo = min(x, na.rm = TRUE),
                  massimo = max(x, na.rm = TRUE),
                  Q1 = quantile(x, 0.25, na.rm = TRUE),
                  Q3 = quantile(x, 0.75, na.rm = TRUE),
                  varianza = var(x, na.rm = TRUE),
                  deviazione_std = sd(x, na.rm = TRUE),
                  coeff_var = sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE),
                  skewness = skewness(x, na.rm = TRUE),
                  kurtosis = kurtosis(x, na.rm = TRUE)
                )
  }
 
}

# funzione per calcolare la frequenza per variabili categoriche
get_frequenze <- function(var) {
  dati %>%
    count({{var}}, name = "Frequenza") %>%
    mutate(Frequenza_percentuale = round(Frequenza / sum(Frequenza) * 100, 2))
}

#indice di GINI
gini.index <- function(x){
  ni = table(x)
  fi = ni/length(x)
  fi2 = fi^2
  j = length(table(x))
  
  gini = 1-sum(fi2)
  gini.normalizzato = gini/((j-1)/j)
  
  return(gini.normalizzato)
  
}

1 Raccolta dei Dati e Struttura del Dataset

L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature.

# carica il dataframe dati da un file csv
dati <- read.csv("neonati.csv")

# conta quanti record ha il dataframe
count_record = nrow(dati)

# conta numero di colonne del dataframe
count_column = length(dati)

# controlliamo se ci sono informazioni mancanti
info_mancanti <- sum(is.na(dati)) 

# tramite la libreria kableExtra mostriamo i primi 5 record del dataset
# per una esplorazione iniziale delle variabili
head(dati) %>%
  kbl(caption = "Prime righe dataset neonati",
      align = "c", 
      booktabs = TRUE,
      row.names = FALSE
      ) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Prime righe dataset neonati
Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
26	0	42	3380	490	325	Nat	osp3	M
21	2	39	3150	490	345	Nat	osp1	F
34	3	38	3640	500	375	Nat	osp2	M
28	1	41	3690	515	365	Nat	osp2	M
20	0	38	3700	480	335	Nat	osp3	F
32	0	40	3200	495	340	Nat	osp2	F

Il dataset è composto da 2500 osservazioni e da 10 variabili,

nel dataset ci sono 0 informazioni mancanti.

Le variabili del dataset sono:

Anni.madre: età della madre in anni (variabile quantitativa continua).
N.gravidanze: quante gravidanze ha avuto la madre (variabile quantitativa discreta).
Fumatrici: indicatore binario del fumo materno [0 = non fumatrice; 1 = fumatrice] (variabile qualitativa nominale dicotomica).
Gestazione: numero di settimane di gestazione (variabile quantitativa continua).
Peso: peso del neonato alla nascita in grammi (variabile quantitativa continua).
Lunghezza: lunghezza del neonato alla nascita in millimetri (variabile quantitativa continua).
Cranio: diametro del cranio del neonato alla nascita in millimetri (variabile quantitativa continua).
Tipo.parto: tipologia del parto [Nat = naturale; Ces = Cesareo] (variabile qualitativa nominale dicotomica).
Ospedale: ospedale di nascita [osp1 = ospedale 1; osp2 = Ospedale 2; osp3 = Ospedale 3] (variabile qualitativa nominale politomica).
Sesso: sesso del neonato [F = femmina; M = maschio] (variabile qualitativa nominale dicotomica).

2 Analisi e Modellizazione

Analisi Preliminare

Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.Inoltre si saggeranno le seguenti ipotesi con i test adatti:

in alcuni ospedali si fanno più parti cesarei.

La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione

Le misure antropometriche sono significativamente diverse tra i due sessi

Variabili Quantitative Discrete

# viene utilizzata la funzione calcola_indici per calcolare gli indici per la # variabile qiantitativa discreta N_gravidanze
riassunto_N_gravidanze <- dati %>%
  select(where(is.numeric)) %>%
  map_dfr(calcola_indici, .id = "variabile") %>%
  filter(variabile == "N.gravidanze")

riassunto_N_gravidanze %>%
  kbl(caption = "Riepilogo statistico della variabile quantitativa discreta (N.gravidanze)",
      align = "c", 
      booktabs = TRUE,
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Riepilogo statistico della variabile quantitativa discreta (N.gravidanze)
variabile	media	mediana	moda	minimo	massimo	Q1	Q3	varianza	deviazione_std	coeff_var	skewness	kurtosis
N.gravidanze	0.9812	1	0	0	12	0	1	1.639903	1.280587	1.305123	2.512746	10.97822

In media, il numero di gravidanze per madre è di 0.98.

La maggior parte delle madri non ha avuto gravidanze in precedenza o ne ha avuta una, il valore più frequente è 0. (donne alla prima gravidanza).

La deviazione standard tuttavia è 1.28 il che indica indica una variabilità piuttosto elevata rispetto alla media; nel dataset infatti sebbene ci siano molte osservazioni con la variabile N.gravidanze = 0 sono presenti anche delle osservazioni con valori di N.gravidanze oltre 10 che aumentano la deviazione standard.

La skewness pari a 2.51 evidenzia una forte asimmetria positiva, con molte osservazioni concentrate su valori bassi e una lunga coda verso destra (ci sono valori molto alti, rari ma presenti).

La kurtosis pari a 10.97 suggerisce che la distribuzione è molto appuntita rispetto alla normale, quindi con valori estremi rispetto alla normale.

Infine, il valore massimo di 12 gravidanze sembra anomalo per tale motivo più avanti (indicata con *) verrà effettuata una verifica (potrebbe rappresentare un outlier o un errore di registrazione).

# Istogramma del numero di gravidanze per madre

ggplot(dati, aes(x = N.gravidanze)) +
  geom_histogram(
    binwidth = 1,              # ampiezza delle classi (1 gravidanza per intervallo)
    fill = "skyblue",          # colore delle barre
    color = "black",           # bordo delle barre
    alpha = 0.7                # trasparenza
  ) +
  labs(
    title = "Distribuzione del numero di gravidanze per madre",
    x = "Numero di gravidanze",
    y = "Frequenza"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"), 
    axis.title = element_text(face = "bold")
  )

Approfondimento sulla variabile N.gravidanze (*)

Nel controllo dei dati sono emerse delle osservazioni con il valore della variabile N.gravidanze uguale a 12, biologicamente è possibile, tuttavia essendo un valore estremo è stata effettuata una verifica: è stata estratta una tabella con le osservazioni che hanno la variabile N.gravidanze maggiore di 6 e creato un boxplot.

# Filtra il dataset per N.gravidanze > 6
dati_gravidanze6 <- dati %>%
  filter(N.gravidanze > 6)

# Visualizza il risultato come tabella formattata
dati_gravidanze6 %>%
  kbl(caption = "Osservazioni con N.gravidanze > 6",
      align = "c",
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 13)

Osservazioni con N.gravidanze > 6
Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
36	8	39	3610	500	351	Ces	osp1	M
35	9	42	3760	540	348	Nat	osp2	F
30	8	40	3850	518	340	Nat	osp3	F
40	8	38	3520	470	341	Nat	osp3	M
30	7	35	2220	470	316	Nat	osp3	M
33	11	43	3400	475	360	Nat	osp1	M
38	12	39	3350	490	344	Nat	osp2	M
36	8	41	3730	480	335	Nat	osp3	M
30	8	39	2860	490	337	Ces	osp2	F
36	8	36	2860	460	334	Nat	osp2	F
35	9	37	3150	490	335	Nat	osp2	M
26	8	40	3250	500	355	Nat	osp2	M
37	8	28	930	355	235	Nat	osp1	F
35	10	39	2950	495	335	Nat	osp1	F
33	10	40	3090	485	353	Nat	osp3	M
34	10	38	2880	470	345	Ces	osp2	M

# Boxplot per il Controllo variabile N.gravidanze
ggplot(dati, aes(x = "", y = N.gravidanze)) +  
  geom_boxplot(fill = "lightblue", color = "darkblue", outlier.color = "red", outlier.shape = 16) +
  labs(
    title = "Controllo della variabile N.gravidanze",
    y = "Numero di gravidanze",
    x = NULL
  ) +
  theme_minimal(base_size = 13)

Il boxplot evidenzia la presenza di alcuni valori outlier nella variabile N.gravidanze; anche la tabella conferma 16 osservazioni con un numero di gravidanze superiore a 6. Tuttavia, osservando la variabile Anni.madre, si nota che tutte le madri con più di 6 gravidanze hanno un’età superiore ai 30 anni, rendendo quindi questi valori plausibili.

Variabili Quantitative Continue

# viene utilizzata la funzione calcola_indici per calcolare gli indici
riassunto_statistico_continue <- dati %>%
  select(where(is.numeric)) %>%
  map_dfr(~calcola_indici(.x, variabile_continua = TRUE), .id = "variabile")


# dal riassunto statistico viene tolta la varibile quantitativa discreta N.gravidanze perchè già calcolate e la variabile quantitativa Fumatrici
riassunto_statistico_continue <- riassunto_statistico_continue %>%
  filter(!variabile %in% c("N.gravidanze","Fumatrici"))

# gli indici vengono arrotondati alla seconda cifra decimale
riassunto_statistico_continue <- riassunto_statistico_continue %>%
  mutate(across(where(is.numeric), ~round(., 2)))

riassunto_statistico_continue %>%
  kbl(caption = "Riepilogo statistico delle variabili quantitative continue",
      align = "c", 
      booktabs = TRUE,
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Riepilogo statistico delle variabili quantitative continue
variabile	media	mediana	minimo	massimo	Q1	Q3	varianza	deviazione_std	coeff_var	skewness	kurtosis
Anni.madre	28.16	28	0	46	25	32	27.81	5.27	0.19	0.04	0.38
Gestazione	38.98	39	25	43	38	40	3.49	1.87	0.05	-2.06	8.25
Peso	3284.08	3300	830	4930	2990	3620	275665.68	525.04	0.16	-0.65	2.03
Lunghezza	494.69	500	310	565	480	510	692.67	26.32	0.05	-1.51	6.48
Cranio	340.03	340	235	390	330	350	269.79	16.43	0.05	-0.78	2.94

Anni.madre:

Media = 28.16, Mediana = 28 -> l’età media delle madri è di circa 28 anni
Skewness è circa 0 -> la distribuzione è simmetrica.
Kurtosis = 0.38 -> una forma vicina alla normale con code leggere.
Minimo = 0 -> valore biologicamente impossibile, verrà eseguita una verifica più avanti (#).

Gestazione (in settimane):

Media = 38.98 -> nel range della durata biologica della gravidanza
Mediana = 39 -> nel range della durata biologica della gravidanza
Minimo = 25 -> valore estremo, verrà eseguita una verifica più avanti (+).
Skewness = 0.05 -> i dati sono distribuiti in modo equilibrato intorno alla media.
Kurtosis = 8.25 -> distribuzione simmetrica, ma con molti outlier e una concentrazione dei dati attorno alla media (appuntita con code pesanti).

Peso:

Media = 3284g -> il peso medio nel range per un neonato sano
Mediana = 3300g -> il peso medio nel range per un neonato sano.
Deviazione standard = 525g -> variabilità moderata
Skewness = -0.65 -> ci sono più neonati che pesano sopra la media, ma altri con peso molto inferiore alla media che spostano la coda verso sinistra.
Kurtosis = 2.03 -> distribuzione leggermente appiattita rispetto alla normale (pochi valori estremi)

Lunghezza:

Media = 49.47cm -> lunghezza nel range di un neonato sano.
Mediana = 50cm -> lunghezza nel range di un neonato sano.
Skewness = -1.51 -> curva asimmetrica a sinistra in quanto la maggior parte dei neonati ha una lunghezza leggermente sopra la media, ma ci sono alcuni neonati con lunghezze significativamente più basse.
Kurtosis = 6.48 -> distribuzione molto appuntita per la presenza di neonati con lunghezze significativamente più basse della media.

Cranio:

Media = 34cm -> circonferenza cranica nella norma.
Mediana = 34cm -> circonferenza cranica nella norma.
Skewness = -0.78 -> leggermente asimettrica verso sinistra (alcuni neonati con circonferenza cranica inferiore alla norma).
Kurtosis = 2,94 -> distribuzione quasi normale.

Approfondimento sulla variabile Anni.madre (#)

ggplot(dati, aes(y = Anni.madre)) +
  geom_boxplot(fill = "lightblue", color = "darkblue") +
  labs(
    title = "Distribuzione dell'età delle madri",
    y = "Età (anni)"
  ) +
  theme_minimal(base_size = 13)

ggplot(dati, aes(x = Anni.madre)) +
  geom_histogram(binwidth = 1, fill = "skyblue", color = "black", alpha = 0.7) +
  labs(
    title = "Istogramma dell'età delle madri",
    x = "Età (anni)",
    y = "Frequenza"
  ) +
  theme_minimal(base_size = 13)

Sia dalla distribuzione che dall’istogramma della variabile Anni.madre si evince che ci sono degli outlier con valori troppo bassi (ad esempio 0 anni).

In questi casi ci sono due possibili approcci:

eliminare le osservazioni errate dal dataset.
Sostituire i valori anomali con un valore plausibile, ad esempio la media o la mediana della variabile.

Si è scelta la seconda opzione (sostituire tutti i valori di Anni.madre minori uguali a 15 con il valore della mediana)

# Crea un nuovo dataset copia dell'originale per non modificare i dati iniziali

dati_corretto <- dati

# Calcola la mediana escludendo valori anomali (<=15)
mediana_valida <- median(dati_corretto$Anni.madre[dati_corretto$Anni.madre > 15], na.rm = TRUE)

# Sostituisci i valori 15 con la mediana
dati_corretto$Anni.madre[dati_corretto$Anni.madre <= 15] <- mediana_valida

# Controlla se ci sono ancora valori di Anni.madre <= 15
somma_anni_madre_minore_15 <- sum(dati_corretto$Anni.madre <= 15)

# viene utilizzata la funzione calcola_indici per calcolare gli indici per la # variabile qiantitativa corretta Anni.madre
riassunto_Anni_madre <- dati_corretto %>%
  select(where(is.numeric)) %>%
  map_dfr(calcola_indici, .id = "variabile") %>%
  filter(variabile == "Anni.madre")

riassunto_Anni_madre %>%
  kbl(caption = "Riepilogo statistico della variabile Anni.madre del dataset dati_corretto",
      align = "c", 
      booktabs = TRUE,
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Riepilogo statistico della variabile Anni.madre del dataset dati_corretto
variabile	media	mediana	moda	minimo	massimo	Q1	Q3	varianza	deviazione_std	coeff_var	skewness	kurtosis
Anni.madre	28.2344	28	30	16	46	25	32	26.52447	5.150191	0.1824084	0.195975	-0.1441597

# estrtazione indici statistici per il dataset originale per un confronto con il nuvovo dataset
riassunto_Anni_madre_dataset_originale <- dati %>%
  select(where(is.numeric)) %>%
  map_dfr(calcola_indici, .id = "variabile") %>%
  filter(variabile == "Anni.madre")

riassunto_Anni_madre_dataset_originale %>%
  kbl(caption = "Riepilogo statistico della variabile Anni.madre del dataset originale",
      align = "c", 
      booktabs = TRUE,
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Riepilogo statistico della variabile Anni.madre del dataset originale
variabile	media	mediana	moda	minimo	massimo	Q1	Q3	varianza	deviazione_std	coeff_var	skewness	kurtosis
Anni.madre	28.164	28	30	0	46	25	32	27.81063	5.273578	0.1872454	0.0427858	0.3777127

Come si può notare dalle tabelle di riepilogo statistico dei due dataset per la variabile Anni.madre gli indici maggiormente differenti sono:

Skewness -> è passato da 0.04 a 0.19.
Kurtosis -> è passato da 0.38 a -0.14.
Minimo -> è passato da 0 a 16.

A conferma delle modifiche apportate è stato rifatto il boxplot che conferma l’assenza di outlier nella parte inferiore della distribuzione.

# Boxplot per la correzione della variabile Anni.madre
ggplot(dati_corretto, aes(y = Anni.madre)) +
  geom_boxplot(fill = "lightblue", color = "darkblue") +
  labs(title = "Età della madre (dopo correzione)", y = "Età (anni)") +
  theme_minimal()

Approfondimento sulla variabile Gestazione (+)

#controllo della varibile Gestazione per valori estremi

Gestazione_uguale_25 <- dati_corretto %>%
  filter(Gestazione == 25)


Gestazione_uguale_25 %>%
  kbl(caption = "Osservazione con la variabile Gestazione uguale a 25",
      align = "c", 
      booktabs = TRUE,
      row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Osservazione con la variabile Gestazione uguale a 25
Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
25	2	0	25	900	325	253	Nat	osp3	F

La tabella mostra l’osservazione con Gestazione = 25 settimane. Nonostante sia al limite estremo, Peso, Lunghezza e circonferenza Cranica risultano coerenti con questa età gestazionale, quindi l’osservazione viene mantenuta nel dataset.

Variabili Qualitative

# usando la funzione get_frequenze calcola la frequenza assoluta e percentuale 
attach(dati)
frequenze_tipo_parto  <- get_frequenze(Tipo.parto)
frequenze_fumatirce   <- get_frequenze(Fumatrici)
frequenze_ospedale  <- get_frequenze(Ospedale)
frequenze_sesso   <- get_frequenze(Sesso)

#tramite le funzioni kbl e kable_styling knitr e kableExtra visualizza i dati un una tabella
frequenze_tipo_parto %>%
  kbl(caption = "Distribuzione di frequenza per la variabile Tipo.parto",
      align = "c", 
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Distribuzione di frequenza per la variabile Tipo.parto
Tipo.parto	Frequenza	Frequenza_percentuale
Ces	728	29.12
Nat	1772	70.88

frequenze_fumatirce %>%
  kbl(caption = "Distribuzione di frequenza per la variabile Fumatrici [0 = Non Fumatrici - 1 = Fumatrici]",
      align = "c", 
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Distribuzione di frequenza per la variabile Fumatrici [0 = Non Fumatrici - 1 = Fumatrici]
Fumatrici	Frequenza	Frequenza_percentuale
0	2396	95.84
1	104	4.16

frequenze_ospedale %>%
  kbl(caption = "Distribuzione di frequenza per la variabile Ospedale",
      align = "c", 
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Distribuzione di frequenza per la variabile Ospedale
Ospedale	Frequenza	Frequenza_percentuale
osp1	816	32.64
osp2	849	33.96
osp3	835	33.40

frequenze_sesso %>%
  kbl(caption = "Distribuzione di frequenza per la variabile Sesso",
      align = "c", 
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Distribuzione di frequenza per la variabile Sesso
Sesso	Frequenza	Frequenza_percentuale
F	1256	50.24
M	1244	49.76

Dai dati estratti dalle frequenze si evince:

il parto Naturale (Nat) ha una frequanza del 70.88%.
le madri Non Fumatrici (valore 0 della variabile) sono il 95.84%; questa variabile potrebbe avere scarsa rilevanza nella costruzione dei modelli di regressione a causa della bassa variabilità.
la variabile Ospedali ha una distribuzione uniforme con percentuali per i 3 ospedali di circa il 33%.
la variabile Sesso ha una distribuzione uniforme con circa il 50% per entrambi i sessi.

In alcuni ospedali si fanno più parti cesarei

-> Falso

# Calcola la percentuale di parti cesarei per ospedale
parti_ospedale <- dati_corretto %>%
  group_by(Ospedale, Tipo.parto) %>%
  summarise(Frequenza = n(), .groups = "drop_last") %>%
  mutate(Totale_ospedale = sum(Frequenza),
         Percentuale = round(Frequenza / Totale_ospedale * 100, 2))

# Mostra la tabella
parti_ospedale %>%
  kbl(caption = "Distribuzione percentuale dei tipi di parto per ospedale",
      align = "c", booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE, font_size = 14)

Distribuzione percentuale dei tipi di parto per ospedale
Ospedale	Tipo.parto	Frequenza	Totale_ospedale	Percentuale
osp1	Ces	242	816	29.66
osp1	Nat	574	816	70.34
osp2	Ces	254	849	29.92
osp2	Nat	595	849	70.08
osp3	Ces	232	835	27.78
osp3	Nat	603	835	72.22

# Grafico a barre
ggplot(parti_ospedale, aes(x = Ospedale, y = Percentuale, fill = Tipo.parto)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Percentuale di parti naturali e cesarei per ospedale",
       x = "Ospedale", y = "Percentuale (%)") +
  theme_minimal() +
  scale_fill_manual(values = c("steelblue", "lightgreen"))

In tutti e 3 gli ospedali la maggior parte dei parti avviene in modo naturale, con una percentuale che si aggira intorno al 70%.

La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione

-> Peso Vero

-> Lunghezza Falso

In Italia il peso medio dei neonati è di circa 3200 grammi - 3400 grammi e la lunghezza media e di circa 495 mm - 505 mm

#controllo dei dati medi del dataset con la media Italiana

media_italia_peso <- 3300     # in grammi
media_italia_lunghezza <- 500  # in mm

# Calcola le medie del dataset
media_campione_peso <- mean(dati_corretto$Peso, na.rm = TRUE)
media_campione_lunghezza <- mean(dati_corretto$Lunghezza, na.rm = TRUE)

# Crea un dataset per il grafico
medie <- data.frame(
  Variabile = c("Peso", "Lunghezza"),
  Campione = c(media_campione_peso, media_campione_lunghezza),
  Italia = c(media_italia_peso, media_italia_lunghezza)
)

# Grafico per  Peso
df_peso <- data.frame(
  Categoria = c("Campione", "Italia"),
  Media = c(media_campione_peso, media_italia_peso)
)

ggplot(df_peso, aes(x = Categoria, y = Media, fill = Categoria)) +
  geom_col(width = 0.5, show.legend = FALSE) +
  scale_fill_manual(values = c("Campione" = "orange", "Italia" = "lightblue")) +
  labs(title = "Confronto tra media Peso neonati (Campione vs Italia)",
       y = "Peso medio (grammi)", x = "") +
  theme_minimal()

# grafico per la lunghezza
df_lunghezza <- data.frame(
  Categoria = c("Campione", "Italia"),
  Media = c(media_campione_lunghezza, media_italia_lunghezza)
)

ggplot(df_lunghezza, aes(x = Categoria, y = Media, fill = Categoria)) +
  geom_col(width = 0.5, show.legend = FALSE) +
  scale_fill_manual(values = c("Campione" = "orange", "Italia" = "lightblue")) +
  labs(title = "Confronto tra media Lunghezza neonati (Campione vs Italia)",
       y = "Lunghezza media (mm)", x = "") +
  theme_minimal()

#  T-test per confrontare il campione con la media Italiana
t_peso <- t.test(dati_corretto$Peso, mu = media_italia_peso)
t_lunghezza <- t.test(dati_corretto$Lunghezza, mu = media_italia_lunghezza)

#stampa risultato t_peso
t_peso

## 
##  One Sample t-test
## 
## data:  dati_corretto$Peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.490 3304.672
## sample estimates:
## mean of x 
##  3284.081

La media del peso del campione è di 3284.081 grammi, il p_value del t-test è 0.1296 percui maggiore di 0.05 questo significa che non ci sono differenze statisticamente significative tra la media del campione e quella della popolazione italiana; l’intervallo di confidenza 95% va da 3263.5 grammi a 3304.7 grammi, che include la media nazionale 3300 grammi, confermando che il campione è compatibile con la popolazione.

#stampa t_lunghezza
t_lunghezza

## 
##  One Sample t-test
## 
## data:  dati_corretto$Lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  493.6598 495.7242
## sample estimates:
## mean of x 
##   494.692

La media della lunghezza del campione è 494.692 millimetri, il p_value del t-test è molto inferiore a 0.05 questo significa che la differenza tra il campione e la media italiana è statisticamente significativa. L’intervallo di confidenza di 95% va da 493.6598 millimetri a 495.7242 millimetri che non include la media italiana. Conclusione la lunghezza media del campione e significativamente più bassa della media italiana.

Le misure antropometriche sono significativamente diverse tra i due sessi

-> Vero

# t-test per le combinazioni Sesso - Peso - Lunghezza - Cranio

# Peso
t_peso_sesso <- t.test(Peso ~ Sesso, data = dati_corretto)
t_peso_sesso

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287.1051 -207.0615
## sample estimates:
## mean in group F mean in group M 
##        3161.132        3408.215

# Lunghezza
t_lunghezza_sesso <- t.test(Lunghezza ~ Sesso, data = dati_corretto)
t_lunghezza_sesso

## 
##  Welch Two Sample t-test
## 
## data:  Lunghezza by Sesso
## t = -9.582, df = 2459.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.929470  -7.876273
## sample estimates:
## mean in group F mean in group M 
##        489.7643        499.6672

# Cranio
t_cranio_sesso <- t.test(Cranio ~ Sesso, data = dati_corretto)
t_cranio_sesso

## 
##  Welch Two Sample t-test
## 
## data:  Cranio by Sesso
## t = -7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.089912 -3.541270
## sample estimates:
## mean in group F mean in group M 
##        337.6330        342.4486

L’analisi del t-test Peso/Sesso indica che la differenza di Peso tra Femmine e Maschi è statisticamente significativa (p-value < 0.05)
L’analisi del t-test Lunghezza/Sesso indica che la differenza di Lunghezza tra Femmine e Maschi è statisticamente significativa (p-value < 0.05)
L’analisi del t-test Cranio/Sesso indica che la differenza di grandezza del Cranio tra Femmine e Maschi è statisticamente significativa (p-value < 0.05)

#tabella per la comparazione delle medie antropometriche per Sesso
medie_sesso <- dati_corretto %>%
  group_by(Sesso) %>%
  summarise(
    Peso = mean(Peso, na.rm = TRUE),
    Lunghezza = mean(Lunghezza, na.rm = TRUE),
    Cranio = mean(Cranio, na.rm = TRUE)
  )

medie_sesso %>%
  kbl(caption = "Medie antropometriche per sesso",
      align = "c",
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Medie antropometriche per sesso
Sesso	Peso	Lunghezza	Cranio
F	3161.132	489.7643	337.6330
M	3408.215	499.6672	342.4486

In conclusione, anche guardando la tabella precedente, si può affermare che tutte e tre le misure antropometriche (Peso, Lunghezza e Cranio) sono significativamente maggiori nei maschi rispetto alle femmine, confermando le differenze fisiologiche tra i sessi nei neonati.

Controllo Relazione tra Fumatrici/Gestazione

# Calcola statistiche per Gestazione in base al fumo materno
gestazione_fumo <- dati_corretto %>%
  group_by(Fumatrici) %>%
  summarise(
    Media = mean(Gestazione, na.rm = TRUE),
    Mediana = median(Gestazione, na.rm = TRUE),
    Deviazione_std = sd(Gestazione, na.rm = TRUE),
    Numerosita = n()
  )

# Mostra la tabella con kbl
gestazione_fumo %>%
  kbl(caption = "Durata media della gravidanza in base al fumo materno (0 = non fumatrici, 1 = fumatrici)",
      align = "c",
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Durata media della gravidanza in base al fumo materno (0 = non fumatrici, 1 = fumatrici)
Fumatrici	Media	Mediana	Deviazione_std	Numerosita
0	38.96786	39	1.884840	2396
1	39.26923	40	1.422638	104

# t-test gestazione/fumo materno
t_test_gestazione <- t.test(Gestazione ~ Fumatrici, data = dati_corretto)
t_test_gestazione

## 
##  Welch Two Sample t-test
## 
## data:  Gestazione by Fumatrici
## t = -2.0824, df = 119.26, p-value = 0.03944
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -0.58791720 -0.01481813
## sample estimates:
## mean in group 0 mean in group 1 
##        38.96786        39.26923

ggplot(dati_corretto, aes(x = factor(Fumatrici), y = Gestazione, fill = factor(Fumatrici))) +
  geom_boxplot() +
  scale_fill_manual(values = c("0" = "lightblue", "1" = "salmon"), labels = c("Non Fumatrici", "Fumatrici")) +
  labs(title = "Durata della gravidanza in base al fumo materno",
       x = "Fumatrici",
       y = "Settimane di Gestazione") +
  theme_minimal()

La media della durata della Gestazione è minore per le non fumatrici:

Non Fumatrici -> 38.97 settimane
Fumatrici -> 39.27 settimane

Il p-value = 0.0394, il che indica una differenza statisticamente significativa. L’intervallo di confidenza al 95% per la differenza tra le medie va da -0.588 a -0.015 settimane, non include 0, confermando la significatività della differenza.

In conclusione nel campione analizzato le madri fumatrici tendono ad avere una gestazione più lunga delle madri non fumatrici.

Controllo relazione Fumatrici/Tipo Parto/Gestazione

#controllo relazione tra Fumatrici- tipo parto -gestazione

# Calcola statistiche e percentuale
gestazione_parto_fumo <- dati_corretto %>%
  group_by(Fumatrici, Tipo.parto) %>%
  summarise(
    Media = mean(Gestazione, na.rm = TRUE),
    Mediana = median(Gestazione, na.rm = TRUE),
    Deviazione_std = sd(Gestazione, na.rm = TRUE),
    Numerosita = n(),
    .groups = "drop"
  ) %>%
  group_by(Fumatrici) %>%
  mutate(Percentuale_gruppo = round(Numerosita / sum(Numerosita) * 100, 2)) %>%
  ungroup()

# Mostra la tabella con kbl
gestazione_parto_fumo %>%
  kbl(caption = "Durata media della gravidanza in base a fumo materno e tipo di parto con percentuale sul totale del gruppo",
      align = "c",
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Durata media della gravidanza in base a fumo materno e tipo di parto con percentuale sul totale del gruppo
Fumatrici	Tipo.parto	Media	Mediana	Deviazione_std	Numerosita	Percentuale_gruppo
0	Ces	39.01994	39.0	1.737281	702	29.3
0	Nat	38.94628	39.0	1.942789	1694	70.7
1	Ces	39.11538	39.5	1.840151	26	25.0
1	Nat	39.32051	40.0	1.263818	78	75.0

Osservazioni generali:

in entrambi i gruppi la numerosità del parto naturale è maggiore del parto cesareo, la durata della gestazione è maggiore nelle madri fumatrici.

Analisi variabile Peso (variabile dipendente)

# creazione del grafico per la distribuzione della densità della variabile Peso
mu_vero <- mean(dati_corretto$Peso, na.rm = TRUE)

ggplot(dati_corretto, aes(x = Peso)) +
  geom_density(fill = "lightblue", alpha = 0.6) +
  geom_vline(xintercept = mu_vero, color = "red", linetype = "dashed", size = 1) +
  labs(title = "Distribuzione di densità del Peso dei neonati",
       x = "Peso (grammi)",
       y = "Densità") +
  theme_minimal()

La distribuzione di densità del peso mostrata nel grafico appare approssimativamente normale con una forma a campana, si nota una leggera asimmetria negativa (coda più lunga a sinistra), coerente con il coefficiente di asimmetria -0.65 riportato nelle statistiche descrittive. Questo suggerisce che la maggior parte dei neonati hanno peso vicino alla media però esistono dei casi con peso inferiore.

Codifica delle variabili Sesso - Tipo.parto - Ospedali

La codifica viene effettuata perché i modelli lineari trattano più facilmente le variabili numeriche rispetto a quelle categoriali.

Sesso -> Femmina = 0; Maschio = 1
Tipo.parto -> Nat = 0; Ces = 1
La variabile Ospedale è stata suddivisa in tre variabili Ospedaleosp1 - Ospedaleosp2 - Ospedaleosp3

# Trasforma Sesso in 0-1 -> F = 0, M = 1
dati_corretto$Sesso <- ifelse(dati_corretto$Sesso == "M", 1, 0)

# Trasforma Tipo.parto in 0-1 -> Nat = 0, Ces = 1
dati_corretto$Tipo.parto <- ifelse(dati_corretto$Tipo.parto == "Nat", 0, 1)

# trasforma la variabile Ospedale in 3 variabili Ospedaleosp1 - Ospedaleosp2 - Ospedaleosp3
dati_dummy <- model.matrix(~ Ospedale - 1, data = dati_corretto)
dati_dummy <- as.data.frame(dati_dummy)
dati_corretto <- cbind(dati_corretto[ , !names(dati_corretto) %in% "Ospedale"], dati_dummy)

head(dati_corretto) %>%
  kbl(caption = "Prime righe del dataset con variabili trasformate",
      align = "c", 
      booktabs = TRUE,
      row.names = FALSE
      ) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Prime righe del dataset con variabili trasformate
Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Sesso	Ospedaleosp1	Ospedaleosp2	Ospedaleosp3
26	0	42	3380	490	325	1	0	0	1
21	2	39	3150	490	345	0	1	0	0
34	3	38	3640	500	375	1	0	1	0
28	1	41	3690	515	365	1	0	1	0
20	0	38	3700	480	335	0	0	0	1
32	0	40	3200	495	340	0	0	1	0

Analisi delle correlazioni tra variabili

# variabile Peso vs altre variabili del dataset per tipologia, ho utilizzato la libreria gridExtra per formattare i grafici nel documento

# Variabili qualitative: boxplot
p1 <- ggplot(dati_corretto, aes(x = Ospedale, y = Peso)) +
  geom_boxplot(fill = "lightblue") +
  labs(title = "Peso per Ospedale") +
  theme_minimal()

p2 <- ggplot(dati_corretto, aes(x = Tipo.parto, y = Peso)) +
  geom_boxplot(fill = "lightgreen") +
  labs(title = "Peso per Tipo.parto") +
  theme_minimal()

p3 <- ggplot(dati_corretto, aes(x = Fumatrici, y = Peso)) +
  geom_boxplot(fill = "lightpink") +
  labs(title = "Peso per Fumatrici") +
  theme_minimal()

p4 <- ggplot(dati_corretto, aes(x = Sesso, y = Peso)) +
  geom_boxplot(fill = "lightyellow") +
  labs(title = "Peso per Sesso") +
  theme_minimal()

# Variabili quantitative: scatterplot con linea di regressione
p5 <- ggplot(dati_corretto, aes(x = Gestazione, y = Peso)) +
  geom_point(color = "steelblue", alpha = 0.6) +
  geom_smooth(method = "lm", col = "red") +
  labs(title = "Peso vs Gestazione") +
  theme_minimal()

p6 <- ggplot(dati_corretto, aes(x = Anni.madre, y = Peso)) +
  geom_point(color = "darkgreen", alpha = 0.6) +
  geom_smooth(method = "lm", col = "red") +
  labs(title = "Peso vs Anni.madre") +
  theme_minimal()

p7 <- ggplot(dati_corretto, aes(x = N.gravidanze, y = Peso)) +
  geom_point(color = "purple", alpha = 0.6) +
  geom_smooth(method = "lm", col = "red") +
  labs(title = "Peso vs N.gravidanze") +
  theme_minimal()

p8 <- ggplot(dati_corretto, aes(x = Lunghezza, y = Peso)) +
  geom_point(color = "orange", alpha = 0.6) +
  geom_smooth(method = "lm", col = "red") +
  labs(title = "Peso vs Lunghezza") +
  theme_minimal()

p9 <- ggplot(dati_corretto, aes(x = Cranio, y = Peso)) +
  geom_point(color = "brown", alpha = 0.6) +
  geom_smooth(method = "lm", col = "red") +
  labs(title = "Peso vs Cranio") +
  theme_minimal()

# Mostra i grafici 3x3

grid.arrange(p1, p2, p3, p4, p5, p6, p7, p8, p9, ncol = 3)

Panoramica generale delle relazioni tra tutte le variabili

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  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 * r)
}

variabili_da_controllare <- dati_corretto %>% 
  select(Peso, Gestazione, Anni.madre, N.gravidanze, Lunghezza, Cranio)

pairs(variabili_da_controllare,upper.panel = panel.smooth, lower.panel = panel.cor)

Commenti grafico:

Peso e Lunghezza -> correlazione molto alta (0.80)
Peso e Cranio -> correlazione alta (0.70)
Peso e Gestazione -> correlazione moderata (0.59)
Peso e N.gravidanze - Anni.madre -> correlazioni deboli (circa 0.1 - 0.2)
Lunghezza e Cranio -> correlazione moderata (0.60)
Gestazione e Lunghezza -> correlazione moderata (0.62)

In conclusione, le variabili fisiche del neonato (Peso/Lunghezza/Cranio) sono fortemente correlate tra loro; la Gestazione ha una correlazione moderata con le variabili fisiche del neonato, più settimane di gestazione corrispondono ad una maggiore grandezza fisica del neonato; le varibili materne (Anni.madre/N.gravidanze) hanno una correlazione debole con la variabile Peso quindi probabilmente non sono buoni predittori diretti del peso.

# Tabella di contingenza
tabella <- table(dati$Fumatrici, dati$Tipo.parto)

# Test Chi-quadrato
test_indipendenza <- chisq.test(tabella)


# Grafico delle proporzioni per visualizzare la relazione
tabella_df <- as.data.frame(tabella)
colnames(tabella_df) <- c("Fumatrici", "TipoParto", "Frequenza")
 
ggplot(tabella_df, aes(x = as.factor(Fumatrici),
                       y = Frequenza,
                       fill = as.factor(TipoParto))) +
  geom_bar(stat = "identity", position = "fill") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(title = "Distribuzione percentuale del Tipo di Parto per Fumatrici",
       x = "Fumatrici (0 = Non Fumatrice, 1 = Fumatrice)",
       y = "Percentuale",
       fill = "Tipo di Parto") +
  theme_minimal(base_size = 14)

# Risultato test indipendenza Fumatrici - Tipo.parto
test_indipendenza

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tabella
## X-squared = 0.69629, df = 1, p-value = 0.404

Il p-value del test di indipendenza è uguale a 0.4 -> non c’è relazione significativa tra le variabili.

Relazione variabile Tipo.parto e Fumatrici con la variabile Peso

Relazione Tipo.parto - Peso

t.test(Peso ~ Tipo.parto, data = dati_corretto)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Tipo.parto
## t = 0.12968, df = 1493, p-value = 0.8968
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -40.54037  46.27992
## sample estimates:
## mean in group 0 mean in group 1 
##        3284.916        3282.047

La differenza media è di circa 2.9 grammi e il p-value di 0.89 -> non c’è relazione significativa tra le variabili.

Relazione Fumatrici - Peso

t.test(Peso ~ Fumatrici, data = dati_corretto)

## 
##  Welch Two Sample t-test
## 
## data:  Peso by Fumatrici
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -45.61354 145.22674
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.153        3236.346

La differenza media è di circa 49 grammi e il p-value di 0.3 -> non c’è relazione significativa tra le variabili.

Creazione del Modello di Regressione

Dalle analisi eseguite per questo campione le variabili correlate al Peso risultano essere:

Lunghezza, Gestazione, circonferenza Cranio, Sesso

Primo modello

variabili del modello -> Gestazione, Lunghezza, Cranio, Sesso, Fumatrici, Tipo.parto

# Modello di regressione lineare con tutte le variabili
modello_primo <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso + Fumatrici + Tipo.parto, data = dati_corretto)

# Riepilogo del modello
summary(modello_primo)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     Fumatrici + Tipo.parto, data = dati_corretto)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1118.88  -185.33   -16.16   161.80  2622.90 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6646.6168   135.4003 -49.089  < 2e-16 ***
## Gestazione     31.4067     3.7882   8.291  < 2e-16 ***
## Lunghezza      10.2276     0.3011  33.969  < 2e-16 ***
## Cranio         10.6379     0.4242  25.075  < 2e-16 ***
## Sesso          79.2477    11.2025   7.074 1.95e-12 ***
## Fumatrici     -27.5245    27.5783  -0.998   0.3184    
## Tipo.parto    -29.3167    12.1132  -2.420   0.0156 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2493 degrees of freedom
## Multiple R-squared:  0.7268, Adjusted R-squared:  0.7262 
## F-statistic:  1106 on 6 and 2493 DF,  p-value: < 2.2e-16

R-squared = 0.7278. Valore ottimale spiega il 72.78% della variabilità del peso del neonato.

p-value < 0.0001: il modello è statisticamente significativo.

Secondo modello

Variabili del modello -> Gestazione, Lunghezza, Cranio, Sesso, Tipo.parto

# Modello di regressione lineare 
modello_secondo <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso + Tipo.parto, data = dati_corretto)

# Riepilogo del modello
summary(modello_secondo)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     Tipo.parto, data = dati_corretto)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1118.43  -185.56   -16.07   161.53  2626.18 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6646.7022   135.4002 -49.089  < 2e-16 ***
## Gestazione     31.1917     3.7821   8.247 2.59e-16 ***
## Lunghezza      10.2412     0.3008  34.049  < 2e-16 ***
## Cranio         10.6397     0.4242  25.080  < 2e-16 ***
## Sesso          79.0670    11.2010   7.059 2.17e-12 ***
## Tipo.parto    -29.1062    12.1113  -2.403   0.0163 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2494 degrees of freedom
## Multiple R-squared:  0.7267, Adjusted R-squared:  0.7262 
## F-statistic:  1326 on 5 and 2494 DF,  p-value: < 2.2e-16

R-squared = 0.7267. Valore ottimale spiega il 72.67% della variabilità del peso del neonato.

p-value < 0.0001: il modello è statisticamente significativo.

Terzo modello

Variabili del modello -> Gestazione, Lunghezza, Cranio, Sesso

# Modello di regressione lineare
modello_terzo <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso, data = dati_corretto)

# Riepilogo del modello
summary(modello_terzo)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso, 
##     data = dati_corretto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.2  -184.3   -17.6   163.3  2627.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6651.1188   135.5172 -49.080  < 2e-16 ***
## Gestazione     31.2737     3.7856   8.261 2.31e-16 ***
## Lunghezza      10.2054     0.3007  33.939  < 2e-16 ***
## Cranio         10.6704     0.4245  25.139  < 2e-16 ***
## Sesso          79.1049    11.2117   7.056 2.22e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared:  0.7261, Adjusted R-squared:  0.7257 
## F-statistic:  1654 on 4 and 2495 DF,  p-value: < 2.2e-16

R-squared = 0.7261. Valore ottimale spiega il 72.61% della variabilità del peso del neonato.

p-value < 0.0001: il modello è statisticamente significativo.

Quarto modello

Variabili del modello -> Gestazione, Lunghezza, Cranio, Fumatrici

# Modello di regressione lineare
modello_quarto <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Fumatrici, data = dati_corretto)

# Riepilogo del modello
summary(modello_quarto)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Fumatrici, 
##     data = dati_corretto)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1106.20  -183.72   -13.46   167.10  2621.03 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6777.3182   135.6503 -49.962   <2e-16 ***
## Gestazione     31.8727     3.8284   8.325   <2e-16 ***
## Lunghezza      10.4140     0.3023  34.449   <2e-16 ***
## Cranio         10.7880     0.4283  25.190   <2e-16 ***
## Fumatrici     -23.2035    27.8670  -0.833    0.405    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 277.7 on 2495 degrees of freedom
## Multiple R-squared:  0.7207, Adjusted R-squared:  0.7203 
## F-statistic:  1610 on 4 and 2495 DF,  p-value: < 2.2e-16

R-squared= 0.7207. Valore ottimale spiega il 72.07% della variabilità del peso del neonato.

p-value < 0.0001: il modello è statisticamente significativo.

Quinto modello

Variabili dell modello -> Gestazione, Lunghezza, N.gravidanze, Sesso

# Modello di regressione lineare
modello_quinto <- lm(Peso ~ Gestazione + Lunghezza + N.gravidanze + Sesso, data = dati_corretto)

# Riepilogo del modello
summary(modello_quinto)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + N.gravidanze + Sesso, 
##     data = dati_corretto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1159.6  -193.5   -17.9   185.5  3619.4 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -5314.7941   138.3023 -38.429  < 2e-16 ***
## Gestazione      46.1998     4.1907  11.024  < 2e-16 ***
## Lunghezza       13.6061     0.2992  45.468  < 2e-16 ***
## N.gravidanze    23.8752     4.8141   4.959 7.55e-07 ***
## Sesso           87.8954    12.4896   7.037 2.52e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 306.4 on 2495 degrees of freedom
## Multiple R-squared:  0.6601, Adjusted R-squared:  0.6595 
## F-statistic:  1211 on 4 and 2495 DF,  p-value: < 2.2e-16

R-squared = 0.6601. Valore non tanto buono spiega il 66.01% della variabilità del peso del neonato.

p-value < 0.0001: il modello è statisticamente significativo.

Sesto modello

Variabili del modello -> Gestazione * Sesso (interazione tra le due variabili), Lunghezza, Cranio, Tipo.parto

# Modello di regressione lineare
# Interazione Gestazione - Sesso
modello_sesto <- lm(Peso ~ Gestazione * Sesso + Lunghezza + Cranio + Tipo.parto, data = dati_corretto)


# Riepilogo del modello
summary(modello_sesto)

## 
## Call:
## lm(formula = Peso ~ Gestazione * Sesso + Lunghezza + Cranio + 
##     Tipo.parto, data = dati_corretto)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1113.2  -184.4   -17.0   163.0  2624.6 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -6528.7502   166.7057 -39.163  < 2e-16 ***
## Gestazione          28.0844     4.5681   6.148 9.12e-10 ***
## Sesso             -205.5068   234.9380  -0.875    0.382    
## Lunghezza           10.2459     0.3008  34.065  < 2e-16 ***
## Cranio              10.6404     0.4242  25.083  < 2e-16 ***
## Tipo.parto         -29.4752    12.1140  -2.433    0.015 *  
## Gestazione:Sesso     7.2921     6.0134   1.213    0.225    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2493 degrees of freedom
## Multiple R-squared:  0.7269, Adjusted R-squared:  0.7262 
## F-statistic:  1106 on 6 and 2493 DF,  p-value: < 2.2e-16

R-squared = 0.7269. Valore ottimale spiega il 72.69% della variabilità del peso del neonato.

p-value < 0.0001: il modello è statisticamente significativo.

Il modello include un’interazione tra Gestazione e Sesso, quindi l’effetto della durata della gestazione sul peso del neonato può variare in base al sesso, e viceversa.

Settimo modello

Variabili del modello -> Gestazione, Lunghezza, Cranio, Sesso, Fumatrici * Tipo.parto (interazione tra le due variabili)

# Modello di regressione lineare
# Interazione Fumatrici - Tipo.parto
modello_settimo <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso + Fumatrici * Tipo.parto, data = dati_corretto)


# Riepilogo del modello
summary(modello_settimo)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     Fumatrici * Tipo.parto, data = dati_corretto)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1115.89  -185.28   -16.71   163.00  2621.06 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -6646.7661   135.3665 -49.102  < 2e-16 ***
## Gestazione              31.5295     3.7882   8.323  < 2e-16 ***
## Lunghezza               10.2228     0.3010  33.960  < 2e-16 ***
## Cranio                  10.6347     0.4241  25.073  < 2e-16 ***
## Sesso                   78.9235    11.2018   7.046 2.38e-12 ***
## Fumatrici              -51.4850    31.8728  -1.615  0.10637    
## Tipo.parto             -32.9052    12.3447  -2.666  0.00774 ** 
## Fumatrici:Tipo.parto    95.0585    63.4385   1.498  0.13415    
## ---
## 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.7271, Adjusted R-squared:  0.7263 
## F-statistic: 948.4 on 7 and 2492 DF,  p-value: < 2.2e-16

R-squared = 0.7271. Valore ottimale spiega il 72.71% della variabilità del peso del neonato.

p-value < 0.0001: il modello è statisticamente significativo.

Il modello include un’interazione tra Fumatrici e Tipo.parto, quindi l’effetto del fumo materno sul peso del neonato può variare a seconda del tipo di parto, e viceversa.

Ottavo modello

Variabili del modello -> Gestazione, Lunghezza * Cranio (interazione tra le due variabili), Sesso, Tipo.parto

# Modello di regressione lineare
# Interazione Lunghezza + Cranio (due variabili morfometriche)
modello_ottavo <- lm(Peso ~ Gestazione + Lunghezza * Cranio + Sesso + Tipo.parto, data = dati_corretto)


# Riepilogo del modello
summary(modello_ottavo)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza * Cranio + Sesso + 
##     Tipo.parto, data = dati_corretto)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1120.29  -184.11   -15.08   159.37  2846.23 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.896e+03  1.019e+03  -1.861   0.0629 .  
## Gestazione        3.677e+01  3.948e+00   9.313  < 2e-16 ***
## Lunghezza        -3.550e-02  2.205e+00  -0.016   0.9872    
## Cranio           -4.247e+00  3.192e+00  -1.330   0.1835    
## Sesso             7.451e+01  1.120e+01   6.655 3.46e-11 ***
## Tipo.parto       -2.761e+01  1.206e+01  -2.288   0.0222 *  
## Lunghezza:Cranio  3.074e-02  6.533e-03   4.705 2.67e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.6 on 2493 degrees of freedom
## Multiple R-squared:  0.7291, Adjusted R-squared:  0.7285 
## F-statistic:  1118 on 6 and 2493 DF,  p-value: < 2.2e-16

R-squared = 0.7291. Valore ottimale spiega il 72.91% della variabilità del peso del neonato.

p-value < 0.0001: il modello è statisticamente significativo.

Il modello include un’interazione tra Lunghezza e Cranio, quindi l’effetto della lunghezza sul peso del neonato può variare a seconda della circonferenza del Cranio, e viceversa.

Nono modello

Variabili del modello -> Gestazione, I(Gestazione^2) [relazione non lineare], Lunghezza, Cranio, Sesso, Tipo.parto

# Modello di regressione lineare
# Effetto non lineare di Gestazione
modello_nono <- lm(Peso ~ Gestazione + I(Gestazione^2) + Lunghezza + Cranio + Sesso + Tipo.parto, data = dati_corretto)


# Riepilogo del modello
summary(modello_nono)

## 
## Call:
## lm(formula = Peso ~ Gestazione + I(Gestazione^2) + Lunghezza + 
##     Cranio + Sesso + Tipo.parto, data = dati_corretto)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1120.27  -184.61   -15.21   164.47  2648.27 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -4676.0064   898.7752  -5.203 2.12e-07 ***
## Gestazione        -78.8252    49.7474  -1.585   0.1132    
## I(Gestazione^2)     1.4686     0.6622   2.218   0.0267 *  
## Lunghezza          10.3420     0.3040  34.024  < 2e-16 ***
## Cranio             10.7344     0.4260  25.195  < 2e-16 ***
## Sesso              76.9383    11.2333   6.849 9.32e-12 ***
## Tipo.parto        -28.6370    12.1037  -2.366   0.0181 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.5 on 2493 degrees of freedom
## Multiple R-squared:  0.7273, Adjusted R-squared:  0.7266 
## F-statistic:  1108 on 6 and 2493 DF,  p-value: < 2.2e-16

R-squared = 0.7273. Valore ottimo che spiega il 72.73% della variabilità del peso del neonato.

p-value < 0.0001: il modello è statisticamente significativo.

Il modello include un termine quadratico di Gestazione, quindi considera che l’effetto della durata della gestazione sul peso del neonato può cambiare in modo non lineare.

Decimo modello

Variabili del modello -> Gestazione, Lunghezza, I(Lunghezza^2) [relazione non lineare], Cranio, Sesso

# Modello di regressione lineare
# Effetto non lineare di Lunghezza
modello_decimo <- lm(Peso ~ Gestazione + Lunghezza + I(Lunghezza^2) + Cranio + Sesso, data = dati_corretto)


# Riepilogo del modello
summary(modello_decimo)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + I(Lunghezza^2) + 
##     Cranio + Sesso, data = dati_corretto)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1156.75  -182.00   -12.52   166.67  1783.83 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    156.673587 724.783055   0.216    0.829    
## Gestazione      41.174410   3.860512  10.666  < 2e-16 ***
## Lunghezza      -19.913124   3.165788  -6.290 3.73e-10 ***
## I(Lunghezza^2)   0.031241   0.003269   9.555  < 2e-16 ***
## Cranio          10.795854   0.417182  25.878  < 2e-16 ***
## Sesso           71.369249  11.043854   6.462 1.24e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 270.2 on 2494 degrees of freedom
## Multiple R-squared:  0.7358, Adjusted R-squared:  0.7352 
## F-statistic:  1389 on 5 and 2494 DF,  p-value: < 2.2e-16

R-squared = 0.7358. Valore ottimo spiega il 73.58% della variabilità del peso del neonato.

p-value < 0.0001: il modello è statisticamente significativo.

Il modello include un termine quadratico di Lunghezza, quindi considera che l’effetto della lunghezzae sul peso del neonato può cambiare in modo non lineare

Selezione del modello ottimale

Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC)

AIC(modello_primo, modello_secondo, modello_terzo, modello_quarto, modello_quinto,modello_sesto,modello_settimo,modello_ottavo,modello_nono,modello_decimo)

##                 df      AIC
## modello_primo    8 35182.82
## modello_secondo  7 35181.82
## modello_terzo    6 35185.60
## modello_quarto   6 35234.30
## modello_quinto   6 35725.52
## modello_sesto    8 35182.34
## modello_settimo  9 35182.57
## modello_ottavo   8 35161.71
## modello_nono     8 35178.89
## modello_decimo   7 35097.71

Il modello con AIC più basso è il modello_decimo con un valore di 35097.71.

BIC(modello_primo, modello_secondo, modello_terzo, modello_quarto, modello_quinto,modello_sesto,modello_settimo,modello_ottavo,modello_nono,modello_decimo)

##                 df      BIC
## modello_primo    8 35229.41
## modello_secondo  7 35222.59
## modello_terzo    6 35220.54
## modello_quarto   6 35269.24
## modello_quinto   6 35760.46
## modello_sesto    8 35228.94
## modello_settimo  9 35234.98
## modello_ottavo   8 35208.31
## modello_nono     8 35225.48
## modello_decimo   7 35138.48

Il modello con BIC più basso è il modello_decimo con un valore di 35138.48.

Conclusione:

In base sia al criterio di informazione di Akaike (AIC) che a quello di Bayes (BIC), il modello_decimo risulta il più parsimonioso, confermandosi come la scelta ottimale per rappresentare la relazione tra il peso del neonato e le variabili esplicative considerate.

Analisi della qualità del modello

Il modello finale presenta un valore di R² pari a 0.7358, indicando che il modello spiega circa 73.58% della variabilità del peso del neonato.
Il valore di R² aggiustato (Adjusted R²), leggermente inferiore, conferma la buona capacità esplicativa anche tenendo conto del numero di variabili incluse.

Analisi del modello VIF

# analisi delle varibili del modello_decimo  VIF
vif(modello_decimo)

##     Gestazione      Lunghezza I(Lunghezza^2)         Cranio          Sesso 
##       1.781859     237.696630     229.653852       1.607724       1.044426

L’analisi dei VIF evidenzia valori elevati per le variabili Lunghezza e I(Lunghezza²) (rispettivamente 237.7 e 229.7), dovuti alla forte correlazione strutturale tra una variabile e il suo termine quadratico. Tale risultato è atteso e non compromette la validità complessiva del modello, poiché il termine quadratico è stato introdotto per modellare un effetto non lineare.

RMSE (Root Mean Squared Error)

Serve a valutare l’errore medio nelle previsioni.

# Calcolo del RMSE
residui <- residuals(modello_decimo)
rmse <- sqrt(mean(residui^2))

Il valore del RMSE è pari a 269.8330401 grammi, indicando che in media il modello commette un errore di circa 270 grammi nella stima del peso del neonato. La media della variabile Peso del dataset è 3300 grammi il margine di errore è dell’8.18%.

Analisi dei Residui (parte erratica)

Per verificare che un modello sia valido i residui devono essere:

distribuiti casualmente
con media circa 0
e con varianza costante (omoschedasticità)

# Grafici diagnostici del modello
par(mfrow = c(2, 2))
plot(modello_decimo)

Grafico Residuals vs Fitted

Mostra la relazione tra i residui e i valori stimati. I punti sono distribuiti in modo abbastanza casuale intorno a 0. Ai bordi si notano alcuni residui più ampi segno di possibili outlier. -> Verificare i casi con residui estremi (osservazione 1551).

Grafico Normal Q–Q

Mostra se e i residui seguono una distribuzione normale. La maggior parte dei punti segue la linea ma alle estremità si discosta. Questo indica una leggera non normalità, probabilmente dovuta a pochi valori molto grandi (osservazione 1551).

Grafico Scale–Location

Mostra la varianza dei residui rispetto ai fitted values. La linea rossa è quasi orizzontale (varianza quasi costante e questo è un bene). Piccole fluttuazioni indicano che l’eteroschedasticità è minima. -> Il modello soddisfa bene l’assunzione di omoschedasticità.

Grafico Residuals vs Leverage

Mostra osservazioni che hanno una forte influenza sul modello. Alcuni punti hanno leve alte (sono influenti). -> Il modello sembra robusto ma si devono controllare le osservazioni evidenziate perchè potrebbero essere o errori di input o casi particolari.

Controllo dei punti influenti

# grafico per i punti influenti (con leva elevata)
lev <- hatvalues(modello_decimo)  
plot(lev)  
p = sum(lev)  
soglia = 2*p/count_record  
abline(h=soglia,col=2)

Il grafico mostra che ci sono diversi punti con leva elevata, corrispondenti a quelli che sono molto al di sopra della linea rossa.

La tabella seguente mostra le osservazioni con leva elevata.

# crea un data frame con indice e valore di leva superiore alla soglia
lev_influenti <- data.frame(
  Osservazione = which(lev > soglia),
  Leva = lev[lev > soglia]
)%>%
  arrange(desc(Leva))  # ordina per leva decrescente

# visualizza la tabella dei punti con leva elevata
kable(lev_influenti, 
      caption = "Osservazioni con leva elevata (ordinate per valore di leva)") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

Osservazioni con leva elevata (ordinate per valore di leva)
	Osservazione	Leva
1551	1551	0.1552797
928	928	0.1139467
2437	2437	0.0894685
1780	1780	0.0792499
1619	1619	0.0551926
2452	2452	0.0501384
2175	2175	0.0405639
2114	2114	0.0361004
312	312	0.0323187
805	805	0.0311605
310	310	0.0287638
2120	2120	0.0267620
1248	1248	0.0249878
2307	2307	0.0243824
1429	1429	0.0241893
2149	2149	0.0189061
1701	1701	0.0175106
492	492	0.0166497
1428	1428	0.0165403
378	378	0.0161125
1014	1014	0.0152629
2040	2040	0.0140453
106	106	0.0138319
2115	2115	0.0128657
748	748	0.0128186
101	101	0.0123723
1385	1385	0.0121959
947	947	0.0120065
2391	2391	0.0113541
2200	2200	0.0112349
151	151	0.0108768
2408	2408	0.0103592
206	206	0.0101433
1091	1091	0.0095473
1686	1686	0.0088416
1610	1610	0.0088046
684	684	0.0087830
587	587	0.0087748
2216	2216	0.0087271
750	750	0.0087214
471	471	0.0086362
1188	1188	0.0086143
956	956	0.0085842
2089	2089	0.0085552
1809	1809	0.0084515
2257	2257	0.0084285
155	155	0.0081963
1273	1273	0.0081902
1513	1513	0.0081603
2458	2458	0.0080766
666	666	0.0079155
1067	1067	0.0078960
445	445	0.0077056
220	220	0.0075408
131	131	0.0074133
964	964	0.0071923
15	15	0.0071186
36	36	0.0070253
1712	1712	0.0069914
1357	1357	0.0069570
328	328	0.0069154
2140	2140	0.0069091
729	729	0.0068476
1323	1323	0.0068476
1553	1553	0.0068429
1283	1283	0.0068359
991	991	0.0067759
2478	2478	0.0067657
1400	1400	0.0065884
1130	1130	0.0065881
383	383	0.0065556
592	592	0.0065415
765	765	0.0064537
1920	1920	0.0063999
638	638	0.0063143
34	34	0.0062821
988	988	0.0061446
1639	1639	0.0061446
1977	1977	0.0061344
1827	1827	0.0061159
61	61	0.0061025
317	317	0.0060716
1628	1628	0.0060590
895	895	0.0060351
958	958	0.0060304
2225	2225	0.0059380
305	305	0.0059187
1181	1181	0.0059013
1556	1556	0.0058960
67	67	0.0058886
1767	1767	0.0058703
697	697	0.0058115
2224	2224	0.0057907
1293	1293	0.0057388
1200	1200	0.0057004
1356	1356	0.0056796
1238	1238	0.0055875
119	119	0.0055316
1222	1222	0.0055316
1370	1370	0.0055213
315	315	0.0054625
190	190	0.0054579
1420	1420	0.0054265
1267	1267	0.0054057
600	600	0.0054012
271	271	0.0053619
205	205	0.0053469
304	304	0.0052879
241	241	0.0051929
1593	1593	0.0051702
1644	1644	0.0051431
1826	1826	0.0051431
702	702	0.0051324
1215	1215	0.0051206
165	165	0.0050984
428	428	0.0050984
1693	1693	0.0050427
1868	1868	0.0050303
486	486	0.0050167
1893	1893	0.0050035
2016	2016	0.0049244
513	513	0.0049069
96	96	0.0048693
2337	2337	0.0048560
656	656	0.0048538
565	565	0.0048537
1692	1692	0.0048504
2215	2215	0.0048452
1375	1375	0.0048272
629	629	0.0048174
1802	1802	0.0048136
1634	1634	0.0048032

Outliers

plot(rstudent(modello_decimo))
abline(h=c(-2,2),col=2)

Dal grafico si può osservare che la maggior parte dei residui è concentrata tra -2 e +2 e questo indica che il modello si adatta ragionevolmente ai dati. Ci sono diversi punti che superano le soglie (sopra alla soglia del +2 sovrastima; al di stotto del -2 sottostima) di poco ma c’è un punto y > 6 che è un outlier estremo. Il modello generalmente è corretto però bisogna indagare sull’outlier molto fuori dalla soglia.

Identifichiamo gli outlier più estremi tramite la seguente tabella:

car::outlierTest(modello_decimo)

##       rstudent unadjusted p-value Bonferroni p
## 1551  7.258271         5.2127e-13   1.3032e-09
## 1306  4.781064         1.8452e-06   4.6130e-03
## 155   4.664817         3.2521e-06   8.1303e-03
## 1399 -4.301603         1.7610e-05   4.4025e-02
## 1694  4.293808         1.8235e-05   4.5588e-02

Sono stati identificati 5 outlier statisticamente significativi:

Osservazione 1551 è la più estrema (è quella del grafico con y>6), Bonferroni p = 1.30e-09: altamente significativo.
Osservazioni 1306, 155, 1694: residui positivi (4.7, 4.6, 4.2), p-value Bonferroni < 0.05 (significativi).
Osservazione 1399 il residuo negativo -4.3, p-value Bonferroni = 0.044 (appena significativo)

Esaminiamo gli outlier

# filtro per il dataframe
outliers <- c(1551, 1306, 155, 1399, 1694)
# dataframe filtrato 
dati_corretto[outliers, ]  %>%
  kbl(caption = "outlier estremi",
      align = "c",
      booktabs = TRUE) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                font_size = 14)

outlier estremi
	Anni.madre	N.gravidanze	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Sesso	Ospedaleosp1	Ospedaleosp2	Ospedaleosp3
1551	35	1	38	4370	315	374	0	0	0	0	1
1306	23	0	41	4900	510	352	0	0	0	1	0
155	30	0	36	3610	410	330	0	1	1	0	0
1399	42	2	38	2560	525	349	1	1	0	1	0
1694	23	1	36	3850	460	334	1	0	0	0	1

Controllo distanza di Cook con il seguente grafico

cook <- cooks.distance(modello_decimo)

plot(cook)

Dal grafico della distanza di Cook emerge chiaramente che c’è un’osservazione estremamente influente che domina tutte le altre che va oltre 1.5.

index_max_cook <- max(cook)

il valore è 1.5812813 ed è un valore altissimo dato che la soglia critica è 1, questa osservazione sembra influire molto sulle capacità del modello_decicmo.

Residui studentizzati: abbiamo identificato 5 outlier, con l’osservazione 1551 che ha un valore estremo (rstudent = 7.26)
Distanza di Cook: max(cook) = 1.5812813 >> 1 (soglia critica), dal grafico si vede un punto che spicca nettamente sopra tutti gli altri

In conclusione: molto probabilmente si tratta della stessa osservazione (1551) che è sia un outlier (molto lontana dalla linea di regressione che un punto influente (modifica sostanzialmente le stime del modello)

Testiamo un modello togliendo l’osservazione 1551

# tramite il which.max(cook) togliamo l'osservazione con l'indice di Cook maggiore

modello_senza_influente <- update(modello_decimo, subset = -which.max(cook))

summary(modello_senza_influente)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + I(Lunghezza^2) + 
##     Cranio + Sesso, data = dati_corretto, subset = -which.max(cook))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1163.72  -182.53   -12.42   164.75  1307.91 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -1.663e+03  7.600e+02  -2.189 0.028702 *  
## Gestazione      3.640e+01  3.877e+00   9.387  < 2e-16 ***
## Lunghezza      -1.137e+01  3.347e+00  -3.398 0.000688 ***
## I(Lunghezza^2)  2.289e-02  3.434e-03   6.665 3.25e-11 ***
## Cranio          1.029e+01  4.187e-01  24.590  < 2e-16 ***
## Sesso           7.359e+01  1.094e+01   6.730 2.10e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 267.4 on 2493 degrees of freedom
## Multiple R-squared:  0.7408, Adjusted R-squared:  0.7403 
## F-statistic:  1425 on 5 and 2493 DF,  p-value: < 2.2e-16

R-squared = 0.7408. Valore ottimo spiega circa il 74.08% della variabilità del peso del neonato.

p-value < 0.0001: il modello è statisticamente significativo.

Utilizziamo la tecnica di Bayes (BIC) per la selezione di un modello.

BIC(modello_decimo,modello_senza_influente)

##                         df      BIC
## modello_decimo           7 35138.48
## modello_senza_influente  7 35073.18

Il modello con BIC più basso è il modello_senza_influente con un valore di 35073.18.

# Calcolo del RMSE
residui2 <- residuals(modello_senza_influente)
rmse2 <- sqrt(mean(residui2^2))

Il valore del RMSE è pari a 267.0797888 grammi, indicando che in media il modello commette un errore di circa 267 grammi nella stima del peso del neonato. La media della variabile Peso del dataset è 3300 grammi il margine di errore è dell’8.09%.

#calcolo del miglioramento,in percentuale, tra i due modelli
miglioramento <- round(((rmse - rmse2) / rmse) * 100,2)

Il miglioramento del nuovo modello è dello 1.02%.

Confronto modelli -> analisi dei residui

Grafici del modello: modello_decimo

# Grafici diagnostici del modello modello_decimo
par(mfrow = c(2, 2))
plot(modello_decimo)

Grafici del modello: modello_senza_influente

# Grafici diagnostici del modello modello_senza_influente
par(mfrow = c(2, 2))
plot(modello_senza_influente)

Residuals vs Fitted -> Renge residui migliorati da -1500 / +1500 a -1000 / +1000 e nessun outlier estremo (miglior modello -> modello_senza_influente).

Normal Q-Q Plot -> I punti seguono molto bene la linea teorica al centro e le code superiori si discostano leggermente, mentre per il modello_decimo le code superiori si discostano nettamente (miglior modello -> modello_senza_influente).

Scale-Location -> Meno punti estremi e distribuzione più omogenea (miglior modello -> modello_senza_influente).

Residuals vs Leverage -> Nel grafico del modello modello_senza_influente nessun punto supera le linee di Cook’s distance e tutti i punti hanno leverage < 0.04 e nessuna osservazione domina (miglior modello -> modello_senza_influente).

Dopo queste osservazioni si è deciso di continuare con il modello -> modello_senza_influente

3 Previsioni e Risultati

Variabili del modello

names(modello_senza_influente$model)

## [1] "Peso"           "Gestazione"     "Lunghezza"      "I(Lunghezza^2)"
## [5] "Cranio"         "Sesso"

La richiesta dell’esercizio è di stimare il peso di una neonata considerando una madre alla terza gravidanza e un parto alla 39ª settimana di gestazione.
Non è possibile effettuare questa stima perché la variabile N.gravidanze non è inclusa nel modello (il modello non tiene conto del numero di gravidanze materne). Sono state tuttavia eseguite delle stime con le variabili presenti nel modello:

Stima 1

# media della lunghezza dei neonati del dataset
Lunghezza.media <- mean(dati_corretto$Lunghezza, na.rm = TRUE)

# media della circonferenza del cranio dei neonati del dataset
Cranio.media <- mean(dati_corretto$Cranio, na.rm = TRUE)

previsione1 <- round(predict(modello_senza_influente
        , newdata = data.frame( Gestazione = 39
                               , Lunghezza = Lunghezza.media
                               , Cranio = Cranio.media
                               ,  Sesso = 0 #femmina
                               )),2)

previsione2 <- round(predict(modello_senza_influente
        , newdata = data.frame( Gestazione = 39
                               , Lunghezza = Lunghezza.media
                               , Cranio = Cranio.media
                               ,  Sesso = 1 #maschio
                               )),2)

Variabili uguali per tutte le stime

Lunghezza = Lunghezza.media
Cranio = Cranio.media
Gestazione = 39 settimane

Risultati:

neonato Femmina -> peso 3231.48 grammi.
neonato Maschio -> peso 3305.07 grammi.

Stima 2

# media della lunghezza dei neonati del dataset
Lunghezza.media <- mean(dati_corretto$Lunghezza, na.rm = TRUE)

# media della circonferenza del cranio dei neonati del dataset
Cranio.media <- mean(dati_corretto$Cranio, na.rm = TRUE)

previsione5 <- round(predict(modello_senza_influente
        , newdata = data.frame( Gestazione = 27
                               , Lunghezza = Lunghezza.media
                               , Cranio = Cranio.media
                               ,  Sesso = 0 #femmina
                               )),2)

previsione6 <- round(predict(modello_senza_influente
        , newdata = data.frame( Gestazione = 27
                               , Lunghezza = Lunghezza.media
                               , Cranio = Cranio.media
                               ,  Sesso = 1 #maschio
                               )),2)

Variabili uguali per tutte le stime

Lunghezza = Lunghezza.media
Cranio = Cranio.media
Gestazione = 27 settimane

Risultati:

neonato Femmina -> peso 2794.73 grammi.
neonato Maschio -> peso 2868.32 grammi.

Stima 3

# media della lunghezza dei neonati del dataset


previsione9 <- round(predict(modello_senza_influente
        , newdata = data.frame( Gestazione = 32
                               , Lunghezza = 465
                               , Cranio = Cranio.media
                               ,  Sesso = 0 #femmina
                               )),2)

previsione10 <- round(predict(modello_senza_influente
        , newdata = data.frame( Gestazione = 32
                               , Lunghezza = 465
                               , Cranio = Cranio.media
                               ,  Sesso = 1 #maschio
                               )),2)


previsione11 <- round(predict(modello_senza_influente
        , newdata = data.frame( Gestazione = 39
                               , Lunghezza = 465
                               , Cranio = Cranio.media
                               ,  Sesso = 1 #maschio
                               )),2)

previsione12 <- round(predict(modello_senza_influente
        , newdata = data.frame( Gestazione = 39
                               , Lunghezza = 465
                               , Cranio = Cranio.media
                               ,  Sesso = 0 #femmina
                               )),2)

Variabili uguali per tutte le stime

Cranio = Cranio.media
Lunghezza = 465

Risultati:

neonato Femmina con Gestazione = 32 settimane -> peso 2662.17 grammi.
neonato Femmina con Gestazione = 39 settimane -> peso 2916.94 grammi.
neonato Maschio con Gestazione = 32 settimane -> peso 2735.76 grammi.
neonato Maschio con Gestazione = 39 settimane -> peso 2990.54 grammi.

4 Visualizzazioni

Dal modello scelto per le previsioni (modello_senza_influente) è stata eliminata l’osservazione 1551, ma nel dataset di riferimento è ancora presente così è stato creato un dataframe copia dell’originale togliendo l’osservazione 1551.

#modello senza la riga 1551
dati_modello_grafici <- dati_corretto[-1551, ]

Relazione tra le variabili Gestazione, Peso e Sesso

ggplot(dati_modello_grafici, aes(x = Gestazione, y = Peso, 
                                  color = factor(Sesso, levels = c(0, 1), labels = c("F", "M")))) +
  geom_point(alpha = 0.6, size = 2) +
  scale_color_manual(
    values = c("F" = "#E74C3C", "M" = "#3498DB"),
    labels = c("F" = "Femmina", "M" = "Maschio")
  ) +
  labs(
    title = "Relazione tra Settimane di Gestazione e Peso alla Nascita",
    x = "Settimane di gestazione",
    y = "Peso (g)",
    color = "Sesso"
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

Il grafico mostra come all’aumentare delle settimane di gestazione aumenta anche il peso; la maggior parte dei dati è tra le 37 e 41 settimane.

Relazione tra Peso medio e Sesso del neonato

ggplot(dati_modello_grafici, aes(x = factor(Sesso, levels = c(0, 1), labels = c("F", "M")), 
                                  y = Peso, 
                                  fill = factor(Sesso, levels = c(0, 1), labels = c("F", "M")))) +
  geom_boxplot(alpha = 0.7) +
  scale_fill_manual(values = c("F" = "#E74C3C", "M" = "#3498DB")) +
  labs(
    title = "Peso medio per sesso",
    x = "Sesso",
    y = "Peso (g)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
        plot.title = element_text(hjust = 0.5, face = "bold")
  )

Dal grafico si evince che i Maschi (box Blue) tendono a pesare di più alla nascita delle Femmine (box rosa); entrambi hanno dei outliers sia in alto che in basso che indicano valori estremi.

Relazione tra le variabili Gestazione, Sesso, Fumatrici

ggplot(dati_modello_grafici, aes(x = Gestazione, y = Peso, 
                                  color = factor(Fumatrici, levels = c(0, 1), 
                                                labels = c("Non fumatrici", "Fumatrici")))) +
  geom_point(alpha = 0.5, size = 2) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.2) +  # Aggiungere intervallo di confidenza
  scale_color_manual(
    values = c("Non fumatrici" = "#2ECC71", "Fumatrici" = "#E74C3C"),
    name = "Stato"
  ) +
  labs(
    title = "Relazione tra Gestazione e Peso alla Nascita",
    subtitle = "Confronto tra madri fumatrici e non fumatrici",
    x = "Settimane di gestazione",
    y = "Peso (g)"
  ) +
  theme_minimal() +
  theme(
    legend.position = "top",
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, color = "gray40")
  )

Il grafico mostra la correlazione positiva tra Gestazione e Peso; la distribuzione è sbilanciata tra Fumatrici e Non Fumatrici (Punti Verdi più numerosi dei Punti Rossi), il numero di fumatrici è nettamente inferiore rispetto al numero delle non fumatrici. La linea verde è più alta della linea rossa e questo suggerisce che la variabile fumo è associato a valori di peso neonatale inferiori.

Confronto tra valori osservati e previsti (valutazione del modello)

#viene aggiunta una colonna nel dataset dati_modello_grafici con nome Peso_previsto e viene valorizzata con il peso previsto dal modello modello_senza_influente, predic calcolerà il valore previsto basandisi sulle variabili che costituiscono il modello Gestazione, Lunghezza, Cranio, Sesso
dati_modello_grafici$Peso_previsto <- predict(modello_senza_influente)

ggplot(dati_modello_grafici, aes(x = Peso, y = Peso_previsto)) +
  geom_point(color = "darkblue", alpha = 0.6) +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
  labs(
    title = "Confronto tra Peso osservato e previsto",
    x = "Peso osservato (g)",
    y = "Peso previsto (g)"
  ) +
  theme_minimal()

Il grafico presenta un diagramma di dispersione che confronta i valori di peso osservati (asse x) con i valori previsti dal modello di regressione lineare multipla (asse y). Il modello utilizzato include come variabili predittive: Gestazione, Lunghezza, Cranio e Sesso.

Gli elementi del grafico sono:

Punti Blu -> ogni punto indica una osservazione del dataset, mostrando la relazione tra Peso effettivo e Peso stimato dal modello.

Linea Rossa tratteggiata -> rappresenta la linea di previsione perfetta (x=y); se il modello prevedesse perfettamente ogni osservazione allora tutte le previsioni sarebbero sulla linea rossa tratteggiata.

Il modello mostra una buona capacità predittiva generale. La maggior parte delle osservazioni si trova in un area concentrata rispetto alla linea rossa tratteggiata. La maggior parte delle osservazioni è collocata nella fascia di peso compresa tra 2500 grammi e 4000 grammi.

Modello Statistico per la Previsione del Peso Neonatale

Fabio Finocchi

2025-10-16

Obiettivo del progetto:

Benefici:

1 Raccolta dei Dati e Struttura del Dataset

L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature.

2 Analisi e Modellizazione

Analisi Preliminare

Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.Inoltre si saggeranno le seguenti ipotesi con i test adatti:

in alcuni ospedali si fanno più parti cesarei.

La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione

Le misure antropometriche sono significativamente diverse tra i due sessi

Analisi variabile Peso (variabile dipendente)

Analisi delle correlazioni tra variabili

Creazione del Modello di Regressione

Selezione del modello ottimale

Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC)

Analisi della qualità del modello

3 Previsioni e Risultati

4 Visualizzazioni