Modello Statistico per la Previsione del Peso Neonatale
DM
2025-05-17
Modello Statistico per la Previsione del Peso Neonatale
1. Struttura del dataset
dati <- read.csv("neonati.csv",
stringsAsFactors = T,
sep = ",")
head(dati)## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1 26 0 0 42 3380 490 325 Nat
## 2 21 2 0 39 3150 490 345 Nat
## 3 34 3 0 38 3640 500 375 Nat
## 4 28 1 0 41 3690 515 365 Nat
## 5 20 0 0 38 3700 480 335 Nat
## 6 32 0 0 40 3200 495 340 Nat
## Ospedale Sesso
## 1 osp3 M
## 2 osp1 F
## 3 osp2 M
## 4 osp2 M
## 5 osp3 F
## 6 osp2 F# converte in variabili numeriche
dati$Peso <- as.numeric(dati$Peso)
dati$Lunghezza <- as.numeric(dati$Lunghezza)
dati$Cranio <- as.numeric(dati$Cranio)attach(dati)
N = nrow(dati)Nel dataset sono presenti 2500 osservazioni e 10 variabili.
2. Analisi e Modellazione
# verifica presenza di valori nulli
colSums(is.na(dati))## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## 0 0 0 0 0 0
## Cranio Tipo.parto Ospedale Sesso
## 0 0 0 0# libreria per indici di forma
library(moments)
# libreria per la funzione kable()
library(knitr)# funzione per Coefficiente di Variazione
CV <-function(x){
return( sd(x)/mean(x) * 100 )
}
# funzione per gli indici di variabilità e forma
stat_index <- function(x) {
as.data.frame(
list(
Range = max(x) - min(x),
IQR = IQR(x),
Var = var(x),
SD = sd(x),
CV = CV(x),
skewness = skewness(x),
kurtosis = kurtosis(x)-3
)
)
}Analisi preliminare
Anni.madre
summary(Anni.madre)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 25.00 28.00 28.16 32.00 46.00stat_index(Anni.madre)## Range IQR Var SD CV skewness kurtosis
## 1 46 7 27.81063 5.273578 18.72454 0.0428115 0.3804165Variabile quasi simmetrica e leggermente leptocurtica.
La variabile Anni.madre ha un minimo di 0, verifichiamo se ci sono altre
anomalie con un boxplot.
boxplot(Anni.madre)
dati[dati$Anni.madre < 10, ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1152 1 1 0 41 3250 490 350
## 1380 0 0 0 39 3060 490 330
## Tipo.parto Ospedale Sesso
## 1152 Nat osp2 F
## 1380 Nat osp3 MEssendo solo 2 osservazioni anomale su 2500 (meno dello 0,1%) optiamo per la rimozione.
df_bkp <- dati
dati <- dati[dati$Anni.madre > 10, ]
N = nrow(dati)N.gravidanze
summary(N.gravidanze)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.9812 1.0000 12.0000freq_ngrav = table(N.gravidanze)
kable(freq_ngrav, caption = "Frequenze N.gravidanze")| N.gravidanze | Freq |
|---|---|
| 0 | 1096 |
| 1 | 818 |
| 2 | 340 |
| 3 | 150 |
| 4 | 48 |
| 5 | 21 |
| 6 | 11 |
| 7 | 1 |
| 8 | 8 |
| 9 | 2 |
| 10 | 3 |
| 11 | 1 |
| 12 | 1 |
La maggior parte delle donne ha avuto 0 o 1 grafidanza; sono presenti
anche cinque osservazioni con 10 o più gravidanze.
Non si rilevano anomalie.
Fumatrici
freq_fuma = table(Fumatrici)
kable(freq_fuma, caption = "Frequenze Fumatrici(1)")| Fumatrici | Freq |
|---|---|
| 0 | 2396 |
| 1 | 104 |
La maggior parte delle osservazioni sono di non fumatrici.
Non si rilevano anomalie
Gestazione
freq_gest = table(Gestazione)
kable(freq_gest, caption = "Frequenze settimane gestazione")| Gestazione | Freq |
|---|---|
| 25 | 1 |
| 26 | 1 |
| 27 | 2 |
| 28 | 4 |
| 29 | 3 |
| 30 | 5 |
| 31 | 8 |
| 32 | 9 |
| 33 | 18 |
| 34 | 16 |
| 35 | 33 |
| 36 | 62 |
| 37 | 192 |
| 38 | 437 |
| 39 | 581 |
| 40 | 741 |
| 41 | 329 |
| 42 | 56 |
| 43 | 2 |
La maggior parte delle osservazioni ha 40 settimane di gestazione che è il valore di riferimento per la gestazione umana.
Non si rilevano anomalie.
Lunghezza
summary(Lunghezza)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 310.0 480.0 500.0 494.7 510.0 565.0stat_index(Lunghezza)## Range IQR Var SD CV skewness kurtosis
## 1 255 30 692.671 26.31864 5.320208 -1.514699 6.487174Variabile con asimmetria negativa e leptocurtica.
La media del campione è quasi 50 cm che è il valore di riferimento per i
neonati.
boxplot(Lunghezza)
Ci sono potenziali outlier, ma non si rilevano anomalie.
Cranio
summary(Cranio)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 235 330 340 340 350 390stat_index(Cranio)## Range IQR Var SD CV skewness kurtosis
## 1 155 20 269.7915 16.42533 4.830565 -0.7850527 2.946206Variabile con asimmetria negativa e leptocurtica.
boxplot(Cranio)
Ci sono potenziali outlier, ma non si rilevano anomalie.
Tipo.parto
freq_tparto = table(Tipo.parto)
kable(freq_tparto, caption = "Frequenze Tipo parto")| Tipo.parto | Freq |
|---|---|
| Ces | 728 |
| Nat | 1772 |
La maggior parte dei parti nel campione è naturale.
Non si rilevano anomalie.
Ospedale
freq_ospe = table(Ospedale)
kable(freq_ospe, caption = "Frequenze Ospedale")| Ospedale | Freq |
|---|---|
| osp1 | 816 |
| osp2 | 849 |
| osp3 | 835 |
Le osservazioni sono quali equamente distribuite tra i tre ospedali.
Sesso del neonato
freq_sex = table(Sesso)
kable(freq_sex, caption = "Frequenze Sesso Neonato")| Sesso | Freq |
|---|---|
| F | 1256 |
| M | 1244 |
Le osservazioni sono equamente distribuite tra i due sessi.
Peso (variabile target)
summary(Peso)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 2990 3300 3284 3620 4930stat_index(Peso)## Range IQR Var SD CV skewness kurtosis
## 1 4100 630 275665.7 525.0387 15.98739 -0.6470308 2.031532boxplot(Peso)
Ci sono potenziali outlier, ma non si rilevano anomalie.
Risulta una leggera asimmetria negativa e una cutosi positiva
(distribuzione leptocurtica).
Verifichiamo la Normalità con il test di Shapiro
shapiro.test(Peso)##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16Il test di Shapiro-Wilk ci porta a rifiutare l’ipotesi nulla: la
variabile Peso non ha una distribuzione Normale.
Non essendo comunque molto distorta non verrano applicate
trasformazioni.
Test di ipotesi
- in alcuni ospedali si fanno più parti cesarei
Per verificare questa ipotesi verifichiamo se le variabili Ospedale e Tipo.parto sono dipendenti con un test χ².
freq_ospe_parto <- table(Ospedale,Tipo.parto)
freq_ospe_parto## Tipo.parto
## Ospedale Ces Nat
## osp1 242 574
## osp2 254 595
## osp3 232 603test.indipendenza<- chisq.test(freq_ospe_parto)
test.indipendenza##
## Pearson's Chi-squared test
##
## data: freq_ospe_parto
## X-squared = 1.0972, df = 2, p-value = 0.5778Non si rifiuta l’ipotesti nulla di indipendenza, quindi le due varibili sono indipendenti, quindi statisticamente non si può affermare che 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
Per verificare questa ipotesti si effettua un T test con i seguenti valori medi della popolazione italiana (dati Ministero della Salute 2022):
media_peso_neonati_popolazione <- 3300
media_lunghezza_neonati_popolazione <- 500t.test(x = dati$Peso,
mu = media_peso_neonati_popolazione,
conf.level = 0.95,
alternative = "two.sided")##
## One Sample t-test
##
## data: dati$Peso
## t = -1.505, df = 2497, p-value = 0.1324
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
## 3263.577 3304.791
## sample estimates:
## mean of x
## 3284.184Il peso risulta significativamente uguale a quello della popolazione.
t.test(x = dati$Lunghezza,
mu = media_lunghezza_neonati_popolazione,
conf.level = 0.95,
alternative = "two.sided")##
## One Sample t-test
##
## data: dati$Lunghezza
## t = -10.069, df = 2497, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 493.6628 495.7287
## sample estimates:
## mean of x
## 494.6958La lunghezza invece NON risulta significativamente uguale alla media.
- Le misure antropometriche sono significativamente diverse tra i due sessi
t.test(Peso ~ Sesso, data = dati)##
## Welch Two Sample t-test
##
## data: Peso by Sesso
## t = -12.115, df = 2488.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.4841 -207.3844
## sample estimates:
## mean in group F mean in group M
## 3161.061 3408.496t.test(Lunghezza ~ Sesso, data = dati)##
## Welch Two Sample t-test
##
## data: Lunghezza by Sesso
## t = -9.5823, df = 2457.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -11.939001 -7.882672
## sample estimates:
## mean in group F mean in group M
## 489.7641 499.6750t.test(Cranio ~ Sesso, data = dati)##
## Welch Two Sample t-test
##
## data: Cranio by Sesso
## t = -7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -6.110504 -3.560417
## sample estimates:
## mean in group F mean in group M
## 337.6231 342.4586Tutte le misure antropometriche sono significativamente diverse tra i due sessi.
Creazione del Modello di Regressione
# Matrice di correlazione variabili quantitative
variabili_qualitative <- c("Tipo.parto", "Ospedale", "Sesso")
dati_solo_numerici <- dati[ , !(names(dati) %in% variabili_qualitative)]
# Visualizzazione grafica
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- (cor(x, y))
txt <- format(c(r, 1), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = 1.5)
}
suppressWarnings({
pairs(dati_solo_numerici, lower.panel = panel.cor, upper.panel = panel.smooth)
})
Dalla matrice di correlazione si evince che le variabili con maggiore correlazione con la variabile target Peso sono: Gestazione, Lunghezza e Cranio, tutte con correlazione positiva.
# Anasili grafica relazione variabile target con varibili categoriche
par(mfrow=c(1,3))
boxplot(Peso~Sesso)
boxplot(Peso~Ospedale)
boxplot(Peso~Tipo.parto)
Graficamente il Peso si distribuisce differentemente tra maschi e femmine (già verificato in precedenza), risulta invece sostanzialmente uguale tra ospedali e parto cesareo e naturale.
Verifichiamo queste ultime due ipotesi con un T test:
pairwise.t.test(Peso, Ospedale,
paired = F,
pool.sd = T,
p.adjust.method = "bonferroni")##
## Pairwise comparisons using t tests with pooled SD
##
## data: Peso and Ospedale
##
## osp1 osp2
## osp2 1.00 -
## osp3 0.33 0.33
##
## P value adjustment method: bonferronit.test(Peso ~ Tipo.parto)##
## 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 Ces and group Nat is not equal to 0
## 95 percent confidence interval:
## -46.27992 40.54037
## sample estimates:
## mean in group Ces mean in group Nat
## 3282.047 3284.916Entrambi i test hanno un valore del p-value superiore al livello di significatività , non ci sono quindi differenze statisticamente significative del Peso medio tra ospedali o tipo parto.
Primo modello con tutte le variabili
# Modello con tutte le variabili
mod1<-lm(Peso~., data=dati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.26 -181.53 -14.45 161.05 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6735.7960 141.4790 -47.610 < 2e-16 ***
## Anni.madre 0.8018 1.1467 0.699 0.4845
## N.gravidanze 11.3812 4.6686 2.438 0.0148 *
## Fumatrici -30.2741 27.5492 -1.099 0.2719
## Gestazione 32.5773 3.8208 8.526 < 2e-16 ***
## Lunghezza 10.2922 0.3009 34.207 < 2e-16 ***
## Cranio 10.4722 0.4263 24.567 < 2e-16 ***
## Tipo.partoNat 29.6335 12.0905 2.451 0.0143 *
## Ospedaleosp2 -11.0912 13.4471 -0.825 0.4096
## Ospedaleosp3 28.2495 13.5054 2.092 0.0366 *
## SessoM 77.5723 11.1865 6.934 5.18e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274 on 2487 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 668.7 on 10 and 2487 DF, p-value: < 2.2e-16Le variabili più significative si confermano Gestazione, Lunghezza, Cranio e Sesso, tutte con un coefficiente positivo.
Selezione del Modello Ottimale
## Procedura Stepwise AIC
mod.stepwise <- MASS::stepAIC(mod1,
direction = "both",
k=log(N))## Start: AIC=28118.61
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 36710 186779904 28111
## - Fumatrici 1 90677 186833870 28112
## - Ospedale 2 687555 187430749 28112
## - N.gravidanze 1 446244 187189438 28117
## - Tipo.parto 1 451073 187194266 28117
## <none> 186743194 28119
## - Sesso 1 3610705 190353899 28159
## - Gestazione 1 5458852 192202046 28183
## - Cranio 1 45318506 232061700 28654
## - Lunghezza 1 87861708 274604902 29074
##
## Step: AIC=28111.28
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 91599 186871503 28105
## - Ospedale 2 693914 187473818 28105
## - Tipo.parto 1 452049 187231953 28110
## <none> 186779904 28111
## - N.gravidanze 1 631082 187410986 28112
## + Anni.madre 1 36710 186743194 28119
## - Sesso 1 3617809 190397713 28151
## - Gestazione 1 5424800 192204704 28175
## - Cranio 1 45569477 232349381 28649
## - Lunghezza 1 87852027 274631931 29066
##
## Step: AIC=28104.68
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 702925 187574428 28098
## - Tipo.parto 1 444404 187315907 28103
## <none> 186871503 28105
## - N.gravidanze 1 608136 187479640 28105
## + Fumatrici 1 91599 186779904 28111
## + Anni.madre 1 37633 186833870 28112
## - Sesso 1 3601860 190473363 28145
## - Gestazione 1 5358199 192229702 28168
## - Cranio 1 45613331 232484834 28642
## - Lunghezza 1 88259386 275130889 29063
##
## Step: AIC=28098.41
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 467626 188042054 28097
## <none> 187574428 28098
## - N.gravidanze 1 648873 188223301 28099
## + Ospedale 2 702925 186871503 28105
## + Fumatrici 1 100610 187473818 28105
## + Anni.madre 1 44184 187530244 28106
## - Sesso 1 3644818 191219246 28139
## - Gestazione 1 5457887 193032315 28162
## - Cranio 1 45747094 233321522 28636
## - Lunghezza 1 87955701 275530129 29051
##
## Step: AIC=28096.81
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 188042054 28097
## - N.gravidanze 1 621053 188663107 28097
## + Tipo.parto 1 467626 187574428 28098
## + Ospedale 2 726146 187315907 28103
## + Fumatrici 1 92548 187949505 28103
## + Anni.madre 1 45366 187996688 28104
## - Sesso 1 3650790 191692844 28137
## - Gestazione 1 5477493 193519547 28161
## - Cranio 1 46098547 234140601 28637
## - Lunghezza 1 87532691 275574744 29044summary(mod.stepwise)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.37 -180.98 -15.57 163.69 2639.09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.7251 135.8036 -49.201 < 2e-16 ***
## N.gravidanze 12.4554 4.3416 2.869 0.00415 **
## Gestazione 32.3827 3.8008 8.520 < 2e-16 ***
## Lunghezza 10.2455 0.3008 34.059 < 2e-16 ***
## Cranio 10.5410 0.4265 24.717 < 2e-16 ***
## SessoM 77.9807 11.2111 6.956 4.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1327 on 5 and 2492 DF, p-value: < 2.2e-16Con la procedure stepwise tipo AIC il modello selezionato comprende
anche la variabile N.Gravidanze, l’R2 aggiustato è sostanzialmente
uguale a quello del modello con tutte le variabili.
Proviamo ora un modello con possibili interazioni tra le variabili
Gestazione e Lunghezza.
mod2<-update(mod.stepwise,~.+Gestazione:Lunghezza)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + Gestazione:Lunghezza, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1133.41 -179.98 -11.52 168.93 2652.65
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.991e+03 9.206e+02 -2.163 0.030631 *
## N.gravidanze 1.303e+01 4.321e+00 3.015 0.002594 **
## Gestazione -9.391e+01 2.481e+01 -3.785 0.000157 ***
## Lunghezza -8.476e-02 2.028e+00 -0.042 0.966661
## Cranio 1.076e+01 4.264e-01 25.234 < 2e-16 ***
## SessoM 7.225e+01 1.121e+01 6.445 1.38e-10 ***
## Gestazione:Lunghezza 2.729e-01 5.298e-02 5.151 2.79e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.3 on 2491 degrees of freedom
## Multiple R-squared: 0.7299, Adjusted R-squared: 0.7292
## F-statistic: 1122 on 6 and 2491 DF, p-value: < 2.2e-16Verifichiamo anche possibili effetti non lineari sulla variabile Gestazione:
mod3<-update(mod.stepwise,~.+I(Gestazione^2))
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2), data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1144.0 -181.5 -12.9 165.8 2661.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4646.7158 898.6322 -5.171 2.52e-07 ***
## N.gravidanze 12.5489 4.3381 2.893 0.00385 **
## Gestazione -81.2309 49.7402 -1.633 0.10257
## Lunghezza 10.3502 0.3040 34.045 < 2e-16 ***
## Cranio 10.6376 0.4282 24.843 < 2e-16 ***
## SessoM 75.7563 11.2435 6.738 1.99e-11 ***
## I(Gestazione^2) 1.5168 0.6621 2.291 0.02206 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.5 on 2491 degrees of freedom
## Multiple R-squared: 0.7276, Adjusted R-squared: 0.7269
## F-statistic: 1109 on 6 and 2491 DF, p-value: < 2.2e-16Per entrambi i modelli testati non si hanno evidenti miglioramenti dell’R2:
mod.stepwise 0.7265
mod2 0.7292 (R2 maggiore)
mod3 0.7269
Analisi della Qualità del Modello
anova(mod.stepwise,mod2,mod3)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## Gestazione:Lunghezza
## Model 3: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## I(Gestazione^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 188042054
## 2 2491 186060295 1 1981759 26.532 2.795e-07 ***
## 3 2491 187646731 0 -1586436
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Con il test ANOVA risulta che il modello 2 spiega una quantità di
varianza maggiore rispetto a mod.stepwise (p-value < 0.05).
Considerando il valore dell’R2 aggiustato, visto in precedenza, però
migliora solo dello 0,003 (0.7265 vs 0.7292).
BIC(mod1,mod2,mod3,mod.stepwise)## df BIC
## mod1 12 35215.45
## mod2 8 35175.01
## mod3 8 35196.21
## mod.stepwise 7 35193.65Secondo il criterio BIC il mod.stepwise risulta il migliore (con valore BIC più basso) in quanto modello più semplice che spiega meglio i dati.
Selezioniamo il mod.stepwise, valutiamo la multicollinearità tra i regressori ed effettuiamo l’analisi dei residui.
car::vif(mod.stepwise)## N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1.023462 1.669779 2.075747 1.624568 1.040184Tutti i valori sono inferiori a 5.
Analisi dei residui:
par(mfrow=c(2,2))
plot(mod.stepwise)
Dall’analisi grafica si vede che i residui tendono a concentrarsi in
una zona, e che i punti sul Q-Q plot agli estremi tendono ad
allontanarsi dalla retta.
Risulta inoltre l’osservazione 1551 con distanza di cook superiore a
0.5.
Verifichiamo con i test:
# test di Breusch-Pagan per omoschedasticitÃ
lmtest::bptest(mod.stepwise)##
## studentized Breusch-Pagan test
##
## data: mod.stepwise
## BP = 90.297, df = 5, p-value < 2.2e-16# test di Durbin-Watson indipendenza degli errori
lmtest::dwtest(mod.stepwise)##
## Durbin-Watson test
##
## data: mod.stepwise
## DW = 1.9532, p-value = 0.1209
## alternative hypothesis: true autocorrelation is greater than 0# test di Shapiro-Wilk distribuzione normale degli errori
shapiro.test(mod.stepwise$residuals)##
## Shapiro-Wilk normality test
##
## data: mod.stepwise$residuals
## W = 0.97414, p-value < 2.2e-16Dai test risulta che i residui soddisfano solo il requisito di
indipendenza.
Dal Q-Q plot si vedeva però che la distribuzione era tutto sommato
normale.
Rimarrebbe quindi da valutare solo l’eteroschedasticità .
Prima di apportare eventuali correzioni, verifichiamo la presenza di leverages o outliers.
# leverage
lev<-hatvalues(mod.stepwise)
plot(lev)
p<-sum(lev)
n<-length(lev)
soglia=2*p/n
abline(h=soglia,col=2)
#lev[lev>soglia]
sum(lev>soglia)## [1] 152# outliers
plot(rstudent(mod.stepwise))
abline(h=c(-2,2))
car::outlierTest(mod.stepwise)## rstudent unadjusted p-value Bonferroni p
## 1551 10.046230 2.6345e-23 6.5810e-20
## 155 5.025345 5.3818e-07 1.3444e-03
## 1306 4.824963 1.4848e-06 3.7092e-03Tra le osservazioni sono presenti 152 leverages e 3 outliers che potrebbero distorcere i residui, valutiamo però i più influenti con la distanza di Cook.
cook<-cooks.distance(mod.stepwise)
plot(cook,ylim = c(0,1))
n <- length(cook)
k <- length(attr(terms(mod.stepwise), "term.labels"))
soglia_cook <- 4 / (n - k - 1)
abline(h = soglia_cook, col = "red")
abline(h = 0.5, col = "red") 
Risulta solo un’osservazione con distanza di Cook superiore a 0,5 (la 1551 rilevata in precedenza).
df_bkp[1551, ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1551 35 1 0 38 4370 315 374
## Tipo.parto Ospedale Sesso
## 1551 Nat osp3 Fdati_femmine <- dati[dati$Sesso == "F", ]
summary(dati_femmine[, c("Peso", "Lunghezza", "Cranio")])## Peso Lunghezza Cranio
## Min. : 830 Min. :310.0 Min. :235.0
## 1st Qu.:2900 1st Qu.:480.0 1st Qu.:330.0
## Median :3160 Median :490.0 Median :340.0
## Mean :3161 Mean :489.8 Mean :337.6
## 3rd Qu.:3470 3rd Qu.:505.0 3rd Qu.:348.0
## Max. :4930 Max. :565.0 Max. :390.0Confrontando i dati con le osservazioni dei neonati di sesso
femminile si evince che la 1551 ha un peso sopra il 3° quartile, una
lunghezza sotto il 1° quartile e un diametro del cranio sopra il 3°
quartile, in pratica è più piccola della media ma con una testa più
grande, non sembra nulla di anomalo.
Prendiamo atto dell’eteroschedasticità dei residui, e confermiamo il
mod.stepwise.
3. Previsioni e risultati
Stima del peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.
# Calcola le medie per Lunghezza e Cranio
media_lunghezza <- mean(Lunghezza)
media_cranio <- mean(Cranio)
nuovi_dati_neonato <- data.frame(
N.gravidanze = 3,
Gestazione = 39,
Lunghezza = media_lunghezza,
Cranio = media_cranio,
Sesso = factor("F", levels = levels(dati$Sesso))
)
# Esegue la previsione del Peso
previsione_peso_neonato <- predict(mod.stepwise, newdata = nuovi_dati_neonato)
# Stampa il risultato
cat("Previsione del Peso alla nascita per la neonata:", round(previsione_peso_neonato, 2), "grammi\n")## Previsione del Peso alla nascita per la neonata: 3271.14 grammiprevisione_intervallo <- predict(mod.stepwise, newdata = nuovi_dati_neonato, interval = "prediction")
cat("\nPrevisione del Peso alla nascita con intervallo di confidelza al 95%:\n")##
## Previsione del Peso alla nascita con intervallo di confidelza al 95%:print(round(previsione_intervallo, 2))## fit lwr upr
## 1 3271.14 2731.98 3810.314. Visualizzazioni
Tra le variabili più siglificative abbiamo le settimane di
Gestazione.
Di seguito il grafico della relazione tra Peso e Gestazione segmentato
per Sesso:
library(ggplot2)
ggplot(data=dati)+
geom_point(aes(x=Gestazione,
y=Peso,
col=Sesso),position = "jitter")+
geom_smooth(aes(x=Gestazione,
y=Peso,
col=Sesso),se=F,method = "lm")## `geom_smooth()` using formula = 'y ~ x'
Tra le variabili non significative abbiamo Fumatrici.
Anche dal grafico si evince la mancanza di correlazione con il Peso:
ggplot(data=dati)+
geom_point(aes(x=Fumatrici,
y=Peso),position = "jitter")+
geom_smooth(aes(x=Fumatrici,
y=Peso),se=F,method = "lm")## `geom_smooth()` using formula = 'y ~ x'