Inferenziale_Peso_Neonati
Luca Pontecorvo
2026-05-12
Importiamo le librerie utili per lo svolgimento del progetto.
library(moments)
library(tibble)
library(dplyr)
library(ggplot2)
library(ggpubr)
library(GGally)
library(car)
library(lmtest)
library(MASS)
library(scatterplot3d)- Introduzione
L’obiettivo di questo lavoro è costruire un modello statistico in grado di stimare il peso alla nascita dei neonati utilizzando le variabili cliniche e demografiche raccolte in tre ospedali nel dataset: neonati.csv Importiamo il dataset:
df <- read.csv("neonati.csv")
head(df,10)## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1 26 0 0 42 3380 490 325
## 2 21 2 0 39 3150 490 345
## 3 34 3 0 38 3640 500 375
## 4 28 1 0 41 3690 515 365
## 5 20 0 0 38 3700 480 335
## 6 32 0 0 40 3200 495 340
## 7 26 1 0 39 3100 480 345
## 8 25 0 0 40 3580 510 349
## 9 22 1 0 40 3670 500 335
## 10 23 0 0 41 3700 510 362
## Tipo.parto Ospedale Sesso
## 1 Nat osp3 M
## 2 Nat osp1 F
## 3 Nat osp2 M
## 4 Nat osp2 M
## 5 Nat osp3 F
## 6 Nat osp2 F
## 7 Nat osp3 F
## 8 Nat osp1 M
## 9 Ces osp2 F
## 10 Ces osp2 FIn particolare, si vuole capire quali fattori influenzano maggiormente il peso alla nascita, con attenzione a variabili come la durata della gestazione, il fumo materno e le caratteristiche fisiche del neonato.
- Descrizione del dataset Il dataset contiene le seguenti variabili: -“Anni.madre”: età della madre, variabile quantitativa continua -“N.gravidenze”: numero di gravidanze precedenti, variabile discreta -“Fumatrici”: variabile binaria 0 = non fumatrice, 1 = fumatrice -“Gestazione”: settimane di gravidanza, variabile continua quantitativa -“Peso”: peso alla nascita (g), variabile quantitativa continua -“Lunghezza”: lunghezza del neonato (mm), variabile quantitativa continua -“Cranio”: diametro cranico (mm), variabile quantitativa continua -“Tipo.parto”: tipo di parto se naturale o cesareo, variabile binaria -“Ospedale”: ospedale di nascita, variabile categorica -“Sesso”: M o F, variabile binaria
df$Fumatrici <- as.factor(df$Fumatrici)
df$Ospedale <- as.factor(df$Ospedale)
summary(df)## Anni.madre N.gravidanze Fumatrici Gestazione Peso
## Min. : 0.00 Min. : 0.0000 0:2396 Min. :25.00 Min. : 830
## 1st Qu.:25.00 1st Qu.: 0.0000 1: 104 1st Qu.:38.00 1st Qu.:2990
## Median :28.00 Median : 1.0000 Median :39.00 Median :3300
## Mean :28.16 Mean : 0.9812 Mean :38.98 Mean :3284
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:40.00 3rd Qu.:3620
## Max. :46.00 Max. :12.0000 Max. :43.00 Max. :4930
## Lunghezza Cranio Tipo.parto Ospedale Sesso
## Min. :310.0 Min. :235 Length:2500 osp1:816 Length:2500
## 1st Qu.:480.0 1st Qu.:330 Class :character osp2:849 Class :character
## Median :500.0 Median :340 Mode :character osp3:835 Mode :character
## Mean :494.7 Mean :340
## 3rd Qu.:510.0 3rd Qu.:350
## Max. :565.0 Max. :390- Analisi esplorativa:
Per prima cosa notiamo dal summary che il valore Minimo di Anni.madre è zero. Applichiamo un età di 12 anni per indagare outlier bilogicamente impossibili.
filter(df, Anni.madre<12)## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1 1 1 0 41 3250 490 350 Nat
## 2 0 0 0 39 3060 490 330 Nat
## Ospedale Sesso
## 1 osp2 F
## 2 osp3 MIl numero di osservazioni con Anni.madre biologicamente impossibili sono 2: età di 0 e 1 anno. Il numero di osservazioni totali è pari a 2500, potremmo pensare di effettuare una imputazione con una sostituzione della mediana di questi due casi biologicamente impossibili, anzichè eliminare totalmente le due osservazioni.
median_eta <- median(df$Anni.madre, na.rm = TRUE)
cat("Anni Madre mediana",median_eta," anni\n")## Anni Madre mediana 28 annidf$Anni.madre[df$Anni.madre < 12] <- median_eta
print("Imputazione avvenuta")## [1] "Imputazione avvenuta"Applichiamo il Test Shapiro-Wilk per determinare se la distribuzione di una variabile è normale o meno.
shapiro.test(df$Anni.madre)##
## Shapiro-Wilk normality test
##
## data: df$Anni.madre
## W = 0.99492, p-value = 1.471e-07Sebbene W è molto vicino ad 1, il che farebbe pensare a una forma quasi normale, il p-value è molto piccolo. Dunque, si rifiuta l’ipotesi nulla che i dati sono normali: La distribuzione degli Anni.madre è non normale.
hist(df$Anni.madre, col = "lightblue",
main = "Distribuzione degli anni madre",
xlab = "Anni")
Effettuiamo alcune ipotesi e saggiamo le seguenti ipotesi Hp: - Hp: In
tutti gli ospedali si fanno lo stesso numero di tipi di parti. Per
questa analisi utilizziamo il Chi quadro test.
hospital_birth_type <- table(df$Tipo.parto,df$Ospedale)
hospital_birth_type##
## osp1 osp2 osp3
## Ces 242 254 232
## Nat 574 595 603chisq.test(hospital_birth_type)##
## Pearson's Chi-squared test
##
## data: hospital_birth_type
## X-squared = 1.0972, df = 2, p-value = 0.5778ggballoonplot(data = as.data.frame(hospital_birth_type),
fill = "value")+
labs(x="Tipo di parto",
y="Sesso",
title = "Sesso vs Tipo di parto",
fill="Frequenza")+
guides(size=F)+
theme(plot.title = element_text(hjust = 0.5))
Il p-value è molto alto, quindi non si rifiuta l’ipotesi nulla. Pertanto
non c’è una preferenza di un ospedale a predilire un tipo di parto. Come
mostra anche il grafico della frequenza sopra, non si nota alcuna
significativa preferenza del parto cesario con l’ospedale.
- Hp: La media del peso e della lunghezza di questo campione sono significativamente uguali a quella della popolazione (al seguente link [https://www.cdc.gov/growthcharts/who_charts.htm]):
La popolazione: - Peso medio maschio: 3.520 g - Peso medio femmina: 3.361 g - Lunghezza media maschio: 51,5 cm - Lunghezza media femmina: 50,1 cm
Applichiamo il T-Student test.
t_weight_F <- t.test(filter(df, Sesso=="F")["Peso"],
mu = 3361,
conf.level = 0.95,
alternative = "two.sided")
t_weight_M <- t.test(filter(df, Sesso=="M")["Peso"],
mu = 3520,
conf.level = 0.95,
alternative = "two.sided")
t_length_F <- t.test(filter(df, Sesso=="F")["Lunghezza"],
mu = 501,
conf.level = 0.95,
alternative = "two.sided")
t_length_M <- t.test(filter(df, Sesso=="M")["Lunghezza"],
mu = 515,
conf.level = 0.95,
alternative = "two.sided")
cat("Pvalue per Peso medio femmina",t_weight_F$p.value,"\n")## Pvalue per Peso medio femmina 1.146491e-38cat("Pvalue per Peso medio maschio",t_weight_M$p.value,"\n")## Pvalue per Peso medio maschio 3.193454e-15cat("Pvalue per Lunghezza media femmina",t_length_F$p.value,"\n")## Pvalue per Lunghezza media femmina 5.864981e-44cat("Pvalue per Lunghezza media maschio",t_length_M$p.value,"\n")## Pvalue per Lunghezza media maschio 2.64661e-94Tutti i p-value relativi a peso e lunghezza (sia per i maschietti che per le femminucce) sono estremamente inferiori a 0.05. I neonati del campione in esame sono significativamente differenti dalla popolazione mondiale.
-Hp: Le misure antropometriche sono significativamente le stesse tra i due sessi.
Per questa analisi si saggia l’ipotesi nulla tramite il test di wilcox applicato alle dimensioni del neonato e al peso.
pairwise.wilcox.test(df$Peso, df$Sesso, p.adjust.method = "none")##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: df$Peso and df$Sesso
##
## F
## M <2e-16
##
## P value adjustment method: nonepairwise.wilcox.test(df$Lunghezza, df$Sesso, p.adjust.method = "none")##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: df$Lunghezza and df$Sesso
##
## F
## M <2e-16
##
## P value adjustment method: nonepairwise.wilcox.test(df$Cranio, df$Sesso, p.adjust.method = "none")##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: df$Cranio and df$Sesso
##
## F
## M 9.6e-15
##
## P value adjustment method: noneIl test di wilcox mostra p value bassissimi sia per la lunghezza, diametro del crano che il peso. Dunque vi è una distribuzione NON omogena con il sesso di queste variabili. Ovvero, le medie dei due gruppi sono significativamente diverse. L’ipotesi nulla è rifiutata e significa che il sesso influenza le dimensioni e il peso del neonato.
Infine, concludiamo l’esplorazione delle variabili con l’investigare l’indipendenza tra variabili. Necessario per definire il set di variabili che farà parte del modello.
- Analisi delle variabili quantitative (continue): calcoliamo la matrice di correlazione che sarà applicata alle sole variabili continue.
num.df <- df[, sapply(df, is.numeric)]
cor(num.df)## Anni.madre N.gravidanze Gestazione Peso Lunghezza
## Anni.madre 1.00000000 0.38328269 -0.1349262 -0.02377346 -0.06495563
## N.gravidanze 0.38328269 1.00000000 -0.1014919 0.00240730 -0.06040371
## Gestazione -0.13492624 -0.10149194 1.0000000 0.59176872 0.61892045
## Peso -0.02377346 0.00240730 0.5917687 1.00000000 0.79603676
## Lunghezza -0.06495563 -0.06040371 0.6189204 0.79603676 1.00000000
## Cranio 0.01620269 0.03900707 0.4608289 0.70480151 0.60334097
## Cranio
## Anni.madre 0.01620269
## N.gravidanze 0.03900707
## Gestazione 0.46082894
## Peso 0.70480151
## Lunghezza 0.60334097
## Cranio 1.00000000La matrice di correlazione mette in evidenza le seguenti correlazioni significative: - Anni della madre, il numero di gravidanze e gestazione. Maggiore è il numero degli anni e maggiore tendenzialmente sarà il numero di gravidanze affrontate. Maggiore è l’età della madre, minore sono le settimane di gestazione. - Il peso, la lunghezza e la dimensione del crano non sono significativamente correlate agli anni della madre, nonchè al numero di gravidanze, bensì al numero di settimane di gestazione. Maggiore sono le settimane di gestazione maggiore saranno le dimensioni.
- Analisi delle variabili categoriche. Creiamo una funzione personalizzata per calcolare il test del Chi-Quadro su tutte le possibili coppie di variabili categoriche. Definiamo una funzione per tale analisi ed eseguiamo.
chisq_calc <- function(df) {
results <- data.frame()
variables <- names(df)
for (i in 1:(length(variables)-1)) {
for (j in (i+1):length(variables)) {
var1 <- df[[variables[i]]]
var2 <- df[[variables[j]]]
contingency_table <- table(var1, var2)
chi_sq_test <- chisq.test(contingency_table)
p_value <- round(chi_sq_test$p.value, 3)
results <- rbind(results, c(variables[i], variables[j], p_value))}}
colnames(results) <- c("Variable1", "Variable2", "P_Value")
return(results)}Eseguiamo la funzione sopra definita. Ma prima trasformomiamo le categoriche in fattori.
df$Sesso <- as.factor(df$Sesso)
df$Ospedale <- as.factor(df$Ospedale)
df$Tipo.parto <- as.factor(df$Tipo.parto)
df$Fumatrici <- as.factor(df$Fumatrici)
categ.df <- df[, sapply(df, is.factor)]
chisq_calc(categ.df)## Variable1 Variable2 P_Value
## 1 Fumatrici Tipo.parto 0.404
## 2 Fumatrici Ospedale 0.698
## 3 Fumatrici Sesso 0.582
## 4 Tipo.parto Ospedale 0.578
## 5 Tipo.parto Sesso 0.843
## 6 Ospedale Sesso 0.737Il pvalue è maggiore di 0.05, quindi non si può rifiutare per tutte le combinazioni l’ipotesi nulla: variabili indipendenti tra loro. Quindi una madre fumatrice non influenza il tipo di parto e non vi sono maggiori concentrazioni di fumatrici in ospedale, cosi come non vi è una maggiore propensione di una donna fumatrice ad avere un sesso piuttosto che l’altro. E cosi via dicendo per tutte le altre combinazioni.
- Costruzione del modello di regressione lineare multipla. Il modello deve includere tutte le variabili rilevanti. Partiamo dal modello totale e togliamo una alla volta le variabili meno significative. Il target della regressione lineare è il peso.
Modello 1: tutte le variabili del dataset.
mod1 <- lm(Peso ~ ., data = df)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.3 -181.2 -14.6 160.7 2612.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6735.1677 141.3977 -47.633 < 2e-16 ***
## Anni.madre 0.7983 1.1463 0.696 0.4862
## N.gravidanze 11.4118 4.6665 2.445 0.0145 *
## Fumatrici1 -30.1567 27.5396 -1.095 0.2736
## Gestazione 32.5265 3.8179 8.520 < 2e-16 ***
## Lunghezza 10.2951 0.3007 34.237 < 2e-16 ***
## Cranio 10.4725 0.4261 24.580 < 2e-16 ***
## Tipo.partoNat 29.5027 12.0848 2.441 0.0147 *
## Ospedaleosp2 -11.2216 13.4388 -0.835 0.4038
## Ospedaleosp3 28.0984 13.4972 2.082 0.0375 *
## SessoM 77.5473 11.1779 6.938 5.07e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 669.1 on 10 and 2489 DF, p-value: < 2.2e-16mod2 <- update(mod1, ~ . - Anni.madre)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1113.93 -180.11 -16.36 160.58 2616.96
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6708.1065 135.9394 -49.346 < 2e-16 ***
## N.gravidanze 12.6085 4.3381 2.906 0.00369 **
## Fumatrici1 -30.3092 27.5359 -1.101 0.27113
## Gestazione 32.2501 3.7968 8.494 < 2e-16 ***
## Lunghezza 10.2944 0.3007 34.239 < 2e-16 ***
## Cranio 10.4876 0.4255 24.651 < 2e-16 ***
## Tipo.partoNat 29.5351 12.0834 2.444 0.01458 *
## Ospedaleosp2 -11.0816 13.4359 -0.825 0.40957
## Ospedaleosp3 28.3660 13.4903 2.103 0.03559 *
## SessoM 77.6205 11.1763 6.945 4.81e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2490 degrees of freedom
## Multiple R-squared: 0.7288, Adjusted R-squared: 0.7278
## F-statistic: 743.6 on 9 and 2490 DF, p-value: < 2.2e-16mod3 <- update(mod2, ~ . - Ospedale)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1130.04 -181.29 -16.31 160.59 2635.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6708.074 135.984 -49.330 < 2e-16 ***
## N.gravidanze 13.012 4.342 2.997 0.00276 **
## Fumatrici1 -31.759 27.570 -1.152 0.24946
## Gestazione 32.541 3.801 8.561 < 2e-16 ***
## Lunghezza 10.272 0.301 34.129 < 2e-16 ***
## Cranio 10.501 0.426 24.648 < 2e-16 ***
## Tipo.partoNat 30.296 12.098 2.504 0.01234 *
## SessoM 78.114 11.191 6.980 3.77e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2492 degrees of freedom
## Multiple R-squared: 0.7278, Adjusted R-squared: 0.7271
## F-statistic: 952 on 7 and 2492 DF, p-value: < 2.2e-16mod4 <- update(mod3, ~ . - Fumatori)
summary(mod4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1130.04 -181.29 -16.31 160.59 2635.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6708.074 135.984 -49.330 < 2e-16 ***
## N.gravidanze 13.012 4.342 2.997 0.00276 **
## Fumatrici1 -31.759 27.570 -1.152 0.24946
## Gestazione 32.541 3.801 8.561 < 2e-16 ***
## Lunghezza 10.272 0.301 34.129 < 2e-16 ***
## Cranio 10.501 0.426 24.648 < 2e-16 ***
## Tipo.partoNat 30.296 12.098 2.504 0.01234 *
## SessoM 78.114 11.191 6.980 3.77e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2492 degrees of freedom
## Multiple R-squared: 0.7278, Adjusted R-squared: 0.7271
## F-statistic: 952 on 7 and 2492 DF, p-value: < 2.2e-16Il modello 5 contiene le sole variabili rilevanti che sono emerse nella fase di esplorazione dati.
mod5 <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso + N.gravidanze, data = df)
summary(mod5)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso +
## N.gravidanze, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.44 -180.81 -15.58 163.64 2639.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.1445 135.7229 -49.226 < 2e-16 ***
## Gestazione 32.3321 3.7980 8.513 < 2e-16 ***
## Lunghezza 10.2486 0.3006 34.090 < 2e-16 ***
## Cranio 10.5402 0.4262 24.728 < 2e-16 ***
## SessoM 77.9927 11.2021 6.962 4.26e-12 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1328 on 5 and 2494 DF, p-value: < 2.2e-16Tutti i modelli hanno un R2 pressochè identico. Il modello 5 però di fatto gestisce le sole variabili indipendenti e questo lo si evince anche dai pvalue che tira fuori il summary, con valori tutti molto piccoli.
- Validazione dei Modelli (AIC e BIC): selezione del modello ottimale. Verifichiamo matematicamente se il modello 5 è davvero il migliore usando i criteri AIC e BIC. Più il valore è basso, migliore è il modello, restituendo un buon equilibrio tra precisione e semplicità.
BIC(mod1, mod2, mod3, mod4, mod5)## df BIC
## mod1 12 35241.97
## mod2 11 35234.64
## mod3 9 35228.24
## mod4 9 35228.24
## mod5 7 35220.10AIC(mod1, mod2, mod3, mod4, mod5)## df AIC
## mod1 12 35172.09
## mod2 11 35170.57
## mod3 9 35175.83
## mod4 9 35175.83
## mod5 7 35179.33Tendiamo a preferire il verdetto del BIC che restituisce il valore minimo al modello 5. Il modello 5 è il modello OTTIMALE.
- Analisi della qualità del modello. Ottenuto il modello finale, valutiamo la sua capacità predittiva utilizzando le metriche come R2 e RMSE.
par(mfrow = c(2, 2))
plot(mod5)
Nel primo grafico in alto a sinistra (Residuals vs Fitted), la linea
rossa di tendenza non è perfettamente piatta e orizzontale, ma mostra
una leggera curvatura. Questo è il tipico segnale della presenza di
effetti non-lineari tra le variabili. Introduciamo un nuovo modello che
tiene conto degli effetti non lineari.Introduciamo gli effetti
quadratici per le variabili Gestazioni, Cranio e Lunghezza.
mod6 <- update(mod5, ~ . + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2))
summary(mod6)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso +
## N.gravidanze + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2),
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1187.03 -182.50 -12.03 162.60 1469.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.217e+03 1.160e+03 -1.049 0.2945
## Gestazione 3.754e+02 6.738e+01 5.572 2.79e-08 ***
## Lunghezza -2.925e+01 4.436e+00 -6.592 5.28e-11 ***
## Cranio -4.979e+00 9.804e+00 -0.508 0.6116
## SessoM 7.232e+01 1.100e+01 6.577 5.83e-11 ***
## N.gravidanze 1.441e+01 4.246e+00 3.394 0.0007 ***
## I(Gestazione^2) -4.374e+00 8.834e-01 -4.951 7.86e-07 ***
## I(Lunghezza^2) 4.075e-02 4.544e-03 8.967 < 2e-16 ***
## I(Cranio^2) 2.276e-02 1.445e-02 1.575 0.1154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.4 on 2491 degrees of freedom
## Multiple R-squared: 0.7395, Adjusted R-squared: 0.7387
## F-statistic: 883.9 on 8 and 2491 DF, p-value: < 2.2e-16par(mfrow = c(2, 2))
plot(mod6)
La linea rossa nei residui ora è decisamente più piatta.Tuttavia, dal
summary si evince che il pvalue delle variabili Cranio e Cranio^2 nel
modello 6 assumono valori maggiori di 0.05, perdendo significatività
statistica. Il modello è tornato ad essere complesso. Creiamo un
ulteriore modello senza l’effetto quadratico di Cranio.
mod7 <- update(mod6, ~ . - I(Cranio^2))
summary(mod7)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso +
## N.gravidanze + I(Gestazione^2) + I(Lunghezza^2), data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1191.07 -182.28 -13.74 163.34 1403.53
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.360e+03 9.053e+02 -2.607 0.00919 **
## Gestazione 3.365e+02 6.270e+01 5.367 8.74e-08 ***
## Lunghezza -3.215e+01 4.037e+00 -7.963 2.53e-15 ***
## Cranio 1.045e+01 4.192e-01 24.922 < 2e-16 ***
## SessoM 7.261e+01 1.100e+01 6.602 4.93e-11 ***
## N.gravidanze 1.448e+01 4.247e+00 3.410 0.00066 ***
## I(Gestazione^2) -3.873e+00 8.245e-01 -4.698 2.77e-06 ***
## I(Lunghezza^2) 4.370e-02 4.140e-03 10.556 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.5 on 2492 degrees of freedom
## Multiple R-squared: 0.7392, Adjusted R-squared: 0.7385
## F-statistic: 1009 on 7 and 2492 DF, p-value: < 2.2e-16par(mfrow = c(2, 2))
plot(mod7)
Per il modello 7, torniamo ad avere tutti i p value bassi < di 0.05,
e il Residuals sono piatti con i Fitted Values. QUest’ultimo modello
sembra il più promettente perchè significativo (verifica del P value) e
considera anche gli effetti quadratici.
- Verifica delle assunzioni matematiche:
- Controllo della multicollinearità (II Test VIF). Ulteriore verifica che le variabili sono indipendenti e che non si stiano influenzando a vicenda
vif(mod7, type = "predictor")## GVIF Df GVIF^(1/(2*Df)) Interacts With
## Gestazione 3.100126 2 1.326920 I(Gestazione^2)
## Lunghezza 3.746905 2 1.391292 I(Lunghezza^2)
## Cranio 1.643288 1 1.281908 --
## Sesso 1.048541 1 1.023983 --
## N.gravidanze 1.025434 1 1.012637 --
## Other Predictors
## Gestazione Lunghezza, Cranio, Sesso, N.gravidanze
## Lunghezza Gestazione, Cranio, Sesso, N.gravidanze
## Cranio Gestazione, Lunghezza, Sesso, N.gravidanze
## Sesso Gestazione, Lunghezza, Cranio, N.gravidanze
## N.gravidanze Gestazione, Lunghezza, Cranio, SessoNon c’è, come atteso, multicollinearità, i valori sono tutti inferiori a 5.
- Controllo della media dei residui. La media degli errori deve essere pari a zero per garantire che il modello non stia sovrastimando o sottostimando sistematicamente.
mean(residuals(mod7))## [1] 5.251133e-15sd(residuals(mod7))## [1] 268.1125Il calcolo restituisce un numero vicinissimo allo zero reale. Assunzione pienamente rispettata.
- Controllo dell’Omoschedasticità (Testi di Breusch-Pagan). Verifichiamo se la varianza degli errori è costante o se cambia a seconda del peso del bambino.
library(lmtest)
bptest(mod7)##
## studentized Breusch-Pagan test
##
## data: mod7
## BP = 97.553, df = 7, p-value < 2.2e-16Il test restituisce un p value inferiore a 0.05. Questo ci costringe a rifiutare l’ipotesi nulla: i nostri residui sono purtroppo eteroschedastici (la loro varianza non è costante). Lo risolviamo con l’introduzione di nuovo modello
- Controllo dell’indipindenza degli errori (Test di Durbin-Watson). Verifichiamo se gli errori commessi su un neonato influenzano l’errore commesso sul neonato successivo.
dwtest(mod7)##
## Durbin-Watson test
##
## data: mod7
## DW = 1.9494, p-value = 0.1028
## alternative hypothesis: true autocorrelation is greater than 0Il pvalue è maggiore di 0.05, non si rifiuta l’ipotesi nulla, dunque gli errori sono indipendenti e correlati.
- Pulizia dei dati e identificazione dell’Outlier (Cook) Applichiamo la distanza di Cook per scovare gli outliers che se presenti spostano anomalamente i coefficienti del modello.
cook <- cooks.distance(mod7)
ggplot() +
geom_point(aes(x = 1:length(cook), y = cook, colour = cook > 1), size = 3) +
geom_hline(aes(yintercept = c(0.5, 1)), linetype = 2, colour = "darkred") +
scale_colour_manual(values = setNames(c("darkred", "black"), c(TRUE, FALSE))) +
labs(title = "Analisi dei residui: Distanza di Cook", x = "Indice", y = "Distanza di Cook") +
theme_minimal() +
theme(plot.title = element_text(size = 22, hjust = 0.5), legend.position = "none")
Dal grafico emerge chiaramente che c’è una osservazione, intorno al
numero 1150 che ha una distanza di Cook superiore a 1 (la soglia di
pericolo). È un valore influente estremo. Eliminiamo tale osservazione
dal dataset e vediamo come risponde il modello senza tale
distorsione.
index <- match(max(cook), cook)
cat("Osservazione da rimuovere è la numero:",index,"\n")## Osservazione da rimuovere è la numero: 1551corrected_nb.df <- df[-index, ]
print("Modello riaddestrato con il nuovo dataframe")## [1] "Modello riaddestrato con il nuovo dataframe"mod8 <- lm(Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
I(Gestazione^2) + I(Lunghezza^2), data = corrected_nb.df)
summary(mod8)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2) + I(Lunghezza^2), data = corrected_nb.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1187.72 -180.63 -13.61 164.12 1320.11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.811e+03 9.019e+02 -3.117 0.001848 **
## N.gravidanze 1.440e+01 4.217e+00 3.414 0.000651 ***
## Gestazione 2.004e+02 6.617e+01 3.028 0.002489 **
## Lunghezza -1.928e+01 4.535e+00 -4.252 2.20e-05 ***
## Cranio 1.009e+01 4.202e-01 24.023 < 2e-16 ***
## SessoM 7.340e+01 1.092e+01 6.721 2.22e-11 ***
## I(Gestazione^2) -2.135e+00 8.673e-01 -2.462 0.013890 *
## I(Lunghezza^2) 3.093e-02 4.618e-03 6.697 2.62e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.6 on 2491 degrees of freedom
## Multiple R-squared: 0.7426, Adjusted R-squared: 0.7419
## F-statistic: 1027 on 7 and 2491 DF, p-value: < 2.2e-16L’osservazione Outlier è precisamente la numero 1551. Tale valore è stato rimosso dal dataset. L’R2 è aumentato leggermente. Quindi bene, tuttavia l’effetto quadratico di Gestazione ha un p value maggiore di 0.05. Rimuoviamolo e vediamo come procede.
mod9 <- update(mod8, ~ . - I(Gestazione^2))
summary(mod9)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2), data = corrected_nb.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1176.79 -178.89 -12.52 164.19 1327.44
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.607e+03 7.586e+02 -2.119 0.034207 *
## N.gravidanze 1.420e+01 4.220e+00 3.364 0.000781 ***
## Gestazione 3.773e+01 3.890e+00 9.700 < 2e-16 ***
## Lunghezza -1.172e+01 3.342e+00 -3.508 0.000459 ***
## Cranio 1.015e+01 4.201e-01 24.157 < 2e-16 ***
## SessoM 7.222e+01 1.092e+01 6.614 4.57e-11 ***
## I(Lunghezza^2) 2.330e-02 3.430e-03 6.795 1.35e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.8 on 2492 degrees of freedom
## Multiple R-squared: 0.742, Adjusted R-squared: 0.7413
## F-statistic: 1194 on 6 and 2492 DF, p-value: < 2.2e-16Il modello 9 ha un R2 alto ed è pulito e significativo per via di tutti i pvalue molto bassi, minori di 0.05.
- Applichiamo il modello. Effettuiamo previsioni del peso di un neonato. Supponiamo una osservazione:
- N.gravidanze=3
- Gestazione=39
- Sesso= F
- Lunghezza e Cranio pari al valore medio della popolazione femminile.
lunghezza_media_F <- mean(df$Lunghezza[df$Sesso == "F"], na.rm = TRUE)
cranio_medio_F <- mean(df$Cranio[df$Sesso == "F"], na.rm = TRUE)
osservazione_pred <- data.frame(
N.gravidanze = 3,
Gestazione = 39,
Lunghezza = lunghezza_media_F,
Cranio = cranio_medio_F,
Sesso = "F")
prediction <- predict(mod9,
newdata = osservazione_pred,
interval = "predict")
prediction## fit lwr upr
## 1 3180.961 2657.185 3704.736Il bambino di sesso F con una gestazione di 39 settimane e con N gravidanze pari a 3, nascerà con un peso di 3,180 Kg.
- Per concludere, si mostra visivamente l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.
corrected_nb.df$Fumatrici <- factor(corrected_nb.df$Fumatrici,
levels = c(0, 1),
labels = c("Non Fumatrice", "Fumatrice"))
ggplot(corrected_nb.df, aes(x = Gestazione, y = Peso, color = Fumatrici)) +
geom_jitter(alpha = 0.3, size = 1.5, width = 0.2) +
geom_smooth(method = "lm", se = TRUE, size = 1.2) +
scale_color_manual(values = c("darkgreen", "red3")) +
labs(
title = "Impatto della Gestazione e del Fumo sul Peso del Neonato",
subtitle = "Analisi degli effetti principali del modello di regressione",
x = "Settimane di Gestazione",
y = "Peso del Neonato (grammi)",
color = "Stato della Madre"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 11, italic = TRUE, hjust = 0.5),
legend.position = "bottom")
Il grafico evidenzia chiaramente la relazione positiva tra la durata
della gestazione e il peso alla nascita: all’aumentare delle settimane,
il peso atteso cresce in modo lineare e costante. Tuttavia, a parità di
settimane di gestazione, la retta associata alle madri fumatrici (in
rosso) si posiziona sistematicamente al di sotto di quella delle non
fumatrici (in verde). Questo gap grafico quantifica visivamente
l’effetto penalizzante del fumo sul peso del neonato, confermando
l’effetto negativo stimato dai coefficienti del modello lineare.