Modello Statistico per la Previsione del Peso Neonatale
- Raccolta dei Dati e Struttura del Dataset
dati <- read.csv("~/Desktop/texas/neonati.csv", stringsAsFactors = T)
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 Fattach(dati)summary(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 Ospedale Sesso
## Min. : 830 Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1256
## 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Nat:1772 osp2:849 M:1244
## Median :3300 Median :500.0 Median :340 osp3:835
## Mean :3284 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :565.0 Max. :390Alcuni dati sembrano errati: in Anni.madre il valore minimo è 0, il che non ha senso, probabilmente si tratta di un errore di registrazione. In N.gravidanze ci sono casi con 12 gravidanze, insolito ma non impossibile. Per correggere i dati, sostituisco i casi in cui Anni.madre = 0 con la mediana.
dati$Anni.madre[dati$Anni.madre == 0] <- median(dati$Anni.madre[dati$Anni.madre > 0])
summary(dati)## Anni.madre N.gravidanze Fumatrici Gestazione
## Min. : 1.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.18 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 Ospedale Sesso
## Min. : 830 Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1256
## 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Nat:1772 osp2:849 M:1244
## Median :3300 Median :500.0 Median :340 osp3:835
## Mean :3284 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :565.0 Max. :390Il valore minimo di Anni.madre ora è 1. Applico nuovamente la correzione.
dati$Anni.madre[dati$Anni.madre == 1] <- median(dati$Anni.madre[dati$Anni.madre > 1])
summary(dati) ## Anni.madre N.gravidanze Fumatrici Gestazione
## Min. :13.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.19 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 Ospedale Sesso
## Min. : 830 Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1256
## 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Nat:1772 osp2:849 M:1244
## Median :3300 Median :500.0 Median :340 osp3:835
## Mean :3284 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :565.0 Max. :390Ora i dati sembrano corretti.
Calcoliamo asimmetria e curtosi per vedere la distribuzione della variabile Peso.
skew <- moments::skewness(Peso)
kurt <- moments::kurtosis(Peso)-3
cat("L'asimmetria è di", round(skew, 3),".\n")## L'asimmetria è di -0.647 .cat("La curtosi è di", round(kurt, 3))## La curtosi è di 2.032plot(density(Peso))
Il valore di -0.647 suggerisce che la distribuzione del peso dei neonati è moderatamente asimmetrica a sinistra. Mentre il valore di 2.032 indica che la distribuzione è leptocurtica, ovvero ha una coda più lunga e appuntita rispetto a una distribuzione normale.
- Analisi e Modellizzazione
Analisi Preliminare
#controlliamo che la variabile risposta(Peso) sia normale
n <- nrow(dati)
shapiro.test(Peso) ##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16Si rifiuta l’ipotesi di normalità(p-value < 0.05)
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.2)
}
#correlazioni
pairs(dati,lower.panel=panel.cor, upper.panel=panel.smooth)
Le variabili più correlate con la variabile risposta sono Lunghezza, Cranio e Gestazione. Si nota, inoltre, una correlazione positiva tra Gestazione e Lunghezza e Cranio e Lunghezza.
par(mfrow=c(1,2))
boxplot(Peso~Fumatrici)
boxplot(Peso~Gestazione) 
Si può notare che i valori medi di Peso tendono ad essere più bassi se la madre è fumatrice. Invece con l’aumentare delle settimane di Gestazione il Peso aumenta.
# 1. In alcuni ospedali si fanno più parti cesarei
table_parto <- table(Ospedale, Tipo.parto)
table_parto## Tipo.parto
## Ospedale Ces Nat
## osp1 242 574
## osp2 254 595
## osp3 232 603chisq.test(table_parto) ##
## Pearson's Chi-squared test
##
## data: table_parto
## X-squared = 1.0972, df = 2, p-value = 0.5778Non si rifiuta l’ipotesi nulla(p-value > 0.05): il risultato non è statisticamente significativo, indicando che il tipo di parto non è associato all’ospedale.
# 2. La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione
# Peso e lunghezza media di un neonato
media_popolazione_peso <- 3300
media_popolazione_lunghezza <- 500
t.test(Peso, mu=media_popolazione_peso)##
## One Sample t-test
##
## data: 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.081t.test(Lunghezza, mu=media_popolazione_lunghezza)##
## 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.692Il test t ha mostrato che la media del peso dei neonati nel campione non differisce significativamente dalla media della popolazione (p-value > 0.05), suggerendo che la media campionaria può essere considerata rappresentativa della popolazione. Invece, la media della lunghezza risulta significativamente inferiore rispetto a quella attesa nella popolazione (p-value < 0.05), quindi c’è una differenza statisticamente significativa tra la media campionaria e quella della popolazione.
# 3. Le misure antropometriche sono significativamente diverse tra i due sessi
t.test(Peso~Sesso)##
## 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.215t.test(Lunghezza~Sesso)##
## 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.6672t.test(Cranio~Sesso)##
## 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.4486In tutti e tre i test, il p-value è estremamente basso, ben al di sotto del valore di soglia di 0.05. Questo significa che la differenza tra i due sessi è statisticamente significativa, come si può osservare anche nel grafico sottostante.
boxplot(Peso~Sesso)
Creazione del Modello di Regressione
# 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.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 *
## Fumatrici -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-16Analisi dei coefficienti calcolati: L’età della madre non ha un effetto significativo sul peso (p-value > 0.05), quindi questa variabile non contribuisce significativamente al modello. Il numero di gravidanze ha un effetto positivo e significativo sul peso. Ogni gravidanza in più aumenta il peso di circa 11.41 grammi. Il fatto che la madre sia fumatrice non sembra avere un effetto significativo sul peso, poiché il p-value è maggiore di 0.05, tuttavia sarebbe opportuno considerare lo stesso questa variabile nel modello, dato che fumare durante la gravidanza può influenzare la salute del neonato. La durata della gestazione ha un effetto positivo e molto significativo sul peso. Ogni settimana aggiuntiva di gestazione comporta un aumento del peso di circa 32.53 grammi. La lunghezza ha un effetto altamente significativo sul peso. Ogni unità aggiuntiva di lunghezza comporta un aumento di circa 10.30 grammi. La dimensione del cranio ha anch’essa un effetto significativo sul peso, con un aumento di 10.47 grammi per ogni unità aggiuntiva di cranio. Il tipo di parto naturale ha un effetto positivo e significativo sul peso, con un aumento medio di 29.50 grammi rispetto al parto cesareo. Il tipo di ospedale non influenza significativamente il peso del neonato. Il sesso maschile ha un effetto significativo e positivo sul peso. I neonati maschi tendono ad avere un peso medio superiore di 77.55 grammi rispetto alle femmine. Questo effetto è altamente significativo, con un p-value estremamente basso.
L’R-quadro aggiustato pari a 0.7278 indica che il modello spiega abbastanza bene i dati, anche se ci sono alcune variabili non significative.
#Proviamo a togliere alcune variabili non significative
mod2 <- update(mod1, ~.-Anni.madre)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso, data = dati)
##
## 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 **
## Fumatrici -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-16La rimozione della variabile “Anni.madre” non ha influito significativamente sulla qualità del modello, poiché l’R-quadro e le statistiche globali sono rimasti pressoché invariati, e come si può osservare dai risultati dei seguenti test(Anova e BIC).
anova(mod2,mod1)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
## Model 2: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2490 186809099
## 2 2489 186772707 1 36392 0.485 0.4862Il p-value è maggiore di 0.05, il che indica che la rimozione della variabile “Anni.madre” non ha un impatto significativo sulla qualità del modello e non causa una perdita significativa di informazione per la previsione del peso.
BIC(mod2,mod1)## df BIC
## mod2 11 35234.64
## mod1 12 35241.97Il modello mod2 ha un BIC più basso rispetto al modello mod1. Un valore di BIC più basso indica che mod2 è preferibile in termini di equilibrio tra la bontà di adattamento e la complessità del modello.
car::vif(mod2)## GVIF Df GVIF^(1/(2*Df))
## N.gravidanze 1.027985 1 1.013896
## Fumatrici 1.007346 1 1.003666
## Gestazione 1.676688 1 1.294870
## Lunghezza 2.085755 1 1.444214
## Cranio 1.626661 1 1.275406
## Tipo.parto 1.004240 1 1.002118
## Ospedale 1.003421 2 1.000854
## Sesso 1.040558 1 1.020077Tutte le variabili hanno un valore minore di 5, questo indica che non c’è pericolo di multicollinearità.
# Procediamo con la rimozione di variabili non significative
mod3 <- update(mod2, ~.-Ospedale)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Sesso, data = dati)
##
## 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 **
## Fumatrici -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, ~.-Tipo.parto)
summary(mod4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1150.3 -181.3 -15.7 163.0 2636.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.6714 135.7178 -49.232 < 2e-16 ***
## N.gravidanze 12.7185 4.3450 2.927 0.00345 **
## Fumatrici -30.4634 27.5948 -1.104 0.26972
## Gestazione 32.5914 3.8051 8.565 < 2e-16 ***
## Lunghezza 10.2341 0.3009 34.011 < 2e-16 ***
## Cranio 10.5359 0.4262 24.718 < 2e-16 ***
## SessoM 78.1713 11.2028 6.978 3.83e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2493 degrees of freedom
## Multiple R-squared: 0.7271, Adjusted R-squared: 0.7265
## F-statistic: 1107 on 6 and 2493 DF, p-value: < 2.2e-16La rimozione delle variabili Ospedale e Tipo.Parto non causa una perdita significativa di informazione per la previsione del peso dei neonati.
mod5 <- update(mod4, ~.-Fumatrici)
summary(mod5)##
## 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-16Selezione del Modello Ottimale
BIC(mod5,mod4,mod3,mod2,mod1)## df BIC
## mod5 7 35220.10
## mod4 8 35226.70
## mod3 9 35228.24
## mod2 11 35234.64
## mod1 12 35241.97Sebbene il BIC è piu basso nel mod5, terrei in considerazione la variabile Fumatrici nel modello perchè, anche se non risulta significativa in questo campione di dati, potrebbe portare a dei risultati considerevoli.
Nei seguenti passaggi, proviamo a inserire l’effetto quadratico della variabile Gestazione, osservato anche dalla curvatura nel grafico delle correlazioni visto in precedenza.
plot(Gestazione, Peso)
mod_q <- lm(Peso ~ N.gravidanze + Fumatrici + Gestazione + I(Gestazione^2) +
Lunghezza + Cranio + Sesso, data = dati)
summary(mod_q)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + I(Gestazione^2) +
## Lunghezza + Cranio + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1145.00 -180.96 -13.12 165.00 2659.10
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4669.1027 898.3246 -5.198 2.18e-07 ***
## N.gravidanze 12.7990 4.3416 2.948 0.00323 **
## Fumatrici -29.2442 27.5772 -1.060 0.28904
## Gestazione -79.7741 49.7260 -1.604 0.10878
## I(Gestazione^2) 1.4999 0.6618 2.266 0.02352 *
## Lunghezza 10.3386 0.3042 33.989 < 2e-16 ***
## Cranio 10.6312 0.4280 24.841 < 2e-16 ***
## SessoM 75.9814 11.2351 6.763 1.68e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.4 on 2492 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.7269
## F-statistic: 951.4 on 7 and 2492 DF, p-value: < 2.2e-16L’effetto lineare della variabile Gestazione è diventato non significativo, mentre l’effetto quadratico risulta poco significativo. Non consideriamo queste variabili nel modello.
Proviamo ad inserire alcuni effetti di interazione tra le variabili.
mod_interazioni <- lm(Peso ~ N.gravidanze + Fumatrici + Gestazione +
Lunghezza + Cranio + Sesso + Gestazione * Fumatrici + Lunghezza * Sesso, data=dati)
summary(mod_interazioni)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso + Gestazione * Fumatrici + Lunghezza * Sesso,
## data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1162.40 -179.33 -16.22 163.06 2571.11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6513.1362 163.9000 -39.738 < 2e-16 ***
## N.gravidanze 13.0268 4.3442 2.999 0.00274 **
## Fumatrici 742.3677 757.2148 0.980 0.32699
## Gestazione 33.4021 3.8385 8.702 < 2e-16 ***
## Lunghezza 9.8486 0.3532 27.885 < 2e-16 ***
## Cranio 10.5011 0.4262 24.637 < 2e-16 ***
## SessoM -358.1206 213.2511 -1.679 0.09321 .
## Fumatrici:Gestazione -19.7174 19.2745 -1.023 0.30642
## Lunghezza:SessoM 0.8818 0.4298 2.051 0.04032 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.4 on 2491 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.7269
## F-statistic: 832.2 on 8 and 2491 DF, p-value: < 2.2e-16Gli effetti tra variabili non sono significativi, anzi annullano anche un predittore importante per il peso, ossia la variabile Sesso.
Proseguiamo con la selezione automatica del modello tramite la procedura stepwise basata sul criterio di informazione di Akaike(AIC).
stepwise.mod <- MASS::stepAIC(mod1,
direction= "both",
k=log(n))## Start: AIC=28139.46
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 36392 186809099 28132
## - Fumatrici 1 89979 186862686 28133
## - Ospedale 2 686397 187459103 28133
## - Tipo.parto 1 447233 187219939 28138
## - N.gravidanze 1 448762 187221469 28138
## <none> 186772707 28140
## - Sesso 1 3611588 190384294 28180
## - Gestazione 1 5446558 192219264 28204
## - Cranio 1 45338051 232110758 28675
## - Lunghezza 1 87959790 274732497 29096
##
## Step: AIC=28132.12
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 90897 186899996 28126
## - Ospedale 2 692738 187501837 28126
## - Tipo.parto 1 448222 187257321 28130
## <none> 186809099 28132
## - N.gravidanze 1 633756 187442855 28133
## + Anni.madre 1 36392 186772707 28140
## - Sesso 1 3618736 190427835 28172
## - Gestazione 1 5412879 192221978 28196
## - Cranio 1 45588236 232397335 28670
## - Lunghezza 1 87950050 274759149 29089
##
## Step: AIC=28125.51
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 701680 187601677 28119
## - Tipo.parto 1 440684 187340680 28124
## <none> 186899996 28126
## - N.gravidanze 1 610840 187510837 28126
## + Fumatrici 1 90897 186809099 28132
## + Anni.madre 1 37311 186862686 28133
## - Sesso 1 3602797 190502794 28165
## - Gestazione 1 5346781 192246777 28188
## - Cranio 1 45632149 232532146 28664
## - Lunghezza 1 88355030 275255027 29086
##
## Step: AIC=28119.23
## 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 2 701680 186899996 28126
## + Fumatrici 1 99840 187501837 28126
## + Anni.madre 1 43845 187557831 28126
## - Sesso 1 3649259 191250936 28160
## - Gestazione 1 5444109 193045786 28183
## - Cranio 1 45758101 233359778 28657
## - Lunghezza 1 88054432 275656108 29074
##
## 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
## + Tipo.parto 1 463870 187601677 28119
## + Ospedale 2 724866 187340680 28124
## + Fumatrici 1 91892 187973654 28124
## + Anni.madre 1 45044 188020502 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-16Anche se dalla procedura stepwise il modello ottimale risulta essere il mod5, ritengo più opportuno scegliere il mod4 con la variabile Fumatrici, con un R^2=0.73 circa.
Analisi della Qualità del Modello
par(mfrow=c(2,2))
plot(mod4)
Nel grafico dei residui si può osservare un pattern leggermente ricurvo e le osservazioni sono concentrate nella parte destra del grafico, tranne l’osservazione 1551 che risulta più spostata. Nel secondo grafico i punti sono disposti correttamente sopra la retta, indicando una distribuzione normale dei residui, a parte per qualche osservazione nella coda superiore. Nel terzo grafico si nota, come nel primo, una leggera curvatura. Nel quarto grafico si visualizzano i potenziali valori influenti sulle stime di regressione: l’osservazione 1551 supera la soglia di avvertimento ed è molto vicina alla soglia di allarme.
#leverage
lev <- hatvalues(mod4)
plot(lev)
p = sum(lev)
soglia = 2*p/n
abline(h=soglia,col=2)
Si osservano parecchi punti con un valore di leverage maggiore della soglia: sono osservazioni che potrebbero avere un’influenza eccessiva sul modello ed essere degli outliers in termini di variabili indipendenti, anche se non necessariamente in termini di variabile dipendente.
#outliers
library(car)
plot(rstudent(mod4))
abline(h=c(-2,2),col=2)
outlierTest(mod4)## rstudent unadjusted p-value Bonferroni p
## 1551 10.039719 2.8060e-23 7.0149e-20
## 155 5.022108 5.4723e-07 1.3681e-03
## 1306 4.823102 1.4986e-06 3.7465e-03Sono presenti diverse osservazioni outliers sia sopra che sotto le soglie: ne spiccano 3 in particolare dopo aver fatto il test di Bonferroni applicato ai residui studentizzati(outlierTest).
#distanza di Cook
cook <- cooks.distance(mod4)
plot(cook)
Si può concludere che l’osservazione 1551 sembra particolaremente influente sulle stime di regressione.
dati[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 Flibrary(lmtest)
bptest(mod4)##
## studentized Breusch-Pagan test
##
## data: mod4
## BP = 89.798, df = 6, p-value < 2.2e-16dwtest(mod4)##
## Durbin-Watson test
##
## data: mod4
## DW = 1.9542, p-value = 0.126
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(residuals(mod4))##
## Shapiro-Wilk normality test
##
## data: residuals(mod4)
## W = 0.9741, p-value < 2.2e-16plot(density(residuals(mod4)))
BP test: p-value < 0.05, rifiuto l’ipotesi nulla(varianza non costante). DW test:p-value > 0.05, i residui non sono autocorrelati. Shapiro test: p-value < 0.05: rifiuto l’ipotesi nulla, i residui non sono normali.
- Previsioni e Risultati
Peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.
media_lunghezzaF <- mean(Lunghezza[Sesso=="F"])
media_CranioF <- mean(Cranio[Sesso=="F"])
#madre non fumatrice
peso_predNF <- predict(mod4, data.frame(N.gravidanze=3, Gestazione=39,
Sesso="F",Lunghezza=media_lunghezzaF,
Cranio=media_CranioF, Fumatrici=0))
cat("Il peso previsto è di", round(peso_predNF, 3), "grammi.")## Il peso previsto è di 3197.088 grammi.# madre fumatrice
peso_predF <- predict(mod4, data.frame(N.gravidanze=3, Gestazione=39,
Sesso="F",Lunghezza=media_lunghezzaF,
Cranio=media_CranioF, Fumatrici=1))
cat("Il peso previsto è di", round(peso_predF, 3), "grammi.")## Il peso previsto è di 3166.625 grammi.Il peso della neonata risulta minore se la madre è fumatrice.
- Visualizzazioni
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")
Come si è osservato dallo studio del modello ottimale, e come si può osservare dal grafico sopra, il peso dei neonati maschi è più elevato rispetto a quello delle neonate. Entrambe le rette, però, hanno la stessa pendenza e il valore del peso aumenta all’aumentare delle settimane di gestazione. Le osservazioni si concentrano tra la 37esima e la 42esima settimana.
ggplot(data=dati) +
geom_point(aes(x=Gestazione,
y=Peso,
col=factor(Fumatrici)),position = "jitter")+
geom_smooth(aes(x=Gestazione,
y=Peso,
col=factor(Fumatrici)),se=F,method = "lm")+
labs(x = "Settimane di Gestazione",
y = "Peso",
color = "Fumatrici") 
Il peso dei neonati tende a diminuire se le madri sono fumatrici: infatti la retta che le rappresenta tende a inclinarsi di pù rispetto a quella delle non-fumatrici.