Progetto Statistica Inferenziale
Saverio
2024-10-30
1. RACCOLTA DEI DATI E STRUTTURA DEL DATASET
Importiamo le librerie necessarie
library(ggplot2)
library(dplyr)##
## Caricamento pacchetto: 'dplyr'## I seguenti oggetti sono mascherati da 'package:stats':
##
## filter, lag## I seguenti oggetti sono mascherati da 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(moments)Importiamo il 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 |
Analisi iniziale del dataset
Il dataset riporta 2500 osservazioni effettuate sulle nascite
avvenute in 3 distinti ospedali.
L’analisi è condotta attraverso l’osservazione e lo studio di 10
variabili:
- Anni.madre -> Età della madre in anni
- N.gravidanze -> Numero di gravidanze avute dalla madre
- Fumatrici -> Fumo materno: Un indicatore binario (0=non fumatrice,
1=fumatrice).
- Gestazione -> Durata della gravidanza in settimane di
gestazione
- Peso -> Peso del neonato in grammi
- Lunghezza; Cranio -> Lunghezza e diametro del cranio del
neonato
- Tipo.parto -> Tipo di parto: Naturale o cesareo.
- Ospedale -> Ospedale di nascita: Ospedale 1, 2 o 3.
- Sesso -> Sesso del neonato: Maschio (M) o femmina (F).
Queste variabili possono essere così categorizzate:
- Quantitative discrete: N.gravidanze.
- Quantitative continue: Peso, Lunghezza, Cranio, Anni.madre,
Gestazione
- Qualitative su scala nominale: Tipo.parto, Ospedale, Sesso,
Fumatrici
L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita: perciò la variabile Peso rappresenta la VARIABILE RISPOSTA, mentre le altre variabili sono VARIABILI ESPLICATIVE.
2. ANALISI E MODELLIZZAZIONE
Analisi preliminare
Analisi della variabile: Anni madre
summary(Anni.madre)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 25.00 28.00 28.16 32.00 46.00table(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 1dati_anni.madre <- data.frame(Asimmetria=round(skewness(Anni.madre),2), Curtosi=round(kurtosis(Anni.madre)-3,2), Dev_standard=round(sd(Anni.madre),2))
knitr::kable((dati_anni.madre), "pipe")| Asimmetria | Curtosi | Dev_standard |
|---|---|---|
| 0.04 | 0.38 | 5.27 |
Dal table riusciamo a notare subito un probabile errore nell’inserimento dei dati, poiché sono state inserite 2 rilevazioni che lasciano presupporre che esistano 2 madri di rispettivamente 0 e 1 anni.
Tuttavia queste rilevazioni non influiscono di molto sulla media avendo valore 28,16 anni non si discosta di molto dalla mediana uguale a 28 anni.
Con una assimmetria uguale a 0.04 abbiamo una distribuzione praticamente simmestrica.
plot(density(Anni.madre))
Analisi della variabile: numero di gravidanze
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 1dati_gravidanze <- data.frame(Asimmetria=round(skewness(N.gravidanze),2),
Curtosi=round(kurtosis(N.gravidanze)-3,2),
Dev_standard=round(sd(N.gravidanze),2))
knitr::kable((dati_gravidanze), "pipe")| Asimmetria | Curtosi | Dev_standard |
|---|---|---|
| 2.51 | 10.99 | 1.28 |
Calcoliamo frequenza, frequenza cumulata e percentuale cumulata gravidanze
frequenza_cumulata <- dati %>%
group_by(N.gravidanze) %>%
summarise(Frequenza = n()) %>%
mutate(
Freq_cumulata = cumsum(Frequenza),
Perc_cumulata = round(cumsum(Frequenza) / sum(Frequenza) * 100, 2)
)
knitr::kable(frequenza_cumulata, "pipe")| N.gravidanze | Frequenza | Freq_cumulata | Perc_cumulata |
|---|---|---|---|
| 0 | 1096 | 1096 | 43.84 |
| 1 | 818 | 1914 | 76.56 |
| 2 | 340 | 2254 | 90.16 |
| 3 | 150 | 2404 | 96.16 |
| 4 | 48 | 2452 | 98.08 |
| 5 | 21 | 2473 | 98.92 |
| 6 | 11 | 2484 | 99.36 |
| 7 | 1 | 2485 | 99.40 |
| 8 | 8 | 2493 | 99.72 |
| 9 | 2 | 2495 | 99.80 |
| 10 | 3 | 2498 | 99.92 |
| 11 | 1 | 2499 | 99.96 |
| 12 | 1 | 2500 | 100.00 |
Dalle frequenze cumulate si capisce chiaramente come la maggior parte
delle donne (il 76.56%) ha massimo un figlio o nessuno.
La media è infatti dello 0.9812, la mediana è pari a 1, la moda è
0.
Il tutto è confermato da un’assimmetria distributiva della variabile
positiva pari a 2.51: il che significa che sono più frequenti dei valori
bassi.
Analisi della 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 2dati_gestazione <- data.frame(Asimmetria=round(skewness(Gestazione),2), Curtosi=round(kurtosis(Gestazione)-3,2), Dev_standard=round(sd(Gestazione),2))
knitr::kable((dati_gestazione), "pipe")| Asimmetria | Curtosi | Dev_standard |
|---|---|---|
| -2.07 | 8.26 | 1.87 |
Le settimane di gestazione osservate vanno un minimo di 25 a un
massimo di 43 settimane.
La media è di 38.98, la mediana è pari a 39.
L’asimmetria negativa di -2.07 indica che la distribuzione è asimmetrica negativamente e che sono quindi più frequenti valori alti.
plot(density(Gestazione))
Analisi della variabile: peso
summary(Peso)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 2990 3300 3284 3620 4930table(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 1dati_peso <- data.frame(Asimmetria=round(skewness(Peso),2),
Curtosi=round(kurtosis(Peso)-3,2),
Dev_standard=round(sd(Peso),2))
knitr::kable((dati_peso), "pipe")| Asimmetria | Curtosi | Dev_standard |
|---|---|---|
| -0.65 | 2.03 | 525.04 |
Le osservazioni della variabile Peso vanno da un minimo di 830g e un
massimo di 4930g.
La media calcolata è di 3284g.
La mediana invece è pari a 3300g.
L’asimmetria negativa è pari a -0.65, ciò signfica che la distribuzione è asimmetrica negativamente, e che quindi sono più frequenti valori di peso alti.
plot(density(Peso))
Analisi della variabile: Lunghezza
summary(Lunghezza)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 310.0 480.0 500.0 494.7 510.0 565.0table(Lunghezza)## Lunghezza
## 310 315 320 325 340 345 355 360 370 380 385 390 400 405 410 420 425 430 435 440
## 1 1 1 1 1 1 2 2 3 3 1 5 4 3 8 8 1 10 2 13
## 445 446 450 455 460 465 470 473 475 480 485 490 492 494 495 497 498 500 502 504
## 13 2 43 16 88 32 126 1 66 195 122 270 1 1 150 1 2 396 1 1
## 505 510 514 515 518 520 523 525 530 535 540 545 550 555 560 565
## 117 273 1 94 1 172 1 57 89 24 32 12 23 3 2 1dati_lunghezza <- data.frame(Asimmetria=round(skewness(Lunghezza),2), Curtosi=round(kurtosis(Lunghezza)-3,2), Dev_standard=round(sd(Lunghezza),2))
knitr::kable((dati_lunghezza), "pipe")| Asimmetria | Curtosi | Dev_standard |
|---|---|---|
| -1.51 | 6.49 | 26.32 |
Le osservazioni della variabile lunghezza vanno da un minimo di 310mm
a un massimo di 565 mm.
La media è pari a 494.7mm, mentre la mediana è pari a 500mm.
L’asimmetria negativa è pari a -1.51, ciò significa che sono più frequenti valori alti.
plot(density(Lunghezza))
Analisi della variabile: tipo di parto
table(Tipo.parto)## Tipo.parto
## Ces Nat
## 728 1772Dalle osservazioni si evidenzia un numero molto maggiore dei parti naturali rispetto ai parti casarei.
Analisi della variabile: ospedale
table(Ospedale)## Ospedale
## osp1 osp2 osp3
## 816 849 835Dalle osservazioni non emergono particolari differenze di numerosità di parti da un ospedale all’altro.
Analisi della variabile: sesso
table(Sesso)## Sesso
## F M
## 1256 1244Le osservazioni di nascite di bambini maschi e femmine non differiscono particolarmente tra di loro.
Analisi della variabile: fumatrici
table(Fumatrici)## Fumatrici
## 0 1
## 2396 104Circa il 4% delle gestanti è fumatrice (ne deriva ovviamente che il 96% non è fumatrice). *Questo dato andrebbe confrontato con statistiche più accurate per accertarsi del fatto che il campione in esame coincida o meno con le statistiche sulla popolazione, e se quindi il 4% di questo campione non stia sottostimando la numerosità delle donne in gestazione fumatrici a livello di popolazione.
Studio della relazione tra le variabili materne (età, fumo, gravidanze precedenti) e il peso del neonato
Relazione età madre - peso neonato
plot(Anni.madre, Peso, pch=20, xlab='Anni della madre', ylab='Peso del neonato')
Tralasciando gli outlier dovuti all’errore di inserimento nel dataset dell’età della madre, possiamo notare visimanete una maggiore concentrazione attorno ai valori alti di peso (da 2500 a 4000), ma non notiamo alcuna rilevante relazione tra le variabili in esame.
Relazione madre fumatrice - peso neonato
dati$Fumatrici <- factor(dati$Fumatrici, labels = c("Non Fumatrice", "Fumatrice"))
library(ggplot2)
ggplot(dati, aes(x = Fumatrici, y = Peso)) +
geom_boxplot(aes(fill = Fumatrici)) +
labs(title = "Impatto della Madre Fumatrice sul Peso del Bambino",
x = "Fumatrici",
y = "Peso del Bambino (g)") +
theme_minimal() +
scale_fill_manual(values = c("Non Fumatrice" = "skyblue", "Fumatrice" = "salmon"))
Per quanto la rilevanza vada lasciata a giudizio di personale medico
o comunque specializzato, si può subito notare come la variabile “Madre
fumatrice” incida sul peso mediano del neonato.
Viene infatti evidenziato come i neonati con madre fumatrice abbiano un
peso mediano leggermente inferiore rispetto ai neonati con madre non
fumatrice.
Relazione Numero Gravidanze precedenti - Peso del neonato
ggplot(dati, aes(x = factor(N.gravidanze), y = Peso)) +
geom_boxplot(aes(fill = factor(N.gravidanze))) +
labs(title = "Distribuzione del Peso del Neonato in base al Numero di Gravidanze Precedenti",
x = "Numero di Gravidanze Precedenti",
y = "Peso del Neonato (g)") +
theme_minimal() +
scale_fill_brewer(palette = "Pastel1")## Warning in RColorBrewer::brewer.pal(n, pal): n too large, allowed maximum for palette Pastel1 is 9
## Returning the palette you asked for with that many colors
Si può notare come a prescindere dal numero di gravidanze il peso
mediano dei neonati non differisca in maniera rilevante.
Sono comunque da evidenziare i casi con almeno 6 gravidanze precedenti,
i quali potrebbero (analizzando visivamente il grafico) incidere sul
peso medio del neonato.
Creazione del modello di regressione
Per prima cosa studiamo la correlazione tra le VARIABILI QUANTITATIVE
quantitative ed il peso dei neonati.
Procediamo con la creazione di una matrice di correlazione.
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
La matrice di regressione ci dice chiaramente che:
- Le settimane di gestazione sono correlate positivamente con il peso e
che all’aumentare delle settimane il peso del neonato aumenta. In
particolare per ogni settimana il peso medio aumenta di 0.59
- La lunghezza e il diamentro del cranio sono correlati positivamente
con il numero di settimane di gestazione e quindi direttamente con il
peso del neonato. In particolare un aumento unitario della lunghezza e
del diametro del cranio porta un aumento medio del peso di,
rispettivamente, 0.8 e 0.7
- Le altre osservazioni non influenzano particolarmente il peso del
neonato
Per le VARIABILI QUALITATIVE invece:
dati$Fumatrici <- factor(dati$Fumatrici, labels = c("Non Fumatrice", "Fumatrice"))
library(ggplot2)
ggplot(dati, aes(x = Fumatrici, y = Peso)) +
geom_boxplot(aes(fill = Fumatrici)) +
labs(title = "Impatto della Madre Fumatrice sul Peso del Bambino",
x = "Fumatrici",
y = "Peso del Bambino (g)") +
theme_minimal() +
scale_fill_manual(values = c("Non Fumatrice" = "skyblue", "Fumatrice" = "salmon"))
dati$Fumatrici_dataframe<-factor(ifelse(Fumatrici==1,'S','N'))
t.test(Peso~dati$Fumatrici_dataframe)##
## Welch Two Sample t-test
##
## data: Peso by dati$Fumatrici_dataframe
## 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.346Dal boxplot elaborato si evidenzia come l’avere o meno una madre
fumatrice incida sul peso del neonato.
Infatti si può notare un peso mediano minore nei neonati con madre
fumatrice.
Tuttavia, un p-value superiore a 0.05 indica che non ci sono prove
sufficienti per rifiutare l’ipotesi nulla, ed in questo caso, con un
p-value di 0.3033, si afferma che non c’è differenza significativa nei
pesi tra i neonati, a prescindere dal fatto le che madri fumino o
meno.
Procediamo con la creazione del modello di regressione lineare multipla e con la selezione del modello ottimale
Si procede con la creazine di un modello di regressione lineare multipla che metta in relazione tutte le variabili con il peso del neonato.
Tuttavia procediamo si da subito ad escludere dal modello la varibile “Ospedale” e la variabile “Tipo di parto”, ritenendo a seguito di opportuna analisi, che non inficino in alcun modo sul peso dei neonati.
boxplot(Peso ~ Ospedale, data = dati,
main = "Peso dei Neonati per Ospedale",
xlab = "Ospedale",
ylab = "Peso (grammi)",
col = "lightblue")
boxplot(Peso ~ Tipo.parto, data = dati,
main = "Peso dei Neonati per Ospedale",
xlab = "Tipo di parto",
ylab = "Peso (grammi)",
col = "lightblue")
Si procede…
mod1<-lm(Peso~ Anni.madre + N.gravidanze + Fumatrici_dataframe + Gestazione + Lunghezza + Cranio + Sesso , data=dati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici_dataframe +
## Gestazione + Lunghezza + Cranio + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1161.56 -181.19 -15.75 163.70 2630.75
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6714.4109 141.1515 -47.569 < 2e-16 ***
## Anni.madre 0.9585 1.1347 0.845 0.3984
## N.gravidanze 11.2756 4.6690 2.415 0.0158 *
## Fumatrici_dataframeS -30.2959 27.5971 -1.098 0.2724
## Gestazione 32.9331 3.8267 8.606 < 2e-16 ***
## Lunghezza 10.2342 0.3009 34.009 < 2e-16 ***
## Cranio 10.5177 0.4268 24.642 < 2e-16 ***
## SessoM 78.0845 11.2039 6.969 4.06e-12 ***
## ---
## 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-16Il modello presenta un “Multiple R-squared” pari a 0.72. Questo è di per se un buon risultato, poiché significa che il 72% della variabilità della variabile Peso dipende dalle variabili sottoposte a modello.
Per migliorarlo possiamo provare a togliere un’altra variabile dal modello. Quella che si ritiene meno incisiva rispetto alle altre è la variabile relativ all’età della madre.
mod2<-update(mod1,~.-Anni.madre)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici_dataframe + 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_dataframeS -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-16Il modello presenta, anche senza la variabile “Anni.madre”, un
“Multiple R-squared” pari a 0.72.
Questo significa che la variabile “Anni.madre” non è significativa nella
variabilità della variabile “Peso”.
Adesso, basandoci sugli studi fatti precedentemente, potremmo pensare di escludere anche la variabile “Fumatrici”.
mod3<-update(mod2,~.-Fumatrici_dataframe)
summary(mod3)##
## 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 senza la variabile “Fumatrici”, il modello presenta un
“Multiple R-squared” pari a 0.72.
Questo significa che la variabile “Fumatrici” non è significativa nella
variabilità della variabile “Peso”.
Rimangono le variabili: Numero di gravidanze, Settimane di gestazione, Lunghezza e diametro del cranio e Sesso.
Procediamo con l’escludere la variabile Numero di gravidanze.
mod4<-update(mod3,~.-N.gravidanze)
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-16Anche senza la variabile “N.Gravidanze”, il modello presenta un
“Multiple R-squared” pari a 0.72.
Questo significa che la variabile “N.Gravidanze” non è significativa
nella variabilità della variabile “Peso”.
Rimangono nel modello di regressione le variabili: Gestazione, Sesso, Lunghezza e diametro del cranio.
Queste 4 variabili determinano il 72% della variabilità della variabile Peso.
Possiamo continuare con la sperimentazione andando a escludere la variabile “Sesso”.
mod5<-update(mod4,~.-Sesso)
summary(mod5)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1105.77 -183.25 -12.83 166.41 2623.80
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6777.1203 135.6417 -49.963 <2e-16 ***
## Gestazione 31.6901 3.8219 8.292 <2e-16 ***
## Lunghezza 10.4252 0.3020 34.522 <2e-16 ***
## Cranio 10.7892 0.4282 25.194 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 277.7 on 2496 degrees of freedom
## Multiple R-squared: 0.7206, Adjusted R-squared: 0.7203
## F-statistic: 2146 on 3 and 2496 DF, p-value: < 2.2e-16Anche adesso (sorprendentemente), il modello presenta un “Multiple
R-squared” pari a 0.72.
Questo significa che la variabile “Sesso” non è significativa nella
variabilità della variabile “Peso”.
Le sole 3 variabili: Gestazione, Lunghezza e Cranio definiscono il 72% della variabilità della variabile Peso.
Per conclusione proviamo delle ultime combinazioni per poi procedere con test definitivo per la scelta del modello migliore.
Escludiamo la variabile Cranio…
mod6<-update(mod5,~.-Cranio)
summary(mod6)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1152.8 -199.6 -20.4 192.0 3627.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5349.8835 138.0018 -38.77 <2e-16 ***
## Gestazione 45.1267 4.2377 10.65 <2e-16 ***
## Lunghezza 13.8973 0.3009 46.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 310.9 on 2497 degrees of freedom
## Multiple R-squared: 0.6496, Adjusted R-squared: 0.6493
## F-statistic: 2314 on 2 and 2497 DF, p-value: < 2.2e-16La “Multiple R-squared” è scesa del 7% (0.65) indicando un significativo apporto della variabile Cranio alla variabilità della variabile Peso.
Proviamo a tornare indietro ed escludere la variabile “Lunghezza”.
mod7<-update(mod5,~.-Lunghezza)
summary(mod7)##
## Call:
## lm(formula = Peso ~ Gestazione + Cranio, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1361.72 -225.99 -8.08 224.47 1469.20
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6391.1327 164.2808 -38.90 <2e-16 ***
## Gestazione 95.2383 4.0705 23.40 <2e-16 ***
## Cranio 17.5361 0.4631 37.87 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 337.5 on 2497 degrees of freedom
## Multiple R-squared: 0.5872, Adjusted R-squared: 0.5869
## F-statistic: 1776 on 2 and 2497 DF, p-value: < 2.2e-16Senza la variabile “Lunghezza” la “Multiple R-squared” scesa addirittura sino al 58.7%, indicando una forte significatività nella variabilità della variabile “Peso”.
In conclusione, proviamo ad escludere la variabile “Gestazione”.
mod8<-update(mod5,~.-Gestazione)
summary(mod8)##
## Call:
## lm(formula = Peso ~ Lunghezza + Cranio, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1295.09 -184.36 -11.95 162.85 2792.61
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6306.9134 124.8791 -50.50 <2e-16 ***
## Lunghezza 11.6312 0.2682 43.37 <2e-16 ***
## Cranio 11.2847 0.4298 26.26 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 281.4 on 2497 degrees of freedom
## Multiple R-squared: 0.7129, Adjusted R-squared: 0.7127
## F-statistic: 3101 on 2 and 2497 DF, p-value: < 2.2e-16Anche la variabile “Gestazione” influisce significativamente sulla variabile “Peso”. In particolare, la “Multiple R-squared” è scesa al 0.71. Molto meno rispetto all’influenza delle variabili Lunghezza e Cranio, le quali da sole influiscono per il 71% sulla variabile “Peso.
Efffettuiamo adesso un test in base al criterio di informazione di Bayes (BIC):
BIC(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8)## df BIC
## mod1 9 35233.81
## mod2 8 35226.70
## mod3 7 35220.10
## mod4 6 35220.54
## mod5 5 35262.11
## mod6 4 35820.73
## mod7 4 36230.13
## mod8 4 35322.22Il modello con il BIC più basso è quello da preferire. Ossia il modello n. 3, il quale tiene conto delle variabili: N.gravidanze, Gestazione, Lunghezza, Cranio, Sesso
Funzione stepAIC
Adesso proviamo a vedere quale modello ci consiglia il software
utilizzando la funzione stepAIC.
Per farlo inseriamo un modello iniziale che tiene conto di tutte le
variabili presenti nel dataset.
n<-nrow(dati)
stepwise.mod<-MASS::stepAIC(lm(Peso~ Anni.madre + N.gravidanze + Fumatrici_dataframe + Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data=dati),
direction = "both",
k=log(n))## Start: AIC=28139.32
## Peso ~ Anni.madre + N.gravidanze + Fumatrici_dataframe + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 46578 186809099 28132
## - Fumatrici_dataframe 1 90019 186852540 28133
## - Ospedale 2 685979 187448501 28133
## - N.gravidanze 1 438452 187200974 28137
## - Tipo.parto 1 447929 187210450 28138
## <none> 186762521 28139
## - Sesso 1 3611021 190373542 28179
## - Gestazione 1 5458403 192220925 28204
## - Cranio 1 45326172 232088693 28675
## - Lunghezza 1 87951062 274713583 29096
##
## Step: AIC=28132.12
## Peso ~ N.gravidanze + Fumatrici_dataframe + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici_dataframe 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 46578 186762521 28139
## - 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_dataframe 1 90897 186809099 28132
## + Anni.madre 1 47456 186852540 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_dataframe 1 99840 187501837 28126
## + Anni.madre 1 54392 187547285 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_dataframe 1 91892 187973654 28124
## + Anni.madre 1 54816 188010731 28125
## - Sesso 1 3655292 191720838 28158
## - Gestazione 1 5464853 193530399 28181
## - Cranio 1 46108583 234174130 28658
## - Lunghezza 1 87632762 275698308 29066La funzione ci suggerisci il modello che tiene conto delle variabili: N.gravidanze, Lunghezza, Cranio, Sesso e Gestazione.
Il modello che ci suggerisce il software è uguale al nostro modello n.3, suggerito anche dal criterio di informazione di Bayes (BIC) applicato in precedenza.
Il MODELLO N.3 è il nostro MODELLO OTTIMALE.
Analisi della Qualità del Modello
Il valore “Multiple R-squared” è già stato precedentemente valutato attraverso la comparazione dei vari modelli.
Un Modello di regressione deve soddisfare le seguenti assunzioni: -
linearità
- normalità
- omoschedasticità
- indipendenza
Analisi dei residui
par(mfrow=c(2,2))
plot(mod3)
#### RESIDUAL VS FITTED: valutare l’ipotesi di linearità
I punti dovrebbero essere distribuiti casualmente attorno alla linea
orizzontale a 0.
Ed infatti, nel nostro caso i residui sembrano essere distribuiti
simmetricamente attorno allo zero.
Q-Q Residuals: i residui seguono una distribuzione normale?
I punti dovrebbero essere allineati lungo la linea diagonale.
Nel nostro caso, per la maggior parte i punti seguono la bisettrice,
eccetto nelle code.
shapiro.test(mod3$residuals)##
## Shapiro-Wilk normality test
##
## data: mod3$residuals
## W = 0.97408, p-value < 2.2e-16Secondo il test di Shapiro-Wilk dovremmo rifiutare l’ipotesi normalità della distribuzione.
Accertiamocene con un grafico visivo.
plot(density(residuals(mod3)))+
lines(density(rnorm(100000,0,250)),col=2)
## integer(0)Visivamente possiamo notare che la distribuzione dei residui segue una una distribuzione normale.
Quindi viene soddisfatta l’assunzione di normalità.
Scale-Location: verificare la costanza della varianza dei residui, ossia la omoschedasticità
I punti dovrebbero essere distribuiti orizzontalmente senza uno schema specifico.
Effettuiamo un test.
lmtest::bptest(mod3)##
## studentized Breusch-Pagan test
##
## data: mod3
## BP = 90.253, df = 5, p-value < 2.2e-16Il test di Breusch-Pagan porta a rifiutare l’ipotesi di omoschedasticità, suggerendo quindi eteroschedasticità.
Residuals vs Leverage: ci sono valori influenti?
C’è un unico valore che supera la soglia di Cook (1551).
Si rendono necessari ulteriori test statistici
Analizziamo gli outlier
plot(rstudent(mod3))
abline(h=c(-2,2),col=2)
car::outlierTest(mod3)## rstudent unadjusted p-value Bonferroni p
## 1551 10.051908 2.4906e-23 6.2265e-20
## 155 5.027798 5.3138e-07 1.3285e-03
## 1306 4.827238 1.4681e-06 3.6702e-03Ci sono 3 outlier: 1551, 155, 1306.
Calcoliamo la distanza di Cook
cook<-cooks.distance(mod3)
plot(cook) 
mydata <- data.frame(Max.dist.di.Cook=round(max(cook),2))
knitr::kable((mydata), "pipe")| Max.dist.di.Cook |
|---|
| 0.83 |
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 | Non Fumatrice | 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
Proviamo a stimare il peso di:
- una neonata (F)
- con madre non fumatrice
- alla sua terza gravidanza
- che partorirà alla 39esima settimana.
Per i dati del modello bisogna inserire anche la lunghezza e il cranio (useremo i valori medi).
previsione=predict(mod3,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 |
|---|---|---|---|---|---|---|
| 3271.08 | F | 3 | 39 | 0 | 494.69 | 340.03 |
Il peso stimato della neonata è pari a 3271.08.
4. VISUALIZZAZIONI
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 del tempo di gestazione e del Sesso')+
theme_classic()## `geom_smooth()` using formula = 'y ~ x'
La nuvola di punti si addensa a partire dalla 36esima settimana di
gestazione circa.
Le rette relative ai due sessi hanno la stessa pendenza, ma altezza
diversa poiché i maschi hanno un peso medio maggiore delle femmine.
Impatto del numero di settimane di gestazione e del fumo della madre sul peso del neonato
ggplot(data=dati)+
geom_point(aes(x=Gestazione,
y=Peso,
col=Fumatrici_dataframe),position='jitter')+
geom_smooth(aes(x=Gestazione,
y=Peso,
col=Fumatrici_dataframe),
se=F, #standardError=False
method = "lm")+
labs(title='Andamento del peso al variare del tempo di gestazione con/senza madri fumatrici',col='Fumatrici')+
theme_classic()## `geom_smooth()` using formula = 'y ~ x'
La scarsa presenza di osservazioni relative a madre fumatrici non
permette di effettuare un confronto efficace. In linea di massima sembra
che il peso medio si aggiri attorno ai 3000g.
Impatto del numero di gravidanze e del sesso sul peso pdel 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 sesso',x='Numero di gravidanze')+
theme_classic()## `geom_smooth()` using formula = 'y ~ x'
Il peso medio dei neonati risulta maggiore di quello delle neonate, e rimane pressochè costante all’aumentare del numero delle gravidanze della gestante.
5. RESOCONTO STUDIO
La variabile Fumatrici, sorprendentemente, non è risultata di importanza significativa nello studio di previsione del peso dei neonati.
Al contrario sono risultate molto più significative le variabili relative al numero di gravidanze precedenti e al tempo di gestazione.
Particolarmente influenti nel modello sono le variabili lunghezza e diametro del cranio.
Il modello costruito presentava un valore R2 del 72% ed è quindi ritenuto un buon modello di previsione per il peso dei neonati.
Inizialmente il modello sembrava non rispettasse le assunzioni necessarie, tuttavia si è constatato (anche visivamente) come in realtà il modello andasse bene.
Nelle analisi è risultato esserci qualche errore di inserimento dei dati, poiché i valori non risultavano coerenti fra loro ma si è riusciti comunque a concludere lo studio efficacemente.