library(ggplot2) library(moments) library(Metrics) library(lmtest) library(car)
neonati <-read.csv("neonati.csv", sep=",")
neonati_puliti <- na.omit(neonati)
Analisi Preliminare a) esploriamo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.
library(moments)
attach(neonati)
neonati$Anni.madre[neonati$Anni.madre < 12] <- NA
summary(neonati$Anni.madre)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 13.00 25.00 28.00 28.19 32.00 46.00 2
IQR(neonati$Anni.madre, na.rm = TRUE)
## [1] 7
range(neonati$Anni.madre, na.rm = TRUE)
## [1] 13 46
var(neonati$Anni.madre, na.rm = TRUE)
## [1] 27.21924
sd(neonati$Anni.madre, na.rm = TRUE)
## [1] 5.217206
skewness(neonati$Anni.madre, na.rm = TRUE)
## [1] 0.1510624
kurtosis(neonati$Anni.madre, na.rm = TRUE) - 3
## [1] -0.1056061
La variabile “Anni.madre” corrisponde all’età della madre espressa in anni. Si tratta di una variabile quantitativa discreta. Dal momento che sono presenti due entries in cui come età viene indicato 0 o 1 , rimuoviamo dal dataset queste due osservazioni in quanto capaci di influenzare i calcoli dei vari indici. La madre più giovane ha 13 anni, la madre più anziana 46. In media le madri interessate dallo studio hanno 28 anni. L’indice di asimmetria è di 0.1510624 (distribuzione asimmetrica positiva), per cui sono più frequenti modalità con valori bassi e abbiamo anche media>mediana>moda. L’’indice di Curtosi è di -0.105661 (distribuzione platicurtica), con una distribuzione dei valori intorno alla media più appiattita rispetto a una distribuzione normale.
summary(N.gravidanze)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.9812 1.0000 12.0000
IQR(N.gravidanze)
## [1] 1
range(N.gravidanze)
## [1] 0 12
var(N.gravidanze)
## [1] 1.639903
sd(N.gravidanze)
## [1] 1.280587
skewness(N.gravidanze)
## [1] 2.514254
kurtosis(N.gravidanze) - 3
## [1] 10.98941
La variabile “N.gravidanze” corrisponde al numero di gravidanze avute
dalla madre. Si tratta di una variabile quantitativa discreta. Il numero
minimo di gravidanze avute trovato è 0, intendendo quindi la gravidanza
al momento dello studio come la prima gravidanza, il numero massimo di
gravidanze trovato è 12, per cui le madri in questione si trovano alla
loro tredicesima gravidanza.
In media le madri interessate dallo studio hanno avuto più o meno una
gravidanza . L’indice di asimmetria è di 2.514254 (distribuzione
asimmetrica positiva), per cui sono più frequenti modalità con valori
bassi (numero di gravidanze basso) e abbiamo anche
media>mediana>moda. L’’indice di Curtosi è di 10.98941
(distribuzione leptocurtica), con una distribuzione dei valori intorno
alla media più assottigliata rispetto a una distribuzione normale.
summary(Fumatrici)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.0416 0.0000 1.0000
IQR(Fumatrici)
## [1] 0
range(Fumatrici)
## [1] 0 1
var(Fumatrici)
## [1] 0.03988539
sd(Fumatrici)
## [1] 0.1997133
skewness(Fumatrici)
## [1] 4.591499
kurtosis(Fumatrici) - 3
## [1] 19.08187
table(Fumatrici)
## Fumatrici
## 0 1
## 2396 104
La variabile Fumatrici è una variabile “dummy” qualitativa codificata con un numero da 0 a 1 a seconda che la madre sia fumatrice o meno.
summary(Gestazione)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 25.00 38.00 39.00 38.98 40.00 43.00
IQR(Gestazione)
## [1] 2
range(Gestazione)
## [1] 25 43
var(Gestazione)
## [1] 3.491813
sd(Gestazione)
## [1] 1.868639
skewness(Gestazione)
## [1] -2.065313
kurtosis(Gestazione) - 3
## [1] 8.25815
table(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 2
La variabile “Gestazione” corrisponde al numero di settimane di
gestazione al momento della nascita del bambino. Si tratta di una
variabile quantitativa discreta. Il numero minimo di settimane di
gestazione è 25, il numero massimo di settimane di gestazione è
43.
In media le madri interessate hanno partorito alla 39esima settimana .
L’indice di asimmetria è di -2.065313(distribuzione asimmetrica
negativa), per cui sono più frequenti modalità con valori alti (elevato
numero di settimane di gestazione) e abbiamo anche
media>mediana>moda. L’’indice di Curtosi è di 8.25815
(distribuzione leptocurtica), con una distribuzione dei valori intorno
alla media più assottigliata rispetto a una distribuzione normale.
summary(Peso)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 2990 3300 3284 3620 4930
IQR(Peso)
## [1] 630
range(Peso)
## [1] 830 4930
var(Peso)
## [1] 275665.7
sd(Peso)
## [1] 525.0387
skewness(Peso)
## [1] -0.6470308
kurtosis(Peso) - 3
## [1] 2.031532
table(Peso)
## Peso
## 830 900 930 980 990 1140 1170 1180 1190 1280 1285 1300 1340 1370 1390 1410
## 1 1 2 1 1 1 2 1 1 3 1 1 1 1 1 1
## 1430 1450 1500 1550 1560 1580 1600 1615 1620 1690 1720 1730 1750 1770 1780 1800
## 1 1 2 1 1 1 1 1 1 1 1 1 3 1 2 1
## 1840 1850 1890 1900 1950 1960 1970 1980 2000 2040 2050 2060 2090 2100 2120 2150
## 1 1 1 1 1 1 2 2 3 2 3 2 1 5 1 2
## 2160 2180 2200 2210 2220 2230 2250 2260 2270 2280 2290 2300 2310 2320 2330 2340
## 2 1 2 1 3 2 2 3 1 1 2 3 1 3 1 3
## 2350 2370 2380 2390 2400 2410 2420 2430 2440 2450 2460 2470 2480 2490 2500 2510
## 2 1 2 1 6 3 2 4 4 6 1 2 2 1 10 3
## 2520 2530 2540 2550 2560 2570 2580 2590 2600 2610 2620 2630 2640 2650 2660 2670
## 5 2 1 7 7 2 6 6 9 2 5 4 3 8 5 6
## 2680 2690 2700 2710 2720 2730 2740 2750 2760 2770 2780 2790 2800 2810 2820 2830
## 11 4 15 7 9 4 16 18 6 5 9 3 17 9 12 12
## 2840 2850 2860 2862 2870 2880 2890 2900 2910 2920 2930 2940 2950 2960 2970 2980
## 10 20 9 1 7 22 11 31 12 13 10 20 23 11 12 17
## 2990 3000 3010 3020 3030 3040 3050 3060 3070 3080 3090 3100 3110 3120 3130 3140
## 14 29 10 11 21 14 29 14 14 20 12 43 17 11 11 23
## 3150 3160 3170 3180 3190 3200 3210 3220 3230 3240 3250 3260 3270 3280 3290 3300
## 33 12 22 26 19 38 9 22 17 18 34 20 17 24 24 56
## 3310 3320 3330 3340 3350 3360 3370 3380 3390 3400 3410 3420 3430 3440 3450 3460
## 18 15 20 21 27 14 19 27 7 35 10 14 13 21 25 11
## 3470 3480 3490 3500 3510 3520 3530 3540 3550 3560 3570 3580 3590 3600 3610 3620
## 9 18 6 45 13 20 20 17 30 13 14 18 13 31 9 15
## 3630 3640 3650 3660 3670 3680 3690 3700 3710 3720 3730 3740 3750 3760 3770 3780
## 14 17 18 10 14 17 8 25 9 12 11 8 22 12 11 13
## 3790 3800 3810 3820 3830 3840 3850 3860 3870 3880 3890 3900 3910 3920 3930 3940
## 8 25 6 18 9 13 17 11 11 4 5 16 11 8 5 7
## 3950 3960 3970 3980 3990 4000 4010 4020 4030 4040 4050 4060 4070 4080 4090 4100
## 20 6 4 3 6 15 3 5 4 2 12 7 3 3 4 6
## 4110 4120 4130 4140 4150 4160 4170 4180 4190 4200 4220 4230 4240 4250 4260 4270
## 3 4 5 8 7 4 1 3 1 7 4 1 4 5 4 2
## 4280 4290 4300 4310 4320 4330 4340 4350 4370 4400 4410 4420 4440 4470 4480 4520
## 2 3 1 3 3 4 1 3 2 4 2 2 2 1 2 1
## 4540 4550 4560 4580 4600 4620 4650 4680 4690 4700 4720 4760 4810 4900 4930
## 1 1 1 1 2 1 1 1 1 1 2 1 1 1 1
La variabile “Peso” corrisponde al peso in grammi del neonato al momento della nascita . Si tratta di una variabile quantitativa continua Il peso più basso è di 830 g, il più alto è di 4930 g. In media il peso dei bambini alla nascita è di 3284 g. L’indice di asimmetria è di -0.6470308(distribuzione asimmetrica negativa), per cui sono più frequenti bambini con peso alla nascita più alto (dato coerente con la distribuzione delle settimane di gestazione) e abbiamo anche media>mediana>moda. L’’indice di Curtosi è di 2.031532 (distribuzione leptocurtica), con una distribuzione dei valori intorno alla media più assottigliata rispetto a una distribuzione normale.
summary(Lunghezza)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 310.0 480.0 500.0 494.7 510.0 565.0
IQR(Lunghezza)
## [1] 30
range(Lunghezza)
## [1] 310 565
var(Lunghezza)
## [1] 692.671
sd(Lunghezza)
## [1] 26.31864
skewness(Lunghezza)
## [1] -1.514699
kurtosis(Lunghezza) - 3
## [1] 6.487174
La variabile “lunghezza”, quantiativa continua si riferisce alla lunghezza dei neonati. Il valore più basso è di 310 mm, quello più alto di 565 mm con un valore medio di 497mm. L’indice di asimmetria è di -1.514699(distribuzione asimmetrica negativa), per cui sono più frequenti bambini con lunghezza craniale maggiore (dato coerente con la distribuzione delle settimane di gestazione) e abbiamo anche media>mediana>moda. L’’indice di Curtosi è di 6.487174 (distribuzione leptocurtica), con una distribuzione dei valori intorno alla media più assottigliata rispetto a una distribuzione normale.
summary(Cranio)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 235 330 340 340 350 390
IQR(Cranio)
## [1] 20
range(Cranio)
## [1] 235 390
var(Cranio)
## [1] 269.7915
sd(Cranio)
## [1] 16.42533
skewness(Cranio)
## [1] -0.7850527
kurtosis(Cranio) - 3
## [1] 2.946206
La variabile “cranio”, quantiativa continua si riferisce al diametro del cranio dei neonati. Il valore più basso è di 235 mm, quello più alto di 390 mm con un valore medio di 340 mm L’indice di asimmetria è di -0.7850527(distribuzione asimmetrica negativa), per cui sono più frequenti bambini con diametro craniale maggiore (dato coerente con la distribuzione delle settimane di gestazione) e abbiamo anche media>mediana>moda. L’’indice di Curtosi è di 2.946206 (distribuzione leptocurtica), con una distribuzione dei valori intorno alla media più assottigliata rispetto a una distribuzione normale.
campione_parto <- ifelse(Tipo.parto=="Nat",1,0)
table(campione_parto)
## campione_parto
## 0 1
## 728 1772
La variabile “Tipo parto” è una variabile qualitativa nominale. Sono presenti 728 donne che hanno avuto un parto cesareo e 1772 che hanno avuto un parto naturale.
table(Ospedale)
## Ospedale
## osp1 osp2 osp3
## 816 849 835
La variabile “Ospedale” è una variabile qualitativa nominale. Vi è una distribuzione abbastanza omogenea delle donne tra i tre ospedali analizzati.
table(Sesso)
## Sesso
## F M
## 1256 1244
Saggiamo l’ipotesi che in alcuni ospedali si facciano più parti cesarei di altri
tabella_parti <- table(Ospedale,Tipo.parto)
attese_tabella <- tabella_parti
margin.table(tabella_parti,1)
## Ospedale
## osp1 osp2 osp3
## 816 849 835
margin.table(tabella_parti,2)
## Tipo.parto
## Ces Nat
## 728 1772
n =margin.table(tabella_parti)
n
## [1] 2500
for(i in 1:nrow(tabella_parti)){
for(j in 1:ncol(tabella_parti)){
attese_tabella[i,j] <-(margin.table(tabella_parti,1)[i]*margin.table(tabella_parti,2)[j])/n
}
}
attese_tabella
## Tipo.parto
## Ospedale Ces Nat
## osp1 237.6192 578.3808
## osp2 247.2288 601.7712
## osp3 243.1520 591.8480
#ORA CALCOLIAMO LA STATISTICA TEST X^2 E CONDUCIAMO TEST DI INDIPENDENZA
# statistica test: sommatoria di (freq. osservate-freq. attese)^2/freq.attese
osservate <- tabella_parti
x_quadro <- sum((osservate- attese_tabella)^2/attese_tabella)
x_quadro
## [1] 1.097202
test.indipendenza <- chisq.test(tabella_parti)
test.indipendenza
##
## Pearson's Chi-squared test
##
## data: tabella_parti
## X-squared = 1.0972, df = 2, p-value = 0.5778
Il valore di X_quadro e il p_value ci portano a non rifiutare H₀, ovvero l’ipotesi nulla secondo cui il tipo di parto non dipende dall’ospedale, quindi non c’è correlazione tra tipo di parto e ospedale.
Saggiamo l’ipotesi che la media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione 1)Media del peso: da ricerche effettuate risulta che la media del peso dei bambini alla nascita sia di circa 3.5 kg, confrontiamola con quella del campione
H₀: media campione = media popolazione H₁: media campione ≠ media popolazione
t.test(Peso,
mu = 3.5,
conf.level = 0.95,
alternative = "two.sided")
##
## One Sample t-test
##
## data: Peso
## t = 312.41, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 3.5
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081
Il p-value ottenuto ha un valore molto basso che ci consente di rifiutare l’ipotesi nulla, per cui possiamo affermare che il peso medio dei neonati del campione osservato è significativamente inferiore a quello della popolazione.
1)Media della lunghezza: da ricerche effettuate risulta che la media del peso dei bambini alla nascita sia di circa 500 mm, confrontiamola con quella del campione
H₀: media campione = media popolazione H₁: media campione ≠ media popolazione
t.test(Lunghezza,
mu = 500,
conf.level = 0.95,
alternative = "two.sided")
##
## One Sample t-test
##
## data: Lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 493.6598 495.7242
## sample estimates:
## mean of x
## 494.692
Il p-value ottenuto ha un valore molto basso che ci consente di rifiutare l’ipotesi nulla, per cui possiamo affermare che la lunghezza media dei neonati del campione osservato è significativamente inferiore a quella della popolazione.
Saggiamo l’ipotesi che le misure antropometriche (peso,lunghezza,cranio) siano significativamente diverse tra i due sessi
Test t per differenza tra medie tra campioni indipendenti
table(Sesso)
## Sesso
## F M
## 1256 1244
#costruiamo un boxplot condizionato
boxplot(Peso~Sesso, data=neonati)
#vediamo le statistiche di sintesi condizionate ai due sessi
summary(neonati$Peso[neonati$Sesso == "M"])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 980 3150 3430 3408 3720 4810
summary(neonati$Peso[neonati$Sesso == "F"])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 2900 3160 3161 3470 4930
t.test(Peso~Sesso,
data=neonati)
##
## 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
Confrontando il peso tra i sessi vediamo che i maschi hanno un peso significativamente maggiore rispetto alle femmine.
table(Sesso)
## Sesso
## F M
## 1256 1244
#costruiamo un boxplot condizionato
boxplot(Lunghezza~Sesso, data=neonati)
#vediamo le statistiche di sintesi condizionate ai due sessi
summary(neonati$Lunghezza[neonati$Sesso == "M"])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 320.0 490.0 500.0 499.7 515.0 560.0
summary(neonati$Lunghezza[neonati$Sesso == "F"])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 310.0 480.0 490.0 489.8 505.0 565.0
t.test(Lunghezza~Sesso,
data=neonati)
##
## 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
Confrontando la lunghezza tra i sessi vediamo che i maschi hanno una lunghezza significativamente maggiore rispetto alle femmine.
table(Sesso)
## Sesso
## F M
## 1256 1244
#costruiamo un boxplot condizionato
boxplot(Cranio~Sesso, data=neonati)
#vediamo le statistiche di sintesi condizionate ai due sessi
summary(neonati$Cranio[neonati$Sesso == "M"])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 265.0 334.0 343.0 342.4 352.0 390.0
summary(neonati$Cranio[neonati$Sesso == "F"])
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 235.0 330.0 340.0 337.6 348.2 390.0
t.test(Cranio~Sesso,
data=neonati)
##
## 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.4486
Confrontando la lunghezza tra i sessi vediamo che i maschi hanno un diametro craniale significativamente maggiore rispetto alle femmine.
Creazione del Modello di Regressione Verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti
# 1) verifichiamo che la variabile risposta sia normale con
# a.indici di forma
# b.test shapiro.wilk
# a. calcoliamo l'indice di asimmetria e curtosi
moments::skewness(Peso)
## [1] -0.6470308
moments::kurtosis(Peso)-3
## [1] 2.031532
# b. shapiro test.
shapiro.test(Peso)
##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16
# 2)Indaghiamo le relazioni tra variabili con la matrice di correlazione
# usiamo la funzione cor su tutto il datase
# selezioniamo solo le variabili numeriche
num_vars <- neonati_puliti[sapply(neonati_puliti, is.numeric)]
# calcoliamo la matrice di correlazione
round(cor(num_vars), 2)
## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## Anni.madre 1.00 0.38 0.01 -0.14 -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.14 -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
## Cranio
## Anni.madre 0.02
## N.gravidanze 0.04
## Fumatrici -0.01
## Gestazione 0.46
## Peso 0.70
## Lunghezza 0.60
## Cranio 1.00
#3)stimiamo il modello di regressione lineare multiple inserendo
# tutte le variabili esplicative (usiamo di nuovo la funzione lm)
mod1 <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Cranio + Lunghezza + Tipo.parto + Ospedale + Sesso, data=neonati_puliti)
summary(mod1)
##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Cranio + Lunghezza + Tipo.parto + Ospedale + Sesso, data = neonati_puliti)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1124.40 -181.66 -14.42 160.91 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6738.4762 141.3087 -47.686 < 2e-16 ***
## Anni.madre 0.8921 1.1323 0.788 0.4308
## N.gravidanze 11.2665 4.6608 2.417 0.0157 *
## Fumatrici -30.1631 27.5386 -1.095 0.2735
## Gestazione 32.5696 3.8187 8.529 < 2e-16 ***
## Cranio 10.4707 0.4260 24.578 < 2e-16 ***
## Lunghezza 10.2945 0.3007 34.236 < 2e-16 ***
## Tipo.partoNat 29.5254 12.0844 2.443 0.0146 *
## Ospedaleosp2 -11.2095 13.4379 -0.834 0.4043
## Ospedaleosp3 28.0958 13.4957 2.082 0.0375 *
## SessoM 77.5409 11.1776 6.937 5.08e-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.2 on 10 and 2489 DF, p-value: < 2.2e-16
Dai test sulla variabile risposta possiamo dedurre che , vista l’asimmetria negativa, i pesi alti sono leggermente più frequenti e che la distribuzione è leptocurtica. Il test di Shapiro Wilk ci porta a rifiutare l’ipotesi nulla di normalità della distribuzione: questo ci porta a prestare più attenzione alla normalità dei residui del modello. Le variabili peso, lunghezza, cranio e gestazione, dall’analisi della matrice di correlazione, sono fortemente correlate tra loro, ad indicare probabile presenza di multicollinearità. Con questo modello vediamo che l’aumento del numero di gravidanze ha un piccolo effetto singificativo sul peso; un aumento di una settimana di gestazione aumenta il peso in maniera significativa , circa 32 g in più; ogni mm in più di diametro del cranio aumenta il peso di circa 10.6 g; ogni mm in più di lunghezza aumenta il peso di circa 10.45 g; i bambini nati con parto naturale pesano in media 29.6 g in più; I maschi pesano in media 77 g in più. Da questo modello risulta non significativo l’effetto di queste variabili: Età della madre, fumo materno e nascita nell’ospedale numero 2.
Proviamo diversi modelli con la selezione stepwise Selezione del Modello Ottimale Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC), selezioneremo il modello più parsimonioso, eliminando le variabili non significative. Verranno considerati anche modelli con interazioni tra le variabili e possibili effetti non lineari.
AIC(mod1)
## [1] 35171.95
BIC(mod1)
## [1] 35241.84
modello_stepwise_AIC <- step(mod1)
## Start: AIC=28075.26
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Cranio +
## Lunghezza + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 46578 186809099 28074
## - Fumatrici 1 90019 186852540 28074
## <none> 186762521 28075
## - N.gravidanze 1 438452 187200974 28079
## - Tipo.parto 1 447929 187210450 28079
## - Ospedale 2 685979 187448501 28080
## - Sesso 1 3611021 190373542 28121
## - Gestazione 1 5458403 192220925 28145
## - Cranio 1 45326172 232088693 28616
## - Lunghezza 1 87951062 274713583 29038
##
## Step: AIC=28073.88
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Cranio + Lunghezza +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 90897 186899996 28073
## <none> 186809099 28074
## - Tipo.parto 1 448222 187257321 28078
## - Ospedale 2 692738 187501837 28079
## - N.gravidanze 1 633756 187442855 28080
## - Sesso 1 3618736 190427835 28120
## - Gestazione 1 5412879 192221978 28143
## - Cranio 1 45588236 232397335 28618
## - Lunghezza 1 87950050 274759149 29036
##
## Step: AIC=28073.1
## Peso ~ N.gravidanze + Gestazione + Cranio + Lunghezza + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 186899996 28073
## - Tipo.parto 1 440684 187340680 28077
## - Ospedale 2 701680 187601677 28078
## - N.gravidanze 1 610840 187510837 28079
## - Sesso 1 3602797 190502794 28119
## - Gestazione 1 5346781 192246777 28142
## - Cranio 1 45632149 232532146 28617
## - Lunghezza 1 88355030 275255027 29039
summary(modello_stepwise_AIC)
##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Cranio + Lunghezza +
## Tipo.parto + Ospedale + Sesso, data = neonati_puliti)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1113.18 -181.16 -16.58 161.01 2620.19
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.4293 135.9438 -49.340 < 2e-16 ***
## N.gravidanze 12.3619 4.3325 2.853 0.00436 **
## Gestazione 31.9909 3.7896 8.442 < 2e-16 ***
## Cranio 10.4922 0.4254 24.661 < 2e-16 ***
## Lunghezza 10.3086 0.3004 34.316 < 2e-16 ***
## Tipo.partoNat 29.2803 12.0817 2.424 0.01544 *
## Ospedaleosp2 -11.0227 13.4363 -0.820 0.41209
## Ospedaleosp3 28.6408 13.4886 2.123 0.03382 *
## SessoM 77.4412 11.1756 6.930 5.36e-12 ***
## ---
## 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.3 on 8 and 2491 DF, p-value: < 2.2e-16
L’analisi di regressione lineare multipla è stata ottimizzata utilizzando una procedura stepwise backward selection, basata sulla minimizzazione del Criterio di Informazione di Akaike (AIC). Questo processo iterativo ha portato alla rimozione delle variabili Anni.madre e Fumatrici, poiché la loro esclusione dal modello ha comportato una riduzione del valore di AIC, indicando una maggiore parsimonia e un migliore equilibrio tra adattamento e complessità.
Il modello finale selezionato presenta un R² di 0.7287 e un R² aggiustato di 0.7278. La prossimità di questi due valori conferma che le variabili rimosse non contribuivano in modo significativo alla capacità predittiva del modello. L’R² aggiustato indica che il modello finale spiega circa il 72.8% della varianza del peso neonatale, un risultato notevole che evidenzia l’ottima capacità esplicativa del modello
Analisi della Qualità del Modello Una volta ottenuto il modello finale, valuteremo la sua capacità predittiva utilizzando metriche come R² e il Root Mean Squared Error (RMSE). Un’attenzione particolare sarà rivolta all’analisi dei residui e alla presenza di valori influenti, che potrebbero distorcere le previsioni, indagando su di essi.
library(Metrics)
#Calcoliamo le previsioni del modello
previsioni <- predict(modello_stepwise_AIC,neonati_puliti)
#calcoliamo l'RMSE
rmse(neonati_puliti$Peso,previsioni)
## [1] 273.4227
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
# Test di normalità dei residui
shapiro.test(residuals(modello_stepwise_AIC))
##
## Shapiro-Wilk normality test
##
## data: residuals(modello_stepwise_AIC)
## W = 0.97408, p-value < 2.2e-16
# Test di omoschedasticità (Breusch-Pagan)
bptest(modello_stepwise_AIC)
##
## studentized Breusch-Pagan test
##
## data: modello_stepwise_AIC
## BP = 91.768, df = 8, p-value < 2.2e-16
# Test di indipendenza (Durbin-Watson)
dwtest(modello_stepwise_AIC)
##
## Durbin-Watson test
##
## data: modello_stepwise_AIC
## DW = 1.9527, p-value = 0.1184
## alternative hypothesis: true autocorrelation is greater than 0
L’analisi dei residui del modello ottimale ha rivelato che due assunzioni chiave della regressione lineare sono state violate. I test hanno indicato la non-normalità dei residui (p-value < 2.2e-16 nel test di Shapiro-Wilk) e la presenza di eteroschedasticità (p-value < 2.2e-16 nel test di Breusch-Pagan).
Previsioni e Risultati Una volta validato il modello, lo useremo per fare previsioni pratiche. Ad esempio, potremo stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.
Per ottenere il probabile valore del peso di una neonata alla nascita considerando una madre alla terza gravidanza che partorisce alla 39esima settimana, analizziamo innanzitutto i regressori di cui disponiamo : abbiamo il valore della variabile “N.gravidanze” (3), “Gestazione”(39) e “Sesso” (femmina); tuttavia nel nostro modello ci sono altre variabili (Cranio, lunghezza, Tipo.parto e Ospedale) di cui non conosciamo il valore, per cui possiamo procedere con una stima del valore medio di queste variabili all’interno del nostro campione . Per quanto riguarda le variabili qualitative “tipo parto” e “ospedale” userei la funzione “table” per vedere qual è la frequenza di ciascuno e sceglierei il tipo di parto e l’ospedale con frequenza più alta”
# Calcolo la media per il Cranio
mean(neonati_puliti$Cranio)
## [1] 340.0292
# Calcolo la media per la Lunghezza
mean(neonati_puliti$Lunghezza)
## [1] 494.692
# Conto le frequenze per Tipo.parto
table(neonati_puliti$Tipo.parto)
##
## Ces Nat
## 728 1772
# Conto le frequenze per Ospedale
table(neonati_puliti$Ospedale)
##
## osp1 osp2 osp3
## 816 849 835
nuovi_dati <- data.frame(
N.gravidanze = 3,
Gestazione = 39,
Cranio = 340.0292,
Lunghezza = 494.692,
Tipo.parto = "Nat",
Ospedale = "osp2",
Sesso = "F"
)
predict(modello_stepwise_AIC,nuovi_dati)
## 1
## 3262.81
Il peso stimato per una neonata con una madre alla terza gravidanza che partorisce alla 39esima settimana (e con valori medi per le altre variabili) è di 3262.81 grammi.
4. Visualizzazioni Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili. Ad esempio, potremmo visualizzare l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.
# Creo il grafico a dispersione con la retta di regressione
library(ggplot2)
ggplot(data = neonati_puliti, aes(x = Gestazione, y = Peso)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE, color = "blue") +
labs(title = "Relazione tra Gestazione e Peso alla Nascita",
subtitle = "con retta di regressione del modello",
x = "Settimane di Gestazione",
y = "Peso (g)") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
Il grafico a dispersione mostra chiaramente una forte relazione positiva
tra le settimane di gestazione e il peso alla nascita.
# Creo un boxplot per confrontare il peso tra i due gruppi
ggplot(data = neonati_puliti, aes(x = as.factor(Fumatrici), y = Peso)) +
geom_boxplot() +
labs(title = "Differenza di Peso tra Neonati di Madri Fumatrici e Non Fumatrici",
x = "Stato Fumatrice (0 = No, 1 = Sì)",
y = "Peso (g)") +
theme_minimal()
Il boxplot conferma visivamente che sebbene ci sia una piccola differenza tra le mediane dei due gruppi, le due scatole si sovrappongono in modo sostanziale. Questo suggerisce che la differenza non è abbastanza marcata da essere considerata significativa in un’analisi che tiene conto anche di altre variabili più influenti.
Conclusioni 1. Introduzione e Obiettivi
Il presente progetto è stato condotto con l’obiettivo di sviluppare uno strumento statistico affidabile per la previsione del peso neonatale. Utilizzando un modello di regressione lineare multipla, l’analisi mira a identificare e quantificare l’influenza di diverse variabili demografiche e cliniche sul peso alla nascita. Il modello predittivo che ne deriva rappresenta una risorsa fondamentale per il personale medico e le strutture sanitarie, consentendo una migliore preparazione e una gestione più efficace di eventuali complicanze legate al peso del neonato.
2. Sviluppo del Modello
Il dataset iniziale, composto da 2500 osservazioni, ha permesso di costruire un modello di regressione lineare multipla completo. Tuttavia, per ottimizzare l’accuratezza e la parsimonia del modello, è stata applicata la selezione delle variabili basata sul criterio di informazione di Akaike (AIC). Questo processo ha portato all’eliminazione di alcune variabili meno influenti, risultando in un modello più robusto, come confermato da un valore di R² elevato e valido.
3. Validazione e Affidabilità
Sebbene l’R² sia elevato l’analisi dei residui abbia evidenziato la presenza di non-normalità e di eteroschedasticità.
4. Previsioni Pratiche
Utilizzando il nostro modello validato e lo strumento
predict in R, è stata effettuata una previsione specifica.
Abbiamo stimato il peso alla nascita di una neonata, considerando una
madre alla terza gravidanza che partorisce alla 39esima settimana. Per
le variabili non specificate (Cranio, Lunghezza, Tipo.parto e Ospedale),
sono stati impiegati i valori medi o più frequenti estratti dal nostro
campione. Il risultato di circa 3262 grammi è pienamente conforme alle
aspettative per un neonato di sesso femminile nato in quella specifica
settimana di gestazione, dimostrando la praticità e l’efficacia del
modello.
5. Visualizzazioni e Risultati Chiave
Infine, le visualizzazioni grafiche hanno reso i risultati immediatamente comprensibili. Un grafico a dispersione ha evidenziato chiaramente la forte relazione positiva tra la durata della gestazione e il peso del neonato, confermando visivamente un risultato atteso e significativo del modello. Un boxplot ha messo a confronto la distribuzione del peso tra i neonati di madri fumatrici e non fumatrici. Il grafico ha mostrato che, sebbene ci sia una leggera differenza, l’impatto del fumo non è statisticamente significativo se si tiene conto degli altri predittori nel modello, un risultato cruciale che conferma le nostre scoperte.
Ho voluto però verificare l’assenza di multicollinearità tra le variabili dei due modelli(modello iniziale e modello ottenuto con la procedura stepwise)
library(car)
## Loading required package: carData
# Calcolo il VIF per il modello completo iniziale
vif(mod1)
## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.187454 1 1.089704
## N.gravidanze 1.186428 1 1.089233
## Fumatrici 1.007392 1 1.003689
## Gestazione 1.695810 1 1.302233
## Cranio 1.630796 1 1.277026
## Lunghezza 2.085755 1 1.444214
## Tipo.parto 1.004242 1 1.002119
## Ospedale 1.004071 2 1.001016
## Sesso 1.040643 1 1.020119
# Calcolo il VIF per il modello finale (modello_stepwise_AIC)
vif(modello_stepwise_AIC)
## GVIF Df GVIF^(1/(2*Df))
## N.gravidanze 1.025244 1 1.012544
## Gestazione 1.670237 1 1.292376
## Cranio 1.626507 1 1.275346
## Lunghezza 2.081921 1 1.442886
## Tipo.parto 1.003872 1 1.001934
## Ospedale 1.003038 2 1.000759
## Sesso 1.040337 1 1.019969
Da questa analisi, noto che sia nel modello iniziale (mod1, con tutte le variabili), sia nel modello finale (dopo la procedura stepwise) i valori di VIF sono tutti inferiori a 5, per cui non riscontro problemi di multicollinearità.
Sebbene questo modello sia stato validato, c’è da considerare che è stato ottenuto a partire da un modello iniziale in cui tutte le variabili erano state prese in considerazione: tra queste variabili vi erano anche variabili intrinsecamente simili alla variabile risposta (lunghezza del neonato e diametro del cranio , ovvero altre misuree corporee) e altre variabili non direttamente coinvolte, a rigor di logica, nel peso del neonato alla nascita, come l’ospedale in cui il parto avviene. Per questi motivi ho deciso di ripetere alcuni passaggi di questo progetto ripartendo questa volta da un modello in cui come regressori scelgo di tenere, nel primo step, solo le cause direttamente correlate al peso del neonato.
Creazione del Modello di Regressione Verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti
#Proviamo un modello alternativo con sole variabili biologiche
mod2 <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Tipo.parto + Sesso, data=neonati_puliti)
summary(mod2)
##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Tipo.parto + Sesso, data = neonati_puliti)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1471.24 -272.64 -13.21 264.17 1896.73
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3306.623 187.157 -17.668 < 2e-16 ***
## Anni.madre 3.794 1.705 2.224 0.02621 *
## N.gravidanze 18.616 7.007 2.657 0.00794 **
## Fumatrici -111.370 41.491 -2.684 0.00732 **
## Gestazione 163.592 4.519 36.201 < 2e-16 ***
## Tipo.partoNat 15.896 18.208 0.873 0.38274
## SessoM 164.945 16.698 9.878 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 413.4 on 2493 degrees of freedom
## Multiple R-squared: 0.3816, Adjusted R-squared: 0.3801
## F-statistic: 256.4 on 6 and 2493 DF, p-value: < 2.2e-16
AIC(mod2)
## [1] 37225.46
BIC(mod2)
## [1] 37272.05
modello2_stepwise_AIC <- step(mod2)
## Start: AIC=30128.76
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Tipo.parto +
## Sesso
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 130235 426131021 30128
## <none> 426000786 30129
## - Anni.madre 1 845501 426846287 30132
## - N.gravidanze 1 1206111 427206897 30134
## - Fumatrici 1 1231187 427231973 30134
## - Sesso 1 16674239 442675025 30223
## - Gestazione 1 223933725 649934511 31183
##
## Step: AIC=30127.53
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 426131021 30128
## - Anni.madre 1 846924 426977945 30130
## - N.gravidanze 1 1188774 427319795 30132
## - Fumatrici 1 1215268 427346289 30133
## - Sesso 1 16667830 442798852 30222
## - Gestazione 1 223815017 649946038 31181
summary(modello2_stepwise_AIC)
##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Sesso, data = neonati_puliti)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1466.6 -271.3 -12.0 261.1 1901.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3292.719 186.469 -17.658 <2e-16 ***
## Anni.madre 3.797 1.705 2.226 0.0261 *
## N.gravidanze 18.477 7.005 2.638 0.0084 **
## Fumatrici -110.625 41.480 -2.667 0.0077 **
## Gestazione 163.525 4.518 36.193 <2e-16 ***
## SessoM 164.913 16.697 9.877 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 413.4 on 2494 degrees of freedom
## Multiple R-squared: 0.3814, Adjusted R-squared: 0.3802
## F-statistic: 307.6 on 5 and 2494 DF, p-value: < 2.2e-16
Il processo stepwise ha creato un modello migliorato eliminando la variabile Tipo.parto (ha migliorato il valore dell’AIC), in quanto non ritenuto significativo per il peso alla nascita. Il modello finale è dunque composto da queste variabili :Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Sesso. Le variabili rimanenti nel modello sono tutte significative e hanno un legame causale con il peso del neonato. - Anni.madre (Età della madre): Per ogni anno in più di età della madre, il peso del neonato aumenta in media di circa 3.8 grammi. - N.gravidanze (Numero di gravidanze): Per ogni gravidanza in più, il peso del neonato aumenta in media di circa 18.5 grammi. * Fumatrici: I neonati di madri fumatrici pesano in media circa 110.6 grammi in meno rispetto ai neonati di madri non fumatrici, a parità di tutte le altre condizioni. - Gestazione: Per ogni settimana in più di gestazione, il peso del neonato aumenta di ben 163.5 grammi. - Sesso: I neonati maschi pesano in media circa 164.9 grammi in più rispetto alle femmine, a parità di tutte le altre condizioni.
Rimuovendo le variabili che misurano direttamente il neonato (Cranio-Lunghezza) il modello ora rileva come significativa la variabile Fumatrici. Analisi della qualità del modello
library(Metrics)
#Calcoliamo le previsioni del modello
previsioni2 <- predict(modello2_stepwise_AIC,neonati_puliti)
#calcoliamo l'RMSE
rmse(neonati_puliti$Peso,previsioni2)
## [1] 412.8588
library(lmtest)
# Test di normalità dei residui
shapiro.test(residuals(modello2_stepwise_AIC))
##
## Shapiro-Wilk normality test
##
## data: residuals(modello2_stepwise_AIC)
## W = 0.99693, p-value = 6.3e-05
# Test di omoschedasticità (Breusch-Pagan)
bptest(modello2_stepwise_AIC)
##
## studentized Breusch-Pagan test
##
## data: modello2_stepwise_AIC
## BP = 8.8749, df = 5, p-value = 0.1142
# Test di indipendenza (Durbin-Watson)
dwtest(modello2_stepwise_AIC)
##
## Durbin-Watson test
##
## data: modello2_stepwise_AIC
## DW = 1.8934, p-value = 0.00383
## alternative hypothesis: true autocorrelation is greater than 0
L’RMSE ha un valore di 412.86g. I residui non sono distribuiti in modo normale ed è presente eteroschedasticità, sono inoltre correlati tra loro.
# Calcolo il VIF per il modello finale alternativo
vif(modello2_stepwise_AIC)
## Anni.madre N.gravidanze Fumatrici Gestazione Sesso
## 1.182927 1.176920 1.003719 1.042552 1.019768
I valori di multicollinearità ottenuti, tutti inferiori a 5, indicano che le variabili sono indipendenti tra loro.
In conclusione, ci troviamo davanti a due modelli ciascuno con i propri punti di forza (un R quadro più alto e un RMSE più basso per il primo modello, dei residui non eteroschedastici ma un R quadro basso nel secondo modello). Non ritengo che uno dei due modelli sia migliore in assoluto. Ho ottenuto due modelli di regressione lineare multipla. Il modello completo spiega circa il 74% della variabilità del peso, mentre con il processo stepwise abbiamo un R² leggermente inferiore ma un RMSE più basso. Le settimane di gestazione e il sesso risultano predatori molto rilevanti mentre il fumo emerge come significativo solo nel secondo modello. Entrambi i modelli mostrano una buona capacità predittiva tuttavia si riscontrano problemi di non normalità e eteroschedasticità.