Progetto Statistica Inferenziale
Diego Amelio
2025-05-31
1. INTRODUZIONE
Scopo del presente lavoro è ricavare una correlazione tra la variabile peso del neonato e le altre variabili nel database.
Carichiamo il database ed esploriamolo
db=read.csv("neonati.csv")
summary(db)## Anni.madre N.gravidanze Fumatrici Gestazione
## Min. : 0.00 Min. : 0.0000 Min. :0.0000 Min. :25.00
## 1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:38.00
## Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00
## Mean :28.16 Mean : 0.9812 Mean :0.0416 Mean :38.98
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:40.00
## Max. :46.00 Max. :12.0000 Max. :1.0000 Max. :43.00
## Peso Lunghezza Cranio Tipo.parto
## Min. : 830 Min. :310.0 Min. :235 Length:2500
## 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Class :character
## Median :3300 Median :500.0 Median :340 Mode :character
## Mean :3284 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :565.0 Max. :390
## Ospedale Sesso
## Length:2500 Length:2500
## Class :character Class :character
## Mode :character Mode :character
##
##
## str(db)## '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 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ 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 : chr "Nat" "Nat" "Nat" "Nat" ...
## $ Ospedale : chr "osp3" "osp1" "osp2" "osp2" ...
## $ Sesso : chr "M" "F" "M" "M" ...n=nrow(db)for (v in names(db)) {
cat("Valori nulli nella variabile", v, "=", sum(is.na(db[[v]])), "\n")
}## Valori nulli nella variabile Anni.madre = 0
## Valori nulli nella variabile N.gravidanze = 0
## Valori nulli nella variabile Fumatrici = 0
## Valori nulli nella variabile Gestazione = 0
## Valori nulli nella variabile Peso = 0
## Valori nulli nella variabile Lunghezza = 0
## Valori nulli nella variabile Cranio = 0
## Valori nulli nella variabile Tipo.parto = 0
## Valori nulli nella variabile Ospedale = 0
## Valori nulli nella variabile Sesso = 0Non ci sono valori N/A nel database
2 STUDIO DELL’INFLUENZA SUL PESO DEI FATTORI
2.1 Il fattore N. di gravidanze
Visualizziamo la distribuzione dei neonati per Numero di gravidanze precedenti della madre.
records_by_N.gravidanze=table(db$N.gravidanze)/n
barplot(records_by_N.gravidanze,
main="Distribuzione dei record del database",
xlab="N. gravidanze precedenti",
ylab="frequenze relative",
ylim = c(0, 0.5), # range asse y
col="lightblue",
border="darkblue"
)
records_by_N.gravidanze##
## 0 1 2 3 4 5 6 7 8 9 10
## 0.4384 0.3272 0.1360 0.0600 0.0192 0.0084 0.0044 0.0004 0.0032 0.0008 0.0012
## 11 12
## 0.0004 0.0004Le donne con un numero di gravidanze >= 4 hanno una frequenza relativa < 5%. Conviene accorpare questi record in un’unica classe “4+”.
db$grav_cat <- ifelse(db$N.gravidanze >= 4, "4+", as.character(db$N.gravidanze))
db$grav_cat <- factor(db$grav_cat, levels = c("0", "1", "2", "3", "4+"))records_by_grav_cat=table(db$grav_cat)/n
barplot(records_by_grav_cat,
main="Distribuzione dei record del database",
xlab="N. gravidanze precedenti",
ylab="frequenze relative",
ylim = c(0, 0.5), # range asse y
col="lightblue",
border="darkblue"
)
Studiamo l’influenza del fattore grav_cat sulla variabile Peso
#boxplot(db$Peso~db$grav_cat)
g=ggplot(data=db)+
labs(title=paste("Distribuzioni di","Peso","in funzione di","Numero di gravidanze precedenti"),
x="Numero di gravidanze precedenti",
y="Peso [g]")
g=g+
geom_half_boxplot(aes(x=grav_cat, y=Peso, fill=grav_cat), side="l" )+
geom_half_violin(aes(x=grav_cat, y=Peso, fill=grav_cat), side="r")+
stat_summary(aes(x = grav_cat, y = Peso, group = grav_cat,linetype = "Media"),
fun = mean,
geom = "crossbar",
width = 0.8,
fatten = 1,
color = "black") +
scale_linetype_manual(
name = "Statistiche", # Titolo della legenda
values = c("Media" = "dashed") # Etichetta + tipo di linea
)#+
#scale_fill_manual(values = set_colori)
g
leveneTest(Peso ~ grav_cat, data = db)#controlliamo che la varianza sia simile tra i gruppi## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 4 2.599 0.03451 *
## 2495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- si rifiuta l’ipotesi che le varianze dei gruppi siano uguali
pairwise.t.test(db$Peso, db$grav_cat, paired = FALSE, pool.sd = FALSE,
p.adjust.method = "bonferroni") #chiamiamo il t test a coppie di gruppi con pool.sd = FALSE##
## Pairwise comparisons using t tests with non-pooled SD
##
## data: db$Peso and db$grav_cat
##
## 0 1 2 3
## 1 1 - - -
## 2 1 1 - -
## 3 1 1 1 -
## 4+ 1 1 1 1
##
## P value adjustment method: bonferroninon si rifiuta l’ipotesi che le medie dei gruppi siano uguali
visivamente i boxplot della distribuzione del peso nei cinque gruppi sembrano equivalenti
2.2 Il fattore Fumatrice
Vediamo come si distribuiscono i record in funzione del fattore Fumatrici.
db$Fumatrici_cat <- ifelse(db$Fumatrici == 0, "NF", "F")
db$Fumatrici_cat <- factor(db$Fumatrici_cat, levels = c("NF", "F"))
records_by_Fumatrici_cat=table(db$Fumatrici_cat)/n
bar_heights=barplot(records_by_Fumatrici_cat,
main="Distribuzione dei record del database",
xlab="Fattore fumo (Non Fumatrice / Fumatrice)",
ylab="frequenze relative",
ylim = c(0, 1.1), # range asse y
col = "lightblue", # colore barre (opzionale)
border = "darkblue", # bordo barre (opzionale)
axes = FALSE # disattiva disegno automatico degli assi
)
# Aggiunge asse Y con tick personalizzati
axis(2, at = seq(0, 1, by = 0.1), las = 1) # las=1 => etichette orizzontali
# Aggiunge griglia orizzontale
abline(h = seq(0, 1, by = 0.1), col = "lightgray", lty = "dotted")
# Aggiunge etichette numeriche sopra ogni barra
text(x = bar_heights, y = records_by_Fumatrici_cat,
labels = round(records_by_Fumatrici_cat, 3),
pos = 3, cex = 0.8, col = "black")
Le donne fumatrici sono il 4.2% dei record.
Studiamo l’influenza del fattore Fumatrici sulla variabile peso
print("statistiche Peso nati con madre Non Fumatrice")## [1] "statistiche Peso nati con madre Non Fumatrice"summary(db$Peso[db$Fumatrici_cat=="NF"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 3000 3300 3286 3620 4900sd(db$Peso[db$Fumatrici_cat=="NF"])## [1] 526.9477print('')## [1] ""print("statistiche Peso nati con madre Fumatrice")## [1] "statistiche Peso nati con madre Fumatrice"summary(db$Peso[db$Fumatrici_cat=="F"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1780 2940 3175 3236 3440 4930sd(db$Peso[db$Fumatrici_cat=="F"])## [1] 478.7972#boxplot(db$Peso~db$Fumatrici_cat)
g=ggplot(data=db)+
labs(title=paste("Distribuzioni di","Peso","in funzione di","Madre Non Fumatrice/ Fumatrice"),
x="Madre Fumatrice / Non Fumatrice",
y="Peso [g]")
g=g+
geom_half_boxplot(aes(x=Fumatrici_cat, y=Peso, fill=Fumatrici_cat), side="l" )+
geom_half_violin(aes(x=Fumatrici_cat, y=Peso, fill=Fumatrici_cat), side="r")+
stat_summary(aes(x = Fumatrici_cat, y = Peso, group = Fumatrici_cat,linetype = "Media"),
fun = mean,
geom = "crossbar",
width = 0.8,
fatten = 1,
color = "black") +
scale_linetype_manual(
name = "Statistiche", # Titolo della legenda
values = c("Media" = "dashed") # Etichetta + tipo di linea
)#+
#scale_fill_manual(values = set_colori)
g
leveneTest(Peso ~ Fumatrici_cat, data = db)#controlliamo che la varianza sia simile tra i gruppi## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.4129 0.2347
## 2498- Non si rifiuta l’ipotesi che le varianze dei gruppi siano uguali
t.test(db$Peso[db$Fumatrici_cat=="NF"],db$Peso[db$Fumatrici_cat=="F"], paired = FALSE, pool.sd = TRUE) #chiamiamo il t test con pool.sd=TRUE##
## Welch Two Sample t-test
##
## data: db$Peso[db$Fumatrici_cat == "NF"] and db$Peso[db$Fumatrici_cat == "F"]
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -45.61354 145.22674
## sample estimates:
## mean of x mean of y
## 3286.153 3236.346non si rifiuta l’ipotesi che le medie dei gruppi siano uguali
visivamente i boxplot della distribuzione del peso nei due gruppi sembrano simili. La distribuzione dei pesi per madri fumatrici è più irregolare e ha media e mediana inferiori. Tuttavia la differenza tra le medie non è statisticamente significativa.
2.3 Il fattore Tipo.parto
Vediamo come si distribuiscono i record in funzione del fattore Tipo.parto.
records_by_Tipo.parto=table(db$Tipo.parto)/n
bar_heights=
barplot(records_by_Tipo.parto,
main="Distribuzione dei record del database",
xlab="Tipo parto (Naturale / Cesareo)",
ylab="frequenze relative",
ylim = c(0, 1.1), # range asse y
col = "lightblue", # colore barre (opzionale)
border = "darkblue", # bordo barre (opzionale)
axes = FALSE # disattiva disegno automatico degli assi
)
# Aggiunge asse Y con tick personalizzati
axis(2, at = seq(0, 1, by = 0.1), las = 1) # las=1 => etichette orizzontali
# Aggiunge griglia orizzontale
abline(h = seq(0, 1, by = 0.1), col = "lightgray", lty = "dotted")
# Aggiunge etichette numeriche sopra ogni barra
text(x = bar_heights, y = records_by_Tipo.parto,
labels = round(records_by_Tipo.parto, 3),
pos = 3, cex = 0.8, col = "black")
I cesarei sono il 29.1% dei parti
Studiamo l’influenza del fattore Tipo.parto sulla variabile peso
print("statistiche Peso nati con parto Naturale")## [1] "statistiche Peso nati con parto Naturale"summary(db$Peso[db$Tipo.parto=="Nat"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 3000 3300 3285 3620 4900sd(db$Peso[db$Tipo.parto=="Nat"])## [1] 540.2184print('')## [1] ""print("statistiche Peso nati con Cesareo")## [1] "statistiche Peso nati con Cesareo"summary(db$Peso[db$Tipo.parto=="Ces"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 930 2980 3290 3282 3610 4930sd(db$Peso[db$Tipo.parto=="Ces"])## [1] 486.4645#boxplot(db$Peso~db$Tipo.parto)
g=ggplot(data=db)+
labs(title=paste("Distribuzioni di","Peso","in funzione di","Tipo di Parto"),
x="Tipo di parto",
y="Peso [g]")
g=g+
geom_half_boxplot(aes(x=Tipo.parto, y=Peso, fill=Tipo.parto), side="l" )+
geom_half_violin(aes(x=Tipo.parto, y=Peso, fill=Tipo.parto), side="r")+
stat_summary(aes(x = Tipo.parto, y = Peso, group = Tipo.parto,linetype = "Media"),
fun = mean,
geom = "crossbar",
width = 0.8,
fatten = 1,
color = "black") +
scale_linetype_manual(
name = "Statistiche", # Titolo della legenda
values = c("Media" = "dashed") # Etichetta + tipo di linea
)#+
#scale_fill_manual(values = set_colori)
g
leveneTest(Peso ~ Tipo.parto, data = db)#controlliamo che la varianza sia simile tra i gruppi## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 4.4896 0.0342 *
## 2498
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- si rifiuta l’ipotesi che le varianze dei gruppi siano uguali
t.test(db$Peso[db$Tipo.parto=="Nat"],db$Peso[db$Tipo.parto=="Ces"], paired = FALSE, pool.sd = FALSE
) #chiamiamo il t test con pool.sd=FALSE##
## Welch Two Sample t-test
##
## data: db$Peso[db$Tipo.parto == "Nat"] and db$Peso[db$Tipo.parto == "Ces"]
## t = 0.12968, df = 1493, p-value = 0.8968
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -40.54037 46.27992
## sample estimates:
## mean of x mean of y
## 3284.916 3282.047- non si rifiuta l’ipotesi che le medie dei gruppi siano uguali
- visivamente i boxplot della distribuzione del peso nei due gruppi sembrano equivalenti.
2.4 Il fattore Ospedale
Vediamo come si distribuiscono i parti tra i vari ospedali
records_by_Ospedale=table(db$Ospedale)/n
bar_heights=
barplot(records_by_Ospedale,
main="Distribuzione dei record del database",
xlab="Ospedale",
ylab="frequenze relative",
ylim = c(0, 1.1), # range asse y
col = "lightblue", # colore barre (opzionale)
border = "darkblue", # bordo barre (opzionale)
axes = FALSE # disattiva disegno automatico degli assi
)
# Aggiunge asse Y con tick personalizzati
axis(2, at = seq(0, 1, by = 0.1), las = 1) # las=1 => etichette orizzontali
# Aggiunge griglia orizzontale
abline(h = seq(0, 1, by = 0.1), col = "lightgray", lty = "dotted")
# Aggiunge etichette numeriche sopra ogni barra
text(x = bar_heights, y = records_by_Ospedale,
labels = round(records_by_Ospedale, 3),
pos = 3, cex = 0.8, col = "black")
I record sono equamente distribuiti tra i tre ospedali.
Studiamo l’influenza del fattore Ospedale sulla variabile peso
#boxplot(db$Peso~db$Ospedale)
leveneTest(Peso ~ Ospedale, data = db)#controlliamo che la varianza sia simile tra i gruppi## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.5892 0.5548
## 2497g=ggplot(data=db)+
labs(title=paste("Distribuzioni di","Peso","in funzione di","Ospedale"),
x="Ospedale",
y="Peso [g]")
g=g+
geom_half_boxplot(aes(x=Ospedale, y=Peso, fill=Ospedale), side="l" )+
geom_half_violin(aes(x=Ospedale, y=Peso, fill=Ospedale), side="r")+
stat_summary(aes(x = Ospedale, y = Peso, group = Ospedale,linetype = "Media"),
fun = mean,
geom = "crossbar",
width = 0.8,
fatten = 1,
color = "black") +
scale_linetype_manual(
name = "Statistiche", # Titolo della legenda
values = c("Media" = "dashed") # Etichetta + tipo di linea
)#+
#scale_fill_manual(values = set_colori)
g
pairwise.t.test(db$Peso, db$Ospedale, paired = FALSE, pool.sd = TRUE,
p.adjust.method = "bonferroni") #chiamiamo il t test a coppie di gruppi con pool.sd=TRUE##
## Pairwise comparisons using t tests with pooled SD
##
## data: db$Peso and db$Ospedale
##
## osp1 osp2
## osp2 1.00 -
## osp3 0.33 0.33
##
## P value adjustment method: bonferroni- non si rifiuta l’ipotesi che le varianze dei gruppi siano uguali
- non si rifiuta l’ipotesi che le medie dei gruppi siano uguali
- visivamente i boxplot della distribuzione del peso nei tre gruppi sembrano equivalenti
2.4.1 I parti cesarei sono distribuiti equamente tra gli ospedali?
Vediamo come si distribuiscono i parti cesarei tra gli ospedali. In particolare siamo interessati a confrontare la frequenza dei cesarei per ciascun ospedale.
Per comodità trasformiamo in una variabile binaria il fattore Tipo.parto.
db$Tipo.parto_num=ifelse(db$Tipo.parto=="Ces",1,0)Calcoliamo la frequenza dei cesarei nei tre ospedali
for (osp in c("osp1", "osp2", "osp3")) {
cat("Frequenza parti cesarei in", osp, ":",
mean(db$Tipo.parto_num[db$Ospedale == osp]), "\n")
}## Frequenza parti cesarei in osp1 : 0.2965686
## Frequenza parti cesarei in osp2 : 0.2991755
## Frequenza parti cesarei in osp3 : 0.2778443leveneTest(Tipo.parto_num ~ Ospedale, data = db)#controlliamo che la varianza sia simile tra i gruppi## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.5482 0.5781
## 2497- non si rifiuta l’ipotesi di varianza uguale tra i gruppi
pairwise.t.test(db$Tipo.parto_num, db$Ospedale, paired = FALSE, pool.sd = TRUE,
p.adjust.method = "bonferroni") #chiamiamo il t test a coppie di gruppi con pool.sd=TRUE##
## Pairwise comparisons using t tests with pooled SD
##
## data: db$Tipo.parto_num and db$Ospedale
##
## osp1 osp2
## osp2 1 -
## osp3 1 1
##
## P value adjustment method: bonferroni- la differenza tra l’incidenza dei cesarei nei vari ospedali non è statisticamente significativa, non si rifiuta l’ipotesi di incidenza uniforme
2.5 Il fattore Sesso
Vediamo come si distribuiscono i parti col sesso del neonato
records_by_Sesso=table(db$Sesso)/n
bar_heights=
barplot(records_by_Sesso,
main="Distribuzione dei record del database",
xlab="Sesso del neonato (Maschio / Femmina)",
ylab="frequenze relative",
ylim = c(0, 1.1), # range asse y
col = "lightblue", # colore barre (opzionale)
border = "darkblue", # bordo barre (opzionale)
axes = FALSE # disattiva disegno automatico degli assi
)
# Aggiunge asse Y con tick personalizzati
axis(2, at = seq(0, 1, by = 0.1), las = 1) # las=1 => etichette orizzontali
# Aggiunge griglia orizzontale
abline(h = seq(0, 1, by = 0.1), col = "lightgray", lty = "dotted")
# Aggiunge etichette numeriche sopra ogni barra
text(x = bar_heights, y = records_by_Sesso,
labels = round(records_by_Sesso, 3),
pos = 3, cex = 0.8, col = "black")
Come prevedibile, il sesso dei neonati è equamente distribuito tra femmine e maschi.
Studiamo l’influenza del fattore Sesso sulla variabile peso
#boxplot(db$Peso~db$Ospedale)
leveneTest(Peso ~ Sesso, data = db)#controlliamo che la varianza sia simile tra i gruppi## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.7931 0.3732
## 2498g=ggplot(data=db)+
labs(title=paste("Distribuzioni di","Peso","in funzione di","Sesso"),
x="Sesso [M/F]",
y="Peso [g]")
g=g+
geom_half_boxplot(aes(x=Sesso, y=Peso, fill=Sesso), side="l" )+
geom_half_violin(aes(x=Sesso, y=Peso, fill=Sesso), side="r")+
stat_summary(aes(x = Sesso, y = Peso, group = Sesso,linetype = "Media"),
fun = mean,
geom = "crossbar",
width = 0.8,
fatten = 1,
color = "black") +
scale_linetype_manual(
name = "Statistiche", # Titolo della legenda
values = c("Media" = "dashed") # Etichetta + tipo di linea
)#+
#scale_fill_manual(values = set_colori)
g
1) Non si rifiuta l’ipotesi di varianza uguale tra i gruppi
t.test(db$Peso[db$Sesso=="F"], db$Peso[db$Sesso=="M"], paired = FALSE, pool.sd = TRUE ) #chiamiamo il t test a coppie di gruppi con pool.sd=TRUE##
## Welch Two Sample t-test
##
## data: db$Peso[db$Sesso == "F"] and db$Peso[db$Sesso == "M"]
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -287.1051 -207.0615
## sample estimates:
## mean of x mean of y
## 3161.132 3408.215si rifiuta, con p-value praticamente nullo, l’ipotesi che le due medie siano uguali
anche graficamente le due distribuzioni sono nettamente sfalsate. La differenza (statisticamente significativa) tra le medie delle due distribuzioni è, col 95% di confidenza, compresa tra 207 e 287 grammi a favore dei maschi.
2.5.1 Circonferenza del cranio Vs Il fattore Sesso
Vediamo come si distribuiscono le circonferenze craniali col sesso del neonato
#boxplot(db$Peso~db$Ospedale)
leveneTest(Cranio ~ Sesso, data = db)#controlliamo che la varianza sia simile tra i gruppi## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.1978 0.2739
## 2498g=ggplot(data=db)+
labs(title=paste("Distribuzioni di","Circonferenza craniale","in funzione di","Sesso"),
x="Sesso [M/F]",
y="Circonferenza Craniale [mm]")
g=g+
geom_half_boxplot(aes(x=Sesso, y=Cranio, fill=Sesso), side="l" )+
geom_half_violin(aes(x=Sesso, y=Cranio, fill=Sesso), side="r")+
stat_summary(aes(x = Sesso, y = Cranio, group = Sesso,linetype = "Media"),
fun = mean,
geom = "crossbar",
width = 0.8,
fatten = 1,
color = "black") +
scale_linetype_manual(
name = "Statistiche", # Titolo della legenda
values = c("Media" = "dashed") # Etichetta + tipo di linea
)#+
#scale_fill_manual(values = set_colori)
g
- non si rifiuta l’ipotesi di varianza omogenea
- graficamente la differenza tra le medie è netta, verifichiamolo con un t.test
t.test(db$Cranio[db$Sesso=="F"], db$Cranio[db$Sesso=="M"], paired = FALSE, pool.sd = TRUE ) #chiamiamo il t test sui due gruppi con pool.sd=TRUE##
## Welch Two Sample t-test
##
## data: db$Cranio[db$Sesso == "F"] and db$Cranio[db$Sesso == "M"]
## t = -7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -6.089912 -3.541270
## sample estimates:
## mean of x mean of y
## 337.6330 342.4486- le due medie sono significativamente diverse. La differenza di circonferenza craniale tra femmine e maschi è compresa, con una confidenza del 95%, tra 3.5 e 6mm a favore dei maschi.
2.5.2 Lunghezza del neonato Vs Il fattore Sesso
Vediamo come si distribuiscono le lunghezze dei neonati col sesso del neonato
#boxplot(db$Peso~db$Ospedale)
leveneTest(Lunghezza ~ Sesso, data = db)#controlliamo che la varianza sia simile tra i gruppi## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 10.496 0.001212 **
## 2498
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1g=ggplot(data=db)+
labs(title=paste("Distribuzioni di","Lunghezza del neonato","in funzione di","Sesso"),
x="Sesso [M/F]",
y="Lunghezza [mm]")
g=g+
geom_half_boxplot(aes(x=Sesso, y=Lunghezza, fill=Sesso), side="l" )+
geom_half_violin(aes(x=Sesso, y=Lunghezza, fill=Sesso), side="r")+
stat_summary(aes(x = Sesso, y = Lunghezza, group = Sesso,linetype = "Media"),
fun = mean,
geom = "crossbar",
width = 0.8,
fatten = 1,
color = "black") +
scale_linetype_manual(
name = "Statistiche", # Titolo della legenda
values = c("Media" = "dashed") # Etichetta + tipo di linea
)#+
#scale_fill_manual(values = set_colori)
g
- Si rifiuta l’ipotesi di varianza omogenea
- graficamente la differenza tra le medie è netta, verifichiamolo con un t.test
t.test(db$Lunghezza[db$Sesso=="F"], db$Lunghezza[db$Sesso=="M"], paired = FALSE, pool.sd = FALSE ) #chiamiamo il t test a coppie di gruppi con pool.sd=FALSE##
## Welch Two Sample t-test
##
## data: db$Lunghezza[db$Sesso == "F"] and db$Lunghezza[db$Sesso == "M"]
## t = -9.582, df = 2459.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -11.929470 -7.876273
## sample estimates:
## mean of x mean of y
## 489.7643 499.6672- le due medie sono significativamente diverse. La differenza media di lunghezza tra femmine e maschi è compresa, con una confidenza del 95%, tra 7.9 e 11.9 mm a favore dei maschi.
Concludendo, peso, circonferenza craniale e lunghezza dei neonati sono significativamente diverse tra femmine e maschi.
2.6 CONCLUSIONI
Il fattore sesso va sicuramente incluso nel modello di predizione del peso. Gli altri fattori non sono così influenti da generare medie significativamente diverse nei sottogruppi che individuano, ma non è da escludere che la regressione lineare individui dei coefficienti significativamente non nulli. Non è nemmeno escluso che i fattori possano avere interazione tra loro o con le variabili quantitative, pertanto, nella ricerca del modello ottimo, si partirà comunque da un modello completo che li include tutti.
3 STUDIO DELLA CORRELAZIONE TRA LE VARIABILI QUANTITATIVE
Dopo aver passato in rassegna gli effetti delle variabili qualitative sul peso, esaminiamo la correlazione tra le variabili quantitative. Utilizziamo la matrice di correlazione, che correla ogni variabile (sulle righe) con ogni altra variabile (sulle colonne).
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
#Definiamo il dataframe sottoinsieme di db, formato dalle sole variabili quantitative
db2= db[, c("Peso", "Anni.madre", "Gestazione", "Lunghezza", "Cranio")]
#usiamo la funzione pairs
pairs(db2, lower.panel = panel.smooth, upper.panel = panel.cor, gap=0, row1attop=FALSE)
Il grafico a matrice ha sulla diagonale secondaria le variabili
quantitative. Nella metà superiore sinistra mostra, nella cella
all’intersezione tra riga e colonna di due diverse variabili, il
coefficiente di correlazione (tra -1 e 1) tra le due variabili, per
comodità ingrandendo il corpo del carattere proporzionalmente al valore
del coefficiente. Nella metà inferiore destra, si mostra invece lo
scatterplot dei record, rappresentati nel piano xy con x = la variabile
corrispondente alla colonna e y = la variabile corrispondente alla
riga.
Possiamo ottenere i coefficienti di correlazione lineare anche usando la funzione cor() sul database con sole variabili quantitative
cor(db2)## Peso Anni.madre Gestazione Lunghezza Cranio
## Peso 1.00000000 -0.02247017 0.5917687 0.79603676 0.7048015
## Anni.madre -0.02247017 1.00000000 -0.1356671 -0.06349157 0.0160767
## Gestazione 0.59176872 -0.13566714 1.0000000 0.61892045 0.4608289
## Lunghezza 0.79603676 -0.06349157 0.6189204 1.00000000 0.6033410
## Cranio 0.70480151 0.01607670 0.4608289 0.60334097 1.0000000La variabile peso è correlata positivamente con Lunghezza (rxy=0.80), Cranio (rxy=0.70), Gestazione (rxy=0.59), in ordine decrescente di correlazione. Non c’è correlazione significativa con Anni.madre (rxy=-0.02). Le curve rosse di regressione appaiono lineari o quasi. Lunghezza, Cranio e Gestazione sono correlate positivamente tra loro, come era prevedibile, infatti, all’aumentare delle settimane di gestazione al momento del parto, tendenzialmente aumenta sia la lunghezza (rxy= 0.62), sia - leggermente meno - la circonferenza craniale (rxy=0.46) (il bambino ha avuto più tempo per crescere). Lunghezza e Cranio sono anch’esse correlate positivamente tra loro (rxy=0.60), sia perché entrambe aumentano con Gestazione, sia perché un bambino che è più lungo di un altro, tendenzialmente avrà anche una circonferenza craniale maggiore, “proporzionata” alla lunghezza maggiore.
Bisognerà valutare eventuali problemi di collinearità tra le variabili, ed eventualemente escluderne qualcuna.
Osservando i grafici si notano alcuni punti “problematici”.
- Nei grafici con Anni.madre ci sono dei punti vicini allo 0, cosa ovviamente impossibile.
Cerchiamo nel database questi punti
db[db$Anni.madre<5,]## 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 grav_cat Fumatrici_cat Tipo.parto_num
## 1152 Nat osp2 F 1 NF 0
## 1380 Nat osp3 M 0 NF 0si tratta delle righe 1152 (Anni.madre = 1) e 1380 (Anni.madre=0). Le righe, a parte la variabile Anni.madre, sembrano valide, quindi non eliminiamole, ma rimpiazziamo l’età della madre con un valore “probabile”.
Proviamo a vedere se l’età della madre è correlata al numero di gravidanze
g=ggplot(data=db)+
labs(title=paste("Distribuzioni di","Anni.madre","in funzione di","grav_cat"),
x="Numero gravidanze precedenti",
y="Età madre [anni]")
g=g+
geom_half_boxplot(aes(x=grav_cat, y=Anni.madre, fill=grav_cat), side="l" )+
geom_half_violin(aes(x=grav_cat, y=Anni.madre, fill=grav_cat), side="r")+
stat_summary(aes(x = grav_cat, y = Anni.madre, group = grav_cat,linetype = "Media"),
fun = mean,
geom = "crossbar",
width = 0.8,
fatten = 1,
color = "black") +
scale_linetype_manual(
name = "Statistiche", # Titolo della legenda
values = c("Media" = "dashed") # Etichetta + tipo di linea
)#+
#scale_fill_manual(values = set_colori)
g
leveneTest(Anni.madre ~ grav_cat, data = db)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 4 0.1063 0.9804
## 2495pairwise.t.test(db$Anni.madre, db$grav_cat, paired = FALSE, pool.sd = TRUE,
p.adjust.method = "bonferroni")##
## Pairwise comparisons using t tests with pooled SD
##
## data: db$Anni.madre and db$grav_cat
##
## 0 1 2 3
## 1 < 2e-16 - - -
## 2 < 2e-16 2.5e-09 - -
## 3 < 2e-16 < 2e-16 0.0015 -
## 4+ < 2e-16 < 2e-16 7.5e-08 0.2239
##
## P value adjustment method: bonferroniCome valore di rimpiazzo per Anni.madre è ragionevole usare l’età media di una madre per quel dato numero di gravidanze precedenti.
db[1152,"Anni.madre"]=mean(db$Anni.madre[db$grav_cat==db[1152,"grav_cat"]])
db[1380,"Anni.madre"]=mean(db$Anni.madre[db$grav_cat==db[1380,"grav_cat"]])Cerchiamo poi altri due outlier nella colonna dei grafici delle varabile lunghezza
db[(db$Lunghezza<350 & db$Gestazione>33)|(db$Lunghezza<350 & db$Peso>3500),]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1551 35 1 0 38 4370 315 374
## Tipo.parto Ospedale Sesso grav_cat Fumatrici_cat Tipo.parto_num
## 1551 Nat osp3 F 1 NF 0Il punto individuato è unico, Gestazione e Peso sono coerenti, è plausibile ci sia stato un errore nella registrazione della lunghezza. Dall’osservazione dei grafici (Cranio - Lunghezza, Peso - Lunghezza), si ipotizza (supponendo una sola cifra errata) un valore originario di 515 mm.
db[1551,"Lunghezza"]=515Individuiamo poi i due punti outlier con Cranio > 330mm e Gestazione < 30, e i due punti con Cranio > 330mm e Peso < 2000g.
db[(db$Cranio>330 & db$Gestazione<30)|(db$Cranio>330 & db$Peso<2000),]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 310 40 3 0 28 1560 420 379
## 1429 24 4 0 29 1280 390 355
## Tipo.parto Ospedale Sesso grav_cat Fumatrici_cat Tipo.parto_num
## 310 Nat osp3 F 3 NF 0
## 1429 Nat osp1 F 4+ NF 0I punti problematici sono gli stessi due, per i quali stavolta Gestazione e Peso sono coerenti e verosimilmente è errata la registrazione della circonferenza craniale. Cranio correla principalmente con Peso. Osservando il grafico Cranio-Peso, e nell’ipotesi di una sola cifra errata, rimpiazzeremo il 379mm con 279mm e il 355mm con 255mm.
db[310,"Cranio"]=279
db[1429,"Cranio"]=255Visualizziamo di nuovo la matrice di correlazione
#Aggiorniamo il database solo quantitativo db2
db2= db[, c("Peso", "Anni.madre", "Gestazione", "Lunghezza", "Cranio")]
#usiamo la funzione pairs
pairs(db2, lower.panel = panel.smooth, upper.panel = panel.cor,
gap=0, row1attop=FALSE)
cor(db2)## Peso Anni.madre Gestazione Lunghezza Cranio
## Peso 1.00000000 -0.02371685 0.5917687 0.80982114 0.71722595
## Anni.madre -0.02371685 1.00000000 -0.1348876 -0.06153161 0.01265594
## Gestazione 0.59176872 -0.13488757 1.0000000 0.62309183 0.48486121
## Lunghezza 0.80982114 -0.06153161 0.6230918 1.00000000 0.62783308
## Cranio 0.71722595 0.01265594 0.4848612 0.62783308 1.00000000La matrice di correlazione mostra dei lievi miglioramenti (coefficienti di correlazione aumentati), dopo le correzioni sugli outlier errati.
4 MODELLO SENZA INTERAZIONI
Definiamo le variabili qualitative Tipo.parto, Ospedale, Sesso come fattori
for (var in c("Tipo.parto", "Ospedale", "Sesso")){
db[[var]]=factor(db[[var]])
}
str(db)## 'data.frame': 2500 obs. of 13 variables:
## $ Anni.madre : num 26 21 34 28 20 32 26 25 22 23 ...
## $ N.gravidanze : int 0 2 3 1 0 0 1 0 1 0 ...
## $ Fumatrici : int 0 0 0 0 0 0 0 0 0 0 ...
## $ 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 : num 490 490 500 515 480 495 480 510 500 510 ...
## $ Cranio : num 325 345 375 365 335 340 345 349 335 362 ...
## $ Tipo.parto : Factor w/ 2 levels "Ces","Nat": 2 2 2 2 2 2 2 2 1 1 ...
## $ Ospedale : Factor w/ 3 levels "osp1","osp2",..: 3 1 2 2 3 2 3 1 2 2 ...
## $ Sesso : Factor w/ 2 levels "F","M": 2 1 2 2 1 1 1 2 1 1 ...
## $ grav_cat : Factor w/ 5 levels "0","1","2","3",..: 1 3 4 2 1 1 2 1 2 1 ...
## $ Fumatrici_cat : Factor w/ 2 levels "NF","F": 1 1 1 1 1 1 1 1 1 1 ...
## $ Tipo.parto_num: num 0 0 0 0 0 0 0 0 1 1 ...Prepariamo il dataframe db3 con tutte le variabili quantitative e i fattori.
db3=db[, c("Peso", "Anni.madre", "Gestazione", "Lunghezza", "Cranio", "Tipo.parto", "Ospedale", "Sesso", "grav_cat", "Fumatrici_cat")]Costruiamo il modello mod1, completo di tutte le variabili e di tutti i fattori, per ora senza interazioni e senza variabili al quadrato.
mod1=lm(Peso ~ Anni.madre + Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso + grav_cat + Fumatrici_cat, data=db3)
summary(mod1)##
## Call:
## lm(formula = Peso ~ Anni.madre + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + grav_cat + Fumatrici_cat,
## data = db3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1156.87 -180.12 -14.23 160.92 1396.67
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6665.1314 136.8932 -48.689 < 2e-16 ***
## Anni.madre 0.1300 1.1305 0.115 0.90844
## Gestazione 26.6798 3.7670 7.082 1.84e-12 ***
## Lunghezza 10.8462 0.3009 36.049 < 2e-16 ***
## Cranio 10.1817 0.4252 23.943 < 2e-16 ***
## Tipo.partoNat 31.1285 11.8068 2.636 0.00843 **
## Ospedaleosp2 -11.7757 13.1336 -0.897 0.37001
## Ospedaleosp3 25.8985 13.1876 1.964 0.04966 *
## SessoM 76.0751 10.9155 6.969 4.06e-12 ***
## grav_cat1 6.5198 12.7278 0.512 0.60852
## grav_cat2 50.7561 17.3412 2.927 0.00345 **
## grav_cat3 43.2616 24.3501 1.777 0.07575 .
## grav_cat4+ 76.0986 29.7713 2.556 0.01064 *
## Fumatrici_catF -28.9548 26.9352 -1.075 0.28249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.5 on 2486 degrees of freedom
## Multiple R-squared: 0.7417, Adjusted R-squared: 0.7403
## F-statistic: 549.1 on 13 and 2486 DF, p-value: < 2.2e-16Il modello mod1, cioè il set di coefficienti di regressione associati a variabili quantitative, all’intercetta, e alle modalità non base dei fattori è significativamente diverso da zero, con p-value praticamente nullo.
I coefficienti significativamente diversi da zero (con livelli di significatività diversi) sono quelli relativi all’intercetta e alle variabili Gestazione, Lunghezza, Cranio, e, tra i fattori, le modalità SessoM, Ospedaleosp3, grav_cat2 e grav_cat4+ (cioè le donne con due o 4+ gravidanze precedenti).
Il coefficiente diverso da zero di Ospedaleosp3 (al limite della soglia di significatività) indicherebbe un’influenza dell’ospedale sul peso del neonato, cosa ovviamente non plausibile. Al limite potrebbe indicare una non omogeneità delle madri che accedono a un dato ospedale (es. se donne con caratteristiche diverse avessero accesso diverso agli ospedali nelle diverse zone della città), ma il fattore andrà eliminato.
I coefficienti diversi da zero per le donne con 2 o 4+ gravidanze, sia pure significativi al test statistico, sono difficili da spiegare (donne con 1 o 3 gravidanze precedenti invece sarebbero non significativamente diverse dalle donne con 0 gravidanze precedenti). Probabilmente il fattore numero di gravidanze andrà anch’esso eliminato.
Il tipo di parto naturale ha un coefficiente positivo, ma per logica è più il parto a essere funzione delle condizioni di feto e madre, pertanto il fattore andrà eliminato.
Anni.madre e fumo non hanno un coefficiente significativamente diverso da zero, ma potrebbe esserci interazione tra loro.
R-squared adjusted è pari al 74%, un buon punto di partenza. Ottimizziamo mod1 aggiungendo o eliminando regressori per minimizzare il parametro BIC.
stepwise.mod1=MASS::stepAIC(mod1, direction="both", k=log(n))## Start: AIC=28042.02
## Peso ~ Anni.madre + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso + grav_cat + Fumatrici_cat
##
## Df Sum of Sq RSS AIC
## - grav_cat 4 1040443 178995428 28025
## - Anni.madre 1 947 177955932 28034
## - Ospedale 2 621141 178576127 28035
## - Fumatrici_cat 1 82720 178037706 28035
## - Tipo.parto 1 497579 178452565 28041
## <none> 177954985 28042
## - Sesso 1 3477015 181432000 28083
## - Gestazione 1 3590642 181545627 28084
## - Cranio 1 41037170 218992155 28553
## - Lunghezza 1 93023675 270978660 29086
##
## Step: AIC=28025.3
## Peso ~ Anni.madre + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso + Fumatrici_cat
##
## Df Sum of Sq RSS AIC
## - Fumatrici_cat 1 55502 179050930 28018
## - Ospedale 2 680426 179675854 28019
## - Anni.madre 1 179163 179174591 28020
## - Tipo.parto 1 450080 179445508 28024
## <none> 178995428 28025
## + grav_cat 4 1040443 177954985 28042
## - Gestazione 1 3477277 182472705 28066
## - Sesso 1 3588668 182584096 28067
## - Cranio 1 42219102 221214530 28547
## - Lunghezza 1 92544997 271540425 29059
##
## Step: AIC=28018.25
## Peso ~ Anni.madre + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 686492 179737422 28012
## - Anni.madre 1 176841 179227771 28013
## - Tipo.parto 1 444605 179495535 28017
## <none> 179050930 28018
## + Fumatrici_cat 1 55502 178995428 28025
## + grav_cat 4 1013225 178037706 28035
## - Gestazione 1 3437641 182488571 28058
## - Sesso 1 3575235 182626165 28060
## - Cranio 1 42230071 221281001 28540
## - Lunghezza 1 92959412 272010342 29056
##
## Step: AIC=28012.17
## Peso ~ Anni.madre + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 197238 179934661 28007
## - Tipo.parto 1 467010 180204432 28011
## <none> 179737422 28012
## + Ospedale 2 686492 179050930 28018
## + Fumatrici_cat 1 61568 179675854 28019
## + grav_cat 4 1070201 178667222 28029
## - Gestazione 1 3509339 183246762 28053
## - Sesso 1 3622751 183360174 28054
## - Cranio 1 42301285 222038708 28533
## - Lunghezza 1 92715469 272452892 29044
##
## Step: AIC=28007.09
## Peso ~ Gestazione + Lunghezza + Cranio + Tipo.parto + Sesso
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 463801 180398462 28006
## <none> 179934661 28007
## + Anni.madre 1 197238 179737422 28012
## + Ospedale 2 706889 179227771 28013
## + Fumatrici_cat 1 59117 179875544 28014
## + grav_cat 4 1265411 178669249 28021
## - Gestazione 1 3348621 183283282 28045
## - Sesso 1 3659258 183593918 28050
## - Cranio 1 43118062 223052723 28536
## - Lunghezza 1 92591578 272526239 29037
##
## Step: AIC=28005.7
## Peso ~ Gestazione + Lunghezza + Cranio + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 180398462 28006
## + Tipo.parto 1 463801 179934661 28007
## + Anni.madre 1 194029 180204432 28011
## + Ospedale 2 729412 179669050 28011
## + Fumatrici_cat 1 53479 180344983 28013
## + grav_cat 4 1214439 179184022 28020
## - Gestazione 1 3369006 183767468 28044
## - Sesso 1 3662330 184060792 28048
## - Cranio 1 43389526 223787988 28537
## - Lunghezza 1 92182868 272581330 29030Visualizziamo le statistiche del modello ottimizzato:
summary(stepwise.mod1)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso,
## data = db3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1157.5 -182.7 -15.9 163.6 1378.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6595.8870 131.1073 -50.309 < 2e-16 ***
## Gestazione 25.4275 3.7251 6.826 1.09e-11 ***
## Lunghezza 10.7626 0.3014 35.706 < 2e-16 ***
## Cranio 10.3690 0.4233 24.497 < 2e-16 ***
## SessoM 78.0044 10.9603 7.117 1.44e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.9 on 2495 degrees of freedom
## Multiple R-squared: 0.7381, Adjusted R-squared: 0.7377
## F-statistic: 1758 on 4 and 2495 DF, p-value: < 2.2e-16I fattori significativi che conviene mantenere per bilanciare la percentuale di varianza spiegata col modello e la complessità del modello sono rimasti i seguenti: - Intercetta - Gestazione - Lunghezza - Cranio - Sesso
Passiamo ai residui:
plot(stepwise.mod1)
Esaminiamo i quattro grafici diagnostici:
- Residuals vs Fitted Obiettivo: valutare linearità
I residui sembrano non perfettamente distribuiti a caso attorno a 0. C’è una leggera curvatura. Il modello sovrastima (residui positivi) i pesi dei (pochi) neonati molto piccoli. Probabilmente andranno aggiunte delle variabili (termini al quadrato e/o interazioni)
- Q-Q residuals Obiettivo: valutare se i residui sono normalmente distribuiti
La maggior parte dei punti è allineata (indicando normalità dei residui), ma le code si discostano dalla retta di riferimento (in alto a destra e in basso a sinistra).
- Scale-Location Obiettivo: verificare omoscedasticità (varianza costante dei residui)
La linea rossa indica una leggera eterogeneità della varianza (prima decrescente poi crescente) con i fitted values.
- Residuals vs Leverage Obiettivo: identificare osservazioni influenti
Vediamo alcuni punti ad alto leverage e residuo elevato, ma nessun punto oltre Cook’s distance (le curve soglia non sono nemmeno visibili in questo grafico), quindi nessuna osservazione è fortemente influente. In altre parole, non ci sono outlier preoccupanti in termini di influenza sul modello.
5 MODELLO CON INTERAZIONI E VARIABILI AL QUADRATO
Aggiungiamo al modello mod1 (Peso ~ Gestazione + Lunghezza + Cranio + Sesso) i termini quadratici delle due variabili quantitative Gestazione e Cranio (che mostrano una leggera curvatura), l’interazione a due vie con effetti principali tra Anni.madre e Fumatrici_cat, l’interazione semplice tra Gestazione e Fumatrici_cat, l’interazione semplice tra Lunghezza e Fumatrici_cat. Ma prima centriamo le variabili quantitative per diminuire la collinearità tra una variabile e il suo quadrato.
db3$Gestazione_cen=db3$Gestazione - mean(db3$Gestazione)
db3$Cranio_cen=db3$Cranio - mean(db3$Cranio)
db3$Lunghezza_cen=db3$Lunghezza - mean(db3$Lunghezza)mod2 =lm(Peso~ Gestazione_cen+ Lunghezza_cen + Cranio_cen+ Sesso + I(Gestazione_cen^2) + I(Lunghezza_cen^2)+I(Cranio_cen^2)+ Anni.madre*Fumatrici_cat +Gestazione_cen:Fumatrici_cat + Lunghezza_cen:Fumatrici_cat,
data = db3)
summary(mod2)##
## Call:
## lm(formula = Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen +
## Sesso + I(Gestazione_cen^2) + I(Lunghezza_cen^2) + I(Cranio_cen^2) +
## Anni.madre * Fumatrici_cat + Gestazione_cen:Fumatrici_cat +
## Lunghezza_cen:Fumatrici_cat, data = db3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1181.81 -181.97 -12.77 163.18 1308.96
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.182e+03 3.076e+01 103.442 < 2e-16 ***
## Gestazione_cen 3.512e+01 4.275e+00 8.215 3.38e-16 ***
## Lunghezza_cen 1.105e+01 3.094e-01 35.718 < 2e-16 ***
## Cranio_cen 1.055e+01 4.244e-01 24.859 < 2e-16 ***
## SessoM 7.181e+01 1.090e+01 6.590 5.34e-11 ***
## I(Gestazione_cen^2) -1.314e+00 9.363e-01 -1.403 0.16072
## I(Lunghezza_cen^2) 2.090e-02 5.071e-03 4.123 3.87e-05 ***
## I(Cranio_cen^2) 3.245e-02 1.461e-02 2.221 0.02647 *
## Anni.madre 1.720e+00 1.058e+00 1.626 0.10397
## Fumatrici_catF -6.343e+01 1.456e+02 -0.436 0.66307
## Anni.madre:Fumatrici_catF 2.449e+00 5.053e+00 0.485 0.62791
## Gestazione_cen:Fumatrici_catF -5.051e+01 2.162e+01 -2.337 0.01954 *
## Lunghezza_cen:Fumatrici_catF 4.083e+00 1.442e+00 2.832 0.00466 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 265.8 on 2487 degrees of freedom
## Multiple R-squared: 0.745, Adjusted R-squared: 0.7438
## F-statistic: 605.4 on 12 and 2487 DF, p-value: < 2.2e-16R quadro è migliorato di un solo punto percentuale (è ora 74%) Risulta significativa l’influenza del quadrato di Lunghezza_cen
Eliminiamo manualmente le variabili con coefficienti di regressione non significativi
mod2.1=update(mod2,~. -I(Gestazione_cen^2)-Anni.madre-Anni.madre:Fumatrici_cat )
summary(mod2.1)##
## Call:
## lm(formula = Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen +
## Sesso + I(Lunghezza_cen^2) + I(Cranio_cen^2) + Fumatrici_cat +
## Gestazione_cen:Fumatrici_cat + Lunghezza_cen:Fumatrici_cat,
## data = db3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1154.92 -180.06 -13.04 164.95 1327.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.230e+03 8.039e+00 401.839 < 2e-16 ***
## Gestazione_cen 3.635e+01 3.960e+00 9.181 < 2e-16 ***
## Lunghezza_cen 1.104e+01 3.095e-01 35.668 < 2e-16 ***
## Cranio_cen 1.067e+01 4.210e-01 25.353 < 2e-16 ***
## SessoM 7.167e+01 1.089e+01 6.581 5.67e-11 ***
## I(Lunghezza_cen^2) 1.792e-02 4.607e-03 3.888 0.000104 ***
## I(Cranio_cen^2) 2.505e-02 1.364e-02 1.836 0.066461 .
## Fumatrici_catF 6.775e+00 2.801e+01 0.242 0.808915
## Gestazione_cen:Fumatrici_catF -5.322e+01 2.152e+01 -2.474 0.013444 *
## Lunghezza_cen:Fumatrici_catF 3.901e+00 1.426e+00 2.736 0.006262 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 265.9 on 2490 degrees of freedom
## Multiple R-squared: 0.7444, Adjusted R-squared: 0.7435
## F-statistic: 805.8 on 9 and 2490 DF, p-value: < 2.2e-16Eliminiamo la variabile Cranio^2
mod2.2=update(mod2.1,~. -I(Cranio_cen^2))
summary(mod2.2)##
## Call:
## lm(formula = Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen +
## Sesso + I(Lunghezza_cen^2) + Fumatrici_cat + Gestazione_cen:Fumatrici_cat +
## Lunghezza_cen:Fumatrici_cat, data = db3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1162.17 -181.22 -14.45 165.40 1293.63
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.233e+03 7.894e+00 409.559 < 2e-16 ***
## Gestazione_cen 3.540e+01 3.927e+00 9.014 < 2e-16 ***
## Lunghezza_cen 1.105e+01 3.095e-01 35.695 < 2e-16 ***
## Cranio_cen 1.063e+01 4.207e-01 25.279 < 2e-16 ***
## SessoM 7.207e+01 1.089e+01 6.616 4.51e-11 ***
## I(Lunghezza_cen^2) 2.358e-02 3.424e-03 6.885 7.29e-12 ***
## Fumatrici_catF 5.589e+00 2.802e+01 0.199 0.84191
## Gestazione_cen:Fumatrici_catF -5.089e+01 2.149e+01 -2.368 0.01795 *
## Lunghezza_cen:Fumatrici_catF 3.903e+00 1.427e+00 2.736 0.00627 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266 on 2491 degrees of freedom
## Multiple R-squared: 0.7441, Adjusted R-squared: 0.7433
## F-statistic: 905.3 on 8 and 2491 DF, p-value: < 2.2e-16library(effects)## lattice theme set by effectsTheme()
## See ?effectsTheme for details.eff_lung_fum = Effect(c("Lunghezza_cen", "Fumatrici_cat"), mod2.2)
plot(eff_lung_fum)
Madri non fumatrici (NF): Il peso parte da circa 2000 g e arriva fino a
4000 g.
La curva è convessa (aumento accelerato), ma con un range più contenuto.
Madri fumatrici (F): Il peso parte da circa 1000 g e raggiunge circa 4300 g.
La curva è anche essa convessa, ma mediamente più ripida.
L’intervallo totale è maggiore (circa 3300 g contro i 2200 dei neonati delle madri non fumatrici).
Il peso nei figli di madri fumatrici è più sensibile alla lunghezza.
Per i neonati più piccoli (più corti e naturalmente con peso minore) il fattore fumo tende ad accentuare il sottopeso, quindi peggiora le condizione di salute dei neonati.
Passiamo ora all’interazione tra fattore fumo e settimane di gestazione
eff_gest_fum = Effect(c("Gestazione_cen", "Fumatrici_cat"), mod2.2)
plot(eff_gest_fum)
Per i neonati di madri non fumatrici la relazione è positiva e lineare:
all’aumentare di Gestazione_cen, il peso del neonato aumenta
chiaramente.
L’intervallo di confidenza (area azzurra) è stretto, il che indica una relazione robusta e stimata con buona precisione.
Per i neonati di madri fumatrici la linea di regressione ha una pendenza leggermente negativa: più giorni di gestazione sembrano associati a una leggera diminuzione del peso (cosa da escludere), e il risultato è altamente incerto.
L’intervallo di confidenza è molto ampio: vuol dire che non possiamo stimare bene l’effetto della gestazione sul peso nei neonati di madri fumatrici, probabilmente per pochi dati in quel gruppo (ricordiamo che le madri fumatrici sono una piccola minoranza del campione).
Nei figli di non fumatrici, più lunga è la gestazione, maggiore è il peso alla nascita, come atteso biologicamente.
Nei figli di fumatrici, l’effetto positivo della gestazione si perde (o addirittura diventa negativo), ma l’incertezza è elevata, per cui concludiamo che l’effetto è statisticamente non significativo.
L’intervallo di confidenza intorno alla retta di regressione per le non fumatrici è così ampio che in effetti ci sarebbe lo spazio per una curva ipotetica monotòna crescente ma meno ripida della curva delle non fumatrici.
Possiamo concludere che l’effetto positivo della durata della gestazione sul peso si osserva (in questo campione) solo nei figli di madri non fumatrici. Nei figli di fumatrici, la relazione è incerta e non significativa, suggerendo che il fumo potrebbe compromettere il normale accrescimento con le settimane di gestazione.
Eliminiamo le variabili superflue:
stepwise.mod2.2=MASS::stepAIC(mod2.2, direction="both", k=log(n))## Start: AIC=27979.61
## Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen + Sesso +
## I(Lunghezza_cen^2) + Fumatrici_cat + Gestazione_cen:Fumatrici_cat +
## Lunghezza_cen:Fumatrici_cat
##
## Df Sum of Sq RSS AIC
## - Gestazione_cen:Fumatrici_cat 1 396926 176701275 27977
## - Lunghezza_cen:Fumatrici_cat 1 529748 176834097 27979
## <none> 176304349 27980
## - Sesso 1 3097752 179402102 28015
## - I(Lunghezza_cen^2) 1 3355076 179659425 28019
## - Cranio_cen 1 45228555 221532904 28543
##
## Step: AIC=27977.41
## Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen + Sesso +
## I(Lunghezza_cen^2) + Fumatrici_cat + Lunghezza_cen:Fumatrici_cat
##
## Df Sum of Sq RSS AIC
## - Lunghezza_cen:Fumatrici_cat 1 229704 176930979 27973
## <none> 176701275 27977
## + Gestazione_cen:Fumatrici_cat 1 396926 176304349 27980
## - Sesso 1 3037452 179738727 28012
## - I(Lunghezza_cen^2) 1 3327501 180028776 28016
## - Gestazione_cen 1 5406622 182107897 28045
## - Cranio_cen 1 45331823 222033098 28541
##
## Step: AIC=27972.83
## Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen + Sesso +
## I(Lunghezza_cen^2) + Fumatrici_cat
##
## Df Sum of Sq RSS AIC
## - Fumatrici_cat 1 34715 176965694 27966
## <none> 176930979 27973
## + Lunghezza_cen:Fumatrici_cat 1 229704 176701275 27977
## + Gestazione_cen:Fumatrici_cat 1 96882 176834097 27979
## - Sesso 1 3113651 180044630 28009
## - I(Lunghezza_cen^2) 1 3414004 180344983 28013
## - Gestazione_cen 1 5362659 182293638 28040
## - Cranio_cen 1 45296523 222227502 28535
## - Lunghezza_cen 1 95140055 272071035 29041
##
## Step: AIC=27965.5
## Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen + Sesso +
## I(Lunghezza_cen^2)
##
## Df Sum of Sq RSS AIC
## <none> 176965694 27966
## + Fumatrici_cat 1 34715 176930979 27973
## - Sesso 1 3102784 180068478 28001
## - I(Lunghezza_cen^2) 1 3432767 180398462 28006
## - Gestazione_cen 1 5333710 182299405 28032
## - Cranio_cen 1 45307983 222273677 28528
## - Lunghezza_cen 1 95582988 272548682 29037summary(stepwise.mod2.2)##
## Call:
## lm(formula = Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen +
## Sesso + I(Lunghezza_cen^2), data = db3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1167.91 -182.63 -13.94 164.35 1295.79
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.232e+03 7.822e+00 413.215 < 2e-16 ***
## Gestazione_cen 3.356e+01 3.871e+00 8.670 < 2e-16 ***
## Lunghezza_cen 1.120e+01 3.052e-01 36.702 < 2e-16 ***
## Cranio_cen 1.064e+01 4.212e-01 25.269 < 2e-16 ***
## SessoM 7.202e+01 1.089e+01 6.613 4.60e-11 ***
## I(Lunghezza_cen^2) 2.381e-02 3.424e-03 6.955 4.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.4 on 2494 degrees of freedom
## Multiple R-squared: 0.7431, Adjusted R-squared: 0.7426
## F-statistic: 1443 on 5 and 2494 DF, p-value: < 2.2e-16L’ottimizzazione BIC ha eliminato gli effetti dell’interazione fumo - lunghezza.
Controlliamo la collinearità delle variabili
vif(stepwise.mod2.2)## Gestazione_cen Lunghezza_cen Cranio_cen Sesso
## 1.842862 2.231165 1.708287 1.044878
## I(Lunghezza_cen^2)
## 1.519015I vif sono tutti inferiori a 5, quindi non ci sono evidenti problemi di collinearità.
BIC(stepwise.mod1, mod2.2, stepwise.mod2.2)## df BIC
## stepwise.mod1 6 35108.22
## mod2.2 10 35082.12
## stepwise.mod2.2 7 35068.01Aver introdotto la variabile Lunghezza_cen^2 ha migliorato (sia pur di poco) il modello.
Controlliamo se la variabile Fumatrici_cat ha un effetto di interazione con Sesso
mod3=update(stepwise.mod2.2, ~. + Fumatrici_cat+ Sesso:Fumatrici_cat)
summary(mod3)##
## Call:
## lm(formula = Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen +
## Sesso + I(Lunghezza_cen^2) + Fumatrici_cat + Sesso:Fumatrici_cat,
## data = db3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1169.04 -182.64 -14.13 164.15 1295.60
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.232e+03 7.970e+00 405.529 < 2e-16 ***
## Gestazione_cen 3.372e+01 3.877e+00 8.698 < 2e-16 ***
## Lunghezza_cen 1.119e+01 3.057e-01 36.614 < 2e-16 ***
## Cranio_cen 1.064e+01 4.212e-01 25.262 < 2e-16 ***
## SessoM 7.359e+01 1.111e+01 6.625 4.24e-11 ***
## I(Lunghezza_cen^2) 2.377e-02 3.426e-03 6.937 5.08e-12 ***
## Fumatrici_catF -7.013e-03 3.888e+01 0.000 1.000
## SessoM:Fumatrici_catF -3.545e+01 5.348e+01 -0.663 0.508
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.4 on 2492 degrees of freedom
## Multiple R-squared: 0.7432, Adjusted R-squared: 0.7425
## F-statistic: 1030 on 7 and 2492 DF, p-value: < 2.2e-16I coefficienti di regressione associati all’interazione Fumatrici_cat*Sesso non sono significativamente diversi da zero, quindi l’interazione verrà eliminata dal modello.
eff_sex_fum = Effect(c("Sesso", "Fumatrici_cat"), mod3)
plot(eff_sex_fum)
eff_sex_fum##
## Sesso*Fumatrici_cat effect
## Fumatrici_cat
## Sesso NF F
## F 3232.089 3232.082
## M 3305.684 3270.230Con madri fumatrici 1) La variabilità aumenta, le due distribuzioni sono più sovrapposte 2) Per i maschi il peso atteso scende leggermente
La variabilità è troppa per considerare questo effetto statisticamente significativo, almeno con il campione oggetto di studio
Testiamo una possibile interazione tra Fumatrici_cat e grav_cat (numero di gravidanze precedenti)
mod4=update(stepwise.mod2.2, ~. + Fumatrici_cat*grav_cat)
summary(mod4)##
## Call:
## lm(formula = Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen +
## Sesso + I(Lunghezza_cen^2) + Fumatrici_cat + grav_cat + Fumatrici_cat:grav_cat,
## data = db3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1206.7 -178.8 -13.8 160.2 1312.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.215e+03 1.002e+01 320.785 < 2e-16 ***
## Gestazione_cen 3.562e+01 3.914e+00 9.101 < 2e-16 ***
## Lunghezza_cen 1.123e+01 3.055e-01 36.758 < 2e-16 ***
## Cranio_cen 1.048e+01 4.233e-01 24.753 < 2e-16 ***
## SessoM 7.072e+01 1.088e+01 6.501 9.60e-11 ***
## I(Lunghezza_cen^2) 2.422e-02 3.424e-03 7.076 1.92e-12 ***
## Fumatrici_catF 1.094e+01 4.842e+01 0.226 0.821236
## grav_cat1 1.251e+01 1.269e+01 0.986 0.324252
## grav_cat2 5.992e+01 1.697e+01 3.530 0.000423 ***
## grav_cat3 5.682e+01 2.397e+01 2.371 0.017830 *
## grav_cat4+ 8.340e+01 2.945e+01 2.832 0.004667 **
## Fumatrici_catF:grav_cat1 -1.949e+01 6.485e+01 -0.301 0.763805
## Fumatrici_catF:grav_cat2 -8.005e+01 8.196e+01 -0.977 0.328866
## Fumatrici_catF:grav_cat3 -1.276e+02 1.035e+02 -1.232 0.217981
## Fumatrici_catF:grav_cat4+ -7.448e+01 1.150e+02 -0.648 0.517269
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 265.7 on 2485 degrees of freedom
## Multiple R-squared: 0.7454, Adjusted R-squared: 0.744
## F-statistic: 519.7 on 14 and 2485 DF, p-value: < 2.2e-16I coefficienti relativi ai gruppi 2, 3 e 4+ gravidanze precedenti risultano diversi da zero (con livelli di significatività diversi), ma non è detto che aggiungere questi termini migliori il modello. I coefficienti relativi all’interazione tra gravidanze e fumo sono tutti negativi (il fumo causerebbe una diminuzione di peso), ma il test non supera la soglia di significatività.
eff_grav_fum = Effect(c("grav_cat", "Fumatrici_cat"), mod4)
plot(eff_grav_fum)
eff_grav_fum##
## grav_cat*Fumatrici_cat effect
## Fumatrici_cat
## grav_cat NF F
## 0 3249.831 3260.773
## 1 3262.342 3253.796
## 2 3309.747 3240.643
## 3 3306.653 3190.005
## 4+ 3333.234 3269.694Per le madri non fumatrici il peso tenderebbe a salire col numero di gravidanze precedenti, mentre l’effetto sembra svanire per le madri fumatrici (anche qui la variabilità è troppa per avere significatività statistica).
Ottimizziamo il modello mod4
stepwise.mod4=MASS::stepAIC(mod4, direction="both", k=log(n))## Start: AIC=28013.56
## Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen + Sesso +
## I(Lunghezza_cen^2) + Fumatrici_cat + grav_cat + Fumatrici_cat:grav_cat
##
## Df Sum of Sq RSS AIC
## - Fumatrici_cat:grav_cat 4 160906 175551348 27985
## <none> 175390442 28014
## - Sesso 1 2983004 178373446 28048
## - I(Lunghezza_cen^2) 1 3533731 178924173 28056
## - Gestazione_cen 1 5846060 181236502 28088
## - Cranio_cen 1 43244603 218635045 28557
## - Lunghezza_cen 1 95362208 270752650 29091
##
## Step: AIC=27984.55
## Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen + Sesso +
## I(Lunghezza_cen^2) + Fumatrici_cat + grav_cat
##
## Df Sum of Sq RSS AIC
## - grav_cat 4 1379631 176930979 27973
## - Fumatrici_cat 1 62452 175613800 27978
## <none> 175551348 27985
## + Fumatrici_cat:grav_cat 4 160906 175390442 28014
## - Sesso 1 2959508 178510857 28019
## - I(Lunghezza_cen^2) 1 3549592 179100940 28027
## - Gestazione_cen 1 5794713 181346062 28058
## - Cranio_cen 1 43283122 218834470 28528
## - Lunghezza_cen 1 95873885 271425234 29066
##
## Step: AIC=27972.83
## Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen + Sesso +
## I(Lunghezza_cen^2) + Fumatrici_cat
##
## Df Sum of Sq RSS AIC
## - Fumatrici_cat 1 34715 176965694 27966
## <none> 176930979 27973
## + grav_cat 4 1379631 175551348 27985
## - Sesso 1 3113651 180044630 28009
## - I(Lunghezza_cen^2) 1 3414004 180344983 28013
## - Gestazione_cen 1 5362659 182293638 28040
## - Cranio_cen 1 45296523 222227502 28535
## - Lunghezza_cen 1 95140055 272071035 29041
##
## Step: AIC=27965.5
## Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen + Sesso +
## I(Lunghezza_cen^2)
##
## Df Sum of Sq RSS AIC
## <none> 176965694 27966
## + Fumatrici_cat 1 34715 176930979 27973
## + grav_cat 4 1351894 175613800 27978
## - Sesso 1 3102784 180068478 28001
## - I(Lunghezza_cen^2) 1 3432767 180398462 28006
## - Gestazione_cen 1 5333710 182299405 28032
## - Cranio_cen 1 45307983 222273677 28528
## - Lunghezza_cen 1 95582988 272548682 29037summary(stepwise.mod4)##
## Call:
## lm(formula = Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen +
## Sesso + I(Lunghezza_cen^2), data = db3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1167.91 -182.63 -13.94 164.35 1295.79
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.232e+03 7.822e+00 413.215 < 2e-16 ***
## Gestazione_cen 3.356e+01 3.871e+00 8.670 < 2e-16 ***
## Lunghezza_cen 1.120e+01 3.052e-01 36.702 < 2e-16 ***
## Cranio_cen 1.064e+01 4.212e-01 25.269 < 2e-16 ***
## SessoM 7.202e+01 1.089e+01 6.613 4.60e-11 ***
## I(Lunghezza_cen^2) 2.381e-02 3.424e-03 6.955 4.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.4 on 2494 degrees of freedom
## Multiple R-squared: 0.7431, Adjusted R-squared: 0.7426
## F-statistic: 1443 on 5 and 2494 DF, p-value: < 2.2e-16La versione ottimizzata di mod4, stepwise.mod4 si è ridotta alla versione stepwise.mod2.
6 SELEZIONE MODELLO
Confrontiamo i modelli fin qui sviluppati
BIC(stepwise.mod1, mod2.2, stepwise.mod2.2, mod3, mod4, stepwise.mod4)## df BIC
## stepwise.mod1 6 35108.22
## mod2.2 10 35082.12
## stepwise.mod2.2 7 35068.01
## mod3 9 35082.73
## mod4 16 35116.08
## stepwise.mod4 7 35068.01models <- list(stepwise.mod1=stepwise.mod1, mod2.2=mod2.2, stepwise.mod2.2=stepwise.mod2.2, mod3=mod3, mod4=mod4, stepwise.mod4=stepwise.mod4)
tibble(
modello = names(models),
R2 = map_dbl(models, ~ summary(.x)$r.squared),
adj_R2 = map_dbl(models, ~ summary(.x)$adj.r.squared),
RMSE = map_dbl(models, ~ sqrt(mean(residuals(.x)^2))),
BIC=map_dbl(models, ~ BIC(.x))
)## # A tibble: 6 × 5
## modello R2 adj_R2 RMSE BIC
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 stepwise.mod1 0.738 0.738 269. 35108.
## 2 mod2.2 0.744 0.743 266. 35082.
## 3 stepwise.mod2.2 0.743 0.743 266. 35068.
## 4 mod3 0.743 0.742 266. 35083.
## 5 mod4 0.745 0.744 265. 35116.
## 6 stepwise.mod4 0.743 0.743 266. 35068.I modelli hanno tutti performance simili. La quota di variabilità spiegata è di circa il 74%. Rimane una quota del 26% di variabilità residua. Scegliamo il modello con BIC minore.
summary(stepwise.mod4)##
## Call:
## lm(formula = Peso ~ Gestazione_cen + Lunghezza_cen + Cranio_cen +
## Sesso + I(Lunghezza_cen^2), data = db3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1167.91 -182.63 -13.94 164.35 1295.79
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.232e+03 7.822e+00 413.215 < 2e-16 ***
## Gestazione_cen 3.356e+01 3.871e+00 8.670 < 2e-16 ***
## Lunghezza_cen 1.120e+01 3.052e-01 36.702 < 2e-16 ***
## Cranio_cen 1.064e+01 4.212e-01 25.269 < 2e-16 ***
## SessoM 7.202e+01 1.089e+01 6.613 4.60e-11 ***
## I(Lunghezza_cen^2) 2.381e-02 3.424e-03 6.955 4.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.4 on 2494 degrees of freedom
## Multiple R-squared: 0.7431, Adjusted R-squared: 0.7426
## F-statistic: 1443 on 5 and 2494 DF, p-value: < 2.2e-16plot(stepwise.mod4)
Aver introdotto Lunghezza_cen^2 ha migliorato la linearità della nuvola dei residui.
La normalità dei residui non è migliorata (le code si discostano ancora dalla diagonale).
L’omoscedasticità è migliorata.
La distanza di Cook dei punti più estremi del modello sembra leggermente peggiorata (ora sono visibili le curve corrispondenti alla soglia di allarme 0.5), ma nessun punto è oltre la soglia 0.5, pertanto rimane invariata la conclusione che non ci sono punti pericolosamente influenti sul modello.
7 PREDIZIONE COL MODELLO SCELTO (VALORE ATTESO)
Possiamo usare il modello selezionato per predire il peso di un neonato di cui si conoscono alcune variabili. Le variabili del modello che non si conoscono verranno poste uguali alla media di quelle variabili.
predict(stepwise.mod4, newdata = data.frame(grav_cat=2, Sesso="F", Gestazione_cen=39-mean(db3$Gestazione), Lunghezza_cen=0, Cranio_cen=0))## 1
## 3232.714Il peso di una neonata alla 39esima settimana, terza figlia, sarebbe 3233 grammi. Tuttavia, assegnando 0 a lunghezza_cen e cranio_cen, stiamo assumendo delle dimensioni attese pari alla media generale dei neonati, mediando da maschi e femmine e nati a una qualunque settimana. Una stima più accurata potrebbe considerare come valore atteso per una femmina la media relativa alle sole neonate femmine di 39 settimane (escludiamo il numero di gravidanze precedenti, variabile non presente nel modello).
LunghezzaF39_cen=mean(db3$Lunghezza_cen[db3$Sesso=="F" & db3$Gestazione==39])
CranioF39_cen=mean(db3$Cranio_cen[db3$Sesso=="F" & db3$Gestazione==39])
predict(stepwise.mod4, newdata = data.frame(grav_cat=2, Sesso="F", Gestazione_cen=39-mean(db3$Gestazione), Lunghezza_cen=LunghezzaF39_cen, Cranio_cen=CranioF39_cen))## 1
## 3198.539Il peso atteso di una neonata a 39 settimane, in assenza di altre informazioni, potrebbe essere 3199 grammi.
8 CONCLUSIONE
Il peso atteso di un neonato dipende significativamente dalle dimensioni (Lunghezza, Lunghezza^2 e circonferenza craniale), dalla settimana di gestazione e dal sesso.
Altri fattori come il numero di gravidanze precedenti o l’età della madre non hanno effetti statisticamente significativi o, se ce li hanno, non migliorano la frazione di variabilità spiegata dal modello a sufficienza per giustificare la propria presenza nel modello.
Il fattore fumo mostra un’interazione sia con la lunghezza (i neonati più corti hanno sottopeso accentuato), le settimane di gestazione (indebolendo o annullando l’effetto di guadagno del peso con le settimane), o il numero di gravidanze (annullando l’effetto di guadagno del peso col numero di gravidanze), ma l’effetto è troppo debole per questo campione con solo il 5% di madri fumatrici per avere dei risultati significativi.
Il sospetto che il fumo danneggi lo sviluppo del feto resta forte ma si consigliano ulteriori studi per confermarlo, magari con un campione più bilanciato.