Modello Statistico per la Previsione del Peso Neonatale
Emanuela Spagnolo
2024-10-23
1. Raccolta dei Dati e Struttura del Dataset
dati <-read.csv("neonati.csv", stringsAsFactors = T)
attach(dati)
knitr::kable((head(dati,5)), "pipe")| 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 |
Il dataset è costituito da 2500 osservazioni e 10 variabili. Esse sono così suddivise:
- Quantitative discrete: N.gravidanze.
- Quantitative continue: Peso, Lunghezza, Cranio, Anni.madre, Gestazione.
- Qualitative su scala nominale: Tipo.parto, Ospedale, Sesso, Fumatrici.
L’obiettivo dello studio è capire se è possibile prevedere il peso del neonato alla nascita basandosi su variabili cliniche raccolte da tre ospedali. Pertanto la variabile Peso è la nostra variabile risposta, tutte le altre sono variabili esplicative.
2. Analisi e Modellizzazione
Analisi preliminare
Quantitative discrete
Sulle variabili quantitative discrete calcoliamo gli indici di posizione e di forma
Analizziamo la variabile N.gravidanze, calcoliamone la tabella delle frequenze assolute, gli indici di posizione e di forma.
summary(N.gravidanze)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.9812 1.0000 12.0000table(N.gravidanze)## N.gravidanze
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1096 818 340 150 48 21 11 1 8 2 3 1 1mydata <- data.frame(Asimmetria=round(skewness(N.gravidanze),2), Curtosi=round(kurtosis(N.gravidanze)-3,2), Dev_standard=round(sd(N.gravidanze),2))
knitr::kable((mydata), "pipe")| Asimmetria | Curtosi | Dev_standard |
|---|---|---|
| 2.51 | 10.99 | 1.28 |
Dalla tabella delle frequenze assolute la maggior parte delle donne ha avuto pochi figli, per lo più nessuno. La media è infatti dello 0.9812, la mediana è pari a 1, la moda è 0. L’asimmetria della distribuzione della variabile è positiva, pari a 2.51, sono quindi più frequenti valori bassi. La curtosi è pari a 10.99, la distribuzione è quindi leptocurtica, cioè risulta più allungata rispetto alla distribuzione normale.
Quantitative continue
Analizziamo le variabili quantitative continue: Peso, Lunghezza, Cranio, Anni.madre e Gestazione. Calcoliamone gli indici di posizione e forma.
Cominciamo dalla variabile Peso.
summary(Peso)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 2990 3300 3284 3620 4930mydata <- data.frame(Asimmetria=round(skewness(Peso),2), Curtosi=round(kurtosis(Peso)-3,2), Dev_standard=round(sd(Peso),2))
knitr::kable((mydata), "pipe")| Asimmetria | Curtosi | Dev_standard |
|---|---|---|
| -0.65 | 2.03 | 525.04 |
Dal summary risulta che la variabile Peso va da un minimo di 830 e un massimo di 4930. La media è pari a 3284, la mediana è pari a 3300.
L’asimmetria è pari a -0.65, pertanto la distribuzione è leggermente asimmetrica negativamente, ossia sono più frequenti valori alti. la curtosi è pari a 2.03, pertanto risulta più allungata rispetto alla distribuzione normale.
normale<-rnorm(100000,3284,525.0387)
plot(density(Peso))+
lines(density(normale), col=2,lwd=3)
## integer(0)Analizziamo ora la variabile Lunghezza. Essa esprime in millimetri la lunghezza del neonato, misurabile anche durante la gravidanza tramite ecografie.
summary(Lunghezza)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 310.0 480.0 500.0 494.7 510.0 565.0mydata <- data.frame(Asimmetria=round(skewness(Lunghezza),2), Curtosi=round(kurtosis(Lunghezza)-3,2), Dev_standard=round(sd(Lunghezza),2))
knitr::kable((mydata), "pipe")| Asimmetria | Curtosi | Dev_standard |
|---|---|---|
| -1.51 | 6.49 | 26.32 |
La variabile lunghezza va da un minimo di 310 a un massimo di 565 mm. La media è pari a 494.7, la mediana è pari a 500.
L’asimmetria è pari a -1.51, ossia sono leggermente più frequenti valori alti. La curtosi è pari a 6.49, ossia la sua distribuzione è leptocurtica, come si vede anche dal grafico sottostante.
normale<-rnorm(100000,494.7,26.31864)
plot(density(Lunghezza))+
lines(density(normale), col=2,lwd=3)
## integer(0)Aanalizziamo la variabile Cranio. Essa rappresenta, in millimetri, la lunghezza del diametro craniale, misurabile anche durante la gravidanza tramite ecografie.
summary(Cranio)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 235 330 340 340 350 390mydata <- data.frame(Asimmetria=round(skewness(Cranio),2), Curtosi=round(kurtosis(Cranio)-3,2), Dev_standard=round(sd(Cranio),2))
knitr::kable((mydata), "pipe")| Asimmetria | Curtosi | Dev_standard |
|---|---|---|
| -0.79 | 2.95 | 16.43 |
La variabile Cranio va da un minimo di 235 a un massimo di 390 mm. La media è pari a 340, la mediana è pari a 340.
L’asimmetria è pari a -0.79, ossia sono leggermente più frequenti valori alti. La curtosi è pari a 2.95, ossia la sua distribuzione è leptocurtica.
normale<-rnorm(100000,340,16.42533)
plot(density(Cranio))+
lines(density(normale), col=2,lwd=3)
## integer(0)summary(Anni.madre)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 25.00 28.00 28.16 32.00 46.00Dai valori individuati dal summary di Anni.madre si nota un valore anomalo per il minimo (0). Vediamone la tabella delle frequenze assolute, media, mediana e moda.
table(Anni.madre)## Anni.madre
## 0 1 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
## 1 1 1 2 6 13 18 24 45 66 74 100 115 131 180 184 197 172 174 200
## 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
## 147 159 110 96 66 64 41 38 27 19 13 8 2 4 1 1mydata <- data.frame(Asimmetria=round(skewness(Anni.madre),2), Curtosi=round(kurtosis(Anni.madre)-3,2), Dev_standard=round(sd(Anni.madre),2))
knitr::kable((mydata), "pipe")| Asimmetria | Curtosi | Dev_standard |
|---|---|---|
| 0.04 | 0.38 | 5.27 |
Dalla tabella delle frequenze assolute anche il valore 1 non può essere un valore reale.
La media è pari a 28.164, la mediana è pari a 28, la moda è pari a 30. L’asimmetria è pari a 0.04 (positiva ma molto vicina allo 0), la curtosi vale 0.38 (positiva ma molto vicina allo 0).
normale<-rnorm(100000,28.164,5.273578)
plot(density(Anni.madre))+
lines(density(normale), col=2,lwd=3)
## integer(0)Dal grafico risulta che la densità della variabile Anni.madre ha una forma molto simile a quella della distribuzione normale con media e deviazione standard pari a quelle di Anni.madre.
Analizziamo la variabile Gestazione.
summary(Gestazione)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 25.00 38.00 39.00 38.98 40.00 43.00table(Gestazione)## Gestazione
## 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
## 1 1 2 4 3 5 8 9 18 16 33 62 192 437 581 741 329 56 2mydata <- data.frame(Asimmetria=round(skewness(Gestazione),2), Curtosi=round(kurtosis(Gestazione)-3,2), Dev_standard=round(sd(Gestazione),2))
knitr::kable((mydata), "pipe")| Asimmetria | Curtosi | Dev_standard |
|---|---|---|
| -2.07 | 8.26 | 1.87 |
Le settimane di gestazione variano da un minimo di 25 a un massimo di 43 settimane. La media è di 38.9804, la mediana è pari a 39, la moda è pari a 40.
Per quanto riguarda l’asimmetria, la distribuzione associata alla variabile Gestazione ha un’asimmetria negativa (-2.07), pertanto sono più frequenti valori alti. La curtosi è pari a 8.26, la distribuzione è quindi leptocurtica.
Quantitative su scala nominale
Vediamo quanti valori diversi hanno le variabili qualitative su scala nominale:
mydata <- data.frame(Tipo.parto=levels(Tipo.parto))
knitr::kable((mydata), "pipe")| Tipo.parto |
|---|
| Ces |
| Nat |
Tipo.parto ha due livelli: “Ces” indica il parto cesareo, “Nat” indica il parto naturale.
table(Tipo.parto)## Tipo.parto
## Ces Nat
## 728 1772Dalla tabella delle frequenze assolute si vede che c’è una netta maggioranza di parti naturali sui parti cesarei.
mydata <- data.frame(Ospedale=levels(Ospedale))
knitr::kable((mydata), "pipe")| Ospedale |
|---|
| osp1 |
| osp2 |
| osp3 |
Ospedale ha tre livelli, ognuno indica l’ospedale di riferimento.
table(Ospedale)## Ospedale
## osp1 osp2 osp3
## 816 849 835Le nascite sono più o meno distribuite in maniera equa nei tre ospedali oggetto dello studio, come mostra la tabella delle frequenze assolute.
mydata <- data.frame(Sesso=levels(Sesso))
knitr::kable((mydata), "pipe")| Sesso |
|---|
| F |
| M |
Sesso ha due livelli, uno riferito al genere femminile, l’altro riferito al genere maschile.
table(Sesso)## Sesso
## F M
## 1256 1244Dalla tabella delle frequenze assolute si nota che è nata qualche bambina in più rispetto ai bambini.
mydata <- data.frame(Fumatrici=levels(Fumatrici))
knitr::kable((mydata), "pipe")Fumatrici non è stata individuata come variabile categorica. Creiamone una variabile dummy tenenendo presente che 0 indica che la madre è non fumatrice, 1 indica che la madre è fumatrice.
dati$Fumatrici_d_f<-factor(ifelse(Fumatrici==1,'S','N'))
mydata <- data.frame(Fumatrici=levels(dati$Fumatrici_d_f))
knitr::kable((mydata), "pipe")| Fumatrici |
|---|
| N |
| S |
Vediamo quante sono le gestanti fumatrici e quante invece non lo sono.
table(dati$Fumatrici_d_f)##
## N S
## 2396 104La quasi totalità delle gestanti non è fumatrice, solo 104 su 2500 donne sono fumatrici.
Focus sull’impatto dell’età della gestante sul peso alla nascita
Prima di procedere con la creazione del modello, creo una nuova variabile per sostituire i valori 0 e 1 di Anni.madre con il valore medio di Anni.madre, che sono risultati essere non coerenti con il resto delle informazioni.
mean_annimadre=round(mean(dati$Anni.madre),0)
dati$Anni.madre <- ifelse(dati$Anni.madre == 0, mean_annimadre,dati$Anni.madre)
dati$Anni.madre <- ifelse(dati$Anni.madre==1,mean_annimadre,dati$Anni.madre)
table(dati$Anni.madre)##
## 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
## 1 2 6 13 18 24 45 66 74 100 115 131 180 184 197 174 174 200 147 159
## 33 34 35 36 37 38 39 40 41 42 43 44 45 46
## 110 96 66 64 41 38 27 19 13 8 2 4 1 1Dalla tabella delle frequenze assolute notiamo che i valori anomali (0 e 1) sono stati sostituiti con il valore medio (28).
plot(dati$Anni.madre, Peso, pch=20, xlab='Anni della madre')
Dal grafico che mette in relazione l’età della madre con il peso del neonato non si vede una netta relazione.
Focus sull’impatto del fumo materno sul peso alla nascita
ggplot(data=dati)+
geom_boxplot(aes(x=Fumatrici_d_f, y=Peso), fill='lightblue')+
theme_classic()+
labs(x='Fumatrice', y='Peso alla nascita')
Dal grafico dei boxplot, che mette il relazione il peso del neonato con la condizione di fumatrice della gestante, non si nota una netta relazione sul peso del bambino, il peso mediano sembra leggermente maggiore nel caso di non fumatrice.
Focus sulle gravidanze precedenti sul peso alla nascita
Vediamo adesso il grafico che mette in relazione il peso con le gravidanze precedenti.
plot(N.gravidanze, Peso, pch=20,xlab='Numero di gravidanze')
Dallo scatterplot non si nota una netta relazione fra il numero delle gravidanze e il peso del neonato.
Creazione del modello di regressione
Visualizziamo la matrice di correlazione. Elimino per il momento le variabili qualitative su scala nominale che verranno analizzate in maniera diversa.
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)
}
dati_select<-dati %>% select(Peso,Anni.madre, N.gravidanze, Gestazione, Lunghezza, Cranio)
pairs(dati_select,lower.panel=panel.cor, upper.panel=panel.smooth)## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
## Warning in par(usr): l'argomento 1 non nomina un parametro grafico
Dalla matrice delle correlazioni risulta che le settimane di gestazione e il peso sono correlate positivamente, questo implica che all’aumentare di una settimana di gestazione il peso aumenta. In particolare ad un aumento di una settimana di gestazione il peso aumenta, in media, di 0.59.
Il peso è anche molto correlato con la lunghezza e il diametro del cranio. Infatti un aumento unitario sulla lunghezza genera un aumento medio sul peso di 0.80. Infine un aumento unitario sul cranio genera un aumento medio sul peso di 0.70.
Gestazione e Lunghezza però sono molto correlate fra loro (0.62), come anche Cranio e Lunghezza (0.60). Questo potrebbe portare a problemi di multicollinearità.
Con le variabili qualitative su scala nominale lo scatterplot e il calcolo delle correlazioni non danno informazioni, si devono usare i boxplot.
Analizziamo ora le gestanti fumatrici.
par(mfrow=c(1,2))
boxplot(Peso)
boxplot(Peso~dati$Fumatrici_d_f, xlab='Fumatrici')
t.test(Peso~dati$Fumatrici_d_f)##
## Welch Two Sample t-test
##
## data: Peso by dati$Fumatrici_d_f
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group N and group S is not equal to 0
## 95 percent confidence interval:
## -45.61354 145.22674
## sample estimates:
## mean in group N mean in group S
## 3286.153 3236.346mydata <- data.frame(
Media.peso.non.fum=round(mean(Peso[dati$Fumatrici_d_f=='N']) ,2),
Media.peso.fum=round(mean(Peso[dati$Fumatrici_d_f=='S']),2)
)
knitr::kable((mydata), "pipe")| Media.peso.non.fum | Media.peso.fum |
|---|---|
| 3286.15 | 3236.35 |
Dal boxplot sembra che il fumo abbia un impatto negativo sul peso del neonato, dato che tutta la scatola delle gestanti fumatrici è più bassa rispetto a quelle non fumatrici, ma il p-value del test t per saggiare l’uguaglianza fra medie per gruppi indipendenti è troppo alto per rifiutare l’ipotesi nulla.
Probabilmente non si riesce ad avere dei risultati particolarmente rilevanti sull’impatto che una gestante fumatrice ha sul peso del bambino poiché nel campione assegnato il numero delle fumatrici è pari a circa il 4% del totale (104/2500). Nella realtà il numero delle gestanti fumatrici è nettamente maggiore nella popolazione, nell’articolo https://www.epicentro.iss.it/fumo/WTD2013Gravidanza si riporta il 23%.
Verifichiamo il sesso.
par(mfrow=c(1,2))
boxplot(Peso)
boxplot(Peso~Sesso)
mydata <- data.frame(
Media.peso.M=round(mean(Peso[Sesso=='M']) ,2),
Media.peso.S=round(mean(Peso[Sesso=='F']),2)
)
knitr::kable((mydata), "pipe")| Media.peso.M | Media.peso.S |
|---|---|
| 3408.22 | 3161.13 |
Si nota che i maschi pesano di più delle femmine. Saggiamo ora l’uguaglianza fra medie per gruppi indipendenti. Per farlo usiamo un t-test.
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.215Il p-value è molto piccolo, con un p-value molto piccolo si rifiuta l’ipotesi nulla, quindi le medie sono diverse.
Le variabili Tipo.parto e Ospedale possono essere escluse a priori dalle analisi poiché il peso del bambino sicuramente non dipende dall’ospedale in cui nasce e dal tipo di parto con cui nasce.
Creazione del modello di regressione
Creiamo un modello che tenga conto di tutte le variabili (per le gestanti fumatrici usiamo la variabile dummy creata).
mod1<-lm(Peso~ Anni.madre + N.gravidanze + Fumatrici_d_f + Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso , data=dati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici_d_f +
## Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale +
## Sesso, 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_d_fS -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-16car::vif(mod1)## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.190207 1 1.090966
## N.gravidanze 1.189252 1 1.090528
## Fumatrici_d_f 1.007410 1 1.003698
## Gestazione 1.695001 1 1.301922
## Lunghezza 2.085773 1 1.444220
## Cranio 1.630914 1 1.277073
## Tipo.parto 1.004255 1 1.002125
## Ospedale 1.004235 2 1.001057
## Sesso 1.040650 1 1.020123I coefficienti individuati dal modello lineare sono Gestazione, Lunghezza, Cranio e Sesso. Non si individua al momento multicollinearità. Il valore di \(R^2\) aggiustato è pari a 0.7278, cerchiamo di migliorare il modello.
Selezione del Modello Ottimale
Proviamo a togliere l’ospedale e il tipo di parto.
mod2<-lm(Peso~ Anni.madre + N.gravidanze + Fumatrici_d_f + Gestazione + Lunghezza + Cranio + Sesso , data=dati)
summary(mod2)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici_d_f +
## Gestazione + Lunghezza + Cranio + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1160.62 -181.17 -15.91 163.47 2631.35
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6711.5440 141.2543 -47.514 < 2e-16 ***
## Anni.madre 0.8772 1.1487 0.764 0.4452
## N.gravidanze 11.4029 4.6745 2.439 0.0148 *
## Fumatrici_d_fS -30.2865 27.5981 -1.097 0.2726
## Gestazione 32.8936 3.8259 8.598 < 2e-16 ***
## Lunghezza 10.2348 0.3009 34.010 < 2e-16 ***
## Cranio 10.5192 0.4268 24.644 < 2e-16 ***
## SessoM 78.0898 11.2042 6.970 4.05e-12 ***
## ---
## 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-16car::vif(mod2)## Anni.madre N.gravidanze Fumatrici_d_f Gestazione Lunghezza
## 1.189220 1.187430 1.006678 1.693694 2.078666
## Cranio Sesso
## 1.628877 1.040366I valori di \(R^2\) aggiustato e del p-value non sono cambiati rispetto al modello 1, pertanto Ospedale e Tipo.parto sono effettivamente variabili non significative per la regressione.
Togliamo gli anni della madre.
mod3<-update(mod2,~.-Anni.madre)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici_d_f + 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_d_fS -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-16car::vif(mod3)## N.gravidanze Fumatrici_d_f Gestazione Lunghezza Cranio
## 1.026120 1.006607 1.675575 2.078644 1.624603
## Sesso
## 1.040271Anche in questo caso i valori di \(R^2\) aggiustato e del p-value non sono cambiati rispetto al modello 2, pertanto Anni.madre è una variabile non significativa per la regressione.
Proviamo a tenere solo Gestazione, Lunghezza, Cranio e Sesso, le variabili maggiormente significative.
mod4<-lm(Peso~ Gestazione + Lunghezza + Cranio + Sesso , data=dati)
summary(mod4)##
## 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-16car::vif(mod4)## Gestazione Lunghezza Cranio Sesso
## 1.653502 2.069517 1.606131 1.038813Il valore di \(R^2\) aggiustato è lievemente diminuito, ma stiamo considerando un modello di regressione con sole quattro variabili esplicative, quindi possiamo comunque considerarlo un buon modello.
Proviamo a togliere Lunghezza e a mettere una dipendenza di Gestazione insieme a Cranio e Sesso.
mod5<-lm(Peso ~ + Gestazione + Cranio + Sesso, data=dati)
summary(mod5)##
## Call:
## lm(formula = Peso ~ +Gestazione + Cranio + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1316.39 -219.02 -10.96 210.38 1532.14
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6214.5285 163.0677 -38.110 <2e-16 ***
## Gestazione 92.6067 4.0208 23.032 <2e-16 ***
## Cranio 17.1449 0.4583 37.408 <2e-16 ***
## SessoM 118.5296 13.4793 8.793 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 332.4 on 2496 degrees of freedom
## Multiple R-squared: 0.5996, Adjusted R-squared: 0.5992
## F-statistic: 1246 on 3 and 2496 DF, p-value: < 2.2e-16car::vif(mod5)## Gestazione Cranio Sesso
## 1.276694 1.281694 1.027662Il valore di \(R^2\) aggiustato è sceso notevolmente, pertanto si deduce che la variabile Lunghezza non può essere eliminata dal modello.
In base al criterio di informazione di Bayes (BIC) il modello 4 è quello da preferire, quello che tiene conto di Gestazione, Lunghezza, Cranio e Sesso.
BIC(mod1, mod2, mod3, mod4, mod5)## df BIC
## mod1 12 35241.97
## mod2 9 35233.94
## mod3 8 35226.70
## mod4 6 35220.54
## mod5 5 36161.68Vediamo quale è il modello selezionato da R tramite la procedura stepwise usando la funzione stepAIC. Inseriamo un modello iniziale che tenga conto di tutte le variabili del dataset.
n<-nrow(dati)
stepwise.mod<-MASS::stepAIC(lm(Peso~ Anni.madre + N.gravidanze + Fumatrici_d_f + Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data=dati),
direction = "both",
k=log(n)) #CRITERIO BIC## Start: AIC=28139.46
## Peso ~ Anni.madre + N.gravidanze + Fumatrici_d_f + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 36392 186809099 28132
## - Fumatrici_d_f 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_d_f + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici_d_f 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_d_f 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_d_f 1 99840 187501837 28126
## + Anni.madre 1 43845 187557831 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
## 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_d_f 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 29066La funzione di R stepAIC, dopo aver fornito un modello con tutte le variabili coinvolte seleziona il modello che tiene conto di N.gravidanze, Lunghezza, Cranio, Sesso e Gestazione. Proviamo a costruire questo modello e a confrontarlo con gli altri.
mod6<-lm(Peso ~ Gestazione + Lunghezza+ Cranio + Sesso + N.gravidanze, data=dati)
summary(mod6)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso +
## N.gravidanze, 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 ***
## Gestazione 32.3321 3.7980 8.513 < 2e-16 ***
## Lunghezza 10.2486 0.3006 34.090 < 2e-16 ***
## Cranio 10.5402 0.4262 24.728 < 2e-16 ***
## SessoM 77.9927 11.2021 6.962 4.26e-12 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1328 on 5 and 2494 DF, p-value: < 2.2e-16Effettivamente dall’analisi si nota un valore di \(R^2\) aggiustato leggermente superiore a quello del modello 4, quindi possiamo dire che il peso del neonato dipende anche dal numero di gravidanze della gestante.
Proviamo a fare un test ANOVA, per vedere se c’è un aumento o diminuzione significativo della varianza.
anova(mod4,mod6)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ Gestazione + Lunghezza + Cranio + Sesso + N.gravidanze
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2495 188688687
## 2 2494 188065546 1 623141 8.2637 0.004079 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Con un p-value di 0.004079 viene preferito il modello 6.
Infine confrontiamoli con il criterio BIC.
BIC(mod4,mod6)## df BIC
## mod4 6 35220.54
## mod6 7 35220.10Il valore di BIC del modello6 è più basso del modello4 (anche se veramente di pochissimo), quindi si preferisce il modello4, dato che con una variabile in meno ha praticamente la stessa accuratezza del modello 6. Per il “Principio di parsimonia” infatti modelli più semplici sono preferiti a modelli complessi, e si devono usare parametri addizionali solo se essi sono strettamente necessari.
Da questo momento in poi il modello4 sarà il nostro modello.
Analisi della Qualità del Modello
Dopo aver valutato il valore di \(R^2\) aggiustato di volta in volta nei modelli selezionati, si procede con la diagnostica sui residui.
Si deve verificare che le seguenti assunzioni del modello non siano violate:
- linearità
- normalità
- omoschedasticità
- indipendenza
e si deve verificare la presenza di valori anomali o influenti.
Il modello selezionato nel precedente paragrafo è il modello 4. Analizziamone i residui.
par(mfrow=c(2,2))
plot(mod4)
Residuals vs Fitted: per essere un buon modello i residui devono posizionarsi attorno allo zero. Nel nostro caso c’è una curvatura per valori molto piccoli, il resto dei residui sembra essere simmetrico rispetto allo zero.
Q-Q Residuals: rappresenta la relazione dei residui con i quantili di una normale. I quantili devono seguire l’andamento della bisettrice; nel nostro caso alle code ci sono dei valori che non seguono l’andamento della bisettrice.
Scale-Location: studia l’omoschedasticità, ossia la varianza costante; i residui devono quindi seguire una linea orizzontale. Nel nostro caso c’è una piccola curvatura per valori bassi.
Residuals vs Leverage: si visualizzano i potenziali valori influenti, cioè i valori molto leverage, molto outliers oppure una combinazione dei due. I valori non devono superare le soglie di Cook (0.5 e 1). Nel nostro caso c’è un unico valore (1551) che ha superato la soglia 0.5.
Dopo l’analisi grafica, procediamo con l’analisi dei residui tramite test statistici.
Durbin-Watson test per l’autocorrelazione dei residui
lmtest::dwtest(mod4)##
## Durbin-Watson test
##
## data: mod4
## DW = 1.9557, p-value = 0.1337
## alternative hypothesis: true autocorrelation is greater than 0Dall’esito del test Durbin-Watson non si rifiuta l’ipotesi nulla, quindi i residui non sono autocorrelati.
Shapiro-Wilk test per la verifica della normalità
shapiro.test(mod4$residuals)##
## Shapiro-Wilk normality test
##
## data: mod4$residuals
## W = 0.9742, p-value < 2.2e-16Considerando il test si rifiuta l’ipotesi di normalità.
Test di Breusch-Pagan, test d’ipotesi di eteroschedasticità
lmtest::bptest(mod4)##
## studentized Breusch-Pagan test
##
## data: mod4
## BP = 89.148, df = 4, p-value < 2.2e-16Dall’esito del Breusch-Pagan test si rifiuta l’ipotesi di omoschedasticità, quindi la varianza non può essere considerata costante.
Guardando il grafico sottostante ci si rende conto che i residui del modello 4 sono molto simili a una distribuzione normale di media 0 e deviazione standard pari a 250, quindi graficamente possiamo confermare che il modello non viola le assunzioni di linearità, normalità, omoschedasticità e indipendenza.
plot(density(residuals(mod4)))+
lines(density(rnorm(100000,0,250)),col=2)
## integer(0)Analizziamo i leverage e gli outliers.
#leverage
lev<-hatvalues(mod4)
p<-sum(lev)
n<-length(lev)
soglia=2*p/n
plot(lev)
abline(h=soglia,col=2)
lev[lev>soglia]## 15 34 42 61 67 96
## 0.006144457 0.006237758 0.004252678 0.004587185 0.005345322 0.004801899
## 101 106 117 131 151 155
## 0.007185969 0.012812305 0.004746840 0.006964953 0.010847336 0.006670119
## 190 205 206 220 249 295
## 0.005288500 0.005297016 0.009402903 0.007376867 0.004655679 0.004008437
## 304 305 310 312 315 348
## 0.004419003 0.005382162 0.028757757 0.013126063 0.005342858 0.004187962
## 378 383 445 471 486 492
## 0.015366923 0.004305772 0.007094099 0.004289364 0.004740595 0.008175223
## 565 587 592 615 629 638
## 0.004635950 0.008375235 0.006313937 0.004575867 0.004041054 0.006216787
## 656 666 684 697 702 726
## 0.004735898 0.004345904 0.008750225 0.005809934 0.004719868 0.004038656
## 748 750 765 805 821 895
## 0.008236046 0.006712718 0.006049647 0.014303716 0.004039952 0.005203974
## 928 947 956 964 968 991
## 0.022095348 0.007803378 0.007670034 0.004618460 0.004075295 0.004259145
## 1014 1067 1091 1096 1130 1157
## 0.008221844 0.007868942 0.008927908 0.004268008 0.006561457 0.004080725
## 1166 1181 1188 1200 1238 1248
## 0.004020427 0.005586402 0.006404596 0.005445729 0.005374725 0.014573080
## 1273 1283 1293 1294 1356 1357
## 0.007080183 0.004055883 0.005555616 0.004735838 0.005285646 0.006526609
## 1361 1385 1395 1400 1402 1420
## 0.004080786 0.012017154 0.004646054 0.005559147 0.004788164 0.005130876
## 1428 1429 1551 1553 1556 1560
## 0.008191659 0.021037040 0.048521821 0.006810651 0.005868384 0.004588261
## 1593 1606 1610 1619 1628 1634
## 0.004855489 0.004746852 0.007761683 0.014504316 0.004663751 0.004527422
## 1686 1692 1693 1701 1712 1735
## 0.008715871 0.004260736 0.005041490 0.010197488 0.006945920 0.004370535
## 1780 1802 1806 1809 1827 1858
## 0.025528862 0.004005956 0.004090967 0.008388773 0.006008547 0.004328741
## 1868 1893 1911 1920 1977 2037
## 0.004746360 0.004245302 0.004481835 0.004131544 0.005384125 0.004234559
## 2040 2066 2089 2114 2115 2120
## 0.011429685 0.004013637 0.006252340 0.012850646 0.011763363 0.018002540
## 2140 2149 2175 2200 2215 2216
## 0.006090376 0.013033782 0.021167345 0.011030260 0.004833388 0.007548562
## 2220 2224 2225 2257 2307 2318
## 0.004023907 0.005780067 0.005408611 0.005767952 0.013898144 0.004770755
## 2337 2359 2391 2408 2437 2452
## 0.004700477 0.004750437 0.004386621 0.009632048 0.023757012 0.023227398
## 2458 2476 2478
## 0.008002531 0.004070348 0.005745434Ci sono molti valori inusuali nello spazio dei regressori.
Analizziamo gli outliers.
#outliers
plot(rstudent(mod4))
abline(h=c(-2,2),col=2)
car::outlierTest(mod4)## rstudent unadjusted p-value Bonferroni p
## 1551 9.986149 4.7193e-23 1.1798e-19
## 155 4.951276 7.8654e-07 1.9663e-03
## 1306 4.781188 1.8440e-06 4.6100e-03Ci sono tre valori outliers (1551, 155 e 1306).
Calcoliamo la distanza di Cook.
#distanza di cook
cook<-cooks.distance(mod4)
plot(cook) 
mydata <- data.frame(Max.dist.di.Cook=round(max(cook),2))
knitr::kable((mydata), "pipe")| Max.dist.di.Cook |
|---|
| 0.98 |
La soglia del 0.5 è stata superata, ed è anzi quasi stata raggiunta la soglia di 1. Il valore è 1551, già notato nelle precedenti diagnostiche. Non superando la soglia di 1 non lo consideriamo come una anomalia vera e propria e lo teniamo nel modello.
knitr::kable((dati[1551,1:10]), "pipe")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1551 | 35 | 1 | 0 | 38 | 4370 | 315 | 374 | Nat | osp3 | F |
Analizzando nel dettaglio la singola osservazione, possiamo notare un peso pari a 4370 contro una lunghezza di 315 e un cranio di 374. Effettivamente la misura della lunghezza e quella del cranio non sono consistenti, poiché la neonata ha un corpo troppo corto per il suo peso e la grandezza del cranio.
3. Previsioni e Risultati
Si potrebbe stimare il peso di una neonata considerando una madre alla terza gravidanza, non fumatrice, che partorirà alla 39esima settimana. Oltre a tali dati forniti per la previsione, per il modello costruito è necessario conoscere anche la lunghezza e il cranio, invece il dato relativo al fumo della gestante e quello relativo al numero delle gravidanze saranno ignorati dal modello, poiché per la creazione del modello di regressione lineare non sono state usate né variabili esplicative associate al fumo né variabili relative al numero delle gravidanze.
Proviamo quindi a prendere come lunghezza la lunghezza media e come cranio il valore medio di cranio.
previsione=predict(mod4,newdata = data.frame(Sesso='F', N.gravidanze=3,Gestazione=39, Fumatrici=0, Lunghezza=round(mean(Lunghezza),2), Cranio=round(mean(Cranio),2))) mydata <- data.frame(Previsione=round(previsione,2), Sesso='F',N.gravidanze=3,Gestazione=39, Fumatrici=0, Mean.Lunghezza=round(mean(Lunghezza),2), Mean.Cranio=round(mean(Cranio),2))
knitr::kable((mydata), "pipe")| Previsione | Sesso | N.gravidanze | Gestazione | Fumatrici | Mean.Lunghezza | Mean.Cranio |
|---|---|---|---|---|---|---|
| 3245.32 | F | 3 | 39 | 0 | 494.69 | 340.03 |
Il peso della neonata, stimato in base al modello di regressione costruito, è pari a 3245.32
4. Visualizzazioni
Visualizziamo l’impatto del numero di settimane di gestazione e del sesso del neonato sul peso previsto.
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, #standardError=False
method = "lm")+
labs(title='Andamento del peso al variare delle settimane di gestazione per genere')+
theme_classic()## `geom_smooth()` using formula = 'y ~ x'
La nuvola di punti è maggiormente concentrata dalle 35 settimane di gestazione in poi, le rette relative ai due generi hanno la stessa pendenza, ma la retta dei maschi è più in alto perché i maschi pesano di più.
Visualizziamo ora l’impatto del numero di settimane di gestazione e del fumo della gestante sul peso previsto per il neonato.
ggplot(data=dati)+
geom_point(aes(x=Gestazione,
y=Peso,
col=Fumatrici_d_f),position='jitter')+
geom_smooth(aes(x=Gestazione,
y=Peso,
col=Fumatrici_d_f),
se=F, #standardError=False
method = "lm")+
labs(title='Andamento del peso al variare delle settimane di gestazione per fumatrici',col='Fumatrici')+
theme_classic()## `geom_smooth()` using formula = 'y ~ x'
Come detto anche in fase di analisi, i dati non ci consentono di fare un confronto fra lo stato di fumatrice e di non fumatrice poiché le gestanti fumatrici sono poche nel campione assegnato.
Visualizziamo ora l’impatto del numero di gravidanze e del genere sul peso previsto per il neonato.
ggplot(data=dati)+
geom_point(aes(x=N.gravidanze,
y=Peso,
col=Sesso),position='jitter')+
geom_smooth(aes(x=N.gravidanze,
y=Peso,
col=Sesso),
se=F, #standardError=False
method = "lm")+
labs(title='Andamento del peso al variare del numero di gravidanze per genere',x='Numero di gravidanze')+
theme_classic()## `geom_smooth()` using formula = 'y ~ x'
Il peso dei neonati risulta essere maggiore di quelle delle neonate, ma rimane costante all’aumentare del numero delle gravidanze della gestante.
Conclusioni
In conclusione, rispetto ai dati a nostra disposizione, siamo riusciti a costruire un modello di regressione lineare per la previsione del peso del neonato utilizzando quattro variabili: Gestazione, Lunghezza, Cranio e Sesso. La variabile Fumatrici invece non ha dato nessun esito utile per la previsione, al contrario di quello che ci aspettavamo. Fra le variabili esplicative non si è riscontrata multicollinearità.
E’ stato individuato un valore di outlier problematico, che superava la soglia di attenzione della distanza di Cook di 0.5, e dalle analisi risultava esserci qualche errore nell’inserimento dei dati, poiché i valori non risultavano coerenti fra loro.
I test statistici sui residui non hanno avuto esito positivo, ma graficamente si notava come essi rientrassero nei parametri e non violassero le assunzioni del modello.
Il modello selezionato ha un valore di \(R^2\) aggiustato di circa \(0.7257\), quindi può essere considerato un buon modello per la previsione del peso dei neonati.