Newborn Prediction Palmisano
Palmisano Giuseppe
2024-10-26
Importo i pacchetti che saranno utili per questa analisi
library(moments)
library(tibble)
library(dplyr)
library(ggplot2)
library(ggpubr)
library(GGally)
library(car)
library(lmtest)
library(MASS)
library(scatterplot3d)Ora siamo pronti a importare il nostro dataset per darci un’occhiata
newborn.df <- read.csv("neonati.csv", stringsAsFactors = T)
newborn.df$Fumatrici <- as.factor(newborn.df$Fumatrici)
newborn.df$Ospedale <- as.factor(newborn.df$Ospedale)
head(newborn.df,10)Fumatrici e Ospedale sono
interpretati come variabili numeriche ma, poiché le vogliamo come
fattori, faremo la conversione usando as.factor.
Le variabili nel dataframe sono:
Anni.madre: età della madre (variabile quantitativa continua);
N.gravidanze: numero di gravidanze che la madre ha già avuto (variabile quantitativa discreta);
Fumatrici: è 0 se la madre non fuma, altrimenti è 1 (variabile qualitativa nominale);
Gestazione: numero di settimane di gestazione (variabile quantitativa continua);
Peso: peso del neonato, in g (variabile quantitativa continua);
Lunghezza: lunghezza del neonato, in mm (variabile quantitativa continua);
Cranio: diametro del cranio del neonato, in mm (variabile quantitativa continua);
Tipo.parto: tipo di parto, Naturale o Cesareo (variabile qualitativa nominale);
Ospedale: ospedale, 1, 2 o 3 (variabile qualitativa nominale);
Sesso: sesso del neonato, Maschio o Femmina (variabile qualitativa nominale).
Il progetto mira a prevedere il peso del neonato, date tutte le altre variabili. Studieremo come le variabili influenzano il peso e vedremo quali di esse svolgono un ruolo rilevante nella sua determinazione. Per raggiungere questo obiettivo, utilizziamo la regressione lineare multipla.
Facciamo un’analisi esplorativa delle variabili.
summary(newborn.df)## Anni.madre N.gravidanze Fumatrici Gestazione Peso
## Min. : 0.00 Min. : 0.0000 0:2396 Min. :25.00 Min. : 830
## 1st Qu.:25.00 1st Qu.: 0.0000 1: 104 1st Qu.:38.00 1st Qu.:2990
## Median :28.00 Median : 1.0000 Median :39.00 Median :3300
## Mean :28.16 Mean : 0.9812 Mean :38.98 Mean :3284
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:40.00 3rd Qu.:3620
## Max. :46.00 Max. :12.0000 Max. :43.00 Max. :4930
## Lunghezza Cranio Tipo.parto Ospedale Sesso
## Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1256
## 1st Qu.:480.0 1st Qu.:330 Nat:1772 osp2:849 M:1244
## Median :500.0 Median :340 osp3:835
## Mean :494.7 Mean :340
## 3rd Qu.:510.0 3rd Qu.:350
## Max. :565.0 Max. :390Anni.madre ha un minimo di 0.0, il che chiaramente non è possibile. Dobbiamo indagare più a fondo.
filter(newborn.df, Anni.madre <13)Con filter vediamo che ci sono due casi con un valore
per Anni.madre che sono biologicamente impossibili.
Poiché le altre variabili per questi due record hanno valori che
sembrano plausibili e abbiamo 2500 osservazioni per ciascuna variabile,
non rimuoviamo questi casi, ma facciamo imputation.
Questo significa che li sostituiamo con la media o la mediana dei valori
rimanenti.
Usiamo il test di normalità di Shapiro-Wilk per determinare se la variabile ha una distribuzione normale (ipotesi nulla) o no. Nel primo caso possiamo usare la media per l’imputation, nel secondo la mediana.
age_verified <- subset(newborn.df,Anni.madre>2)$Anni.madre
shapiro.test(age_verified)##
## Shapiro-Wilk normality test
##
## data: age_verified
## W = 0.99491, p-value = 1.477e-07I risultati ci portano a rifiutare l’ipotesi nulla, quindi usiamo la mediana.
median_age <- median(age_verified)
newborn.df$Anni.madre <- replace(newborn.df$Anni.madre,newborn.df$Anni.madre<2,median_age)
summary(newborn.df)## Anni.madre N.gravidanze Fumatrici Gestazione Peso
## Min. :13.00 Min. : 0.0000 0:2396 Min. :25.00 Min. : 830
## 1st Qu.:25.00 1st Qu.: 0.0000 1: 104 1st Qu.:38.00 1st Qu.:2990
## Median :28.00 Median : 1.0000 Median :39.00 Median :3300
## Mean :28.19 Mean : 0.9812 Mean :38.98 Mean :3284
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:40.00 3rd Qu.:3620
## Max. :46.00 Max. :12.0000 Max. :43.00 Max. :4930
## Lunghezza Cranio Tipo.parto Ospedale Sesso
## Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1256
## 1st Qu.:480.0 1st Qu.:330 Nat:1772 osp2:849 M:1244
## Median :500.0 Median :340 osp3:835
## Mean :494.7 Mean :340
## 3rd Qu.:510.0 3rd Qu.:350
## Max. :565.0 Max. :390Passiamo a confrontare i valori di Peso e Lunghezza con quelli della popolazione. I dati della popolazione sono stati presi da https://www.cdc.gov/growthcharts/who_charts.htm. In particolare:
- Peso Maschio: 3.521 kg
- Peso Femmina: 3.362 kg
- Lunghezza Maschio: 51,4 cm
- Lunghezza Femmina: 50,3 cm
Usiamo il Test T di Student per questo confronto.
t_Peso_F <- t.test(filter(newborn.df, Sesso=="F")["Peso"],
mu = 3362,
conf.level = 0.95,
alternative = "two.sided")
t_Peso_M <- t.test(filter(newborn.df, Sesso=="M")["Peso"],
mu = 3521,
conf.level = 0.95,
alternative = "two.sided")
t_Lunghezza_F <- t.test(filter(newborn.df, Sesso=="F")["Lunghezza"],
mu = 503,
conf.level = 0.95,
alternative = "two.sided")
t_Lunghezza_M <- t.test(filter(newborn.df, Sesso=="M")["Lunghezza"],
mu = 514,
conf.level = 0.95,
alternative = "two.sided")
t_Peso_F$p.value## [1] 5.163174e-39t_Peso_M$p.value## [1] 1.837195e-15t_Lunghezza_F$p.value## [1] 1.048128e-58t_Lunghezza_M$p.value## [1] 3.045549e-84Da questi risultati dobbiamo rifiutare l’ipotesi nulla e concludere che questo dataset non appartiene alla popolazione. È importante notare che non sappiamo quando questi dati siano stati raccolti e, citando il sito, grandi cambiamenti sono avvenuti negli ultimi 10 anni.
Possiamo usare un Test T a due campioni per capire se Sesso è statisticamente significativo in relazione a Peso e Lunghezza e visualizzarlo attraverso un boxplot.
t_test_Peso <-t.test(data = newborn.df,
Peso ~ Sesso)
t_test_Lunghezza <-t.test(data = newborn.df,
Lunghezza ~ Sesso)
t_test_Peso$p.value## [1] 8.029743e-33t_test_Lunghezza$p.value## [1] 2.237243e-21Sesso_colors <- c("pink", "lightblue")
box_Peso <- ggplot(newborn.df, aes(x=Sesso,y=Peso))+
geom_boxplot(aes(color = Sesso)
)+
scale_color_manual(values = Sesso_colors)+
scale_y_continuous(breaks = seq(500,5500,500))+
labs(x="Sesso",
y="Peso (g)",
title = "Nuova nascite Peso per Sesso")+
theme_minimal()
box_Lunghezza <- ggplot(newborn.df, aes(x=Sesso,y=Lunghezza))+
geom_boxplot(aes(color = Sesso))+
scale_color_manual(values = Sesso_colors)+
scale_y_continuous(breaks = seq(300,600,50))+
labs(x="Sesso",
y="Lunghezza (mm)",
title = "Newborn's Lunghezza per Sesso")+
theme_minimal()
ggarrange(box_Peso,box_Lunghezza,nrow = 1)
In entrambi i casi, il p-value è molto piccolo, quindi rifiutiamo l’ipotesi nulla, concludendo che la differenza tra i due valori medi è statisticamente significativa.
L’ultima analisi che andremo a eseguire riguarda la correlazione tra il tipo di parto e l’Ospedale. Per farla, creiamo una tabella di contingenza e poi eseguiamo un Test Chi Quadrato.
Ospedale_birth_type <- table(newborn.df$Tipo.parto,newborn.df$Ospedale)
Ospedale_birth_type##
## osp1 osp2 osp3
## Ces 242 254 232
## Nat 574 595 603chisq.test(Ospedale_birth_type)##
## Pearson's Chi-squared test
##
## data: Ospedale_birth_type
## X-squared = 1.0972, df = 2, p-value = 0.5778Visualizziamo i risultati.
ggballoonplot(data = as.data.frame(Ospedale_birth_type),
fill = "value")+
labs(x="Birth type",
y="Sesso",
title = "Genere vs tipo di parto",
fill="Frequenza")+
guides(size=F)+
theme(plot.title = element_text(hjust = 0.5))
Ora indaghiamo la relazione tra ciascuna variabile e le altre. Per semplificare questo processo, dividiamo il dataframe in due sotto-dataframe: uno per le variabili continue e uno per quelle categoriche.
num_newborn.df <- newborn.df[, sapply(newborn.df, is.numeric)]
fact_newborn.df <- newborn.df[, sapply(newborn.df, is.factor)]Per capire la relazione tra le variabili continue utilizziamo il Test di Correlazione e lo visualizziamo graficamente.
ggcorr(num_newborn.df, label = TRUE, size= 2.5)
Troviamo che la variabile Peso è altamente correlata con le variabili Gestazione, Lunghezza e Cranio. Anni.madre e N.gravidanze non hanno alcuna correlazione con Peso.
Per le variabili categoriche utilizziamo il Test Chi Quadrato. Creiamo una funzione per iterare attraverso il dataframe che abbiamo creato e ottenere come output il p-value di ciascun test.
chisq_calc <- function(df) {
results <- data.frame()
variables <- names(df)
for (i in 1:(length(variables) - 1)) {
for (j in (i + 1):length(variables)) {
var1 <- df[[variables[i]]]
var2 <- df[[variables[j]]]
contingency_table <- table(var1, var2)
chi_sq_test <- chisq.test(contingency_table)
p_value <- round(chi_sq_test$p.value, 3)
results <- rbind(results, c(variables[i], variables[j], p_value))
}
}
colnames(results) <- c("Variable1", "Variable2", "P_Value")
return(results)
}
chisq_calc(fact_newborn.df)Nessun p-value è inferiore a 0,05, quindi non possiamo rifiutare l’ipotesi nulla e possiamo affermare che non c’è dipendenza tra le variabili.
L’ultimo passo è indagare la relazione tra variabili continue e variabili categoriche. Questa volta utilizziamo il Test T di Student (o il Test T a coppie quando una variabile categorica ha 3 livelli, come Ospedale).
t_test_calc <- function(df1, df2) {
variable1 <- names(df1)
variable2 <- names(df2)
res <- data.frame()
for (i in variable1) {
for (j in variable2){
var1 <- df1[[i]]
var2 <- df2[[j]]
if (nlevels(var2) == 2) {
test <- t.test(var1 ~ var2)
pval <- round(test$p.value, 3)
res <- bind_rows(res, data.frame(Variable1 = i,
Variable2 = j,
P_Value = pval))
} else {
test <- pairwise.t.test(var1, var2,
pool.sd = TRUE,
p.adjust.method = "holm")
pval12 <- round(test$p.value[1,1], 3)
pval13 <- round(test$p.value[2,1], 3)
pval23 <- round(test$p.value[2,2], 3)
lab12 <- paste(j, "1-2")
lab13 <- paste(j, "1-3")
lab23 <- paste(j, "2-3")
res <- bind_rows(res, data.frame(Variable1 = i,
Variable2 = lab12,
P_Value = pval12))
res <- bind_rows(res, data.frame(Variable1 = i,
Variable2 = lab13,
P_Value = pval13))
res <- bind_rows(res, data.frame(Variable1 = i,
Variable2 = lab23,
P_Value = pval23))
}
}
}
return(res)
}
# Esecuzione della funzione
t_test_calc(num_newborn.df, fact_newborn.df)Si scopre che Peso (la variabile di risposta) è fortemente influenzata da Sesso (lo sapevamo già), mentre non dipende in modo significativo dalle variabili Fumatrici, Tipo.parto e Ospedale.
È ora di concentrarsi sull’obiettivo principale di questo studio. Vogliamo fare previsioni sul Peso dei neonati, quindi dobbiamo creare un modello di regressione lineare multipla per raggiungere questo scopo. Partiamo dal modello che contiene tutte le variabili.
mod1 <- lm(Peso ~ .,data= newborn.df)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = newborn.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.3 -181.2 -14.6 160.7 2612.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6735.1677 141.3977 -47.633 < 2e-16 ***
## Anni.madre 0.7983 1.1463 0.696 0.4862
## N.gravidanze 11.4118 4.6665 2.445 0.0145 *
## Fumatrici1 -30.1567 27.5396 -1.095 0.2736
## Gestazione 32.5265 3.8179 8.520 < 2e-16 ***
## Lunghezza 10.2951 0.3007 34.237 < 2e-16 ***
## Cranio 10.4725 0.4261 24.580 < 2e-16 ***
## Tipo.partoNat 29.5027 12.0848 2.441 0.0147 *
## Ospedaleosp2 -11.2216 13.4388 -0.835 0.4038
## Ospedaleosp3 28.0984 13.4972 2.082 0.0375 *
## SessoM 77.5473 11.1779 6.938 5.07e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 669.1 on 10 and 2489 DF, p-value: < 2.2e-16Da questo modello possiamo vedere che il valore di R^2 è abbastanza buono, ma alcune variabili hanno coefficienti con un valore elevato (>0,05), il che significa che la loro significatività è bassa.
Per trovare il miglior modello rimuoviamo, volta per volta, le variabili con il valore più alto dei coefficienti e aggiorniamo il modello. La prima variabile da eliminare è Anni.madre.
mod2 <- update(mod1,~.-Anni.madre)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso, data = newborn.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1113.93 -180.11 -16.36 160.58 2616.96
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6708.1065 135.9394 -49.346 < 2e-16 ***
## N.gravidanze 12.6085 4.3381 2.906 0.00369 **
## Fumatrici1 -30.3092 27.5359 -1.101 0.27113
## Gestazione 32.2501 3.7968 8.494 < 2e-16 ***
## Lunghezza 10.2944 0.3007 34.239 < 2e-16 ***
## Cranio 10.4876 0.4255 24.651 < 2e-16 ***
## Tipo.partoNat 29.5351 12.0834 2.444 0.01458 *
## Ospedaleosp2 -11.0816 13.4359 -0.825 0.40957
## Ospedaleosp3 28.3660 13.4903 2.103 0.03559 *
## SessoM 77.6205 11.1763 6.945 4.81e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2490 degrees of freedom
## Multiple R-squared: 0.7288, Adjusted R-squared: 0.7278
## F-statistic: 743.6 on 9 and 2490 DF, p-value: < 2.2e-16Ora è Ospedale ad avere il valore più alto. Rimuoviamolo.
mod3 <- update(mod2,~.-Ospedale)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Sesso, data = newborn.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1130.04 -181.29 -16.31 160.59 2635.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6708.074 135.984 -49.330 < 2e-16 ***
## N.gravidanze 13.012 4.342 2.997 0.00276 **
## Fumatrici1 -31.759 27.570 -1.152 0.24946
## Gestazione 32.541 3.801 8.561 < 2e-16 ***
## Lunghezza 10.272 0.301 34.129 < 2e-16 ***
## Cranio 10.501 0.426 24.648 < 2e-16 ***
## Tipo.partoNat 30.296 12.098 2.504 0.01234 *
## SessoM 78.114 11.191 6.980 3.77e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2492 degrees of freedom
## Multiple R-squared: 0.7278, Adjusted R-squared: 0.7271
## F-statistic: 952 on 7 and 2492 DF, p-value: < 2.2e-16Lo stesso vale per Fumatrici. Rimuoviamola.
mod4 <- update(mod3,~.-Fumatrici)
summary(mod4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso, data = newborn.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1129.31 -181.70 -16.31 161.07 2638.85
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.2971 135.9911 -49.322 < 2e-16 ***
## N.gravidanze 12.7558 4.3366 2.941 0.0033 **
## Gestazione 32.2713 3.7941 8.506 < 2e-16 ***
## Lunghezza 10.2864 0.3007 34.207 < 2e-16 ***
## Cranio 10.5057 0.4260 24.659 < 2e-16 ***
## Tipo.partoNat 30.0342 12.0969 2.483 0.0131 *
## SessoM 77.9285 11.1905 6.964 4.22e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2493 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.727
## F-statistic: 1110 on 6 and 2493 DF, p-value: < 2.2e-16Ora è il turno di Tipo.parto. Rimuoviamola.
mod5 <- update(mod4,~.-Tipo.parto)
summary(mod5)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = newborn.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.44 -180.81 -15.58 163.64 2639.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.1445 135.7229 -49.226 < 2e-16 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## Gestazione 32.3321 3.7980 8.513 < 2e-16 ***
## Lunghezza 10.2486 0.3006 34.090 < 2e-16 ***
## Cranio 10.5402 0.4262 24.728 < 2e-16 ***
## SessoM 77.9927 11.2021 6.962 4.26e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1328 on 5 and 2494 DF, p-value: < 2.2e-16Il valore di R^2 non è cambiato molto durante questi aggiornamenti,
ma ora tutte le variabili incluse risultano essere rilevanti.
mod5 sembra essere un buon candidato. Vediamo se
otteniamo lo stesso risultato utilizzando le funzioni BIC e
AIC.
BIC(mod1,mod2,mod3,mod4,mod5)AIC(mod1,mod2,mod3,mod4,mod5)Per entrambi (AIC e BIC), il modello migliore è quello con il valore più basso associato. Quindi è confermato che mod5 è il miglior modello finora.
Diamo un’occhiata ai residui di questo modello.
par(mfrow=c(2,2))
plot(mod5)
Dai grafici sembra esserci un effetto non lineare tra le variabili. Aggiorniamo il modello aggiungendo i termini quadratici delle variabili incluse in mod5.
mod6 <- update(mod5,~.+I(Gestazione^2)+I(Lunghezza^2)+I(Cranio^2))
summary(mod6)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2), data = newborn.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1187.03 -182.50 -12.03 162.60 1469.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.217e+03 1.160e+03 -1.049 0.2945
## N.gravidanze 1.441e+01 4.246e+00 3.394 0.0007 ***
## Gestazione 3.754e+02 6.738e+01 5.572 2.79e-08 ***
## Lunghezza -2.925e+01 4.436e+00 -6.592 5.28e-11 ***
## Cranio -4.979e+00 9.804e+00 -0.508 0.6116
## SessoM 7.232e+01 1.100e+01 6.577 5.83e-11 ***
## I(Gestazione^2) -4.374e+00 8.834e-01 -4.951 7.86e-07 ***
## I(Lunghezza^2) 4.075e-02 4.544e-03 8.967 < 2e-16 ***
## I(Cranio^2) 2.276e-02 1.445e-02 1.575 0.1154
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.4 on 2491 degrees of freedom
## Multiple R-squared: 0.7395, Adjusted R-squared: 0.7387
## F-statistic: 883.9 on 8 and 2491 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(mod6)
I residui (grafico in alto a sinistra) mostrano una linea più orizzontale rispetto a quella del modello precedente. Tuttavia, Cranio ha perso significatività, e anche il suo termine quadratico ha una bassa significatività. Possiamo provare a rimuovere il termine quadratico di Cranio e vedere cosa succede al modello.
mod7 <- update(mod6,~.-I(Cranio^2))
summary(mod7)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2) + I(Lunghezza^2), data = newborn.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1191.07 -182.28 -13.74 163.34 1403.53
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.360e+03 9.053e+02 -2.607 0.00919 **
## N.gravidanze 1.448e+01 4.247e+00 3.410 0.00066 ***
## Gestazione 3.365e+02 6.270e+01 5.367 8.74e-08 ***
## Lunghezza -3.215e+01 4.037e+00 -7.963 2.53e-15 ***
## Cranio 1.045e+01 4.192e-01 24.922 < 2e-16 ***
## SessoM 7.261e+01 1.100e+01 6.602 4.93e-11 ***
## I(Gestazione^2) -3.873e+00 8.245e-01 -4.698 2.77e-06 ***
## I(Lunghezza^2) 4.370e-02 4.140e-03 10.556 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.5 on 2492 degrees of freedom
## Multiple R-squared: 0.7392, Adjusted R-squared: 0.7385
## F-statistic: 1009 on 7 and 2492 DF, p-value: < 2.2e-16par(mfrow=c(2,2))
plot(mod7)
Con questo aggiornamento Cranio ha recuperato un
buon livello di significatività, mentre i grafici sono cambiati solo
leggermente. A questo punto, mod7 è il principale
candidato. Usiamo nuovamente le funzioni BIC e
AIC per verificare se possiamo confermarlo.
BIC(mod5,mod6,mod7)AIC(mod5,mod6,mod7)BIC ci conferma la scelta del modello, mentre
AIC preferisce leggermente mod6. Questo
accade perché il BIC applica una penalizzazione maggiore ai
modelli con molte variabili. Poiché la differenza tra
mod6 e mod7 è trascurabile, se
confrontata, ad esempio, con la differenza tra mod5 e
mod6, scegliamo mod7 come miglior
modello, in quanto più semplice (meno parametri).
Calcoliamo i variance inflation factors (VIF) di mod7 per valutare se è presente multicollinearità.
vif(mod7, type="predictor")I risultati ci mostrano che nessun parametro supera 5, quindi possiamo concludere che non c’è multicollinearità.
Possiamo verificare la media dei residui di mod7 per vedere se è 0, il che indicherebbe una distribuzione normale.
mean(residuals(mod7))## [1] 1.461354e-14sd(residuals(mod7))## [1] 268.1125Il risultato indica che è circa 0, quindi possiamo confermare questo punto.
Verifichiamo ora l’ipotesi di omoschedasticità attraverso il Test di Breusch-Pagan.
bptest(mod7)##
## studentized Breusch-Pagan test
##
## data: mod7
## BP = 97.553, df = 7, p-value < 2.2e-16Dobbiamo rifiutare l’ipotesi nulla e concludere che i residui sono eteroschedastici, indicando che non hanno una varianza costante.
Ora eseguiamo il Test di Durbin-Watson per verificare l’ipotesi di indipendenza dei residui:
dwtest(mod7)##
## Durbin-Watson test
##
## data: mod7
## DW = 1.9494, p-value = 0.1028
## alternative hypothesis: true autocorrelation is greater than 0Questi risultati ci portano a non rifiutare l’ipotesi nulla, quindi possiamo affermare che gli errori sono non correlati.
Infine, diamo un’occhiata agli outlier, cioè quei punti in cui la variabile di risposta (Peso in questo caso) è distante dal valore predetto dal modello, e ai leverages, ovvero osservazioni con valori insoliti nelle variabili. Punti con outlier o alto leverage possono distorcere il risultato e l’accuratezza di un’analisi di regressione. La Distanza di Cook è la statistica che possiamo utilizzare per determinare l’influenza di questi punti.
cook <- cooks.distance(mod7)
ggplot() +
geom_point(aes(x = 1:length(cook),
y = cook,
colour = cook > 1),
size = 3) +
geom_hline(aes(yintercept = c(0.5, 1)),
linetype = 2,
colour = "darkred") +
scale_colour_manual(values = setNames(c("darkred", "black"), c(TRUE, FALSE))) +
labs(title = "Analysis of residuals: Cook's distance",
x = "Index",
y = "Cook's distance") +
theme_minimal() +
theme(plot.title = element_text(size = 22, hjust = 0.5),
axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14),
axis.title = element_text(size = 16),
legend.position = "none")
Sembra che l’unico valore con una Distanza di Cook superiore al valore critico di 1 sia l’osservazione 1551. Possiamo provare a rimuoverla dal nostro dataframe e ricreare il modello.
index <- match(max(cook),cook)
corrected_nb.df <- newborn.df[-index,]
mod8 <- lm(Peso ~ N.gravidanze + Gestazione +
Lunghezza + Cranio + Sesso + I(Gestazione^2) + I(Lunghezza^2),
data = corrected_nb.df)
summary(mod8)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2) + I(Lunghezza^2), data = corrected_nb.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1187.72 -180.63 -13.61 164.12 1320.11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.811e+03 9.019e+02 -3.117 0.001848 **
## N.gravidanze 1.440e+01 4.217e+00 3.414 0.000651 ***
## Gestazione 2.004e+02 6.617e+01 3.028 0.002489 **
## Lunghezza -1.928e+01 4.535e+00 -4.252 2.20e-05 ***
## Cranio 1.009e+01 4.202e-01 24.023 < 2e-16 ***
## SessoM 7.340e+01 1.092e+01 6.721 2.22e-11 ***
## I(Gestazione^2) -2.135e+00 8.673e-01 -2.462 0.013890 *
## I(Lunghezza^2) 3.093e-02 4.618e-03 6.697 2.62e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.6 on 2491 degrees of freedom
## Multiple R-squared: 0.7426, Adjusted R-squared: 0.7419
## F-statistic: 1027 on 7 and 2491 DF, p-value: < 2.2e-16Dopo questi passaggi, R^2 è aumentato, ma lo stesso è accaduto per il p-value del termine quadratico di Gestazione. Anche se è ancora significativo, possiamo provare a rimuoverlo e osservare i risultati del modello.
mod9 <- update(mod8,~.-I(Gestazione^2))
summary(mod9)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2), data = corrected_nb.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1176.79 -178.89 -12.52 164.19 1327.44
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.607e+03 7.586e+02 -2.119 0.034207 *
## N.gravidanze 1.420e+01 4.220e+00 3.364 0.000781 ***
## Gestazione 3.773e+01 3.890e+00 9.700 < 2e-16 ***
## Lunghezza -1.172e+01 3.342e+00 -3.508 0.000459 ***
## Cranio 1.015e+01 4.201e-01 24.157 < 2e-16 ***
## SessoM 7.222e+01 1.092e+01 6.614 4.57e-11 ***
## I(Lunghezza^2) 2.330e-02 3.430e-03 6.795 1.35e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.8 on 2492 degrees of freedom
## Multiple R-squared: 0.742, Adjusted R-squared: 0.7413
## F-statistic: 1194 on 6 and 2492 DF, p-value: < 2.2e-16Ora tutti i parametri hanno lo stesso alto livello di significatività, sebbene R^2 sia leggermente diminuito. Possiamo fare un confronto tra gli ultimi tre modelli.
BIC(mod7,mod8,mod9)AIC(mod7,mod8,mod9)Come accaduto in precedenza, il BIC evidenzia mod9 come il miglior modello, poiché ha meno parametri, mentre l’AIC preferisce mod8. Come prima, scegliamo mod9 poiché è più semplice.
Ora possiamo indagare i residui come fatto per mod7.
par(mfrow=c(2,2))
plot(mod9)
dwtest(mod9)##
## Durbin-Watson test
##
## data: mod9
## DW = 1.9498, p-value = 0.1046
## alternative hypothesis: true autocorrelation is greater than 0bptest(mod9)##
## studentized Breusch-Pagan test
##
## data: mod9
## BP = 15.47, df = 6, p-value = 0.0169Questa volta possiamo assumere l’omoschedasticità dei residui. Infatti, il Test di Breusch-Pagan ha un p-value più alto e possiamo impostare il livello di significatività a 0,01 (invece di 0,05). L’assenza di correlazione è garantita come prima. Possiamo concludere che mod9 è il miglior modello che possiamo ottenere e può essere affidabile per prevedere il Peso dei neonati.
È il momento di testare il nostro modello. Facciamo una previsione per il Peso di un neonato. Ad esempio, consideriamo una madre che:
ha già avuto 3 gravidanze;
partorirà il bambino durante la 39ª settimana.
Supponiamo anche di non avere informazioni sulla Lunghezza e sul diametro del Cranio. In questo caso, possiamo utilizzare i valori medi relativi alle femmine per fornire una stima di questi parametri.
prediction = predict(mod9,
newdata = data.frame(Sesso="F",Gestazione=39,N.gravidanze=3,
Lunghezza=mean(newborn.df$Lunghezza[newborn.df$Sesso=="F"]),
Cranio=mean(newborn.df$Cranio[newborn.df$Sesso=="F"])),
interval = "predict")
prediction## fit lwr upr
## 1 3180.961 2657.185 3704.736I due valori nell’output rappresentano i valori estremi del corrispondente intervallo di previsione al 95%.
prediction[3] - prediction[1]## [1] 523.7752prediction[1] - prediction[2]## [1] 523.7752Possiamo affermare che questo bambino avrà un Peso di 3180,96 \(\pm\) 523,77 g.
Infine, dopo aver effettuato la nostra previsione, possiamo provare a visualizzare una versione semplificata dei nostri dati in un grafico a dispersione 3D.
par(mfrow = c(1, 1))
colors <- c("lightblue", "coral")
colors <- colors[as.numeric(newborn.df$Sesso)]
s3d <- scatterplot3d(newborn.df$N.gravidanze, newborn.df$Lunghezza, newborn.df$Peso,
color = colors,
pch = 16,
angle = 50,
main = "3D Scatter plot",
xlab = "N° gravidanze",
ylab = "Lunghezza (mm)",
zlab = "Peso (g)",
grid = TRUE,
box = FALSE)
legend("right", legend = levels(newborn.df$Sesso),
col = c("lightblue", "coral"),
pch = 16,
xpd = TRUE,
xjust = 0,
cex = 0.8,
box.lty = 0,
bg = "transparent")
F_data <- subset(newborn.df, Ospedale == "osp1")
M_data <- subset(newborn.df, Ospedale == "osp2")
F_lm <- lm(Peso ~ N.gravidanze + Lunghezza, data = F_data)
M_lm <- lm(Peso ~ N.gravidanze + Lunghezza, data = M_data)
print(summary(F_lm))##
## Call:
## lm(formula = Peso ~ N.gravidanze + Lunghezza, data = F_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1007.94 -219.85 -21.27 206.93 1683.05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4531.2354 206.0524 -21.991 <2e-16 ***
## N.gravidanze 18.3707 8.7092 2.109 0.0352 *
## Lunghezza 15.7517 0.4145 38.005 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 320.8 on 813 degrees of freedom
## Multiple R-squared: 0.6401, Adjusted R-squared: 0.6392
## F-statistic: 723 on 2 and 813 DF, p-value: < 2.2e-16print(summary(M_lm))##
## Call:
## lm(formula = Peso ~ N.gravidanze + Lunghezza, data = M_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1459.1 -191.6 -18.2 179.2 1399.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4969.4271 204.9905 -24.242 <2e-16 ***
## N.gravidanze 14.4663 8.1340 1.779 0.0757 .
## Lunghezza 16.6081 0.4117 40.337 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 306.4 on 846 degrees of freedom
## Multiple R-squared: 0.6583, Adjusted R-squared: 0.6575
## F-statistic: 814.9 on 2 and 846 DF, p-value: < 2.2e-16s3d$plane3d(F_lm, plane_color = "coral1", lwd = 1.5)
s3d$plane3d(M_lm, plane_color = "lightblue", lwd = 1.5)
Abbiamo scelto Lunghezza e N.gravidanze rispetto a Peso e le abbiamo correlate a Sesso. Dal grafico 3D possiamo dedurre che, mantenendo fisso il numero di gravidanze, il Peso aumenterà di più per una femmina rispetto a un maschio.