Modello Statistico per la Previsione del Peso Neonatale (Progetto finale modulo statistica inferenziale)

### Caricamento delle librerie necessarie

# Lettura dei dati
data <- read.csv("neonati.csv")

Vengono mostrate le prime righe del dataset per visualizzare la tipologia di variabili.

nomi_nuovi <- c(
  "N.gravidanze" = "Numero.gravidanze",
  "Cranio" = "Circonferenza.cranica",
  "SessoM" = "Sesso.(Maschio=1)"
)


names(data) <- ifelse(
  names(data) %in% names(nomi_nuovi),
  nomi_nuovi[names(data)],
  names(data) 
)

kable(head(data), caption = "Prime 6 righe del dataset")

Prime 6 righe del dataset
Anni.madre	Numero.gravidanze	Gestazione	Peso	Lunghezza	Circonferenza.cranica	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

Mappiamo le variabili categoriche in interi per eventuali successive analisi nella fase di creazione del modello di regressione

data$Tipo.parto <- as.factor(data$Tipo.parto)
data$Ospedale <- as.factor(data$Ospedale)
data$Sesso <- as.factor(data$Sesso)

data$Tipo.parto_num <- as.numeric(data$Tipo.parto)
data$osp_num <- as.numeric(data$Ospedale)
data$sesso_num <- as.numeric(data$Sesso)

print(levels(data$Tipo.parto))

## [1] "Ces" "Nat"

print(levels(data$Ospedale))

## [1] "osp1" "osp2" "osp3"

print(levels(data$Sesso))

## [1] "F" "M"

kable(head(data), caption = "Prime 6 righe del dataset con le nuove variabili")

Prime 6 righe del dataset con le nuove variabili
Anni.madre	Numero.gravidanze	Gestazione	Peso	Lunghezza	Circonferenza.cranica	Tipo.parto	Ospedale	Sesso	Tipo.parto_num	osp_num	sesso_num
26	0	42	3380	490	325	Nat	osp3	M	2	3	2
21	2	39	3150	490	345	Nat	osp1	F	2	1	1
34	3	38	3640	500	375	Nat	osp2	M	2	2	2
28	1	41	3690	515	365	Nat	osp2	M	2	2	2
20	0	38	3700	480	335	Nat	osp3	F	2	3	1
32	0	40	3200	495	340	Nat	osp2	F	2	2	1

# Caricamento del pacchetto necessario
library(knitr)

# Supponiamo che 'data' sia il tuo data frame e 'Peso' la variabile di interesse
# Calcolo delle statistiche riassuntive
statistiche <- c(summary(data$Peso), "Deviazione standard" = sd(data$Peso),
  "Range (max-min)" = diff(range(data$Peso)))

# Conversione delle statistiche in un data frame
statistiche_peso <- data.frame(Statistiche = names(statistiche), Valore = as.numeric(statistiche))

# Visualizzazione della tabella con kable
kable(statistiche_peso, col.names = c("Statistiche", "Valore"), caption = "Statistiche Riassuntive del Peso")

Statistiche Riassuntive del Peso
Statistiche	Valore
Min.	830.0000
1st Qu.	2990.0000
Median	3300.0000
Mean	3284.0808
3rd Qu.	3620.0000
Max.	4930.0000
Deviazione standard	525.0387
Range (max-min)	4100.0000

Suddividiamo in classi la variabile peso per visualizzarne la distribuzione.

breaks <- seq(min(data$Peso),
              max(data$Peso),
              length.out=6)

data$Peso_class<-cut(x = data$Peso,breaks = breaks,dig.lab = 10)

Visualizziamo la nuova variabile Peso suddivisa in classi

