Modello Statistico per la Previsione del Peso Neonatale
Simone Paghera x Neonatal Health Solutions
2026-06-07
Librerie e Caricamento Dati
library(moments) # per skewness e kurtosis
library(lmtest) # per bptest e dwtest
library(car) # per vif e outlierTest
library(MASS) # per stepAIC
library(ggplot2) # per grafici
library(knitr) # per visualizzazione tabellare
dati = read.csv("neonati.csv", na.strings = T)
str(dati)## 'data.frame': 2500 obs. of 10 variables:
## $ Anni.madre : int 26 21 34 28 20 32 26 25 22 23 ...
## $ N.gravidanze: int 0 2 3 1 0 0 1 0 1 0 ...
## $ Fumatrici : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Gestazione : int 42 39 38 41 38 40 39 40 40 41 ...
## $ Peso : int 3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
## $ Lunghezza : int 490 490 500 515 480 495 480 510 500 510 ...
## $ Cranio : int 325 345 375 365 335 340 345 349 335 362 ...
## $ Tipo.parto : chr "Nat" "Nat" "Nat" "Nat" ...
## $ Ospedale : chr "osp3" "osp1" "osp2" "osp2" ...
## $ Sesso : chr "M" "F" "M" "M" ...kable(head(dati))| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso |
|---|---|---|---|---|---|---|---|---|---|
| 26 | 0 | 0 | 42 | 3380 | 490 | 325 | Nat | osp3 | M |
| 21 | 2 | 0 | 39 | 3150 | 490 | 345 | Nat | osp1 | F |
| 34 | 3 | 0 | 38 | 3640 | 500 | 375 | Nat | osp2 | M |
| 28 | 1 | 0 | 41 | 3690 | 515 | 365 | Nat | osp2 | M |
| 20 | 0 | 0 | 38 | 3700 | 480 | 335 | Nat | osp3 | F |
| 32 | 0 | 0 | 40 | 3200 | 495 | 340 | Nat | osp2 | F |
attach(dati)Analisi delle Variabili
Le variabili presenti nel dataset sono:
- Anni.madre: Variabile QUANTITATIVA CONTINUA. Indica gli anni della madre.
- N.gravidanze: Variabile QUANTITATIVA DISCRETA. Indica un conteggio di quante gravidanze ha avuto la madre.
- Fumatrici: Variabile QUALITATIVA NOMINALE (dicotomica). Indicatore binario che indica se la madre fuma o no.
- Gestazione: Variabile QUANTITATIVA CONTINUA. Indica il numero di settimane di gestazione.
- Peso: Variabile QUANTITATIVA CONTINUA. Indica il peso alla nascita in grammi.
- Lunghezza: Variabile QUANTITATIVA CONTINUA. Indica la lunghezza del neonato in millimetri.
- Cranio: Variabile QUANTITATIVA CONTINUA. Indica il diametro del cranio del neonato in millimetri.
- Tipo.parto: Variabile QUALITATIVA NOMINALE (binaria). Indica il tipo di parto della madre.
- Ospedale: Variabile QUALITATIVA NOMINALE. Indica l’ospedale di appartenenza.
- Sesso: Variabile QUALITATIVA NOMINALE (binaria). Indica il sesso del neonato.
Indici di Posizione, Variabilità e Forma
Le variabili su cui ha senso calcolare gli indici sono:
Anni.madre, N.gravidanze,
Gestazione, Peso, Lunghezza,
Cranio. Per Fumatrici, Tipo.parto
e Sesso creo distribuzioni di frequenza.
# Funzione per indici
indici_stat = function(x, nome_var){
cat("Variabile:", nome_var, "\n")
cat("Minimo:", min(x, na.rm = T), "\n")
cat("Q1:", quantile(x, 0.25, na.rm = T), "\n")
cat("Mediana:", median(x, na.rm = T), "\n")
cat("Media:", mean(x, na.rm = T), "\n")
cat("Q3:", quantile(x, 0.75, na.rm = T), "\n")
cat("Massimo:", max(x, na.rm = T), "\n")
cat("Range:", max(x, na.rm = T)-min(x, na.rm = T), "\n")
cat("IQR:", IQR(x, na.rm = T), "\n")
cat("Dev. Standard:", sd(x, na.rm = T), "\n")
cat("Coeff. Variazione:", sd(x, na.rm = T)/mean(x, na.rm = T)*100, "\n")
cat("Assimmetria:", skewness(x, na.rm = T), "\n")
cat("Curtosi:", kurtosis(x, na.rm = T)-3)
}
indici_stat(Anni.madre, "Anni delle madri")## Variabile: Anni delle madri
## Minimo: 0
## Q1: 25
## Mediana: 28
## Media: 28.164
## Q3: 32
## Massimo: 46
## Range: 46
## IQR: 7
## Dev. Standard: 5.273578
## Coeff. Variazione: 18.72454
## Assimmetria: 0.0428115
## Curtosi: 0.3804165indici_stat(N.gravidanze, "Numero di gravidanze/madre")## Variabile: Numero di gravidanze/madre
## Minimo: 0
## Q1: 0
## Mediana: 1
## Media: 0.9812
## Q3: 1
## Massimo: 12
## Range: 12
## IQR: 1
## Dev. Standard: 1.280587
## Coeff. Variazione: 130.5123
## Assimmetria: 2.514254
## Curtosi: 10.98941indici_stat(Gestazione, "Settimane di gestazione")## Variabile: Settimane di gestazione
## Minimo: 25
## Q1: 38
## Mediana: 39
## Media: 38.9804
## Q3: 40
## Massimo: 43
## Range: 18
## IQR: 2
## Dev. Standard: 1.868639
## Coeff. Variazione: 4.793792
## Assimmetria: -2.065313
## Curtosi: 8.25815indici_stat(Peso, "Peso neonato (grammi)")## Variabile: Peso neonato (grammi)
## Minimo: 830
## Q1: 2990
## Mediana: 3300
## Media: 3284.081
## Q3: 3620
## Massimo: 4930
## Range: 4100
## IQR: 630
## Dev. Standard: 525.0387
## Coeff. Variazione: 15.98739
## Assimmetria: -0.6470308
## Curtosi: 2.031532indici_stat(Lunghezza, "Lunghezza neonato (millimetri)")## Variabile: Lunghezza neonato (millimetri)
## Minimo: 310
## Q1: 480
## Mediana: 500
## Media: 494.692
## Q3: 510
## Massimo: 565
## Range: 255
## IQR: 30
## Dev. Standard: 26.31864
## Coeff. Variazione: 5.320208
## Assimmetria: -1.514699
## Curtosi: 6.487174indici_stat(Cranio, "Diametro cranio neonato (millimetri)")## Variabile: Diametro cranio neonato (millimetri)
## Minimo: 235
## Q1: 330
## Mediana: 340
## Media: 340.0292
## Q3: 350
## Massimo: 390
## Range: 155
## IQR: 20
## Dev. Standard: 16.42533
## Coeff. Variazione: 4.830565
## Assimmetria: -0.7850527
## Curtosi: 2.946206# Distribuzioni di frequenza
freq_fumo = table(Fumatrici)
kable(as.data.frame(freq_fumo), col.names = c("Fumatrici", "Frequenza"))| Fumatrici | Frequenza |
|---|---|
| 0 | 2396 |
| 1 | 104 |
freq_parto = table(Tipo.parto)
kable(as.data.frame(freq_parto), col.names = c("Tipo di parto", "Frequenza"))| Tipo di parto | Frequenza |
|---|---|
| Ces | 728 |
| Nat | 1772 |
freq_sesso = table(Sesso)
kable(as.data.frame(freq_sesso), col.names = c("Sesso", "Frequenza"))| Sesso | Frequenza |
|---|---|
| F | 1256 |
| M | 1244 |
Commento: La variabile
Anni.madrepresenta un’anomalia in quanto ci sono due valori (0 e 1) non compatibili con una gravidanza. Vi sono due donne che hanno avuto rispettivamente 11 e 12 gravidanze: per quanto anomali, questi valori sono compatibili con la realtà. Tutte e tre le variabiliPeso,LunghezzaeCraniosono leggermente asimmetriche negative, con il peso che mostra la maggior variabilità rispetto alle altre due (CV del 16%). Le variabili con variabilità più bassa sonoGestazioneeDiametro del cranio, con un CV pari a 4.8% per entrambe.
# Boxplot per visualizzare le variabili quantitative
par(mfrow = c(2, 3))
boxplot(Anni.madre, main = "Età madre", col = "pink")
boxplot(N.gravidanze, main = "N. gravidanze", col = "pink2")
boxplot(Gestazione, main = "Gestazione (sett)", col = "pink3")
boxplot(Peso, main = "Peso (g)", col = "red")
boxplot(Lunghezza, main = "Lunghezza (mm)", col = "blue")
boxplot(Cranio, main = "Cranio (mm)", col = "green")
Test di Ipotesi
Test 1 — Differenze nei parti cesarei tra ospedali
H₀: La proporzione di parti cesarei è uguale tra gli
ospedali.
H₁: La proporzione di parti cesarei è diversa tra gli
ospedali.
freq_parto_osp = table(Ospedale, Tipo.parto)
kable(as.data.frame(freq_parto_osp), col.names = c("Ospedale", "Tipo di parto", "Frequenza"))| Ospedale | Tipo di parto | Frequenza |
|---|---|---|
| osp1 | Ces | 242 |
| osp2 | Ces | 254 |
| osp3 | Ces | 232 |
| osp1 | Nat | 574 |
| osp2 | Nat | 595 |
| osp3 | Nat | 603 |
chisq.test(freq_parto_osp)##
## Pearson's Chi-squared test
##
## data: freq_parto_osp
## X-squared = 1.0972, df = 2, p-value = 0.5778Il p-value (0.577) è maggiore di 0.05, quindi non si rifiuta l’ipotesi nulla: non ci sono evidenze statistiche di differenza tra gli ospedali.
Test 2 — Confronto del peso e della lunghezza con la popolazione
H₀: La media del campione è uguale a quella della
popolazione.
H₁: La media del campione è diversa da quella della
popolazione.
Sulla base dei dati raccolti da Stanford Medicine, il peso medio di un neonato a termine (nato tra la 37ª e 41ª settimana) è di 3200 g, mentre la lunghezza media è di 500 mm. Questi valori vengono usati come riferimento.
# Test t per il peso
t.test(Peso,
mu = 3200,
conf.level = 0.95,
alternative = "two.sided")##
## One Sample t-test
##
## data: Peso
## t = 8.0071, df = 2499, p-value = 1.782e-15
## alternative hypothesis: true mean is not equal to 3200
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081Il t-test ha restituito un p-value < 0.05. Questo risultato permette di rifiutare l’ipotesi nulla, suggerendo che i dati del campione non sono compatibili con la media della popolazione di riferimento.
# Test t per la lunghezza
t.test(Lunghezza,
mu = 500,
conf.level = 0.95,
alternative = "two.sided")##
## One Sample t-test
##
## data: 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.692Anche in questo caso il t-test ha restituito un p-value < 0.05. Pertanto si rifiuta l’ipotesi nulla, indicando come la lunghezza dei neonati del campione sia differente da quella della popolazione di riferimento.
Test 3 — Misure antropometriche per sesso
H₀: Il cranio dei maschi è uguale al cranio delle femmine H₁: Il cranio dei maschi è diverso dal cranio delle femmine
t.test(data = dati,
Cranio ~ Sesso,
alternative = "two.sided")##
## 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.4486Il t test ha restituito un p value < 0.05. Questo consente di rifiutare l’ipotesi nulla, affermando che il diametro del cranio differisce tra i sessi.
H₀: La lunghezza dei maschi è uguale alla lunghezza delle femmine H₁: La lunghezzadei maschi è diversa dalla lunghezza delle femmine
t.test(data=dati,
Lunghezza~Sesso,
alternative ="two.sided")##
## 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.6672Il t test ha restituito un p value < 0.05. Questo consente di rifiutare l’ipotesi nulla, affermando che la lunghezza differisce tra i sessi.
H₀: Il peso dei maschi è uguale al peso delle femmine H₁: Il peso dei maschi è diverso dal peso delle femmine
t.test(data=dati,
Peso~Sesso,
alternative ="two.sided")##
## 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.215Anche in questo ultimo caso, il t test ha restituito un p value < 0.05. Questo consente di rifiutare l’ipotesi nulla, affermando che il peso differisce tra i sessi.
par(mfrow = c(1, 3))
boxplot(Cranio~Sesso, data=dati, col="lightblue", main="Diametro del cranio (mm)")
boxplot(Lunghezza~Sesso, data=dati, col="lightpink", main="Lunghezza (mm)")
boxplot(Peso~Sesso, data=dati, col="lightgreen", main="Peso (g)")
Modello di Regressione Lineare Multipla
Come prima cosa verifichiamo che la variabile risposta
Peso sia distribuita normalmente.
moments::skewness(Peso)## [1] -0.6470308moments::kurtosis(Peso) - 3## [1] 2.031532shapiro.test(Peso)##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16Sulla base del test di Shapiro, la variabile
Pesonon si distribuisce normalmente. Si prosegue comunque con il modello, dato che è importante che siano i residui a distribuirsi normalmente.
Matrice di Correlazione
# Seleziono solo le variabili numeriche
dati_num = dati[, c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")]
kable(round(cor(dati_num), 2), caption = "Matrice di correlazione")| Anni.madre | N.gravidanze | Gestazione | Peso | Lunghezza | Cranio | |
|---|---|---|---|---|---|---|
| Anni.madre | 1.00 | 0.38 | -0.14 | -0.02 | -0.06 | 0.02 |
| N.gravidanze | 0.38 | 1.00 | -0.10 | 0.00 | -0.06 | 0.04 |
| Gestazione | -0.14 | -0.10 | 1.00 | 0.59 | 0.62 | 0.46 |
| Peso | -0.02 | 0.00 | 0.59 | 1.00 | 0.80 | 0.70 |
| Lunghezza | -0.06 | -0.06 | 0.62 | 0.80 | 1.00 | 0.60 |
| Cranio | 0.02 | 0.04 | 0.46 | 0.70 | 0.60 | 1.00 |
# Utilizzo la formula presente nel help della funzione par
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)
}
pairs(dati_num, upper.panel = panel.smooth, lower.panel = panel.cor)
Guardando la matrice risulta evidente che le variabili
Gestazione,Peso,LunghezzaeCraniorisultano ovviamente correlate: all’aumentare dei mesi di gestazione aumentano di conseguenza il peso, la lunghezza e la dimensione del cranio. Viceversa, gli anni della madre e il numero di gravidanze non hanno un effetto evidente su peso, lunghezza e dimensione del cranio.
Costruzione dei Modelli
Si escludono a priori Tipo.parto e Ospedale
in quanto non collegate all’obiettivo del modello. La variabile
Fumatrici viene resa dicotomica.
dati$Fumatrici = factor(dati$Fumatrici, levels = c(0, 1))
# Modello 1: tutte le variabili
mod1 = lm(Peso ~ Anni.madre + N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso + Fumatrici, data = dati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione +
## Lunghezza + Cranio + Sesso + Fumatrici, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1161.56 -181.19 -15.75 163.70 2630.75
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6714.4109 141.1515 -47.569 < 2e-16 ***
## Anni.madre 0.9585 1.1347 0.845 0.3984
## N.gravidanze 11.2756 4.6690 2.415 0.0158 *
## Gestazione 32.9331 3.8267 8.606 < 2e-16 ***
## Lunghezza 10.2342 0.3009 34.009 < 2e-16 ***
## Cranio 10.5177 0.4268 24.642 < 2e-16 ***
## SessoM 78.0845 11.2039 6.969 4.06e-12 ***
## Fumatrici1 -30.2959 27.5971 -1.098 0.2724
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2492 degrees of freedom
## Multiple R-squared: 0.7272, Adjusted R-squared: 0.7264
## F-statistic: 949 on 7 and 2492 DF, p-value: < 2.2e-16Risultati del modello 1: 1. R² aggiustato = 0.7264 2.
Anni.madrenon ha nessun impatto significativo sul peso 3.N.gravidanzeè leggermente significativo (): ogni gravidanza il peso aumenta di 11 g 4.Gestazione,LunghezzaeCraniosono molto significativi (**), confermando l’aspettativa logica 5. C’è una differenza significativa nei maschi rispetto alle femmine di 78 g 6. Le madri fumatrici partoriscono neonati con un peso medio inferiore di 30 g, senza significatività statistica
# Modello 2: rimozione di Anni.madre e Fumatrici
mod2 = update(mod1, ~. - Anni.madre - Fumatrici)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## 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 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## 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 ***
## ---
## 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-16# Modello 3: rimozione anche di N.gravidanze
mod3 = update(mod2, ~. - N.gravidanze)
summary(mod3)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso,
## data = dati)
##
## 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 ***
## SessoM 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-16R² rimane lo stesso in
mod2(0.7265), mentre cala a 0.7257 inmod3. Sulla base del principio del Rasoio di Occam, a parità di R², si preferiscono modelli più semplici: ad oramod2è il migliore.
kable(data.frame(AIC = AIC(mod1, mod2, mod3)$AIC,
BIC =BIC(mod1, mod2, mod3)$BIC,
row.names = c("mod1","mod2","mod3")),
caption = "Confronto AIC e BIC")| AIC | BIC | |
|---|---|---|
| mod1 | 35181.39 | 35233.81 |
| mod2 | 35179.33 | 35220.10 |
| mod3 | 35185.60 | 35220.54 |
mod2conferma i valori di AIC e BIC più bassi (35179.33 e 35220.10 rispettivamente).
# Modello 4: interazione Gestazione*Lunghezza
mod4 = lm(Peso ~ N.gravidanze + Gestazione*Lunghezza + Cranio + Sesso, data = dati)
summary(mod4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione * Lunghezza + Cranio +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1133.50 -179.45 -11.74 168.77 2653.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.991e+03 9.202e+02 -2.163 0.030619 *
## N.gravidanze 1.304e+01 4.319e+00 3.020 0.002551 **
## Gestazione -9.396e+01 2.480e+01 -3.789 0.000155 ***
## Lunghezza -8.111e-02 2.027e+00 -0.040 0.968082
## Cranio 1.076e+01 4.262e-01 25.245 < 2e-16 ***
## SessoM 7.228e+01 1.120e+01 6.454 1.31e-10 ***
## Gestazione:Lunghezza 2.729e-01 5.295e-02 5.153 2.76e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.2 on 2493 degrees of freedom
## Multiple R-squared: 0.7299, Adjusted R-squared: 0.7292
## F-statistic: 1123 on 6 and 2493 DF, p-value: < 2.2e-16# Modello 5: interazione Gestazione*Cranio
mod5 = lm(Peso ~ N.gravidanze + Lunghezza + Gestazione*Cranio + Sesso, data = dati)
summary(mod5)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Lunghezza + Gestazione * Cranio +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1137.16 -180.86 -12.44 167.43 2696.03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -193.11866 1106.42819 -0.175 0.86145
## N.gravidanze 13.14181 4.31191 3.048 0.00233 **
## Lunghezza 10.47042 0.30096 34.790 < 2e-16 ***
## Gestazione -140.67784 29.52621 -4.765 2.00e-06 ***
## Cranio -9.83681 3.47497 -2.831 0.00468 **
## SessoM 72.05680 11.17200 6.450 1.34e-10 ***
## Gestazione:Cranio 0.53339 0.09028 5.908 3.94e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 272.8 on 2493 degrees of freedom
## Multiple R-squared: 0.7308, Adjusted R-squared: 0.7301
## F-statistic: 1128 on 6 and 2493 DF, p-value: < 2.2e-16kable(data.frame(AIC = AIC(mod1, mod2, mod3, mod4, mod5)$AIC,
BIC =BIC(mod1, mod2, mod3, mod4, mod5)$BIC,
row.names = c("mod1","mod2","mod3", "mod4", "mod5")),
caption = "Confronto AIC e BIC")| AIC | BIC | |
|---|---|---|
| mod1 | 35181.39 | 35233.81 |
| mod2 | 35179.33 | 35220.10 |
| mod3 | 35185.60 | 35220.54 |
| mod4 | 35154.84 | 35201.44 |
| mod5 | 35146.57 | 35193.16 |
I valori di AIC e BIC indicano
mod5come il migliore tra quelli proposti.
Valutazione della Multicollinearità
vif(mod5)## N.gravidanze Lunghezza Gestazione Cranio
## 1.024177 2.107498 102.254853 109.432847
## Sesso Gestazione:Cranio
## 1.048535 301.445595I valori di VIF per
GestazioneeCraniosono maggiori di 5, indicando alta correlazione tra le variabili e difficoltà nello stimare separatamente il loro effetto. Il modello migliore resta quindimod2.
Selezione Stepwise con AIC
n = nrow(dati)
stepwise.mod = MASS::stepAIC(mod1,
direction = "both",
k = log(n))## Start: AIC=28131.29
## Peso ~ Anni.madre + N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + Fumatrici
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 53803 187973654 28124
## - Fumatrici 1 90879 188010731 28125
## - N.gravidanze 1 439797 188359648 28129
## <none> 187919851 28131
## - Sesso 1 3662833 191582684 28172
## - Gestazione 1 5585184 193505035 28197
## - Cranio 1 45791026 233710877 28669
## - Lunghezza 1 87220887 275140738 29077
##
## Step: AIC=28124.18
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## Fumatrici
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 91892 188065546 28118
## <none> 187973654 28124
## - N.gravidanze 1 646039 188619694 28125
## + Anni.madre 1 53803 187919851 28131
## - Sesso 1 3671289 191644943 28165
## - Gestazione 1 5531705 193505359 28189
## - Cranio 1 46066755 234040410 28664
## - Lunghezza 1 87218857 275192512 29069
##
## Step: AIC=28117.58
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 188065546 28118
## - N.gravidanze 1 623141 188688687 28118
## + Fumatrici 1 91892 187973654 28124
## + Anni.madre 1 54816 188010731 28125
## - Sesso 1 3655292 191720838 28158
## - Gestazione 1 5464853 193530399 28181
## - Cranio 1 46108583 234174130 28658
## - Lunghezza 1 87632762 275698308 29066summary(stepwise.mod)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## 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 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## 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 ***
## ---
## 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-16La funzione stepAIC conferma effettivamente
mod2come il modello ottimale.
Analisi dei Residui
par(mfrow = c(2, 2))
plot(mod2)
Il Q-Q plot dei residui evidenzia come le code tendano a perdere linearità, con un valore molto lontano dalla retta (osservazione 1551). Anche il grafico della distanza di Cook evidenzia come questo valore sia al di sopra della soglia di 0.5.
Valori di Leva e Outlier
lev = hatvalues(mod2)
par(mfrow = c(1, 1))
plot(lev)
p = sum(lev)
soglia = 2 * p / n
abline(h = soglia, col = 2)
lev[lev > soglia]## 13 15 34 67 89 96
## 0.005630918 0.007050422 0.006744143 0.005891590 0.012816470 0.005351586
## 101 106 131 134 151 155
## 0.007526751 0.014479871 0.007229339 0.007552819 0.010883937 0.007207682
## 161 189 190 204 205 206
## 0.020335483 0.004893297 0.005366557 0.014490604 0.005351634 0.009476652
## 220 294 305 310 312 315
## 0.007393997 0.005912765 0.005442061 0.028812123 0.013169272 0.005385800
## 378 440 442 445 486 492
## 0.015934061 0.005404736 0.007723662 0.007509382 0.005164446 0.008274018
## 497 516 582 587 592 614
## 0.005166306 0.013079851 0.011665555 0.008412325 0.006384116 0.005299262
## 638 656 657 684 697 702
## 0.006688287 0.005927777 0.005322685 0.008818987 0.005863826 0.005202259
## 729 748 750 757 765 805
## 0.005023115 0.008565543 0.006942097 0.008145491 0.006070298 0.014356657
## 828 893 895 913 928 946
## 0.007179817 0.005075205 0.005295896 0.005571144 0.022742332 0.006909044
## 947 956 985 1008 1014 1049
## 0.008409465 0.007784123 0.007039416 0.005343037 0.008470133 0.004956169
## 1067 1091 1106 1130 1166 1181
## 0.008465430 0.008933360 0.005967317 0.031728597 0.005513559 0.005677676
## 1188 1200 1219 1238 1248 1273
## 0.006477203 0.005492370 0.030694311 0.005908078 0.014622914 0.007085831
## 1291 1293 1311 1321 1325 1356
## 0.006117497 0.006073639 0.009625908 0.009293111 0.004857169 0.005303442
## 1357 1385 1395 1400 1402 1411
## 0.006965051 0.012636943 0.005126697 0.005925069 0.004811441 0.008048184
## 1420 1428 1429 1450 1505 1551
## 0.005155654 0.008192811 0.021757172 0.015104831 0.013330439 0.048769569
## 1553 1556 1573 1593 1606 1610
## 0.008504889 0.005919673 0.005047204 0.005623758 0.005001812 0.008722184
## 1617 1619 1628 1686 1693 1701
## 0.004866796 0.015067498 0.005069731 0.009349313 0.005077858 0.010842957
## 1712 1718 1727 1735 1780 1781
## 0.006992084 0.006958857 0.013300523 0.004884846 0.025538678 0.016832361
## 1809 1827 1868 1892 1962 1967
## 0.008707504 0.006065698 0.005205637 0.005332812 0.005540442 0.005337356
## 1977 2037 2040 2046 2086 2089
## 0.006927281 0.004889127 0.011494872 0.005471670 0.013193090 0.006293550
## 2098 2114 2115 2120 2140 2146
## 0.005094455 0.013316875 0.011772090 0.018659995 0.006244232 0.005802168
## 2148 2149 2157 2175 2200 2215
## 0.007926839 0.013583436 0.005907225 0.032527273 0.011670024 0.004892265
## 2216 2220 2221 2224 2225 2244
## 0.008117864 0.005414040 0.021628717 0.005838076 0.005591261 0.006929217
## 2257 2307 2317 2318 2337 2359
## 0.006170254 0.013965608 0.007673614 0.004831118 0.005230450 0.010067364
## 2408 2422 2436 2437 2452 2458
## 0.009696691 0.021532808 0.004986522 0.023943328 0.023838489 0.008506087
## 2471 2478
## 0.020903740 0.005775173outlierTest(mod2)## rstudent unadjusted p-value Bonferroni p
## 1551 10.051908 2.4906e-23 6.2265e-20
## 155 5.027798 5.3138e-07 1.3285e-03
## 1306 4.827238 1.4681e-06 3.6702e-03I valori outlier identificati sono: 1551, 155 e 1306.
Distanza di Cook
cook = cooks.distance(mod2)
plot(cook)
max(cook)## [1] 0.8300965Il valore in corrispondenza dell’osservazione 1551 presenta una distanza di Cook pari a 0.83, l’unica al di sopra della soglia di 0.5: questo singolo valore è in grado di influenzare il modello.
Test sui Residui
bptest(mod2)##
## studentized Breusch-Pagan test
##
## data: mod2
## BP = 90.253, df = 5, p-value < 2.2e-16dwtest(mod2)##
## Durbin-Watson test
##
## data: mod2
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(residuals(mod2))##
## Shapiro-Wilk normality test
##
## data: residuals(mod2)
## W = 0.97408, p-value < 2.2e-16
- Il test di Breusch–Pagan evidenzia la presenza di eteroschedasticità (p < 0.001): la varianza degli errori non è costante nel modello.
- Il test di Durbin–Watson (DW = 1.95, p = 0.12) non evidenzia autocorrelazione, indicando indipendenza degli errori.
- Il test di Shapiro–Wilk suggerisce una deviazione dalla normalità (p < 0.001), verosimilmente influenzata dalla numerosità campionaria.
Previsioni con il Modello
# Caso 1: donna alla terza gravidanza, 39 settimane, neonata femmina (L=475, C=322)
newdati = data.frame(N.gravidanze = 3,
Gestazione = 39,
Lunghezza = 475,
Cranio = 322,
Sesso = "F")
predict(mod2, newdata = newdati)## 1
## 2879.244Il peso predetto è di circa 2879 grammi.
# Caso 2: neonato maschio, 500 mm, cranio 374, prima gravidanza, 40 settimane
newdati2 = data.frame(N.gravidanze = 1,
Gestazione = 40,
Lunghezza = 500,
Cranio = 374,
Sesso = "M")
predict(mod2, newdata = newdati2)## 1
## 3768.923Il peso previsto è di circa 3768 grammi.
Visualizzazioni
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")
È evidente come, per entrambi i sessi, l’aumentare delle settimane di gestazione porti ad un aumento del peso del neonato.
ggplot(data = dati) +
geom_boxplot(aes(x = Fumatrici, y = Peso, fill = Fumatrici))
Mediante il boxplot è possibile apprezzare come le madri non fumatrici concepiscano in media neonati con un peso maggiore.