Modello Statistico per la Previsione del Peso Neonatale
Vincenzo Avellino
# Import dataframe
dati = read.csv('neonati.csv')Esplorazione del dataframe
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" ...n = nrow(dati) # memorizzo il numero delle osservazioni nella variabile nIl data frame contiene 2500 osservazioni di 10 variabili.
# controllo valori mancanti
anyNA(dati)## [1] FALSEsummary(dati)## Anni.madre N.gravidanze Fumatrici Gestazione
## Min. : 0.00 Min. : 0.0000 Min. :0.0000 Min. :25.00
## 1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:38.00
## Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00
## Mean :28.16 Mean : 0.9812 Mean :0.0416 Mean :38.98
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:40.00
## Max. :46.00 Max. :12.0000 Max. :1.0000 Max. :43.00
## Peso Lunghezza Cranio Tipo.parto
## Min. : 830 Min. :310.0 Min. :235 Length:2500
## 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Class :character
## Median :3300 Median :500.0 Median :340 Mode :character
## Mean :3284 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :565.0 Max. :390
## Ospedale Sesso
## Length:2500 Length:2500
## Class :character Class :character
## Mode :character Mode :character
##
##
## Descrizione delle variabili:
| Variabile | Tipo | Descrizione |
|---|---|---|
| Anni.madre | Quantitativa Discreta | Età in anni (conteggio intero) |
| N.gravidanze | Quantitativa Discreta | Numero di gravidanze precedenti |
| Fumatrici | Categoriale Dummy | 0 = non fumatrice, 1 = fumatrice |
| Gestazione | Quantitativa Discreta | Settimane di gestazione |
| Peso | Quantitativa Continua | Peso del neonato in grammi |
| Lunghezza | Quantitativa Continua | Lunghezza del neonato in millimetri |
| Cranio | Quantitativa Continua | Circonferenza cranica in millimetri |
| Tipo.parto | Categoriale | Tipo di parto (es. naturale) |
| Ospedale | Categoriale | Ospedale di nascita |
| Sesso | Categoriale | Sesso del neonato (M/F) |
Il peso medio dei neonati è circa 3284g, con un range da 830g a 4930g.
Lunghezza media: circa 495 mm.
Le madri hanno in media 28 anni e 1 gravidanza.
Dagli indici di posizione si nota il valore minimo della variabile ‘Anni.madre’ anomalo.
Outlier o Anomalie
par(mfrow = c(1, 5))
boxplot(dati$Anni.madre, main= "Anni Madre")
boxplot(dati$Gestazione , main= "Gestazione")
boxplot(dati$Peso , main= "Peso")
boxplot(dati$Lunghezza , main= "Lunghezza")
boxplot(dati$Cranio , main= "Cranio")
Nel primo boxplot “Anni.madre” sono evidentemente isolati i due valori vicino allo 0.
# Ordine crescente e visualizzo i primi valori della variabile Anni.madre
head(sort(dati$Anni.madre))## [1] 0 1 13 14 14 15Le due osservazioni con valori 0 e 1 sono anomale. Potrebbe trattarsi di rilevazione errata. Sostituiamo i valori con la media della stessa variabile-
Anche le altre variabili presentano outliers ma a differenza di Anni.madre non è possibile definirli errori di relavazione.
# media escludendo i due valori anomali
media_anni = mean(dati$Anni.madre[dati$Anni.madre>1])
# sostituzione con il valore medio
dati$Anni.madre[dati$Anni.madre <= 1] = media_anni
boxplot(dati$Anni.madre, main= "Anni Madre")
Statistiche Test
Test 1:
PARTI CESAREI
In alcuni ospedali si fanno più parti cesarei
Ipotesi nulla (H₀) :La proporzione di parti cesarei è uguale in tutti gli ospedali
Ipotesi alternativa (H₁):La proporzione di parti cesarei varia tra ospedali.
chisq.test(table(dati$Ospedale, dati$Tipo.parto))##
## Pearson's Chi-squared test
##
## data: table(dati$Ospedale, dati$Tipo.parto)
## X-squared = 1.0972, df = 2, p-value = 0.5778Osservazioni:
P-value signifativamente maggiore del valore soglia del 5%: mon ci sono differenze significative nella frequenza di parti cesarei tra ospedali.
Ipotesi nulla non rifiutabile.
Test 2:
La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione?
PESO
Ipotesi nulla (H₀): la media del peso del campione è uguale dalla media della popolazione
\[ H_0: \mu(peso) = \mu_0 \]Ipotesi alternativa (H₁): la media del peso del campione è diversa alla media della popolazione
\[ H_1: \mu(peso) \ne \mu_0 \]
# fonte dati "ospedale bambino gesu"
P0_Peso = 3300
# media del campione
P_cap_peso = mean(dati$Peso)
# Analisi Normalità
shapiro.test(dati$Peso) ##
## Shapiro-Wilk normality test
##
## data: dati$Peso
## W = 0.97066, p-value < 2.2e-16# T - test
t.test(dati$Peso, mu = P0_Peso, conf.level = 0.95)##
## One Sample t-test
##
## data: dati$Peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081plot(density(dati$Peso), main = "Peso",
xlab = "Peso (g)", col = "darkblue", lwd = 2)
abline(v = P0_Peso, col = "orange", lwd = 1) # media popolazione
abline(v = P_cap_peso, col = "green", lwd = 1) # media campione
abline(v = c(3263.490, 3304.672), col = "grey", lty = 2, lwd = 1)
legend("topleft", legend = c("Media Popolazione", "Media Campione", "Int. Conf."),
col = c("orange", "green","grey"), lwd = 2)
Osservazioni:
Il p-value è maggiore di 0.05, non possiamo rifiutare l’ipotesi nulla a un livello di significatività del 5%. Non ci sono prove sufficienti per affermare che il peso medio del campione sia diverso da 3300. P-value = 0.1296.
La media del campione rappresenta quella della popolazione.
LUNGHEZZA
Ipotesi nulla (H₀): la media della lunghezza del campione è uguale dalla media della popolazione
\[ H_0: \mu(lungheza) = \mu_0 \]Ipotesi alternativa (H₁): la media della lunghezza del campione è diversa alla media della popolazione
\[ H_1: \mu(lungheza) \ne \mu_0 \]
# fonte dati "ospedale bambino gesu"
P0_Lungh = 500
# media del campione
P_cap_Lungh = mean(dati$Lunghezza)
# Analisi Normalità
shapiro.test(dati$Lunghezza)##
## Shapiro-Wilk normality test
##
## data: dati$Lunghezza
## W = 0.90941, p-value < 2.2e-16# T - test
t.test(dati$Lunghezza, mu = P0_Lungh, conf.level = 0.95)##
## One Sample t-test
##
## data: dati$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.692plot(density(dati$Lunghezza), col = "darkblue", main = "Lunghezza",
xlab = "Lunghezza (mm)", lwd = 2)
abline(v = P0_Lungh, col = "orange", lwd = 1) # media popolazione
abline(v = P_cap_Lungh, col = "red", lwd = 1) # media campione
abline(v = c(493.6598, 495.7242), col = "grey", lty = 2, lwd = 1)
legend("topleft", legend = c("Media Popolazione", "Media Campione","Int. Conf."),
col = c("orange", "red", "grey"), lwd = 2)
Osservazioni:
Il p-value è significativamente minore di 0.05, possiamo rifiutare l’ipotesi nulla. La media del campione è diversa da 500. p-value < 2.2e-16
La media del campione non rappresenta la popolazione
Test 3:
Le misure antropometriche sono significativamente diverse tra i due sessi?
par(mfrow = c(1, 3))
boxplot(Peso~Sesso, data = dati, main = "Peso")
boxplot(Lunghezza~Sesso, data = dati, main = "Lunghezza")
boxplot(Cranio~Sesso, data = dati, main = "Cranio")
Ipotesi nulla (H₀): la media delle misure antropometriche di Maschi è uguale dalla media delle misure antropometrich delle Femmine \[ H_0: \mu(M) = \mu(F) \]
Ipotesi alternativa (H₁): la media delle misure antropometriche di Maschi è diversa dalla media delle misure antropometrich delle Femmine \[ H_1: \mu(M) \ne \mu(F) \]
# Test peso
t.test(Peso ~ Sesso, data = dati)##
## Welch Two Sample t-test
##
## data: Peso by Sesso
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -287.1051 -207.0615
## sample estimates:
## mean in group F mean in group M
## 3161.132 3408.215# Test lunghezza
t.test(Lunghezza ~ Sesso, data = dati)##
## Welch Two Sample t-test
##
## data: Lunghezza by Sesso
## t = -9.582, df = 2459.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -11.929470 -7.876273
## sample estimates:
## mean in group F mean in group M
## 489.7643 499.6672# Test Cranio
t.test(Cranio ~ Sesso, data = dati)##
## 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.4486Osservazioni:
| Misura | Femmine (mu) | Maschi (mu) |
|---|---|---|
| Peso | 3161.132 | 3408.215 |
| Lunghezza | 489.80 | 499.70 |
| Cranio | 337.60 | 345.40 |
Tutte le misure antropometriche (peso, lunghezza, cranio) risultano significativamente maggiori nei maschi rispetto alle femmine, con p-value estremamente bassi.
Rifiutiamo l’ipotesi nulla a favore dell’ipotesi alternativa.
Modello di Regressione
Trasformazione delle variabili categoriche in variabili dummy (Tipo.Parto, Ospedale e Sesso)
# ces = 1 Nat = 0
dati$Tipo.parto = ifelse(dati$Tipo.parto == "Ces", 1, 0)
# M = 1 F = 0
dati$Sesso = ifelse(dati$Sesso == "M", 1, 0)
dati$Ospedale = as.numeric(as.factor(dati$Ospedale))Variabile risposta: Peso
Analisi della distribuzione
# Asimmetria e Curtosi
moments::skewness(dati$Peso)## [1] -0.6470308moments::kurtosis(dati$Peso)-3## [1] 2.031532# test di normalità
shapiro.test(dati$Peso)##
## Shapiro-Wilk normality test
##
## data: dati$Peso
## W = 0.97066, p-value < 2.2e-16Osservazioni:
Asimmetrica negatiava e leptocurtica Si rifiuta l’ipotesi di normalità: P-value minore del livello di significatività fissato al 5%
In questo caso procediamo comunque con la costruzione del modello in quanto si tratta di un modello con un numero elevato di osservazioni
Matrice di correlazione
round(cor(dati), 2)## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## Anni.madre 1.00 0.38 0.01 -0.13 -0.02 -0.06
## N.gravidanze 0.38 1.00 0.05 -0.10 0.00 -0.06
## Fumatrici 0.01 0.05 1.00 0.03 -0.02 -0.02
## Gestazione -0.13 -0.10 0.03 1.00 0.59 0.62
## Peso -0.02 0.00 -0.02 0.59 1.00 0.80
## Lunghezza -0.06 -0.06 -0.02 0.62 0.80 1.00
## Cranio 0.02 0.04 -0.01 0.46 0.70 0.60
## Tipo.parto 0.00 0.02 -0.02 0.01 0.00 0.04
## Ospedale 0.03 0.01 -0.02 0.02 0.03 0.01
## Sesso 0.01 0.02 0.01 0.13 0.24 0.19
## Cranio Tipo.parto Ospedale Sesso
## Anni.madre 0.02 0.00 0.03 0.01
## N.gravidanze 0.04 0.02 0.01 0.02
## Fumatrici -0.01 -0.02 -0.02 0.01
## Gestazione 0.46 0.01 0.02 0.13
## Peso 0.70 0.00 0.03 0.24
## Lunghezza 0.60 0.04 0.01 0.19
## Cranio 1.00 0.00 0.01 0.15
## Tipo.parto 0.00 1.00 -0.02 0.00
## Ospedale 0.01 -0.02 1.00 0.01
## Sesso 0.15 0.00 0.01 1.00panel.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, upper.panel = panel.smooth, lower.panel = panel.cor)
Osservazioni:
| Variabile | Correlazione con il peso |
|---|---|
| Gestazione | +0.59 forte positiva |
| Lunghezza | +0.80 molto forte |
| Cranio | +0.70 forte positiva |
| Sesso | +0.24 leggera positiva |
| Altre (Anni.madre, Fumatrici…) | trascurabili |
Il peso del neonato è fortemente associato a variabili biometriche e di gestazione rispetto a fattori materni come età, numero di gravidanze o fumo.
Creazione del modello di regressione
Modello 1
Iniziamo con il modello base contenente tutte le variabili
mod1 = lm(dati$Peso ~ ., data=dati)
summary(mod1)##
## Call:
## lm(formula = dati$Peso ~ ., data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1139.63 -182.97 -15.91 163.33 2617.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6729.1513 141.4096 -47.586 < 2e-16 ***
## Anni.madre 0.7956 1.1472 0.693 0.4881
## N.gravidanze 11.7043 4.6681 2.507 0.0122 *
## Fumatrici -30.5268 27.5598 -1.108 0.2681
## Gestazione 32.6751 3.8201 8.553 < 2e-16 ***
## Lunghezza 10.2754 0.3008 34.161 < 2e-16 ***
## Cranio 10.4852 0.4263 24.594 < 2e-16 ***
## Tipo.parto -29.8162 12.0931 -2.466 0.0137 *
## Ospedale 14.1345 6.7536 2.093 0.0365 *
## Sesso 77.9335 11.1850 6.968 4.11e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.1 on 2490 degrees of freedom
## Multiple R-squared: 0.7284, Adjusted R-squared: 0.7274
## F-statistic: 741.8 on 9 and 2490 DF, p-value: < 2.2e-16Osservazioni:
| Variabile | Coeff. | p-valore | Interpretazione |
|---|---|---|---|
| Gestazione | +32.7 | < 2e-16 | Ogni settimana in più +32.7g |
| Lunghezza | +10.3 | < 2e-16 | Ogni mm in più +10.3g |
| Cranio | +10.5 | < 2e-16 | Ogni mm in più +10.5g |
| Sesso | +77.9 | 4.12e-12 | Differenza di peso tra sessi: mediamente +78g |
| Anni madre | +0.89 | 0.43 | Non significativo |
| N. gravidanze | +11.6 | 0.013 | Ogni gravidanza +11.6g |
| Fumatrici | -30.5 | 0.268 | Non significativo |
| Tipo parto | -29.8 | 0.014 | L’effetto (es. parto cesareo) riduce peso di ~30g |
| Ospedale | +14.1 | 0.037 | Alcuni ospedali associati a peso maggiore |
R² = 0.7284: Il modello spiega il 72.8% della variabilità del peso neonatale.
R² aggiustato = 0.7274: Molto vicino al R², quindi il modello non è sovradimensionato.
Modello 2
Rimuoviamo le variabili indipendenti Anni.madre e Fumatrici data la bassa significatività
mod2 = update(mod1, ~. -Anni.madre -Fumatrici)
summary(mod2)##
## Call:
## lm(formula = dati$Peso ~ N.gravidanze + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1129.55 -183.36 -16.91 163.67 2625.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6702.0240 135.9906 -49.283 < 2e-16 ***
## N.gravidanze 12.6481 4.3338 2.918 0.00355 **
## Gestazione 32.1378 3.7919 8.475 < 2e-16 ***
## Lunghezza 10.2891 0.3005 34.240 < 2e-16 ***
## Cranio 10.5050 0.4257 24.675 < 2e-16 ***
## Tipo.parto -29.5917 12.0901 -2.448 0.01445 *
## Ospedale 14.4072 6.7493 2.135 0.03289 *
## Sesso 77.8261 11.1827 6.960 4.35e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.1 on 2492 degrees of freedom
## Multiple R-squared: 0.7282, Adjusted R-squared: 0.7274
## F-statistic: 953.7 on 7 and 2492 DF, p-value: < 2.2e-16Osservazioni:
Aumenta la significatività per la variabile N.gravidanze con un p-value più vicino allo 0.
Anche l’R² aggiustato non subisce variazioni quindi possiamo dire che le variabile Anni.madre e Fumatrici sono inutile alla spiegazione della variabile risposta.
Il p-value per tutte le variabili è inferiore al 5%, tutte significative.
Modello 3
Viusualizziamo la correalzione tra Peso e Gestazione.
plot(dati$Peso, dati$Gestazione, pch = 20)
Osservazioni:
Dallo scatter plot si nota una leggera curvatura nella relazione tra Peso e Gestazione che potrebbe suggerire un effetto non lineare.
Proviamo ad aggiungere l’effetto quadratico della variabile Gestazione.
mod3 = update(mod2, ~ . + I(Gestazione^2))
summary(mod3)##
## Call:
## lm(formula = dati$Peso ~ N.gravidanze + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso + I(Gestazione^2),
## data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1141.6 -182.9 -13.4 163.0 2647.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4730.1032 897.0269 -5.273 1.46e-07 ***
## N.gravidanze 12.7335 4.3306 2.940 0.00331 **
## Gestazione -77.9332 49.6383 -1.570 0.11654
## Lunghezza 10.3902 0.3037 34.214 < 2e-16 ***
## Cranio 10.5988 0.4275 24.794 < 2e-16 ***
## Tipo.parto -29.1337 12.0823 -2.411 0.01597 *
## Ospedale 14.2323 6.7444 2.110 0.03494 *
## Sesso 75.6899 11.2150 6.749 1.85e-11 ***
## I(Gestazione^2) 1.4695 0.6607 2.224 0.02624 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared: 0.7287, Adjusted R-squared: 0.7278
## F-statistic: 836.4 on 8 and 2491 DF, p-value: < 2.2e-16Osservazioni:
Si inverte il valore della variabile Gestazione (possibile collinearità). P-value sotto il 5% di significatività, spiega che esiste una leggera curvatura. R² subisce una variazione trascurabile, quindi l’aggiunta risulta inutile.
Modello 4
Costruiamo un quarto modello, partendo dal modello 2 poichè al momento risulta il più performante e aggiungiamo l’iterazione tra varibil Lunghezza - Sesso e Cranio - Sesso.
mod4 = update(mod2, ~ . +Lunghezza*Sesso+Lunghezza*Sesso)
summary(mod4)##
## Call:
## lm(formula = dati$Peso ~ N.gravidanze + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso + Lunghezza:Sesso,
## data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1142.22 -180.65 -15.76 161.17 2557.15
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6510.2766 162.8600 -39.975 <2e-16 ***
## N.gravidanze 12.9155 4.3325 2.981 0.0029 **
## Gestazione 32.4244 3.7916 8.552 <2e-16 ***
## Lunghezza 9.8952 0.3524 28.081 <2e-16 ***
## Cranio 10.4738 0.4257 24.605 <2e-16 ***
## Tipo.parto -29.9728 12.0828 -2.481 0.0132 *
## Ospedale 14.6332 6.7453 2.169 0.0301 *
## Sesso -376.0459 212.7516 -1.768 0.0773 .
## Lunghezza:Sesso 0.9161 0.4288 2.136 0.0328 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared: 0.7287, Adjusted R-squared: 0.7278
## F-statistic: 836.2 on 8 and 2491 DF, p-value: < 2.2e-16Osservazioni:
Anche in questo caso l’iterazione tra varibili non porta vantaggi.
Il Modello 2 rimane il più semplice e performante.
Scelta del modello
Eseguiamo i test per la scelta del modello predittivo
AIC(mod1, mod2, mod3, mod4)## df AIC
## mod1 11 35174.86
## mod2 9 35172.59
## mod3 10 35169.63
## mod4 10 35170.01BIC(mod1, mod2, mod3, mod4)## df BIC
## mod1 11 35238.93
## mod2 9 35225.01
## mod3 10 35227.87
## mod4 10 35228.25RIepilogo modelli
| Modello | Adj R² | AIC | BIC |
|---|---|---|---|
| mod1 | 0.7274 | 35174.9 | 35238.9 |
| mod2 | 0.7274 | 35172.6 | 35225.0 |
| mod3 | 0.7278 | 35169.6 | 35227.9 |
| mod4 | 0.7277 | 35172.0 | 35236.0 |
Per AIC è preferibile il modello 3, per BIC il modello 2. BIC favorisce il modello più semplice.
Test ANOVA
Confrontiamo i modelli 2 e 3.
anova(mod2, mod3)## Analysis of Variance Table
##
## Model 1: dati$Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
## Model 2: dati$Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + I(Gestazione^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 187259272
## 2 2491 186888201 1 371072 4.946 0.02624 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Il p-value = 0.02624 è < 0.05, quindi l’aggiunta del termine Gestazione^2 migliora il modello rispetto a quello lineare. Questo suggerisce che l’effetto della gestazione sul peso del neonato non è perfettamente lineare, ma segue una leggera curvatura (quadratica).
Preferiamo comunque il modello più semplice (rasoio di Occam), modello 2.
Pacchetto MAAS
Proviamo ad elaborare il modello utilizzando direttamente la funzione stepAIC.
library(car)## Caricamento del pacchetto richiesto: carDatastepwise.mod = MASS::stepAIC(mod1,
direction = "both",
k = log(n))## Start: AIC=28136.41
## dati$Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 36143 187166141 28129
## - Fumatrici 1 92205 187222202 28130
## - Ospedale 1 329182 187459180 28133
## - Tipo.parto 1 456854 187586851 28135
## - N.gravidanze 1 472455 187602452 28135
## <none> 187129998 28136
## - Sesso 1 3648566 190778564 28177
## - Gestazione 1 5498236 192628234 28201
## - Cranio 1 45457162 232587160 28672
## - Lunghezza 1 87702588 274832585 29090
##
## Step: AIC=28129.07
## dati$Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 93132 187259272 28123
## - Ospedale 1 335696 187501837 28126
## - Tipo.parto 1 457853 187623993 28127
## <none> 187166141 28129
## - N.gravidanze 1 663688 187829828 28130
## + Anni.madre 1 36143 187129998 28136
## - Sesso 1 3655721 190821862 28170
## - Gestazione 1 5465000 192631141 28193
## - Cranio 1 45707341 232873481 28668
## - Lunghezza 1 87692908 274859048 29082
##
## Step: AIC=28122.49
## dati$Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Ospedale 1 342404 187601677 28119
## - Tipo.parto 1 450169 187709441 28121
## <none> 187259272 28123
## - N.gravidanze 1 640036 187899309 28123
## + Fumatrici 1 93132 187166141 28129
## + Anni.madre 1 37070 187222202 28130
## - Sesso 1 3639601 190898873 28163
## - Gestazione 1 5397695 192656967 28186
## - Cranio 1 45752230 233011503 28661
## - Lunghezza 1 88097837 275357110 29079
##
## Step: AIC=28119.23
## dati$Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 463870 188065546 28118
## <none> 187601677 28119
## - N.gravidanze 1 651066 188252743 28120
## + Ospedale 1 342404 187259272 28123
## + Fumatrici 1 99840 187501837 28126
## + Anni.madre 1 43769 187557908 28127
## - Sesso 1 3649259 191250936 28160
## - Gestazione 1 5444109 193045786 28183
## - Cranio 1 45758101 233359778 28657
## - Lunghezza 1 88054432 275656108 29074
##
## Step: AIC=28117.58
## dati$Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso
##
## Df Sum of Sq RSS AIC
## <none> 188065546 28118
## - N.gravidanze 1 623141 188688687 28118
## + Tipo.parto 1 463870 187601677 28119
## + Ospedale 1 356105 187709441 28121
## + Fumatrici 1 91892 187973654 28124
## + Anni.madre 1 44972 188020574 28125
## - Sesso 1 3655292 191720838 28158
## - Gestazione 1 5464853 193530399 28181
## - Cranio 1 46108583 234174130 28658
## - Lunghezza 1 87632762 275698308 29066Osservazioni:
Emerge che possiamo semplificare maggiormente il modello selezionato in precedenze, cioè il modello 2, escludendo anche le variabili TIpo.parto e Ospedale
# aggiornamento modello 2
mod2 = update(mod2, ~ . - Ospedale - Tipo.parto)
summary(mod2)##
## Call:
## lm(formula = dati$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 ***
## Sesso 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-16vif(mod2)## N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1.023475 1.669189 2.074689 1.624465 1.040054Osservazioni:
Tutti i VIF sono ben al di sotto del livello critico (di solito 5 o 10), quindi non c’è evidenza di collinearità preoccupante tra le variabili indipendenti.
Diagnostica sui residui
par(mfrow = c(2,2))
plot(mod2)
Osservazioni di leva
# leverage
lev = hatvalues(mod2)
plot(lev)
p = sum(lev)
soglia = 2 * p/n
abline(h = soglia, col ="red")
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.005775173Outliers
plot(rstudent(mod2))
abline(h=c(-2,2), col = "red")
outlierTest(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-03Distanza di Cook
cook = cooks.distance(mod2)
plot(cook)
max(cook)## [1] 0.8300965Osservazioni:
Sono state individuate numerose osservazioni sopra soglia
Osservazioni 1551, 155 e 1306 hanno residui significativamente fuori dall’intervallo [−2,2] e p-value Bonferroni corretti molto bassi. Potenzialmente influenti per la stima dei coefficienti.
Osservazione 1551 sotto la soglia critica ma comunque meritevole di attenzione. Alcune osservazioni possono avere impatto significativo sulla stima del modello.
Test residui
library(lmtest)## Caricamento del pacchetto richiesto: zoo##
## Caricamento pacchetto: 'zoo'## I seguenti oggetti sono mascherati da 'package:base':
##
## as.Date, as.Date.numericbptest(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-16plot(density(residuals(mod2)))
qqnorm(residuals(mod2))
qqline(residuals(mod2), col = "red")
Osservazioni:
Il modello viola i presupposti di omoscedasticità e normalità dei residui ma non quello dell’indipendenza. Questo potrebbe influenzare la validità dei test inferenziali.
La maggior parte dei punti segue bene la linea rossa, il che indica che i residui sono quasi normali. Le deviazioni visibili alle estremità(code) suggeriscono una leggera asimmetria o presenza di outlier, ma non grave.
Trattandosi di un campione che contiene un numero elevato di osservazioni, gli esiti negativi dei test non invalidano del tutto il modello di regressione.
Root Mean Squared Error
rmse = sqrt(mean(residuals(mod2)^2))
rmse## [1] 274.274Errore medio tra i valori osservati e i valori previsti di circa 274g
Previsioni
Peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.
# dati per la previsione
caso_prev = data.frame(
N.gravidanze = 3,
Sesso = 0,
Gestazione = 39,
Lunghezza = mean(dati$Lunghezza[dati$Sesso == 0]),
Cranio = mean(dati$Cranio[dati$Sesso == 0])
)
prev = predict(mod2, newdata = caso_prev, interval = "confidence")
prev## fit lwr upr
## 1 3195.331 3171.967 3218.696Stimato il peso per un neonato Femmina (F) a 39 settimane di gestazione, con valori medi di lunghezza e circonferenza cranica (rispetto al sesso). Il modello prevede un peso di 3190,7 g, con un intervallo di confidenza al 95% compreso tra 3167,0 g e 3214,4 g
Visualizzazioni risultati
Risultato Previsione
plot(density(dati$Peso[dati$Sesso == 0]), main ="Prev Peso Neonata, Sesso = F")
abline(v=mean(dati$Peso[dati$Sesso == 0]), col="red")
abline(v=prev[1], col="green")
legend("topleft", legend = c("Media Peso Campione", "Peso Stimato"),
col = c("red", "green"), lwd = 2)
mean(dati$Peso[dati$Sesso == 0])-prev[1]## [1] -34.19928mean(dati$Peso[dati$Sesso == 0])- P0_Peso## [1] -138.8678Rispetto alla media del campione il peso previsto è inferiore di 34g. Differenza di circa 139g rispetto alla media della popolazione.
Confronto Peso e Gestazione per Fumatrici
par(mfrow = c(1, 2))
boxplot(dati$Peso~dati$Fumatrici, main= "Peso ~ Fumo", xlab = "Fumo", ylab = "Peso")
boxplot(dati$Gestazione~dati$Fumatrici, main= "Gestazione ~ Fumo", xlab = "Fumo", ylab = "Gestazione")
I neonati di madri fumatrici hanno un peso medio più basso rispetto a quelli di madri non fumatrici. Il grafico supporta l’ipotesi che il fumo materno abbia un effetto negativo sul peso neonatale, ma non sulla durata della gestazione.
mean(dati$Peso[dati$Fumatrici == 0]) - mean(dati$Peso[dati$Fumatrici == 1])## [1] 49.8066Per le Fumatrici si registra un peso alla nascita inferiore di circa 50g.
Differeze misure antropometriche tra sessi
library(ggplot2)
Sesso_2 = ifelse(dati$Sesso == 1, "M", "F")
ggplot(data=dati)+
geom_point(aes(x= Lunghezza,
y= Peso,
col = Sesso_2), position = "jitter")+
geom_smooth(aes(x= Lunghezza,
y= Peso,
col = Sesso_2), se=F, method = "lm")+
labs(title = "Correlazione tra Peso e Lunghezza")## `geom_smooth()` using formula = 'y ~ x'
Per Lunghezza maggiore di 400 il peso dei Maschi supera il peso delle Femmine
ggplot(data=dati)+
geom_point(aes(x= Cranio,
y= Peso,
col = Sesso_2), position = "jitter")+
geom_smooth(aes(x= Cranio,
y= Peso,
col = Sesso_2), se=F, method = "lm")+
labs(title = "Correlazione tra Peso e Cranio")## `geom_smooth()` using formula = 'y ~ x'
In relazione alla circoferenza del cranio, il peso dei Maschi è sempre superiore a quello delle Femmine
ggplot(dati, aes(x = Lunghezza, y = Peso, size = Cranio, color = Sesso)) +
geom_point(alpha = 0.4)+
labs(title = "Correlazione tra Peso, Lunghezza e Cranio")
Nella parte più alta e a destra del grafico abbiamo una concentrazione maggiore di azzurro, a dimostrazione che nei Maschi Peso, Lunghezza e circonferenza del cranio registrano valori più alti.
Conclusioni:
L’analisi condotta ha evidenziato come il peso neonatale sia significativamente influenzato da lunghezza e circonferenza cranica alla nascita, durata della gestazione e sesso del neonato che infatti risultano essere predittori altamente significativi. Il numero di gravidanze precedenti mostra un impatto minore.
L’introduzione di un termine quadratico per la gestazione ha migliorato lievemente il modello, suggerendo una relazione non lineare tra la durata della gravidanza e il peso alla nascita.
L’interazione tra sesso e lunghezza ha confermato la presenza di differenze di crescita tra maschi e femmine.
I modelli sviluppati spiegano circa il 73% della variabilità del peso neonatale, indicando una buona capacità predittiva con un erroe medio di 274 g.