ggplot(data = data) +
  geom_bar(aes(x = Peso_class)) +
  labs(x = "Peso(6 intervalli)", y = "Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Distribuzione del peso dei neonati su Grafico Boxplot

# Calcolo della mediana per la variabile peso( per evidenziarla nel grafico)


minimo <- statistiche["Min."]
primo_quartile <- statistiche["1st Qu."]
mediana <- statistiche["Median"]
media <- statistiche["Mean"]
terzo_quartile <- statistiche["3rd Qu."]
massimo <- statistiche["Max."]


library(ggplot2)
# Grafico boxplot
ggplot(data, aes(y = Peso)) +
  geom_boxplot(
    outlier.color = "red",      
    outlier.shape = 16,          
    outlier.size = 0.5            
  ) +
  
  scale_y_continuous(limits = c(0, 10000))+
  geom_hline(
    yintercept = mediana, 
    color = "red")+
  annotate("text", x = -0.5, y = mediana, label = paste( round(mediana, 2)),
           color = "red", vjust = -1, hjust = 0) +
  geom_hline(
    yintercept = minimo, 
    color = "red")+
  annotate("text", x = -0.5, y = minimo, label = paste( round(minimo, 2)),
           color = "red", vjust = -1, hjust = 0)+
  geom_hline(
    yintercept = massimo, 
    color = "red")+
  annotate("text", x = -0.5, y = massimo, label = paste( round(massimo, 2)),
           color = "red", vjust = -1, hjust = 0) +
  theme_classic() +
   scale_y_log10() +
  labs(title = "Distribuzione del peso dei neonati su Grafico Boxplot")

## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

Ho applicato una scala logaritmica in modo da evidenziare gli outliers per i pesi piu bassi poichè sono quelli che meritano maggiore attenzione.

Gli outliers hanno valori comunque in linea con le dimensioni del fenomeno oggetto di studio(830g, per esempio,è un peso plausibile per un neonato alla nascita sebbene estremamente basso.), pertanto non si vede la necessita di trattarli in modo particolare.

Visualizziamo un istogramma per verificare graficamente se il campione ha una distribuzione Normale.

breaks <- seq(min(data$Peso),
              max(data$Peso),
              length.out=120)

ggplot(data = data)+
  labs(x= "Peso(120 intervalli)", y = "Count" )+
  geom_histogram(aes(x=cut(x = data$Peso,breaks = breaks,dig.lab = 10)),stat="count")

Dal grafico a istogrammi si puo constatare che la distribuzione ha una distribuzione molto simile ad una normale.

Verifica analitica sulla normalita della distribuzione.

Eseguo un test per verificare se i dati del campione si distribuiscono come una normale.

shapiro.test(data$Peso)

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

Il test di Shapiro-Wilk ha un valore W molto vicino ad 1, pertanto non si rifiuta l’ipotesi nulla, ovvero il campione si distribuisce normalmente.

Si effettua un test t sulla variabile peso per saggiare l’ipotesi che la media del campione sia uguale alla media della popolazione, che in questo caso non è nota.

mean(data$Peso)

## [1] 3284.081

t.test(data$Peso)$conf.int

## [1] 3263.490 3304.672
## attr(,"conf.level")
## [1] 0.95

Con il 95% di confidenza, la media vera del peso nella popolazione è tra 3263.490 e 3304.672.La media del campione appartiene all’intervallo individuato dal t test, pertanto è signicativamente uguale alla media della popolazione.

cesareo_table <- table(data$Tipo.parto, data$Ospedale)

cesareo_df <- as.data.frame(cesareo_table)
names(cesareo_df) <- c("Tipo_parto", "Ospedale", "Conteggio")


ggplot(cesareo_df, aes(x = Ospedale, y = Conteggio, fill = Tipo_parto)) +
  geom_bar(stat = "identity", position = "dodge") +  
  labs(
    title = "Conteggio dei tipi di parto per ospedale",
    x = "Ospedale",
    y = "Numero di parti",
    fill = "Tipo parto"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

A livello grafico si può osservare che non ci sono significative differenze fra il tipo di parto per ospedale.

Si effettua il test Chi-quadro per verificare analiticamente una differenza significativa nei tipi di parto fra ospedali.

tabella_contingenza <- table(data$Tipo.parto, data$Ospedale)

print(tabella_contingenza)

##      
##       osp1 osp2 osp3
##   Ces  242  254  232
##   Nat  574  595  603

risultato_chi_quadro <- chisq.test(tabella_contingenza)

risultati <- data.frame(
  Statistica = c("Chi-quadro", "Gradi di libertà", "p-value"),
  Valore = c(
    round(risultato_chi_quadro$statistic, 3),
    risultato_chi_quadro$parameter,
    format.pval(risultato_chi_quadro$p.value, digits = 3, eps = 0.001)
  )
)

kable(risultati, col.names = c("Parametro", "Valore"), align = c('l', 'c'))

	Parametro	Valore
X-squared	Chi-quadro	1.097
df	Gradi di libertà	2
	p-value	0.578

Interpretazione del test

Il p-value è leggermente sopra la soglia per la quale rifiuteremmo l’ipotesi di non differenza significativa del numero dei parti cesarei per ospedale. Pertanto non si rifiuta l’ipotesi nulla.

Eseguiamo un t test a due campioni indipendenti per saggiare l’ipotesi che il peso medio sia significativamente diverso rispetto al sesso.

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

df_results <- data.frame(
  "Metrica" = c("Statistica t", "Gradi di libertà", "p-value", 
                "Confidenza (low)", "Confidenza (high)",
                "Media gruppo F", "Media gruppo M",
                "Ipotesi alternativa", "Metodo"),
  "Valore" = c(test$statistic, test$parameter, test$p.value,
               test$conf.int[1], test$conf.int[2],
               test$estimate[1], test$estimate[2],
               test$alternative, test$method)
)

kable(df_results, align = 'l')

Metrica	Valore
Statistica t	-12.1061453419536
Gradi di libertà	2490.71605821777
p-value	8.02974271144157e-33
Confidenza (low)	-287.105070589075
Confidenza (high)	-207.061466367937
Media gruppo F	3161.1321656051
Media gruppo M	3408.2154340836
Ipotesi alternativa	two.sided
Metodo	Welch Two Sample t-test

Il p-value osservato è molto inferiore alla soglia di 0.05,pertanto si rifiuta l’ipotesi nulla di uguaglianza della media del Peso fra i due gruppi(Maschio,Femmina)

Confronto della lunghezza per sesso

ggplot(data, aes(x = Sesso, y = data$`Lunghezza`, fill = Sesso)) +
  geom_boxplot() +
  labs(
    title = "Confronto della lunghezza per sesso",
    x = "Sesso",
    y = "Lunghezza"
  ) +
  theme_minimal()

t.test(data$`Lunghezza` ~ data$Sesso, data = data)

## 
##  Welch Two Sample t-test
## 
## data:  data$Lunghezza by data$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

È stato effettuato un t-test a due campioni per saggiare l’ipotesi di differenza in media della lunghezza fra i due gruppi di sesso.Il p-value di < 2.2e-16 (< 0.05) indica che ci sono differenze significative del peso in media tra i sessi. I maschi hanno una lunghezza media di 499.67 cm, mentre le femmine hanno una lunghezza media di 489.76 cm.

Confronto della dimensione del cranio per sesso.

ggplot(data, aes(x = Sesso, y = data$`Circonferenza.cranica`, fill = Sesso)) +
  geom_boxplot() +
  labs(
    title = "Confronto del diametro del cranio per Sesso",
    x = "Sesso",
    y = "Cranio (mm)"
  ) +
  theme_minimal()

data_numeric <- data[ ,c("Anni.madre","Numero.gravidanze","Fumatrici","Gestazione","Peso","Lunghezza","Circonferenza.cranica","sesso_num","Tipo.parto_num","osp_num")]
cor_matrix<-round(cor(data_numeric),2)
kable(cor_matrix, digits = 2, caption = "Matrice di Correlazione")

Matrice di Correlazione
	Anni.madre	Numero.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Circonferenza.cranica	sesso_num	Tipo.parto_num	osp_num
Anni.madre	1.00	0.38	0.01	-0.14	-0.02	-0.06	0.02	0.01	0.00	0.03
Numero.gravidanze	0.38	1.00	0.05	-0.10	0.00	-0.06	0.04	0.02	-0.02	0.01
Fumatrici	0.01	0.05	1.00	0.03	-0.02	-0.02	-0.01	0.01	0.02	-0.02
Gestazione	-0.14	-0.10	0.03	1.00	0.59	0.62	0.46	0.13	-0.01	0.02
Peso	-0.02	0.00	-0.02	0.59	1.00	0.80	0.70	0.24	0.00	0.03
Lunghezza	-0.06	-0.06	-0.02	0.62	0.80	1.00	0.60	0.19	-0.04	0.01
Circonferenza.cranica	0.02	0.04	-0.01	0.46	0.70	0.60	1.00	0.15	0.00	0.01
sesso_num	0.01	0.02	0.01	0.13	0.24	0.19	0.15	1.00	0.00	0.01
Tipo.parto_num	0.00	-0.02	0.02	-0.01	0.00	-0.04	0.00	0.00	1.00	0.02
osp_num	0.03	0.01	-0.02	0.02	0.03	0.01	0.01	0.01	0.02	1.00

Dalla tabella delle correlazioni si può vedere che le variabili ospedale e tipo parto non hanno una incidenza significativa con nessuna delle altre variabili, pertanto non verranno considerate a priori nella costruzione del modello di regressione.

Heatmap delle Correlazioni

cor_matrix <- cor(data %>% dplyr::select(Anni.madre, `Numero.gravidanze`, Fumatrici, Gestazione, Peso, `Lunghezza`, `Circonferenza.cranica`))
cor_melted <- melt(cor_matrix)

ggplot(cor_melted, aes(x = Var1, y = Var2, fill = value)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, 
                       limits = c(-1, 1), name = "Correlazione") +
  labs(title = "Heatmap delle correlazioni", x = "", y = "") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Per quanto riguarda la correlazione fra durata della gravidanza e se la madre sia una fumatrice, questa è trascurabile.

Invece c’è una correlazione negativa significativa, sebbene non elevata,fra durata della gravidanza ,eta della madre e anche con il numero di gravidanze.

model <- lm(Peso ~ Anni.madre + `Numero.gravidanze` + Fumatrici + 
            Gestazione + `Lunghezza` + `Circonferenza.cranica` + 
            Sesso, 
            data = data)
summary_model <- summary(model)

coeff_df <- as.data.frame(summary_model$coefficients)

kable(coeff_df, digits = 3, caption = "Risultati del Modello di Regressione")

Risultati del Modello di Regressione
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	-6714.411	141.151	-47.569	0.000
Anni.madre	0.958	1.135	0.845	0.398
Numero.gravidanze	11.276	4.669	2.415	0.016
Fumatrici	-30.296	27.597	-1.098	0.272
Gestazione	32.933	3.827	8.606	0.000
Lunghezza	10.234	0.301	34.009	0.000
Circonferenza.cranica	10.518	0.427	24.642	0.000
SessoM	78.085	11.204	6.969	0.000

Il modello iniziale include tutte le variabili predittive. Il summary mostra i coefficienti stimati, i p-value e l’R². Variabili con p-value < 0.05 sono considerate significative.

Calcoliamo la Variance Inflation Factor per verificare una situazione di potenziale multicollinearita fra le features.

vif_results <- vif(model)

kable(
  data.frame(Variable = names(vif_results), VIF = vif_results),
  caption = "Analisi multicollinearità (VIF)",
  col.names = c("Variabile", "VIF"),
  digits = 2
)

Analisi multicollinearità (VIF)
	Variabile	VIF
Anni.madre	Anni.madre	1.19
Numero.gravidanze	Numero.gravidanze	1.18
Fumatrici	Fumatrici	1.01
Gestazione	Gestazione	1.69
Lunghezza	Lunghezza	2.08
Circonferenza.cranica	Circonferenza.cranica	1.63
Sesso	Sesso	1.04

Osserviamo che i valori sono tutti relativamente bassi (< 10, valori oltre questa soglia indicano una multicollinearita elevata.)

Selezioniamo le variabili in modo ottimale con il metodo Step-wise utilizzando il criterio AIC (Akaike Information Criterion), che, a parità di “bontà”, adotta il modello con il minor numero di features.

stepwise_model <- stepAIC(model, direction = "both")

## Start:  AIC=28084.7
## Peso ~ Anni.madre + Numero.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Circonferenza.cranica + Sesso
## 
##                         Df Sum of Sq       RSS   AIC
## - Anni.madre             1     53803 187973654 28083
## - Fumatrici              1     90879 188010731 28084
## <none>                               187919851 28085
## - Numero.gravidanze      1    439797 188359648 28088
## - Sesso                  1   3662833 191582684 28131
## - Gestazione             1   5585184 193505035 28156
## - Circonferenza.cranica  1  45791026 233710877 28628
## - Lunghezza              1  87220887 275140738 29036
## 
## Step:  AIC=28083.42
## Peso ~ Numero.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Circonferenza.cranica + Sesso
## 
##                         Df Sum of Sq       RSS   AIC
## - Fumatrici              1     91892 188065546 28083
## <none>                               187973654 28083
## + Anni.madre             1     53803 187919851 28085
## - Numero.gravidanze      1    646039 188619694 28090
## - Sesso                  1   3671289 191644943 28130
## - Gestazione             1   5531705 193505359 28154
## - Circonferenza.cranica  1  46066755 234040410 28629
## - Lunghezza              1  87218857 275192512 29034
## 
## Step:  AIC=28082.64
## Peso ~ Numero.gravidanze + Gestazione + Lunghezza + Circonferenza.cranica + 
##     Sesso
## 
##                         Df Sum of Sq       RSS   AIC
## <none>                               188065546 28083
## + Fumatrici              1     91892 187973654 28083
## + Anni.madre             1     54816 188010731 28084
## - Numero.gravidanze      1    623141 188688687 28089
## - Sesso                  1   3655292 191720838 28129
## - Gestazione             1   5464853 193530399 28152
## - Circonferenza.cranica  1  46108583 234174130 28629
## - Lunghezza              1  87632762 275698308 29037

summary_model <- summary(stepwise_model)

coeff_df <- as.data.frame(summary_model$coefficients)

kable(coeff_df, digits = 3, caption = "Risultati del Modello di Regressione dopo la selezione delle variabili con metodo stepwise AIC")

Risultati del Modello di Regressione dopo la selezione delle variabili con metodo stepwise AIC
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	-6681.144	135.723	-49.226	0.000
Numero.gravidanze	12.475	4.340	2.875	0.004
Gestazione	32.332	3.798	8.513	0.000
Lunghezza	10.249	0.301	34.090	0.000
Circonferenza.cranica	10.540	0.426	24.728	0.000
SessoM	77.993	11.202	6.962	0.000

interpretazioni <- character()

coeff_df <- coeff_df[-1, , drop = FALSE]
coeff_df <- coeff_df[rownames(coeff_df) != "SessoM", , drop = FALSE]
coeff_df <- coeff_df[rownames(coeff_df) != "I(Lunghezza^2)", , drop = FALSE]


for (row in 1:nrow(coeff_df)) {
  var_name <- rownames(coeff_df)[row]
  estimate <- round(coeff_df[row, "Estimate"], 2)
  p_value <- coeff_df[row, "Pr(>|t|)"]
  

  if (p_value < 0.001) {
    significance <- "altamente significativo"
  } else if (p_value < 0.05) {
    significance <- "significativo"
  } else {
    significance <- "non significativo"
  }
  

  interpretazione <- paste0(
    "### Interpretazione per ", var_name, "\n",
    "Il coefficiente della Feature ", var_name, " è ", estimate, 
    "(p-value = ", format.pval(p_value), "), ",
    "indicando che per ogni incremento unitario in ", var_name, ", ",
    "il peso del neonato varia ",
    ifelse(estimate > 0, "in aumento", "in diminuzione")," di ", abs(estimate), " grammi. "
  )
  
  interpretazioni <- c(interpretazioni, interpretazione)
  
  
}

kable(
  data.frame(Interpretazioni = interpretazioni),
  col.names = "",  
  escape = FALSE,  
  caption = "Risultati dell'analisi di regressione"
)

Risultati dell’analisi di regressione
### Interpretazione per Numero.gravidanze
Il coefficiente della Feature Numero.gravidanze è 12.47(p-value = 0.0040789), indicando che per ogni incremento unitario in Numero.gravidanze, il peso del neonato varia in aumento di 12.47 grammi.
### Interpretazione per Gestazione
Il coefficiente della Feature Gestazione è 32.33(p-value = < 2.22e-16), indicando che per ogni incremento unitario in Gestazione, il peso del neonato varia in aumento di 32.33 grammi.
### Interpretazione per Lunghezza
Il coefficiente della Feature Lunghezza è 10.25(p-value = < 2.22e-16), indicando che per ogni incremento unitario in Lunghezza, il peso del neonato varia in aumento di 10.25 grammi.
### Interpretazione per Circonferenza.cranica
Il coefficiente della Feature Circonferenza.cranica è 10.54(p-value = < 2.22e-16), indicando che per ogni incremento unitario in Circonferenza.cranica, il peso del neonato varia in aumento di 10.54 grammi.

La selezione stepwise basata sull’AIC ha rimosso le variabili non significative. Il modello finale include solo le variabili più rilevanti per prevedere il peso dei neonati.

par(mfrow = c(2,2))
plot(stepwise_model)

Il grafico dei Residui vs Valori Predetti mostra che i valori si distribuiscono in modo casuale solo dopo i 2000 grammi circa. Per cercare di migliorare la diagnostica si verifichera successivamente la presenza di una relazione non lineare di alcune variabili.

# Selezione delle variabili numeriche
variabili_numeriche <- data[, c("Anni.madre", "Numero.gravidanze", "Gestazione", "Peso", "Lunghezza", "Circonferenza.cranica")]

# Creazione della matrice di scatter plot
pairs(variabili_numeriche, main = "Matrice di Scatter Plot delle Variabili Numeriche")

Introduciamo nel modello variabile Lunghezza al quadrato perchè nella matrice di scatterplot sembra avere una relazione quadreatica con il peso.

model2 <- update(stepwise_model,~. +I(`Lunghezza`^2))
summary_model <- summary(model2)

coeff_df <- as.data.frame(summary_model$coefficients)

kable(coeff_df, digits = 3, caption = "Risultati del Modello di Regressione con inclusa la relazione quadratica Lunghezza")

Risultati del Modello di Regressione con inclusa la relazione quadratica Lunghezza
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	215.091	723.560	0.297	0.766
Numero.gravidanze	14.098	4.264	3.306	0.001
Gestazione	42.505	3.874	10.972	0.000
Lunghezza	-20.273	3.161	-6.413	0.000
Circonferenza.cranica	10.650	0.419	25.439	0.000
SessoM	70.007	11.030	6.347	0.000
I(Lunghezza^2)	0.032	0.003	9.697	0.000

La relazione quadratica Lunghezza^2 ha un p-value < 0.05 pertanto è significativa.

par(mfrow = c(2,2))
plot(model2)

Il grafico dei Residui VS Valori Predetti ha una distribuzione piu randomica rispetto alla precedente diagnostica, pertanto il modello si puo considerare migliorato.

Validazione del modello con train-test split

set.seed(123)
train_index <- createDataPartition(data$Peso, p = 0.7, list = FALSE)
train_data <- data[train_index, ]
test_data <- data[-train_index, ]

Addestramento del modello sul training set

stepwise_model_train <- stepAIC(stepwise_model, direction = "both", data = train_data)

## Start:  AIC=28082.64
## Peso ~ Numero.gravidanze + Gestazione + Lunghezza + Circonferenza.cranica + 
##     Sesso
## 
##                         Df Sum of Sq       RSS   AIC
## <none>                               188065546 28083
## - Numero.gravidanze      1    623141 188688687 28089
## - Sesso                  1   3655292 191720838 28129
## - Gestazione             1   5464853 193530399 28152
## - Circonferenza.cranica  1  46108583 234174130 28629
## - Lunghezza              1  87632762 275698308 29037

summary_model <- summary(stepwise_model_train)

coeff_df <- as.data.frame(summary_model$coefficients)

kable(coeff_df, digits = 3, caption = "Risultati del Modello di Regressione con inclusa la relazione quadratica Lunghezza")

Risultati del Modello di Regressione con inclusa la relazione quadratica Lunghezza
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	-6681.144	135.723	-49.226	0.000
Numero.gravidanze	12.475	4.340	2.875	0.004
Gestazione	32.332	3.798	8.513	0.000
Lunghezza	10.249	0.301	34.090	0.000
Circonferenza.cranica	10.540	0.426	24.728	0.000
SessoM	77.993	11.202	6.962	0.000

Previsioni sul test set (Primo modello)

predictions_test <- predict(stepwise_model_train, newdata = test_data)

Metriche di performance

rmse_test <- sqrt(mean((test_data$Peso - predictions_test)^2))
r2_test <- cor(test_data$Peso, predictions_test)^2
cat("RMSE (Test):", rmse_test, "\nR-squared (Test):", r2_test, "\n")

## RMSE (Test): 278.1909 
## R-squared (Test): 0.730429

Addestramento del modello sul training set(Secondo Modello)

stepwise_model_train <- stepAIC(model2, direction = "both", data = train_data)

## Start:  AIC=27992.08
## Peso ~ Numero.gravidanze + Gestazione + Lunghezza + Circonferenza.cranica + 
##     Sesso + I(Lunghezza^2)
## 
##                         Df Sum of Sq       RSS   AIC
## <none>                               181230058 27992
## - Numero.gravidanze      1    794616 182024674 28001
## - Sesso                  1   2928693 184158750 28030
## - Lunghezza              1   2989429 184219486 28031
## - I(Lunghezza^2)         1   6835489 188065546 28083
## - Gestazione             1   8752028 189982085 28108
## - Circonferenza.cranica  1  47043457 228273515 28567

summary_model <- summary(stepwise_model_train)

coeff_df <- as.data.frame(summary_model$coefficients)

kable(coeff_df, digits = 3, caption = "Risultati del Modello di Regressione con inclusa la relazione quadratica Lunghezza")

Risultati del Modello di Regressione con inclusa la relazione quadratica Lunghezza
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	215.091	723.560	0.297	0.766
Numero.gravidanze	14.098	4.264	3.306	0.001
Gestazione	42.505	3.874	10.972	0.000
Lunghezza	-20.273	3.161	-6.413	0.000
Circonferenza.cranica	10.650	0.419	25.439	0.000
SessoM	70.007	11.030	6.347	0.000
I(Lunghezza^2)	0.032	0.003	9.697	0.000

Previsioni sul test set (Secondo modello)

predictions_test <- predict(stepwise_model_train, newdata = test_data)

Metriche di performance

rmse_test <- sqrt(mean((test_data$Peso - predictions_test)^2))
r2_test <- cor(test_data$Peso, predictions_test)^2
cat("RMSE (Test):", rmse_test, "\nR-squared (Test):", r2_test, "\n")

## RMSE (Test): 267.6717 
## R-squared (Test): 0.7503816

Il secondo modello, che considera la lunghezza al quadrato, ha un indice R quadro piu alto rispetto al primo, pertanto si puo dire che è piu affidabile nel effettuare stime predittive.

L’RMSE di 267.67 grammi indica che, in media, le previsioni del modello si discostano dai valori reali di circa 267.67 grammi. L’R² di 0.75 suggerisce che il modello spiega circa il 75% della variabilità del peso dei neonati.

Analisi comparativa tra ospedali

Anova

anova_result <- aov(Peso ~ Ospedale, data = data)

kable(
  summary(anova_result)[[1]], 
  caption = "Risultati ANOVA per il peso dei neonati per ospedale",
  digits = 3,
  align = 'c'
)

Risultati ANOVA per il peso dei neonati per ospedale
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Ospedale	2	936236.8	468118.4	1.699	0.183
Residuals	2497	687952304.9	275511.5	NA	NA

L’ANOVA confronta il peso medio dei neonati tra i tre ospedali. Il p-value di 0.183 (> 0.05) indica che non ci sono differenze significative.

Visualizzazioni Grafiche

Relazione tra durata della gravidanza e peso del neonato

  ggplot(data, aes(x = Gestazione, y = Peso)) + 
  geom_point(alpha = 0.5) + 
  geom_smooth(method = "lm", col = "blue") + 
  ggtitle("Relazione tra durata della gravidanza e peso del neonato")

## `geom_smooth()` using formula = 'y ~ x'

Effetto del fumo materno sul neonato

ggplot(data, aes(x = as.factor(Fumatrici), y = Peso)) + 
  geom_boxplot(fill = "lightblue") + 
  scale_x_discrete(labels = c("0" = "Non fumatrice", "1" = "Fumatrice")) + 
  ggtitle("Effetto del fumo materno sul peso del neonato") +
  xlab("Fumo materno") + 
  theme_minimal()

I neonati di madri fumatrici tendono ad avere un peso leggermente inferiore rispetto a quelli di madri non fumatrici.

# Creazione di un nuovo caso da predire
new_case <- data.frame(
  Eta_Madre = mean(data$Anni.madre),
  Numero.gravidanze = 3,  
  Fumatrici = 0,
  Gestazione = 39,  # 
  Lunghezza = mean(data$`Lunghezza`), 
  Circonferenza.cranica = mean(data$`Circonferenza.cranica`), 
  Tipo.parto = "Nat",  
  Ospedale = "osp1",  
  Sesso = "F"  
)

predicted_weight <- predict(stepwise_model, newdata = new_case)


risultato_text_output <- paste0(
  "Il peso previsto della neonata per una madre non fumatrice, alla terza gravidanza, ",
  "con un tempo di gestazione di 39 settimane, con parto naturale e con sesso Femmina è: ",
  round(predicted_weight, 2), " grammi"
)

# Visualizziamo in una tabella con kable
kable(risultato_text_output, col.names = "Risultato")

Risultato
Il peso previsto della neonata per una madre non fumatrice, alla terza gravidanza, con un tempo di gestazione di 39 settimane, con parto naturale e con sesso Femmina è: 3271.09 grammi

Modello Statistico per la Previsione del Peso Neonatale (Progetto finale modulo statistica inferenziale)

Feliziano Copponi

2025-03-04

Vengono mostrate le prime righe del dataset per visualizzare la tipologia di variabili.

Mappiamo le variabili categoriche in interi per eventuali successive analisi nella fase di creazione del modello di regressione

Suddividiamo in classi la variabile peso per visualizzarne la distribuzione.

Visualizziamo la nuova variabile Peso suddivisa in classi

Distribuzione del peso dei neonati su Grafico Boxplot

Ho applicato una scala logaritmica in modo da evidenziare gli outliers per i pesi piu bassi poichè sono quelli che meritano maggiore attenzione.

Gli outliers hanno valori comunque in linea con le dimensioni del fenomeno oggetto di studio(830g, per esempio,è un peso plausibile per un neonato alla nascita sebbene estremamente basso.), pertanto non si vede la necessita di trattarli in modo particolare.

Visualizziamo un istogramma per verificare graficamente se il campione ha una distribuzione Normale.

Dal grafico a istogrammi si puo constatare che la distribuzione ha una distribuzione molto simile ad una normale.

Verifica analitica sulla normalita della distribuzione.

Eseguo un test per verificare se i dati del campione si distribuiscono come una normale.

Il test di Shapiro-Wilk ha un valore W molto vicino ad 1, pertanto non si rifiuta l’ipotesi nulla, ovvero il campione si distribuisce normalmente.

Si effettua un test t sulla variabile peso per saggiare l’ipotesi che la media del campione sia uguale alla media della popolazione, che in questo caso non è nota.

Con il 95% di confidenza, la media vera del peso nella popolazione è tra 3263.490 e 3304.672.La media del campione appartiene all’intervallo individuato dal t test, pertanto è signicativamente uguale alla media della popolazione.

A livello grafico si può osservare che non ci sono significative differenze fra il tipo di parto per ospedale.

Si effettua il test Chi-quadro per verificare analiticamente una differenza significativa nei tipi di parto fra ospedali.

Interpretazione del test

Il p-value è leggermente sopra la soglia per la quale rifiuteremmo l’ipotesi di non differenza significativa del numero dei parti cesarei per ospedale. Pertanto non si rifiuta l’ipotesi nulla.

Eseguiamo un t test a due campioni indipendenti per saggiare l’ipotesi che il peso medio sia significativamente diverso rispetto al sesso.

Il p-value osservato è molto inferiore alla soglia di 0.05,pertanto si rifiuta l’ipotesi nulla di uguaglianza della media del Peso fra i due gruppi(Maschio,Femmina)

Confronto della lunghezza per sesso

Confronto della dimensione del cranio per sesso.

Dalla tabella delle correlazioni si può vedere che le variabili ospedale e tipo parto non hanno una incidenza significativa con nessuna delle altre variabili, pertanto non verranno considerate a priori nella costruzione del modello di regressione.

Heatmap delle Correlazioni

Per quanto riguarda la correlazione fra durata della gravidanza e se la madre sia una fumatrice, questa è trascurabile.

Invece c’è una correlazione negativa significativa, sebbene non elevata,fra durata della gravidanza ,eta della madre e anche con il numero di gravidanze.

Il modello iniziale include tutte le variabili predittive. Il summary mostra i coefficienti stimati, i p-value e l’R². Variabili con p-value < 0.05 sono considerate significative.

Calcoliamo la Variance Inflation Factor per verificare una situazione di potenziale multicollinearita fra le features.

Osserviamo che i valori sono tutti relativamente bassi (< 10, valori oltre questa soglia indicano una multicollinearita elevata.)

Selezioniamo le variabili in modo ottimale con il metodo Step-wise utilizzando il criterio AIC (Akaike Information Criterion), che, a parità di “bontà”, adotta il modello con il minor numero di features.

La selezione stepwise basata sull’AIC ha rimosso le variabili non significative. Il modello finale include solo le variabili più rilevanti per prevedere il peso dei neonati.

Il grafico dei Residui vs Valori Predetti mostra che i valori si distribuiscono in modo casuale solo dopo i 2000 grammi circa. Per cercare di migliorare la diagnostica si verifichera successivamente la presenza di una relazione non lineare di alcune variabili.

Introduciamo nel modello variabile Lunghezza al quadrato perchè nella matrice di scatterplot sembra avere una relazione quadreatica con il peso.

La relazione quadratica Lunghezza^2 ha un p-value < 0.05 pertanto è significativa.

Il grafico dei Residui VS Valori Predetti ha una distribuzione piu randomica rispetto alla precedente diagnostica, pertanto il modello si puo considerare migliorato.

Validazione del modello con train-test split

Addestramento del modello sul training set

Previsioni sul test set (Primo modello)

Metriche di performance

Addestramento del modello sul training set(Secondo Modello)

Previsioni sul test set (Secondo modello)

Metriche di performance

Il secondo modello, che considera la lunghezza al quadrato, ha un indice R quadro piu alto rispetto al primo, pertanto si puo dire che è piu affidabile nel effettuare stime predittive.

L’RMSE di 267.67 grammi indica che, in media, le previsioni del modello si discostano dai valori reali di circa 267.67 grammi. L’R² di 0.75 suggerisce che il modello spiega circa il 75% della variabilità del peso dei neonati.

Analisi comparativa tra ospedali

Anova

L’ANOVA confronta il peso medio dei neonati tra i tre ospedali. Il p-value di 0.183 (> 0.05) indica che non ci sono differenze significative.

Visualizzazioni Grafiche

Relazione tra durata della gravidanza e peso del neonato

Effetto del fumo materno sul neonato

I neonati di madri fumatrici tendono ad avere un peso leggermente inferiore rispetto a quelli di madri non fumatrici.