Modello Statistico per la Previsione del Peso Neonatale
Eduardo De Filippis
0 - Preparazione dell’ambiente di lavoro
In questa sezione vengono richiamati i packages utili all’analisi da
svolgere.
options(repos = c(CRAN = "https://cran.rstudio.com/")) #per creazione del documento con knit
library(ggplot2)
library(carData)
library(car) #per diagnostica del modello
library(moments) #per calcolo di asimmetria e curtosi
library(gridExtra) #per personalizzazione della finestra output
library(knitr) #per formattare meglio alcune tabelle
library(kableExtra) #per formattare meglio alcune tabelle1 - Preparazione dei dati
In questa sezione si procede dapprima all’importazione del dataset.
Poi vengono manipolate alcune variabili qualitative per una codifica più
appropriata, vengono quindi estratte le singole colonne dal dataset,
nonchè mostrate le prime 10 osservazioni e la struttura generale dei
dati. Infine vengono descritte le singole variabili esplicitando
l’eventuale codifica.
neonati <- read.csv("neonati.csv", header = TRUE, stringsAsFactors = FALSE) #import dati#codifica e conversione variabili
neonati$Fumatrici <- factor(neonati$Fumatrici) #già codificata
neonati$Tipo.parto[neonati$Tipo.parto == "Nat"] <- "0"
neonati$Tipo.parto[neonati$Tipo.parto == "Ces"] <- "1"
neonati$Tipo.parto <- factor(neonati$Tipo.parto)
neonati$Ospedale[neonati$Ospedale == "osp1"] <- "1"
neonati$Ospedale[neonati$Ospedale == "osp2"] <- "2"
neonati$Ospedale[neonati$Ospedale == "osp3"] <- "3"
neonati$Ospedale <- factor(neonati$Ospedale)
neonati$Sesso[neonati$Sesso == "M"] <- "0"
neonati$Sesso[neonati$Sesso == "F"] <- "1"
neonati$Sesso <- factor(neonati$Sesso)
attach(neonati) #estrazione colonne#visualizzazione prime 10 righe
kable(head(neonati,10), align = rep('c', 10)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso |
|---|---|---|---|---|---|---|---|---|---|
| 26 | 0 | 0 | 42 | 3380 | 490 | 325 | 0 | 3 | 0 |
| 21 | 2 | 0 | 39 | 3150 | 490 | 345 | 0 | 1 | 1 |
| 34 | 3 | 0 | 38 | 3640 | 500 | 375 | 0 | 2 | 0 |
| 28 | 1 | 0 | 41 | 3690 | 515 | 365 | 0 | 2 | 0 |
| 20 | 0 | 0 | 38 | 3700 | 480 | 335 | 0 | 3 | 1 |
| 32 | 0 | 0 | 40 | 3200 | 495 | 340 | 0 | 2 | 1 |
| 26 | 1 | 0 | 39 | 3100 | 480 | 345 | 0 | 3 | 1 |
| 25 | 0 | 0 | 40 | 3580 | 510 | 349 | 0 | 1 | 0 |
| 22 | 1 | 0 | 40 | 3670 | 500 | 335 | 1 | 2 | 1 |
| 23 | 0 | 0 | 41 | 3700 | 510 | 362 | 1 | 2 | 1 |
str(neonati) #struttura generale## 'data.frame': 2500 obs. of 10 variables:
## $ Anni.madre : int 26 21 34 28 20 32 26 25 22 23 ...
## $ N.gravidanze: int 0 2 3 1 0 0 1 0 1 0 ...
## $ Fumatrici : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ Gestazione : int 42 39 38 41 38 40 39 40 40 41 ...
## $ Peso : int 3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
## $ Lunghezza : int 490 490 500 515 480 495 480 510 500 510 ...
## $ Cranio : int 325 345 375 365 335 340 345 349 335 362 ...
## $ Tipo.parto : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 2 2 ...
## $ Ospedale : Factor w/ 3 levels "1","2","3": 3 1 2 2 3 2 3 1 2 2 ...
## $ Sesso : Factor w/ 2 levels "0","1": 1 2 1 1 2 2 2 1 2 2 ...Alcune variabili sono a livelli che partono da 1, invece che da 0,
per impostazione di default quando sono state fatte conversioni in
fattori (ciò non influenza l’analisi e le variabili binarie saranno
comunque con codifica effettiva 1 e 0).
A dimostrazione della
considerazione precedente viene mostrata la prima riga del dataset in
cui, per esempio, secondo il comando str, il sesso è pari a 1 e invece
effettivamente è valorizzato a 0.
kable(neonati[1,], align = rep('c', 10)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso |
|---|---|---|---|---|---|---|---|---|---|
| 26 | 0 | 0 | 42 | 3380 | 490 | 325 | 0 | 3 | 0 |
Il dataset è composto da 2500 osservazioni e 10 variabili, le quali
sono:
- Anni.madre: Età in anni della madre. Variabile quantitativa
continua
- N.gravidanze: Numedi di gravidanze avuto dalla madre.
Variabile quantitativa discreta.
- Fumatrici: Qualifica di fumatrice
o non fumatrice per la madre (0=non fumatrice, 1=fumatrice). Variabile
qualitativa binaria.
- Gestazione: Numero di settimane di
gestazione, quindi durata della gravidanza. Variabile quantitativa
continua
- Peso: Peso alla nascita, in grammi, del neonato.
Variabile quantitativa continua.
- Lunghezza: Lunghezza del neonato.
Variabile quantitativa continua.
- Cranio: Diametro craniale.
Variabile quantitativa continua.
- Tipo.parto: Tipo di parto
(0=naturale, 1=cesareo). Variabile qualitativa binaria.
- Ospedale:
Ospedale di nascita (ospedale1=1, ospedale2=2, ospedale3=3). Variabile
qualitativa categoriale.
- Sesso: Sesso del neonato (maschio=0,
femmina=1). Variabile qualitativa binaria.
2 - Analisi esplorativa e descrittiva
In questa sezione viene proposta un’analisi esplorativa sui dati: statistiche descrittive, distribuzioni e boxplots.
2.1 - Statistiche descrittive generali
options(knitr.kable.NA = "") #celle vuote
kable(summary(neonati), align = rep('c', 10)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| Min. : 0.00 | Min. : 0.0000 | 0:2396 | Min. :25.00 | Min. : 830 | Min. :310.0 | Min. :235 | 0:1772 | 1:816 | 0:1244 | |
| 1st Qu.:25.00 | 1st Qu.: 0.0000 | 1: 104 | 1st Qu.:38.00 | 1st Qu.:2990 | 1st Qu.:480.0 | 1st Qu.:330 | 1: 728 | 2:849 | 1:1256 | |
| Median :28.00 | Median : 1.0000 | Median :39.00 | Median :3300 | Median :500.0 | Median :340 | 3:835 | ||||
| Mean :28.16 | Mean : 0.9812 | Mean :38.98 | Mean :3284 | Mean :494.7 | Mean :340 | |||||
| 3rd Qu.:32.00 | 3rd Qu.: 1.0000 | 3rd Qu.:40.00 | 3rd Qu.:3620 | 3rd Qu.:510.0 | 3rd Qu.:350 | |||||
| Max. :46.00 | Max. :12.0000 | Max. :43.00 | Max. :4930 | Max. :565.0 | Max. :390 |
Si nota un problema: Anni.madre ha come valore minimo 0. Ciò è molto strano quindi i valori della colonna vengono ordinati per individuare eventuali errori e vengono stampati i primi 20.
tmp<-sort(neonati$Anni.madre)
print(head(tmp, 20))## [1] 0 1 13 14 14 15 15 15 15 15 15 16 16 16 16 16 16 16 16 16Sono presenti due valori, impossibili da verificarsi, che sono pari a 0 e 1. Vengono eliminati dal vettore per calcolare la media delle restanti valorizzazioni, poi tale media viene imputata alle osservazioni del dataset con Anni.madre=0 e 1.
tmp<-tmp[3:length(tmp)]
mean(tmp)## [1] 28.18615neonati[neonati$Anni.madre %in% c(0, 1), ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1152 1 1 0 41 3250 490 350
## 1380 0 0 0 39 3060 490 330
## Tipo.parto Ospedale Sesso
## 1152 0 2 1
## 1380 0 3 0#imputazione con media nel dataset
neonati$Anni.madre[1380]=28.18615
neonati$Anni.madre[1152]=28.18615Viene riproposto un summary corretto del dataset.
kable(summary(neonati), align = rep('c', 10)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| Min. :13.00 | Min. : 0.0000 | 0:2396 | Min. :25.00 | Min. : 830 | Min. :310.0 | Min. :235 | 0:1772 | 1:816 | 0:1244 | |
| 1st Qu.:25.00 | 1st Qu.: 0.0000 | 1: 104 | 1st Qu.:38.00 | 1st Qu.:2990 | 1st Qu.:480.0 | 1st Qu.:330 | 1: 728 | 2:849 | 1:1256 | |
| Median :28.00 | Median : 1.0000 | Median :39.00 | Median :3300 | Median :500.0 | Median :340 | 3:835 | ||||
| Mean :28.19 | Mean : 0.9812 | Mean :38.98 | Mean :3284 | Mean :494.7 | Mean :340 | |||||
| 3rd Qu.:32.00 | 3rd Qu.: 1.0000 | 3rd Qu.:40.00 | 3rd Qu.:3620 | 3rd Qu.:510.0 | 3rd Qu.:350 | |||||
| Max. :46.00 | Max. :12.0000 | Max. :43.00 | Max. :4930 | Max. :565.0 | Max. :390 |
Si potrebbero discutere tanti aspetti delle statistiche descrittive.
Ciò che però si evidenzia in questa sede è che le variabili quantitative
presentano tutte delle medie che sono prossime alla mediana, quindi
probabilmente sono distribuzioni simmetriche. Ciononostante sarà
presentata un’analisi grafica e degli indici di forma per valutare
ciò.
Indicatori di variabilità sono omessi poichè saranno deducibili
dai modelli presentati in seguito.
Di seguito vengono presentate le
distribuzioni di frequenza per le variabili qualitative.
freq_ass <- table(Fumatrici)
freq_rel <- table(Fumatrici)/length(Fumatrici) # Tabella distribuzione frequenze relative
distr_freq_Fumatrici<-data.frame(cbind(freq_ass,freq_rel))
distr_freq_Fumatrici## freq_ass freq_rel
## 0 2396 0.9584
## 1 104 0.0416Se si considera la caratteristica della madre “fumatrice o non fumatrice” il dataset è totalmente sbilanciato, infatti oltre il 95% dei neonati ha una madre che non è fumatrice.
freq_ass <- table(Tipo.parto)
freq_rel <- table(Tipo.parto)/length(Tipo.parto)
distr_freq_Tipo.parto<-data.frame(cbind(freq_ass,freq_rel))
distr_freq_Tipo.parto## freq_ass freq_rel
## 0 1772 0.7088
## 1 728 0.2912Se si considera il tipo di parto il dataset risulta, anche in questo caso, altamente sbilanciato, infatti meno del 30% dei neonati è nato tramite parto cesareo.
freq_ass <- table(Ospedale)
freq_rel <- table(Ospedale)/length(Ospedale)
distr_freq_Ospedale<-data.frame(cbind(freq_ass,freq_rel))
distr_freq_Ospedale## freq_ass freq_rel
## 1 816 0.3264
## 2 849 0.3396
## 3 835 0.3340Se si considera l’ospedale di nascita a cui è associata la singola osservazione, il dataset risulta più o meno bilaciato. Infatti circa il 33% dei neonati è nato in ognuno degli ospedali considerati.
freq_ass <- table(Sesso)
freq_rel <- table(Sesso)/length(Sesso)
distr_freq_Sesso<-data.frame(cbind(freq_ass,freq_rel))
distr_freq_Sesso## freq_ass freq_rel
## 0 1244 0.4976
## 1 1256 0.5024Se si considera il sesso del neonato, il dataset è pressocchè bilanciato. Infatti quasi la metà è maschio e poco più della metà è femmina.
# La variabile numero gravidanze (quantitativa discreta) presenta poche modalità quindi per fini di esplorazione si stampano le frequenze
freq_ass <- table(N.gravidanze)
freq_rel <- table(N.gravidanze)/length(N.gravidanze)
distr_freq_ng<-data.frame(cbind(freq_ass,freq_rel))
kable(distr_freq_ng, align = rep('c', 3)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))| freq_ass | freq_rel | |
|---|---|---|
| 0 | 1096 | 0.4384 |
| 1 | 818 | 0.3272 |
| 2 | 340 | 0.1360 |
| 3 | 150 | 0.0600 |
| 4 | 48 | 0.0192 |
| 5 | 21 | 0.0084 |
| 6 | 11 | 0.0044 |
| 7 | 1 | 0.0004 |
| 8 | 8 | 0.0032 |
| 9 | 2 | 0.0008 |
| 10 | 3 | 0.0012 |
| 11 | 1 | 0.0004 |
| 12 | 1 | 0.0004 |
Oltre l’80% dei neonati ha una madre che ha avuto al più 3 gravidanze.
2.2 - Distribuzioni empiriche, stime del kernel e indici di asimmetria e curtosi
# Distribuzione Anni della madre
g1<-ggplot(neonati, aes(x = Anni.madre)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 1, fill = "steelblue", color = "black") +
geom_density(color = "red", size = 0.7) +
labs(x = "Anni.madre", y = "Frequenza")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# Distribuzione del Numero di gravidanze associato alla madre
g2<-ggplot(neonati, aes(x = N.gravidanze)) +
geom_histogram(aes(y = after_stat(count / sum(count))), binwidth = 0.3, fill = "steelblue", color = "black") +
labs(x = "N.gravidanzee", y = "Frequenza")
# Distribuzione delle settimane di gestazione
g3<-ggplot(neonati, aes(x = Gestazione)) +
geom_histogram(aes(y = after_stat(count / sum(count))), binwidth = 0.3, fill = "steelblue", color = "black") +
labs(x = "Gestazione", y = "Frequenza")
#Distribuzione del Peso del neonato (variabile risposta)
g4<-ggplot(neonati, aes(x = Peso)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 70, fill = "steelblue", color = "black") +
geom_density(color = "red", size = 0.7) +
labs(x = "Peso", y = "Frequenza")
#Distribuzione della lunghezza del neonato
g5<-ggplot(neonati, aes(x = Lunghezza)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 10, fill = "steelblue", color = "black") +
geom_density(color = "red", size = 0.7) +
labs(x = "Lunghezza", y = "Frequenza")
# Distribuzione della grandezza del cranio del neonato
g6<-ggplot(neonati, aes(x = Cranio)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 5, fill = "steelblue", color = "black") +
geom_density(color = "red", size = 0.7) +
labs(x = "Cranio", y = "Frequenza")
grid.arrange(g1, g2, g3, g4, g5, g6, ncol = 3)
La distribuzione di Anni.madre sembra approssimativamente una
normale.
La distribuzione di N.gravidanze segue un andamento tipico
della distribuzione di Poisson oppure della distribuzione Zeta (o legge
di Zipf).
La distribuzione di Gestazione è asimmetrica negativa.
Segue un andamento tipico della distribuzione binomiale negativa.
La
distribuzione del peso sembra essere leggermente asimmetrica negativa e
l’andamento è campanulare.
La distribuzione di Lunghezza è
asimmetrica negativa e l’andamento è campanulare.
La distribuzione
della dimensione del cranio è asimmetrica negativa e l’andamento è
campanulare anche in questo caso.
skewness_values <- c(
skewness(Anni.madre),
skewness(N.gravidanze),
NA,# non ha senso calcolare questi indicatori su variabili come fumatrici, tipo parto, ospedale e sesso
skewness(Gestazione),
skewness(Peso),
skewness(Lunghezza),
skewness(Cranio),
NA,
NA,
NA
)
kurtosis_values <- c(
kurtosis(Anni.madre)-3,
kurtosis(N.gravidanze)-3,
NA,
kurtosis(Gestazione)-3,
kurtosis(Peso)-3,
kurtosis(Lunghezza)-3,
kurtosis(Cranio)-3,
NA,
NA,
NA
)
# Creazione del dataframe
df <- data.frame(
Variabile = c("Anni.madre","N.gravidanze", "Fumatrici", "Gestazione", "Peso", "Lunghezza", "Cranio", "Tipo.parto", "Ospedale", "Sesso"),
Asimmetria = skewness_values,
Curtosi = kurtosis_values
)
kable(df, align = rep('c', 3)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))| Variabile | Asimmetria | Curtosi |
|---|---|---|
| Anni.madre | 0.0428115 | 0.3804165 |
| N.gravidanze | 2.5142541 | 10.9894064 |
| Fumatrici | ||
| Gestazione | -2.0653133 | 8.2581496 |
| Peso | -0.6470308 | 2.0315318 |
| Lunghezza | -1.5146991 | 6.4871739 |
| Cranio | -0.7850527 | 2.9462063 |
| Tipo.parto | ||
| Ospedale | ||
| Sesso |
Quanto ipotizzato dall’analisi grafica si rivela corretto. L’unica
eccezione è nel caso del numero di gravidanze, che è un caso estremo di
distribuzione dove non esiste una moda.
Per quanto riguarda la
curtosi, tutte le distribuzioni raggiungono dei massimi più elevati di
quello di una normale (leptocurtosi). Anni della madre si rivela
leggermente leptocurtica (e altrettanto leggermente asimmetrica
positiva).
2.3 - Boxplots
p1 <- ggplot(neonati, aes(y = Anni.madre)) +
geom_boxplot(fill = "lightblue", color = "darkblue") +
labs(title = "Anni.madre", y = "Anni") +
theme_minimal()
p2 <- ggplot(neonati, aes(y = N.gravidanze)) +
geom_boxplot(fill = "lightblue", color = "darkblue") +
labs(title = "N.gravidanze", y = "Numero") +
theme_minimal()
p3 <- ggplot(neonati, aes(y = Gestazione)) +
geom_boxplot(fill = "lightblue", color = "darkblue") +
labs(title = "Gestazione", y = "Settimane") +
theme_minimal()
p4 <- ggplot(neonati, aes(y = Peso)) +
geom_boxplot(fill = "lightblue", color = "darkblue") +
labs(title = "Peso", y = "Peso (grammi)") +
theme_minimal()
p5 <- ggplot(neonati, aes(y = Lunghezza)) +
geom_boxplot(fill = "lightblue", color = "darkblue") +
labs(title = "Lunghezza", y = "Lunghezza (cm)") +
theme_minimal()
p6 <- ggplot(neonati, aes(y = Cranio)) +
geom_boxplot(fill = "lightblue", color = "darkblue") +
labs(title = "Cranio", y = "Circonferenza (cm)") +
theme_minimal()
#boxplots relativi a Fumatrici, Tipo.parto, Ospedale e Sesso sono omessi poichè variabili binarie o con tre modalità (Ospedale)
grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 3)
Dei diversi boxplots, l’elemento più interessante è che sono presenti
diversi valori “anomali” e che si discostano dalla maggior parte degli
altri valori, talvolta verso l’alto e talvolta verso il basso.
L’eventuale influenza sarà discussa in seguito.
2.4 - Boxplot del Peso in base ad alcune features
h1<-ggplot(neonati, aes(x = Fumatrici, y = Peso, fill = Fumatrici)) +
geom_boxplot() +
labs(title = "Peso|Fumatrice", x = "fumatrici", y = "Peso (grammi)")
h2<-ggplot(neonati, aes(x = Tipo.parto, y = Peso, fill = Tipo.parto)) +
geom_boxplot() +
labs(title = "Peso|Tipo.parto", x = "tipo parto", y = "Peso (grammi)")
h3<-ggplot(neonati, aes(x = Ospedale, y = Peso, fill = Ospedale)) +
geom_boxplot() +
labs(title = "Peso|Ospedale", x = "ospedale", y = "Peso (grammi)")
h4<-ggplot(neonati, aes(x = Sesso, y = Peso, fill = Sesso)) +
geom_boxplot() +
labs(title = "Peso|Sesso", x = "Sesso", y = "Peso (grammi)")
grid.arrange(h1, h2, h3, h4, ncol = 2)
Per quanto riguarda i boxplots di Peso in base a Fumatrici,
Tipo.parto, Ospedale e Sesso, l’elemento più interessante è che le
differenze in media tra i gruppi non sembrano particolarmente grandi
(ciò comunque andrà valutato con le variabili dummy e categoriali in
fase di stima e tramite tests). Inoltre sono presenti molti valori
anomali.
3 - Verifica delle ipotesi
In questa sezione vengono testate alcune ipotesi.
3.1 - Ipotesi: In alcuni ospedali si fanno più parti cesarei
# Tabella di contingenza per Ospedale e Tipo di Parto
tab_parto <- table(Ospedale, Tipo.parto)
print(tab_parto)## Tipo.parto
## Ospedale 0 1
## 1 574 242
## 2 595 254
## 3 603 232# Test chi-quadro per verificare l'associazione tra variabili categoriali (in distribuzioni doppie)
test_parto <- chisq.test(tab_parto)
print(test_parto)##
## Pearson's Chi-squared test
##
## data: tab_parto
## X-squared = 1.0972, df = 2, p-value = 0.5778Il test del chi quadrato confronta queste frequenze osservate con le
frequenze attese, che rappresentano i valori che ci si aspetterebbe se
non ci fosse alcuna associazione tra le due variabili.
In questo
caso il test è il seguente
H0: Assenza di associazione VS H1:
Associazione
Il p-value è molto elevato quindi non si rifiuta
l’ipotesi nulla e non ci sono evidenze statistiche per concludere che
esista un’associazione tra tipo di parto e ospedale.
3.2 - Ipotesi: La media del peso e della lunghezza sono significativamente uguali a quelle della popolazione
# Definiamo dei valori ipotetici di riferimento:
pop_mean_peso <- 3250 # Secondo la letteratura, il peso della popolazione di riferimento (neonati) si aggira intorno a 3284 grammi
pop_mean_lunghezza <- 49.5 # Secondo la letteratura, la lunghezza della popolazione di riferimento (neonati) si aggira intorno a 49.5 centimetri
#fonte: https://www.who.int/publications/i/item/924154693X
# Test t a campione singolo per il Peso
t_test_peso <- t.test(Peso, mu = pop_mean_peso)
print(t_test_peso)##
## One Sample t-test
##
## data: Peso
## t = 3.2456, df = 2499, p-value = 0.001188
## alternative hypothesis: true mean is not equal to 3250
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081# Supponendo che la lunghezza cranio
t_test_lunghezza <- t.test(Lunghezza, mu = pop_mean_lunghezza)
print(t_test_lunghezza)##
## One Sample t-test
##
## data: Lunghezza
## t = 845.77, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 49.5
## 95 percent confidence interval:
## 493.6598 495.7242
## sample estimates:
## mean of x
## 494.692# Test aggiuntivo peso-----1.1
t_test_lunghezza <- t.test(Peso, mu = pop_mean_peso, alternative = "greater")
print(t_test_lunghezza)##
## One Sample t-test
##
## data: Peso
## t = 3.2456, df = 2499, p-value = 0.0005939
## alternative hypothesis: true mean is greater than 3250
## 95 percent confidence interval:
## 3266.802 Inf
## sample estimates:
## mean of x
## 3284.081# Test aggiuntivo peso------1.2
t_test_lunghezza <- t.test(Peso, mu = pop_mean_peso, alternative = "less")
print(t_test_lunghezza)##
## One Sample t-test
##
## data: Peso
## t = 3.2456, df = 2499, p-value = 0.9994
## alternative hypothesis: true mean is less than 3250
## 95 percent confidence interval:
## -Inf 3301.359
## sample estimates:
## mean of x
## 3284.081# Test aggiuntivo lunghezza cranio---------2.1
t_test_lunghezza <- t.test(Lunghezza, mu = pop_mean_lunghezza, alternative = "greater")
print(t_test_lunghezza)##
## One Sample t-test
##
## data: Lunghezza
## t = 845.77, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is greater than 49.5
## 95 percent confidence interval:
## 493.8259 Inf
## sample estimates:
## mean of x
## 494.692# Test aggiuntivo lunghezza cranio-------2.2
t_test_lunghezza <- t.test(Lunghezza, mu = pop_mean_lunghezza, alternative = "less")
print(t_test_lunghezza)##
## One Sample t-test
##
## data: Lunghezza
## t = 845.77, df = 2499, p-value = 1
## alternative hypothesis: true mean is less than 49.5
## 95 percent confidence interval:
## -Inf 495.5581
## sample estimates:
## mean of x
## 494.692I test svolti sono i seguenti:
1) H0: Peso medio (nel campione) =
3250 VS H1: Peso medio (nel campione) ≠ 3250
2) H0: Lunghezza media
del neonato (nel campione) = 49.5 VS H1: Lunghezza media del neonato
(nel campione) ≠ 49.5
1.1) H0: Peso medio del neonato (nel campione) = 3250 VS H1: Peso
medio del neonato (nel campione) > 3250
1.2) H0: Peso medio del
neonato (nel campione) = 3250 VS H1: Peso medio del neonato (nel
campione) < 3250
2.1) H0: Lunghezza media del neonato (nel campione) = 49.5 VS H1:
Lunghezza media del neonato (nel campione) > 49.5
2.2) H0:
Lunghezza media del neonato (nel campione) = 49.5 VS H1: Lunghezza media
del neonato (nel campione) < 49.5
Nel primo test si rifiuta
H0, quindi il vero peso del campione è statisticamente diverso da 3250
grammi.
Nel secondo test si rifiuta H0, quindi la vera lunghezza è
statisticamente diversa da 49.5 centimetri.
I test 1.1 e 1.2 mirano a vedere il verso della disuguaglianza del
primo test. Più precisamente viene mostrato che il peso medio del
neonato, nel campione considerato è statisticamente maggiore di 3250
grammi (infatti nel test 1.1 si rifiuta H0 ma nel test 1.2 no).
I
test 2.1 e 2.2 mirano a vedere il verso della disuguaglianza del secondo
test. Più precisamente viene mostrato che la lunghezza media del
neonato, nel campione considerato, è statisticamente maggiore di 49.5
centrimetri (infatti nel test 2.1 si rifiuta H0 ma nel test 2.2
no).
3.3 - Ipotesi: Le misure antropometriche sono significativamente diverse tra i due sessi
# Confronto per il Peso
t_test_sex_peso <- t.test(Peso ~ Sesso, data = neonati)
print(t_test_sex_peso)##
## 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 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## 207.0615 287.1051
## sample estimates:
## mean in group 0 mean in group 1
## 3408.215 3161.132# Confronto per la Lunghezza
t_test_sex_lunghezza <- t.test(Lunghezza ~ Sesso, data = neonati)
print(t_test_sex_lunghezza)##
## 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 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## 7.876273 11.929470
## sample estimates:
## mean in group 0 mean in group 1
## 499.6672 489.7643# Confronto per Cranio
t_test_sex_cranio <- t.test(Cranio ~ Sesso, data = neonati)
print(t_test_sex_cranio)##
## 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 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## 3.541270 6.089912
## sample estimates:
## mean in group 0 mean in group 1
## 342.4486 337.6330I test svolti sono i seguenti:
1) H0: Differenza nel peso medio
dei neonati maschi e femmine = 0 VS H1: Peso medio dei neonati maschi è
diverso da quello delle femmine (quindi differenza ≠ 0)
2) H0:
Differenza nella lunghezza media dei neonati maschi e femmine = 0 VS H1:
Lunghezza media dei neonati maschi è diversa da quella delle femmine
(quindi differenza ≠ 0)
3) H0: Differenza nella grandezza media del
cranio dei neonati maschi e femmine = 0 VS H1: Grandezza media del
cranio dei neonati maschi è diversa da quella delle femmine (quindi
differenza ≠ 0)
In tutti i test si rifiuta l’ipotesi nulla H0
quindi c’è una differenza statisticamente significativa nelle misure
antropometriche (peso, grandezza cranio e lunghezza del bambino).
4 - Modelli
In questa sezione vengono mostrati gli scatterplot per dare una idea orientativa sul verso e intensità della potenziale relazione tra le variabili e Peso, inoltre viene mostrato il grado di correlazione lineare tra le variabili. In seguito vengono costruiti e valutati alcuni modelli.
4.1 - Grafici (scatterplot) preliminari alla costruzione modelli regressivi
w1 <- ggplot(neonati, aes(x = Anni.madre, y = Peso)) +
geom_point(color = "blue") +
labs(title = "Scatterplot Anni.madre-Peso", x = "Anni.madre", y = "Peso")
w2 <- ggplot(neonati, aes(x = N.gravidanze, y = Peso)) +
geom_point(color = "blue") +
labs(title = "Scatterplot N.gravidanze-Peso", x = "N.gravidanze", y = "Peso")
w3 <- ggplot(neonati, aes(x = Gestazione, y = Peso)) +
geom_point(color = "blue") +
labs(title = "Scatterplot Gestazione-Peso", x = "Settimane di gestazione", y = "Peso")
w4 <- ggplot(neonati, aes(x = Lunghezza, y = Peso)) +
geom_point(color = "blue") +
labs(title = "Scatterplot Lunghezza-Peso", x = "Lunghezza", y = "Peso")
w5 <- ggplot(neonati, aes(x = Cranio, y = Peso)) +
geom_point(color = "blue") +
labs(title = "Scatterplot Cranio-Peso", x = "Grandezza cranio", y = "Peso")
# Disposizione dei grafici in una griglia
grid.arrange(w1, w2, w3, w4, w5, ncol = 3)
Le variabili relative alla grandezza cranio, lunghezza neonato e
settimane gestazione hanno una relazione positiva con il peso
(l’intensità e l’eventuale non linearità sono da valutare in
seguito).
Per quanto riguarda il numero gravidanze e gli anni della
madre, sembrano avere una relazione ambigua con il peso.
4.2 - Correlazioni
# grafico 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)
}
# sotto le correlazioni, sopra i grafici associati
suppressWarnings({ # soppressione "Warning: argument 1 does not name a graphical parameter"
pairs(neonati, lower.panel = panel.cor, upper.panel = panel.smooth)
})
Gli unici casi di correlazione elevata, a parte quella di alcune variabili con l’output, sono: lunghezza-gestazione e cranio-lunghezza.
4.3 Modelli regressivi
Dapprima è bene evidenziare che i tipo di parto e ospedale non
verranno considerate in alcun modello poichè non utili a spiegare e
prevedere il peso di un neonato.
Si tenga conto che per i dati a
disposizione sorgono certamente almeno due problemi:
- Qualsiasi
tipo di stima e qualsiasi modello si usi, genereranno stime più o meno
distorte a causa di variabili omesse (esempio caratteristiche del padre
o altre caratteristiche rilevanti). - È molto probabile che si
riscontrino problemi di doppia causalità (lunghezza, cranio e sesso
incidono su peso ma vale anche il contrario, oppure il tipo di peso
incide sul peso ma anche viceversa). Dunque la stima sarà sicuramente
distorta. Una soluzione appropriata è quella di adottare una
modellistica basata sulle variabili strumentali e approfondire la
questione della esogenità ed endogenità. La modellistica presentata,
regressione lineare, NON è appropriata.
m0.1 <- lm(Peso ~ Anni.madre)
summary(m0.1)##
## Call:
## lm(formula = Peso ~ Anni.madre)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2461.16 -289.97 8.84 332.49 1641.08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3347.087 57.063 58.656 <2e-16 ***
## Anni.madre -2.237 1.991 -1.123 0.261
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 525 on 2498 degrees of freedom
## Multiple R-squared: 0.0005049, Adjusted R-squared: 0.0001048
## F-statistic: 1.262 on 1 and 2498 DF, p-value: 0.2614m0.2 <- lm(Peso ~ I(Anni.madre))
summary(m0.2)##
## Call:
## lm(formula = Peso ~ I(Anni.madre))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2461.16 -289.97 8.84 332.49 1641.08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3347.087 57.063 58.656 <2e-16 ***
## I(Anni.madre) -2.237 1.991 -1.123 0.261
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 525 on 2498 degrees of freedom
## Multiple R-squared: 0.0005049, Adjusted R-squared: 0.0001048
## F-statistic: 1.262 on 1 and 2498 DF, p-value: 0.2614m0.3<- lm(Peso ~ N.gravidanze)
summary(m0.3)##
## Call:
## lm(formula = Peso ~ N.gravidanze)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2453.11 -293.11 14.91 336.15 1646.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3283.112 13.232 248.11 <2e-16 ***
## N.gravidanze 0.987 8.203 0.12 0.904
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 525.1 on 2498 degrees of freedom
## Multiple R-squared: 5.795e-06, Adjusted R-squared: -0.0003945
## F-statistic: 0.01448 on 1 and 2498 DF, p-value: 0.9042m0.4 <- lm(Peso ~ Fumatrici)
summary(m0.4)##
## Call:
## lm(formula = Peso ~ Fumatrici)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2456.15 -286.15 13.85 333.85 1693.65
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3286.15 10.73 306.359 <2e-16 ***
## Fumatrici1 -49.81 52.59 -0.947 0.344
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 525 on 2498 degrees of freedom
## Multiple R-squared: 0.0003589, Adjusted R-squared: -4.125e-05
## F-statistic: 0.8969 on 1 and 2498 DF, p-value: 0.3437m0.5 <- lm(Peso ~ Gestazione)
summary(m0.5)##
## Call:
## lm(formula = Peso ~ Gestazione)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1603.61 -287.34 -11.07 280.46 1897.75
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3197.250 176.851 -18.08 <2e-16 ***
## Gestazione 166.272 4.532 36.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 423.3 on 2498 degrees of freedom
## Multiple R-squared: 0.3502, Adjusted R-squared: 0.3499
## F-statistic: 1346 on 1 and 2498 DF, p-value: < 2.2e-16m0.6 <- lm(Peso ~ I(Gestazione))
summary(m0.6)##
## Call:
## lm(formula = Peso ~ I(Gestazione))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1603.61 -287.34 -11.07 280.46 1897.75
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3197.250 176.851 -18.08 <2e-16 ***
## I(Gestazione) 166.272 4.532 36.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 423.3 on 2498 degrees of freedom
## Multiple R-squared: 0.3502, Adjusted R-squared: 0.3499
## F-statistic: 1346 on 1 and 2498 DF, p-value: < 2.2e-16m0.7 <- lm(Peso ~ Lunghezza)
summary(m0.7)##
## Call:
## lm(formula = Peso ~ Lunghezza)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1478.4 -206.0 -22.4 190.4 3939.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4571.8176 119.6778 -38.20 <2e-16 ***
## Lunghezza 15.8804 0.2416 65.73 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 317.8 on 2498 degrees of freedom
## Multiple R-squared: 0.6337, Adjusted R-squared: 0.6335
## F-statistic: 4321 on 1 and 2498 DF, p-value: < 2.2e-16m0.8 <- lm(Peso ~ I(Lunghezza))
summary(m0.8)##
## Call:
## lm(formula = Peso ~ I(Lunghezza))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1478.4 -206.0 -22.4 190.4 3939.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4571.8176 119.6778 -38.20 <2e-16 ***
## I(Lunghezza) 15.8804 0.2416 65.73 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 317.8 on 2498 degrees of freedom
## Multiple R-squared: 0.6337, Adjusted R-squared: 0.6335
## F-statistic: 4321 on 1 and 2498 DF, p-value: < 2.2e-16m0.9 <- lm(Peso ~ Cranio)
summary(m0.9)##
## Call:
## lm(formula = Peso ~ Cranio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2602.06 -241.01 -8.19 241.69 1421.29
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4376.4751 154.4533 -28.34 <2e-16 ***
## Cranio 22.5291 0.4537 49.66 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 372.5 on 2498 degrees of freedom
## Multiple R-squared: 0.4967, Adjusted R-squared: 0.4965
## F-statistic: 2466 on 1 and 2498 DF, p-value: < 2.2e-16m0.10 <- lm(Peso ~ I(Cranio))
summary(m0.10)##
## Call:
## lm(formula = Peso ~ I(Cranio))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2602.06 -241.01 -8.19 241.69 1421.29
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4376.4751 154.4533 -28.34 <2e-16 ***
## I(Cranio) 22.5291 0.4537 49.66 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 372.5 on 2498 degrees of freedom
## Multiple R-squared: 0.4967, Adjusted R-squared: 0.4965
## F-statistic: 2466 on 1 and 2498 DF, p-value: < 2.2e-16m0.11 <- lm(Peso ~ Sesso)
summary(m0.11)##
## Call:
## lm(formula = Peso ~ Sesso)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2428.22 -261.13 8.87 311.78 1768.87
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3408.22 14.47 235.5 <2e-16 ***
## Sesso1 -247.08 20.42 -12.1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 510.4 on 2498 degrees of freedom
## Multiple R-squared: 0.05539, Adjusted R-squared: 0.05501
## F-statistic: 146.5 on 1 and 2498 DF, p-value: < 2.2e-16I modelli m0.1-m0.11 mostrano che:
- Anni.madre, N.gravidanze e
Fumatrici, non sembrano avere un impatto significativo sul peso di un
neonato.
- Gestazione, Lunghezza, Cranio e Sesso sembrano avere un
effetto significativo sul peso.
- I termini non lineari delle
variabili che sembrano significative non apportano alcun miglioramento a
livello di modello (nè in termini di r2adj, ne in termini di residui),
quindi la modellazione verrà semplificata con l’utilizzo di termini
lineari.
m1 <- lm(Peso ~ Anni.madre+N.gravidanze+Fumatrici+Gestazione+Lunghezza+Cranio+Sesso)
summary(m1)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Sesso)
##
## 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) -6636.3263 143.0540 -46.390 < 2e-16 ***
## Anni.madre 0.9585 1.1347 0.845 0.3984
## N.gravidanze 11.2756 4.6690 2.415 0.0158 *
## Fumatrici1 -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 ***
## Sesso1 -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-16Nel modello pieno (m1) ci sono variabili non significative quindi si procede alla rimozione prima di Anni.madre (m2) e poi della variabile dummy Fumatrici (m3).
m2<- lm(Peso ~ N.gravidanze+Fumatrici+Gestazione+Lunghezza+Cranio+Sesso)
summary(m2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso)
##
## 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) -6603.5001 137.6662 -47.967 < 2e-16 ***
## N.gravidanze 12.7185 4.3450 2.927 0.00345 **
## Fumatrici1 -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 ***
## Sesso1 -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-16m3 <- lm(Peso ~ N.gravidanze+Gestazione+Lunghezza+Cranio+Sesso)
summary(m3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso)
##
## 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) -6603.1518 137.6719 -47.963 < 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 ***
## Sesso1 -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-16- (Intercept): Il valore di -6603.1518 rappresenta l’intercetta del
modello, ovvero il peso previsto del neonato quando tutte le altre
variabili sono zero (inoltre si intende che anche Sesso=0 quindi la
stima è riferita ai neonati maschi). Tuttavia, poiché valori nulli per
le variabili come lunghezza o cranio non sono realistici, l’intercetta
ha un significato limitato nel contesto pratico.
- N.gravidanze: Il coefficiente di 12.4750 indica che, a parità delle
altre variabili, ogni gravidanza è associata a un aumento medio di circa
12.48 grammi nel peso del neonato.
- Gestazione: Con un coefficiente di 32.3321, ogni settimana
aggiuntiva di gestazione è associata a un aumento medio di circa 32.33
grammi nel peso del neonato.
- Lunghezza: Il coefficiente di 10.2486 suggerisce che ogni centimetro
in più nella lunghezza del neonato è associato a un aumento medio di
circa 10.25 grammi nel peso.
- Cranio: Con un coefficiente di 10.5402, ogni centimetro in più nella
circonferenza cranica è associato a un aumento medio di circa 10.54
grammi nel peso del neonato.
- Sesso1: Il coefficiente di -77.9927 indica che, a parità delle altre
variabili, i neonati di sesso femminile (codificati come 1) pesano in
media circa 78 grammi in meno rispetto ai maschi. Quindi è
statisticamente significativa la differenza di peso alla nascita tra i
sessi. Questo risultato era osservabile già nel confronto dei boxplot
del peso rispetto a maschi e femmine.
In sintesi, il modello evidenzia che la durata della gestazione, la lunghezza e la circonferenza cranica del neonato sono fortemente associate al peso alla nascita. Inoltre, il numero di gravidanze precedenti e il sesso del neonato influenzano significativamente il peso, con le femmine che tendono a pesare meno dei maschi (in questo campione).
anova(m2,m3)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2493 187973654
## 2 2494 188065546 -1 -91892 1.2187 0.2697Il test ANOVA, sopra, vuole verificare se la rimozione della variabile dummy Fumatrici è una scelta opportuna poichè Fumatrici=1 non è significativa ma l’intercetta del modello si (Fumatrici=0). Il test rivela che l’aggiunta della variabile Fumatrice non apporta miglioramenti significativi al modello.
vif(m3)## N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1.023475 1.669189 2.074689 1.624465 1.040054Il calcolo del vif non mostra possibili problemi di collinearità.
Il modello m4, sfrutta le variabili del modello m3 tenendo anche
conto di possibili interazioni tra le variabili dummy Sesso e fumatrici
(altre interazioni con le altre variabili generano altissima
multicollinearità). Nella proposta m4 la dummy Fumatrici e l’interazione
sono non significative, in m5 Fumatrici resta non significativa
(p-value>0.05) e inoltre genera la non significatività di
N.gravidanze.
#interazioni
m4 <- lm(Peso ~ N.gravidanze+Gestazione+Lunghezza+Cranio+Sesso+Sesso*Fumatrici)
summary(m4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + Sesso * Fumatrici)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1151.1 -182.1 -16.2 163.1 2637.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6605.0163 137.7028 -47.966 < 2e-16 ***
## N.gravidanze 12.7379 4.3457 2.931 0.00341 **
## Gestazione 32.6165 3.8057 8.570 < 2e-16 ***
## Lunghezza 10.2357 0.3010 34.011 < 2e-16 ***
## Cranio 10.5371 0.4263 24.717 < 2e-16 ***
## Sesso1 -79.6007 11.4238 -6.968 4.1e-12 ***
## Fumatrici1 -47.1677 37.9551 -1.243 0.21409
## Sesso1:Fumatrici1 35.3381 55.1226 0.641 0.52153
## ---
## 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: 948.9 on 7 and 2492 DF, p-value: < 2.2e-16vif(m4)## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif## N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1.026170 1.675753 2.078800 1.624634 1.081469
## Fumatrici Sesso:Fumatrici
## 1.903891 1.935450m5 <- lm(Peso ~ Gestazione+Cranio+N.gravidanze+Sesso*Fumatrici)
summary(m5)##
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + N.gravidanze + Sesso *
## Fumatrici)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1317.12 -219.17 -12.71 210.62 1597.14
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6111.1127 165.6631 -36.889 <2e-16 ***
## Gestazione 93.5695 4.0617 23.037 <2e-16 ***
## Cranio 17.0723 0.4604 37.083 <2e-16 ***
## N.gravidanze 5.6941 5.2514 1.084 0.2783
## Sesso1 -119.2045 13.7485 -8.670 <2e-16 ***
## Fumatrici1 -80.4223 45.9026 -1.752 0.0799 .
## Sesso1:Fumatrici1 19.0737 66.6845 0.286 0.7749
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 332.2 on 2493 degrees of freedom
## Multiple R-squared: 0.6005, Adjusted R-squared: 0.5996
## F-statistic: 624.7 on 6 and 2493 DF, p-value: < 2.2e-16vif(m5)## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif## Gestazione Cranio N.gravidanze Sesso Fumatrici
## 1.304143 1.294569 1.023839 1.070233 1.902628
## Sesso:Fumatrici
## 1.935304In conclusione, basandoci sulla significatività, errore standard dei residui e r2adj, m3 è il miglior modello a cui si può giungere secondo la metodologia seguita.
kable(cbind(AIC(m1,m2,m3,m4,m5),BIC(m1,m2,m3,m4,m5)), align = rep('c', 5)) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))| df | AIC | df | BIC | |
|---|---|---|---|---|
| m1 | 9 | 35181.39 | 9 | 35233.81 |
| m2 | 8 | 35180.11 | 8 | 35226.70 |
| m3 | 7 | 35179.33 | 7 | 35220.10 |
| m4 | 9 | 35181.70 | 9 | 35234.11 |
| m5 | 8 | 36132.95 | 8 | 36179.54 |
La considerazione fatta sopra resta valida anche seguendo i criteri
AIC e BIC (si ricorda che m4 e m5 presentano problemi di collinearità e
distorsioni che fanno ritornare al modello m3).
Per completezza
vengono riportati di seguito due indicatori già discussi sopra.
r_squared <- summary(m1)$adj.r.squared
rmse <- sqrt(mean(m1$residuals^2))
cbind(paste("m1:::","R² =", r_squared),paste("RMSE =", rmse))## [,1] [,2]
## [1,] "m1::: R² = 0.726446725977325" "RMSE = 274.167723309726"r_squared <- summary(m2)$adj.r.squared
rmse <- sqrt(mean(m2$residuals^2))
cbind(paste("m2:::","R² =", r_squared),paste("RMSE =", rmse))## [,1] [,2]
## [1,] "m2::: R² = 0.726478165258627" "RMSE = 274.206968804178"r_squared <- summary(m3)$adj.r.squared
rmse <- sqrt(mean(m3$residuals^2))
cbind(paste("m3:::","R² =", r_squared),paste("RMSE =", rmse))## [,1] [,2]
## [1,] "m3::: R² = 0.726454178140911" "RMSE = 274.273984419762"r_squared <- summary(m4)$adj.r.squared
rmse <- sqrt(mean(m4$residuals^2))
cbind(paste("m4:::","R² =", r_squared),paste("RMSE =", rmse))## [,1] [,2]
## [1,] "m4::: R² = 0.726413525764128" "RMSE = 274.18436019834"r_squared <- summary(m5)$adj.r.squared
rmse <- sqrt(mean(m5$residuals^2))
cbind(paste("m5:::","R² =", r_squared),paste("RMSE =", rmse))## [,1] [,2]
## [1,] "m5::: R² = 0.599578510820831" "RMSE = 331.773103766806"#il migliore è m3 secondo r2adj, rmse aic, bic e vif5 - Analisi dei residui e diagnostica del modello m3
In questa sezione viene mostrata una valutazione più approfondita del modello m3 che risulta il migliore in un contesto di regressione lineare multipla.
# Creazione dell'istogramma con 80 barre e scala di densità
hist(residuals(m3), breaks = 80, freq = FALSE,
main = "Istogramma dei residui con stima kernel",
xlab = "Residui", col = "lightgray", border = "white")
# Sovrapposizione della curva di densità kernel
lines(density(residuals(m3)), col = "red", lwd = 1)
La stima della funzione di densità dei residui sembra avere la forma di una normale, però la coda destra è notevolmente più allungata
# Plot diagnostici
par(mfrow=c(2,2))
plot(m3)
par(mfrow=c(1,1)) # reset finestra- I residui rispetto ai valori adattati dal modello (residuals vs
fitted) sembrano essere sparsi attorno a zero, però si nota anche un
possibile trend a “U” (quindi una possibile eteroschedasticità
latente).
- Il valori che rappresentano la radice quadrata dei valori assoluti
dei residui standardizzati rispetto ai valori adattati (scale-location)
sembrano essere addensati attorno ad un particolare punto ed inoltre
sembra esserci un pattern con un trend a “U” oppure sinusoidale (dunque
possibile eteroschedasticità).
- La distribuzione dei residui sembra essere in linea con la
distribuzione teorica della normale (Q-Q residuals) e solamente i valori
più estremi sembrano discostarsi dalle stime teoriche.
- Il grafico residuals vs leverage mostra i residui standardizzati
rispetto alla leva per ciascun punto e l’unico valore potenzialmente
influente nel modello è l’osservazione 1551.
shapiro.test(m3$residuals)##
## Shapiro-Wilk normality test
##
## data: m3$residuals
## W = 0.97408, p-value < 2.2e-16lmtest::bptest(m3)##
## studentized Breusch-Pagan test
##
## data: m3
## BP = 90.253, df = 5, p-value < 2.2e-16lmtest::dwtest(m3)##
## Durbin-Watson test
##
## data: m3
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0Nei primi 2 test si rifiuta h0 quindi la distribuzione dei residui non è normale e non ha varianza costante (quindi c’è eteroschedasticità, come già ipotizzato). Si conclude che questo modello è inappropriato e distrorto! Probabilmente sarebbe un buon modello se si inserisse la robustezza all’eteroschedasticità (che genera stime più attendibili sia nel caso sia di eteroschedasticità che di omoschedasticità).
# leverage
lev<-hatvalues(m3)
plot(lev)
p<-sum(lev)
n<-length(lev)
soglia=2*p/n
abline(h=soglia,col=2)
lev[lev>soglia]## 13 15 34 67 89 96
## 0.005630918 0.007050422 0.006744143 0.005891590 0.012816470 0.005351586
## 101 106 131 134 151 155
## 0.007526751 0.014479871 0.007229339 0.007552819 0.010883937 0.007207682
## 161 189 190 204 205 206
## 0.020335483 0.004893297 0.005366557 0.014490604 0.005351634 0.009476652
## 220 294 305 310 312 315
## 0.007393997 0.005912765 0.005442061 0.028812123 0.013169272 0.005385800
## 378 440 442 445 486 492
## 0.015934061 0.005404736 0.007723662 0.007509382 0.005164446 0.008274018
## 497 516 582 587 592 614
## 0.005166306 0.013079851 0.011665555 0.008412325 0.006384116 0.005299262
## 638 656 657 684 697 702
## 0.006688287 0.005927777 0.005322685 0.008818987 0.005863826 0.005202259
## 729 748 750 757 765 805
## 0.005023115 0.008565543 0.006942097 0.008145491 0.006070298 0.014356657
## 828 893 895 913 928 946
## 0.007179817 0.005075205 0.005295896 0.005571144 0.022742332 0.006909044
## 947 956 985 1008 1014 1049
## 0.008409465 0.007784123 0.007039416 0.005343037 0.008470133 0.004956169
## 1067 1091 1106 1130 1166 1181
## 0.008465430 0.008933360 0.005967317 0.031728597 0.005513559 0.005677676
## 1188 1200 1219 1238 1248 1273
## 0.006477203 0.005492370 0.030694311 0.005908078 0.014622914 0.007085831
## 1291 1293 1311 1321 1325 1356
## 0.006117497 0.006073639 0.009625908 0.009293111 0.004857169 0.005303442
## 1357 1385 1395 1400 1402 1411
## 0.006965051 0.012636943 0.005126697 0.005925069 0.004811441 0.008048184
## 1420 1428 1429 1450 1505 1551
## 0.005155654 0.008192811 0.021757172 0.015104831 0.013330439 0.048769569
## 1553 1556 1573 1593 1606 1610
## 0.008504889 0.005919673 0.005047204 0.005623758 0.005001812 0.008722184
## 1617 1619 1628 1686 1693 1701
## 0.004866796 0.015067498 0.005069731 0.009349313 0.005077858 0.010842957
## 1712 1718 1727 1735 1780 1781
## 0.006992084 0.006958857 0.013300523 0.004884846 0.025538678 0.016832361
## 1809 1827 1868 1892 1962 1967
## 0.008707504 0.006065698 0.005205637 0.005332812 0.005540442 0.005337356
## 1977 2037 2040 2046 2086 2089
## 0.006927281 0.004889127 0.011494872 0.005471670 0.013193090 0.006293550
## 2098 2114 2115 2120 2140 2146
## 0.005094455 0.013316875 0.011772090 0.018659995 0.006244232 0.005802168
## 2148 2149 2157 2175 2200 2215
## 0.007926839 0.013583436 0.005907225 0.032527273 0.011670024 0.004892265
## 2216 2220 2221 2224 2225 2244
## 0.008117864 0.005414040 0.021628717 0.005838076 0.005591261 0.006929217
## 2257 2307 2317 2318 2337 2359
## 0.006170254 0.013965608 0.007673614 0.004831118 0.005230450 0.010067364
## 2408 2422 2436 2437 2452 2458
## 0.009696691 0.021532808 0.004986522 0.023943328 0.023838489 0.008506087
## 2471 2478
## 0.020903740 0.005775173lev[lev>0.04]## 1551
## 0.04876957Le osservazioni con valori di leva superiori alla soglia 2 * p / n
sono considerate potenzialmente influenti nello spazio dei regressori.
Questi punti hanno combinazioni di valori delle variabili indipendenti
che si discostano significativamente dalla media delle altre
osservazioni, conferendo loro un peso maggiore nella determinazione dei
coefficienti del modello. Sembrano essere molti.
L’osservazione
associata al valore di leva più alto (0.04876957) è quella con indice
1551 (osservazione che sembrava problematica già nell’analisi dei
residui).
# outliers
plot(rstudent(m3))
abline(h=c(-2,2))
Sono presenti diversi valori influenti, ovvero outliers, nello spazio della variabile di output. Di seguito un test per valutarne l’effettiva influenza.
car::outlierTest(m3)## 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-03Il p-value non aggiustato (unadjusted p-value) rappresenta la
probabilità associata al residuo studentizzato ed indica la possibilità
di osservare un residuo così estremo sotto l’ipotesi nulla che tutti i
dati seguano il modello senza outlier. Un valore p basso suggerisce che
il residuo è significativamente diverso dagli altri.
Poiché testare
ogni osservazione individualmente aumenta la probabilità di falsi
positivi (identificare erroneamente un outlier), si applica la
correzione di Bonferroni. Questa correzione moltiplica il p-value non
aggiustato per il numero totale di osservazioni, adattando il livello di
significatività per il numero di confronti effettuati.
Nei tre test
proposti si rifiuta H0 quindi le osservazioni 155,1551 e 1306 sono
outlier.
# distanza di cook
cook<-cooks.distance(m3)
plot(cook,ylim = c(0,0.9)) 
max(cook)## [1] 0.8300965# cook ----output omesso poichè troppo lungounname(which.max(cook)) #indice osservazione## [1] 1551L’osservazione 1551 risulta un vero e proprio problema perchè sia
nell’ambito dello spazio dei regressori, che nello spazio dell’output
genera influenza notevole.
L’osservazione 1551 non sembra avere
valori delle variabili particolarmente diverse da altre. Però,
evidentemente, la combinazione di valori che assume congiuntamente è
notevolmente diversa dalle altre osservazioni. È improbabile un errore
nell’immissione delle cifre, quindi bisognerebbe (a) verificare il dato
nella fonte primaria che ha generato il dato stesso e stimare un modello
dove questa osservazione viene completamente eliminata.
6 - Previsione
In questa sezione è presentato un esempio di previsione (con il
modello m3) del peso di un neonato nel caso in cui la madre sia alla
terza gravidanza, con 39 settimane di gestazione.
Come mostrato di
seguito, in queste condizioni si stima un peso del neonato pari a
3271.09 (stima puntuale), ovvero un peso compreso tra 2732.11 e 3810.07
(stima intervallare).
# con valori medi per le altre variabili (modifica in base alla distribuzione dei dati).
nuovo_caso <- data.frame(
N.gravidanze = 3,
Gestazione = 39,
Lunghezza = mean(Lunghezza),
Cranio = mean(Cranio),
Sesso = as.factor(1)
)
predizione <- predict(m3, newdata = nuovo_caso, interval = "prediction")
print("Previsione modello 3: ")## [1] "Previsione modello 3: "print(predizione)## fit lwr upr
## 1 3271.09 2732.11 3810.07I grafici di seguito raffigurano dei modelli di regressione semplice in base al sesso ed a specifiche variabili (quelle che compongono il modello m3). Ovviamente i modello m3, più preciso di modelli lineari semplici con un sono regressore (come le raffigurazioni di seguito), non è rappresentabile poichè presenta una dimensionalità maggiore di 3.
# Grafico 1: Peso vs Numero di Gravidanze
f1 <- ggplot(neonati, aes(x = N.gravidanze, y = Peso)) +
geom_point(aes(color = Sesso), alpha = 0.3) +
geom_smooth(aes(group = Sesso), method = "lm", formula = y ~ x, se = FALSE, color = "black", size = 1) +
geom_smooth(aes(color = Sesso), method = "lm", formula = y ~ x, se = FALSE, size = 0.5) +
labs(
x = "N.gravidanze",
y = "Peso")
# Grafico 2: Peso vs Settimane di Gestazione
f2 <- ggplot(neonati, aes(x = Gestazione, y = Peso)) +
geom_point(aes(color = Sesso), alpha = 0.3) +
geom_smooth(aes(group = Sesso), method = "lm", formula = y ~ x, se = FALSE, color = "black", size = 1) +
geom_smooth(aes(color = Sesso), method = "lm", formula = y ~ x, se = FALSE, size = 0.5) +
labs(
x = "Gestazione",
y = "Peso")
# Grafico 3: Peso vs Lunghezza
f3 <- ggplot(neonati, aes(x = Lunghezza, y = Peso)) +
geom_point(aes(color = Sesso), alpha = 0.3) +
geom_smooth(aes(group = Sesso), method = "lm", formula = y ~ x, se = FALSE, color = "black", size = 1) +
geom_smooth(aes(color = Sesso), method = "lm", formula = y ~ x, se = FALSE, size = 0.5) +
labs(
x = "Lunghezza",
y = "Peso")
# Grafico 4: Peso vs Circonferenza Cranica
f4 <- ggplot(neonati, aes(x = Cranio, y = Peso)) +
geom_point(aes(color = Sesso), alpha = 0.3) +
geom_smooth(aes(group = Sesso), method = "lm", formula = y ~ x, se = FALSE, color = "black", size = 1) +
geom_smooth(aes(color = Sesso), method = "lm", formula = y ~ x, se = FALSE, size = 0.5) +
labs(
x = "Cranio",
y = "Peso")
# Disposizione dei grafici in una griglia
grid.arrange(f1, f2, f3, f4, ncol = 2)
Di seguito vengono, invece rappresentate le relazioni variabile-Peso (a parità di altre condizioni) utilizzando il modello m3. Ovviamente i grafici seguenti differiscono dai precedenti poichè i coefficienti (effetti marginali delle singole variabili) sono stimati col modello di regressione multipla (le figure sopra erano stimati con modelli di regressione semplice come i modelli m0.1-m0.10).
# Caricamento del pacchetto ggeffects
library(ggeffects)# Ottenere le previsioni marginali con interazione tra la variabile di interesse e Sesso
effetti_ngravidanze <- ggpredict(m3, terms = c("N.gravidanze", "Sesso"))
effetti_gestazione <- ggpredict(m3, terms = c("Gestazione", "Sesso"))
effetti_cranio <- ggpredict(m3, terms = c("Cranio", "Sesso"))
effetti_lunghezza <- ggpredict(m3, terms = c("Lunghezza", "Sesso"))
# Visualizzare le previsioni con le osservazioni originali
b1 <- plot(effetti_ngravidanze, show_data = TRUE)+ labs(title = NULL)## Data points may overlap. Use the `jitter` argument to add some amount of
## random variation to the location of data points and avoid overplotting.b2 <- plot(effetti_gestazione, show_data = TRUE)+ labs(title = NULL)## Data points may overlap. Use the `jitter` argument to add some amount of
## random variation to the location of data points and avoid overplotting.b3 <- plot(effetti_lunghezza, show_data = TRUE)+ labs(title = NULL)## Data points may overlap. Use the `jitter` argument to add some amount of
## random variation to the location of data points and avoid overplotting.b4 <- plot(effetti_cranio, show_data = TRUE)+ labs(title = NULL)## Data points may overlap. Use the `jitter` argument to add some amount of
## random variation to the location of data points and avoid overplotting.# Disposizione dei grafici in una griglia
grid.arrange(b1, b2, b3, b4, ncol = 2)