Modulo_3_Project
Jacopo Carretta
2025-11-30
INTRODUZIONE
Lo scopo di questo progetto è lo svuluppo di un modello statistico che possa prevedere con precisione il peso dei neonati alla nascita. La base dati è stata raccolta su 2500 neonati di 3 diversi ospedali.
library(ggplot2)
library(moments)
library(dplyr)##
## Caricamento pacchetto: 'dplyr'## I seguenti oggetti sono mascherati da 'package:stats':
##
## filter, lag## I seguenti oggetti sono mascherati da 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(zoo)##
## Caricamento pacchetto: 'zoo'## I seguenti oggetti sono mascherati da 'package:base':
##
## as.Date, as.Date.numericlibrary(MASS)## Warning: il pacchetto 'MASS' è stato creato con R versione 4.5.2##
## Caricamento pacchetto: 'MASS'## Il seguente oggetto è mascherato da 'package:dplyr':
##
## selectlibrary(car)## Warning: il pacchetto 'car' è stato creato con R versione 4.5.2## Caricamento del pacchetto richiesto: carData## Warning: il pacchetto 'carData' è stato creato con R versione 4.5.2##
## Caricamento pacchetto: 'car'## Il seguente oggetto è mascherato da 'package:dplyr':
##
## recodesetwd("C:/Users/Utente/OneDrive/Desktop")
# Carico il dataset in R
dati <- read.csv("neonati.csv")ANALISI PRELIMINARE
Di seguito verrà svolta una veloce analisi descrittiva delle variabili prese in considerazione.
Anni madre: quantitativa discreta su scala di rapporti. N. Gravidanze: quantitativa discreta su scala di rapporti. Fumatrici: qualitativa binaria dummy. Gestazione: quantitativa discreta su scala di rapporti. Peso: quantitativa continua su scala di rapporti. Lunghezza: quantitativa continua su scala di rapporti. Cranio: quantitativa continua su scala di rapporti. Parto: qualitativa binaria. Ospedale: qualitativa su scala nominale. Sesso: qualitativa binaria.
summary(dati)## 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
##
##
## Dai risultati del summary probabilmente ci saranno degli outlier nella variabile “Anni.madre” (è impossibile abbia 0 anni una madre al momento del parto) e anche nella variabile “N.gravidanze” (probabilmente il dato 12 gravidanze sarà un outlier).
Al fine di approfondire la distribuzione e la presenza di outlier per le variabili “Anni.madre”, “N.gravidanze”, “Gestazione”, “Peso”, “Lunghezza” e “Cranio” saranno utilizzati dei boxplots:
boxplot_var <- c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")
for (i in boxplot_var) {
graphs <- ggplot(dati, aes_string(y = i)) +
geom_boxplot() +
labs(x = i, title = paste("Boxplot di", i))
print(graphs)
}## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Nei boxplot precedenti è riassunta la distribuzione delle variabili e
sottolineano la presenza di molti outlier.
Prima di procedere, riassumiamo anche la distribuzione delle variabili qualitative utilizzando anche dei grafici a barre per dare un indicazione visiva:
barplot_var <- c("Fumatrici", "Tipo.parto", "Ospedale", "Sesso")
for (i in barplot_var) {
cat("\nTabella di", i, "\n")
print(table(dati[[i]]))
cat("Percentuali:\n")
print(round(prop.table(table(dati[[i]])) * 100, 2))
graphs<- ggplot(dati, aes(x = .data[[i]])) +
geom_bar(fill = "steelblue") +
labs(
x = i,
y = "Frequenza",
title = paste("Distribuzione di", i)
) +
theme_minimal()
print(graphs)
}##
## Tabella di Fumatrici
##
## 0 1
## 2396 104
## Percentuali:
##
## 0 1
## 95.84 4.16
##
## Tabella di Tipo.parto
##
## Ces Nat
## 728 1772
## Percentuali:
##
## Ces Nat
## 29.12 70.88
##
## Tabella di Ospedale
##
## osp1 osp2 osp3
## 816 849 835
## Percentuali:
##
## osp1 osp2 osp3
## 32.64 33.96 33.40
##
## Tabella di Sesso
##
## F M
## 1256 1244
## Percentuali:
##
## F M
## 50.24 49.76
Per saggiare l’ipotesi che in alcuni ospedali si fanno più parti cesarei
verrà utilizzato il test chi-quadrato:
chisq.test(dati$Ospedale,dati$Tipo.parto)##
## Pearson's Chi-squared test
##
## data: dati$Ospedale and dati$Tipo.parto
## X-squared = 1.0972, df = 2, p-value = 0.5778Dal test chi-quadrato si ottiene un p-value maggiore di 0.05 dunque NON si rifiuta l’ipotesi nulla e perciò non c’è evidenza statistica che in alcuni ospedali ci siano più parti cesarei. Al fine di visualizzare la situazione graficamente in seguito un barplot con la % di parti cesarei rispetto all’ospedale.
perc <- prop.table(table(dati$Ospedale, dati$Tipo.parto),
margin = 1)[, "Ces"] * 100
df_plot <- data.frame(
Ospedale = names(perc),
Perc_cesarei = as.numeric(perc)
)
ggplot(df_plot, aes(x = Ospedale, y = Perc_cesarei)) +
geom_col(fill = "steelblue") +
geom_text(aes(label = paste0(round(Perc_cesarei, 1), "%")),
vjust = -0.3) +
ylim(0, 100) +
labs(
x = "Ospedale",
y = "Percentuale di parti cesarei",
title = "Percentuale di parti cesarei per ospedale"
) +
theme_minimal()
Il grafico conferma il risultato del test chi-quadrato, infatti la percentale di parti cesarei per ospedale è molto simile.
Dalla letteratura, il peso medio e la lunghezza media di un neonato in Italia è di rispettivamente 3255 g e 496 mm (National, longitudinal NASCITA birth cohort study to investigate the health of Italian children and potential influencing factors; C. Pandolfini, A. Clavenna, M. Cartabia, R. Campi, M. Bonati).
Prima di tutto si controlla se i dati di “Peso”, “Lunghezza” e “Cranio” si distribuiscono secondo una distribuzione normale tramite lo Shapiro test:
shapiro.test(dati$Peso)##
## Shapiro-Wilk normality test
##
## data: dati$Peso
## W = 0.97066, p-value < 2.2e-16shapiro.test(dati$Lunghezza)##
## Shapiro-Wilk normality test
##
## data: dati$Lunghezza
## W = 0.90941, p-value < 2.2e-16shapiro.test(dati$Cranio)##
## Shapiro-Wilk normality test
##
## data: dati$Cranio
## W = 0.96357, p-value < 2.2e-16Al fine di saggiare l’ipotesi che la media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione utilizzerò il t-test poichè la distribuzione dei dati non è normale come si può vedere dai risultati dello Shapiro-Wilk test.
t.test(dati$Peso,
mu = 3255,
alternative = "two.sided",
conf.level = 0.95)##
## One Sample t-test
##
## data: dati$Peso
## t = 2.7694, df = 2499, p-value = 0.005658
## alternative hypothesis: true mean is not equal to 3255
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081t.test(dati$Lunghezza,
mu = 496,
alternative = "two.sided",
conf.level = 0.95)##
## One Sample t-test
##
## data: dati$Lunghezza
## t = -2.4849, df = 2499, p-value = 0.01302
## alternative hypothesis: true mean is not equal to 496
## 95 percent confidence interval:
## 493.6598 495.7242
## sample estimates:
## mean of x
## 494.692Il p-value risultante da entrambi i t-test riportano un p-value minore di 0.05 perciò si rifiuta l’ipotesi nulla e quindi c’è evidenza statistica che la media del nostro campione è diversa da quella della popolazione (valore da letteratura). In particolare il peso medio della popolazione risulta minore all’estremo inferiore dell’intervallo di confidenza dei pesi del nostro campione mentre la lunghezza media della popolazione risulta maggiore all’estremo superiore dell’intervallo di confidenza del nostro campione.
Al fine di saggiare l’ipotesi che le misure antropometriche sono significativamente diverse tra i due sessi si utilizzerà il t-test per confronto tra gruppi indipendenti; anche se i campioni non seguono una distribuzione normale si può procedere considerando la numerosità del campione e dunque l’applicabilità del test è assicurata dal teorema dei limiti centrali.
# Peso alla nascita
t.test(Peso ~ Sesso, data = dati)##
## Welch Two Sample t-test
##
## data: Peso by Sesso
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -287.1051 -207.0615
## sample estimates:
## mean in group F mean in group M
## 3161.132 3408.215# Lunghezza
t.test(Lunghezza ~ Sesso, data = dati)##
## Welch Two Sample t-test
##
## data: Lunghezza by Sesso
## t = -9.582, df = 2459.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -11.929470 -7.876273
## sample estimates:
## mean in group F mean in group M
## 489.7643 499.6672# Diametro cranico
t.test(Cranio ~ Sesso, data = dati)##
## Welch Two Sample t-test
##
## data: Cranio by Sesso
## t = -7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -6.089912 -3.541270
## sample estimates:
## mean in group F mean in group M
## 337.6330 342.4486Tutti i t-test per campioni indipendenti riportano un p-value minore di 0.05 dunque si rifiuta l’ipotesi nulla e perciò abbiamo la conferma che le misure antropometriche sono significativamente diverse tra i due sessi.
Prima di procedere con la creazione del modello vediamo graficamente tramite alcuni scatterplot se c’è correlazione tra le variabili numeriche e la variabile “Peso”, per le stesse variabili viene fatta anche la matrice di correlazione.
plot(dati$Peso,dati$Anni.madre)
plot(dati$Peso,dati$N.gravidanze)
plot(dati$Peso,dati$Gestazione)
plot(dati$Peso,dati$Lunghezza)
plot(dati$Peso,dati$Cranio)
idx <- which(names(dati) %in% barplot_var)
dati_num <- dati[, -idx, drop = FALSE]
cor_mat <- cor(dati_num, use = "complete.obs", method = "pearson")
round(cor_mat, 2)## Anni.madre N.gravidanze Gestazione Peso Lunghezza Cranio
## Anni.madre 1.00 0.38 -0.14 -0.02 -0.06 0.02
## N.gravidanze 0.38 1.00 -0.10 0.00 -0.06 0.04
## Gestazione -0.14 -0.10 1.00 0.59 0.62 0.46
## Peso -0.02 0.00 0.59 1.00 0.80 0.70
## Lunghezza -0.06 -0.06 0.62 0.80 1.00 0.60
## Cranio 0.02 0.04 0.46 0.70 0.60 1.00Come si può notare sia dagli scatterplot che dalla matrice di correlazione le variabile “Anni.madre” e “N.gravidanze” sono indipendenti dalla variabile “Peso”; mentre per le altre variabili (“Gestazione”, “Lunghezza” e “Cranio”) si possono notare trend positivi e perciò sono correlate positivamente con coefficienti di correlazione >0.5.
Di seguito, per la valutazione delle variabili qualitative, si vedranno dei boxplot e dei test di significatività tra gruppi.
dati$Fumatrici <- factor(dati$Fumatrici, levels = c(0, 1), labels = c("No", "Sì"))
vars <- c("Fumatrici", "Tipo.parto", "Ospedale", "Sesso")
for (v in vars) {
# Se è numerica con soli 2 valori (0/1), convertila in factor
if (is.numeric(dati[[v]]) && length(unique(dati[[v]])) == 2) {
dati[[v]] <- factor(dati[[v]])
}
p <- ggplot(dati, aes(x = .data[[v]], y = Peso)) +
geom_boxplot() +
labs(
title = paste("Peso vs", v),
x = v,
y = "Peso"
) +
theme_minimal()
if (v == "Ospedale") {
p <- p + theme(axis.text.x = element_text(angle = 45, hjust = 1))
}
print(p)
}
dati$Fumatrici <- factor(dati$Fumatrici, levels = c("No", "Sì"), labels = c(0, 1))
summary(aov(Peso ~ Ospedale, data = dati))## Df Sum Sq Mean Sq F value Pr(>F)
## Ospedale 2 936237 468118 1.699 0.183
## Residuals 2497 687952305 275512t.test(Peso ~ Fumatrici, data = dati)##
## Welch Two Sample t-test
##
## data: Peso by Fumatrici
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -45.61354 145.22674
## sample estimates:
## mean in group 0 mean in group 1
## 3286.153 3236.346t.test(Peso ~ Sesso, data = dati)##
## Welch Two Sample t-test
##
## data: Peso by Sesso
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -287.1051 -207.0615
## sample estimates:
## mean in group F mean in group M
## 3161.132 3408.215t.test(Peso ~ Tipo.parto, data = dati)##
## Welch Two Sample t-test
##
## data: Peso by Tipo.parto
## t = -0.12968, df = 1493, p-value = 0.8968
## alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
## 95 percent confidence interval:
## -46.27992 40.54037
## sample estimates:
## mean in group Ces mean in group Nat
## 3282.047 3284.916cat(
"Risultati test per Peso:\n",
"Ospedale (ANOVA): p =", 0.183, "\n",
"Fumatrici (t-test): p =", 0.3033, "\n",
"Sesso (t-test): p <", 2.2e-16, "\n",
"Tipo.parto (t-test): p =", 0.8968, "\n"
)## Risultati test per Peso:
## Ospedale (ANOVA): p = 0.183
## Fumatrici (t-test): p = 0.3033
## Sesso (t-test): p < 2.2e-16
## Tipo.parto (t-test): p = 0.8968Si noti che la variabile “Sesso” è l’unica in cui la differenza di peso tra i gruppi è statisticamente significativa dunque sarà l’unica, tra queste, che ci si aspetta di trovare nel modello finale.
Al fine di sviluppare un modello, si partirà dalla crezione di questo inserento tutte le variabile anche se, come abbiamo visto, alcune feature sono ininfluenti sul nostro target. Fatto ciò si procederà con l’ottimizzazione del modello (si proverà ad utilizzare anche la procedura automatica del pacchetto “mass”):
mod1 <- lm(Peso ~ ., data=dati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1124.40 -181.66 -14.42 160.91 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6738.4762 141.3087 -47.686 < 2e-16 ***
## Anni.madre 0.8921 1.1323 0.788 0.4308
## N.gravidanze 11.2665 4.6608 2.417 0.0157 *
## Fumatrici1 -30.1631 27.5386 -1.095 0.2735
## Gestazione 32.5696 3.8187 8.529 < 2e-16 ***
## Lunghezza 10.2945 0.3007 34.236 < 2e-16 ***
## Cranio 10.4707 0.4260 24.578 < 2e-16 ***
## Tipo.partoNat 29.5254 12.0844 2.443 0.0146 *
## Ospedaleosp2 -11.2095 13.4379 -0.834 0.4043
## Ospedaleosp3 28.0958 13.4957 2.082 0.0375 *
## SessoM 77.5409 11.1776 6.937 5.08e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 669.2 on 10 and 2489 DF, p-value: < 2.2e-16rmse <- sqrt(mean(residuals(mod1)^2))
rmse## [1] 273.3222Dal summary del modello generato possiamo vedere che le variabili quantitative significative sono “Lunghezza”, “Cranio” e “Gestazione” (come si era già visto dagli scatterplot e dalla matrice di correlazione) e la variabile qualitativa “Sesso” (come visto dal boxplot e dal t-test). Si ottiene un R2 complessivo del 73% ed un RMSE di 273.3, che non è male, ma nella fase di miglioramento del modello si procederà col rimuovere le variabili superflue al fine di ridurre la complessità del modello mantenendone alta la performance.
Si procede perciò col togliere la variabile “Anni.madre”, “Fumatrici” e “Ospedale”. Si lasciano momentaneamente le variabili “N.Gravidanze” e “Tipo.parto”.
mod2 <- update(mod1,.~.-Fumatrici - Anni.madre -Ospedale)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso, data = dati)
##
## 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-16rmse <- sqrt(mean(residuals(mod2)^2))
rmse## [1] 273.9355Il nuovo modello mantiene un R2 del 73% e un RMSE di 274 però con meno feature e perciò risultando meno complesso; è stato deciso di mantenere la variabile “Tipo.parto” per vedere il suo comportamente con la rimozione delle altre variabili ma in realtà è una variabile correlata a Gestazione poichè è ragionevole supporre che un parto cesareo avvenga con un numero di settimane di gestazione minore in molti casi mentre un parto naturale in condizioni di settimane di gestazione maggiore. Per procedere con metodo si studieranno ora le perfomance del modello ottenuto.
AIC(mod2,mod1)## df AIC
## mod2 8 35175.16
## mod1 12 35171.95BIC(mod2,mod1)## df BIC
## mod2 8 35221.75
## mod1 12 35241.84car::vif(mod2)## N.gravidanze Gestazione Lunghezza Cranio Tipo.parto Sesso
## 1.024171 1.669258 2.080039 1.626199 1.003438 1.040060In seguito alle valutazioni fatte non si nota una grossa differenza tra il modello 1 e il modello 2, infatti nel caso dell’AIC è più performante il modello 1 mentre nel caso del BIC è il contrario, in questo caso si tende a preferire il modello meno complesso perciò il modello 2. Tramite la valutazione della multicollinearità si nota che nessun valore supera il valore soglia di 5 perciò non c’è multicollinearità (sfata la supposizione di dipendenza tra “Gestazione” e “Tipo.parto”).
Si continua con l’ottimizzazione: prima togliendo la variabile “Tipo.parto” e valutando il modello e poi togliendo la variabile “N.gravidanze” e valutando il modello.
mod3 <- update(mod2,.~.-Tipo.parto)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.44 -180.81 -15.58 163.64 2639.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.1445 135.7229 -49.226 < 2e-16 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## Gestazione 32.3321 3.7980 8.513 < 2e-16 ***
## Lunghezza 10.2486 0.3006 34.090 < 2e-16 ***
## Cranio 10.5402 0.4262 24.728 < 2e-16 ***
## SessoM 77.9927 11.2021 6.962 4.26e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1328 on 5 and 2494 DF, p-value: < 2.2e-16rmse <- sqrt(mean(residuals(mod3)^2))
rmse## [1] 274.274R2 continua a rimanere vicino al 73% (72.65%) e l’RMSE abbastanza stabile (274.3) ma con una variabile in meno.
mod4 <- update(mod3,.~.-N.gravidanze)
summary(mod4)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso,
## data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1138.2 -184.3 -17.6 163.3 2627.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6651.1188 135.5172 -49.080 < 2e-16 ***
## Gestazione 31.2737 3.7856 8.261 2.31e-16 ***
## Lunghezza 10.2054 0.3007 33.939 < 2e-16 ***
## Cranio 10.6704 0.4245 25.139 < 2e-16 ***
## SessoM 79.1049 11.2117 7.056 2.22e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared: 0.7261, Adjusted R-squared: 0.7257
## F-statistic: 1654 on 4 and 2495 DF, p-value: < 2.2e-16rmse <- sqrt(mean(residuals(mod4)^2))
rmse## [1] 274.728Anche il modello 4 continua a manenere un R2 intorno al 73% (72.6%) e un RMSE di 274.7.
Si mettono a confronto ora le perfomance dei 4 modelli.
AIC(mod4,mod3,mod2,mod1)## df AIC
## mod4 6 35185.60
## mod3 7 35179.33
## mod2 8 35175.16
## mod1 12 35171.95BIC(mod4,mod3,mod2,mod1)## df BIC
## mod4 6 35220.54
## mod3 7 35220.10
## mod2 8 35221.75
## mod1 12 35241.84car::vif(mod4)## Gestazione Lunghezza Cranio Sesso
## 1.653502 2.069517 1.606131 1.038813Si riconferma che i valori di AIC e BIC tra i modelli sono simili e che il modello AIC tende a preferire il modello 1 mentre il BIC il modello 3/4 ma a parità di questi valori (sono molto simili tra loro) si preferisce il modello meno complesso perciò il modello 4.
Al fine di confermare viene provata la procedura stepwise automatica del pacchetto “MASS” sia con metodo AIC (k=2), sia con metodo BIC (k=loh(2500)).
mod_aic <- MASS::stepAIC(mod1, direction = "both", k = 2)## Start: AIC=28075.26
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 46578 186809099 28074
## - Fumatrici 1 90019 186852540 28075
## <none> 186762521 28075
## - N.gravidanze 1 438452 187200974 28079
## - Tipo.parto 1 447929 187210450 28079
## - Ospedale 2 685979 187448501 28080
## - Sesso 1 3611021 190373542 28121
## - Gestazione 1 5458403 192220925 28145
## - Cranio 1 45326172 232088693 28617
## - Lunghezza 1 87951062 274713583 29038
##
## Step: AIC=28073.88
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 90897 186899996 28073
## <none> 186809099 28074
## + Anni.madre 1 46578 186762521 28075
## - Tipo.parto 1 448222 187257321 28078
## - Ospedale 2 692738 187501837 28079
## - N.gravidanze 1 633756 187442855 28080
## - Sesso 1 3618736 190427835 28120
## - Gestazione 1 5412879 192221978 28143
## - Cranio 1 45588236 232397335 28618
## - Lunghezza 1 87950050 274759149 29036
##
## Step: AIC=28073.1
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 186899996 28073
## + Fumatrici 1 90897 186809099 28074
## + Anni.madre 1 47456 186852540 28075
## - Tipo.parto 1 440684 187340680 28077
## - Ospedale 2 701680 187601677 28079
## - N.gravidanze 1 610840 187510837 28079
## - Sesso 1 3602797 190502794 28119
## - Gestazione 1 5346781 192246777 28142
## - Cranio 1 45632149 232532146 28617
## - Lunghezza 1 88355030 275255027 29039summary(mod_aic)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1113.18 -181.16 -16.58 161.01 2620.19
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.4293 135.9438 -49.340 < 2e-16 ***
## N.gravidanze 12.3619 4.3325 2.853 0.00436 **
## Gestazione 31.9909 3.7896 8.442 < 2e-16 ***
## Lunghezza 10.3086 0.3004 34.316 < 2e-16 ***
## Cranio 10.4922 0.4254 24.661 < 2e-16 ***
## Tipo.partoNat 29.2803 12.0817 2.424 0.01544 *
## Ospedaleosp2 -11.0227 13.4363 -0.820 0.41209
## Ospedaleosp3 28.6408 13.4886 2.123 0.03382 *
## SessoM 77.4412 11.1756 6.930 5.36e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared: 0.7287, Adjusted R-squared: 0.7278
## F-statistic: 836.3 on 8 and 2491 DF, p-value: < 2.2e-16mod_bic <- MASS::stepAIC(mod1, direction = "both", k = log(2500))## Start: AIC=28139.32
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 46578 186809099 28132
## - Fumatrici 1 90019 186852540 28133
## - Ospedale 2 685979 187448501 28133
## - N.gravidanze 1 438452 187200974 28137
## - Tipo.parto 1 447929 187210450 28138
## <none> 186762521 28139
## - Sesso 1 3611021 190373542 28179
## - Gestazione 1 5458403 192220925 28204
## - Cranio 1 45326172 232088693 28675
## - Lunghezza 1 87951062 274713583 29096
##
## Step: AIC=28132.12
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 90897 186899996 28126
## - Ospedale 2 692738 187501837 28126
## - Tipo.parto 1 448222 187257321 28130
## <none> 186809099 28132
## - N.gravidanze 1 633756 187442855 28133
## + Anni.madre 1 46578 186762521 28139
## - Sesso 1 3618736 190427835 28172
## - Gestazione 1 5412879 192221978 28196
## - Cranio 1 45588236 232397335 28670
## - Lunghezza 1 87950050 274759149 29089
##
## Step: AIC=28125.51
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 701680 187601677 28119
## - Tipo.parto 1 440684 187340680 28124
## <none> 186899996 28126
## - N.gravidanze 1 610840 187510837 28126
## + Fumatrici 1 90897 186809099 28132
## + Anni.madre 1 47456 186852540 28133
## - Sesso 1 3602797 190502794 28165
## - Gestazione 1 5346781 192246777 28188
## - Cranio 1 45632149 232532146 28664
## - Lunghezza 1 88355030 275255027 29086
##
## Step: AIC=28119.23
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 463870 188065546 28118
## <none> 187601677 28119
## - N.gravidanze 1 651066 188252743 28120
## + Ospedale 2 701680 186899996 28126
## + Fumatrici 1 99840 187501837 28126
## + Anni.madre 1 54392 187547285 28126
## - Sesso 1 3649259 191250936 28160
## - Gestazione 1 5444109 193045786 28183
## - Cranio 1 45758101 233359778 28657
## - Lunghezza 1 88054432 275656108 29074
##
## Step: AIC=28117.58
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 188065546 28118
## - N.gravidanze 1 623141 188688687 28118
## + Tipo.parto 1 463870 187601677 28119
## + Ospedale 2 724866 187340680 28124
## + Fumatrici 1 91892 187973654 28124
## + Anni.madre 1 54816 188010731 28125
## - Sesso 1 3655292 191720838 28158
## - Gestazione 1 5464853 193530399 28181
## - Cranio 1 46108583 234174130 28658
## - Lunghezza 1 87632762 275698308 29066summary(mod_bic)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.44 -180.81 -15.58 163.64 2639.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.1445 135.7229 -49.226 < 2e-16 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## Gestazione 32.3321 3.7980 8.513 < 2e-16 ***
## Lunghezza 10.2486 0.3006 34.090 < 2e-16 ***
## Cranio 10.5402 0.4262 24.728 < 2e-16 ***
## SessoM 77.9927 11.2021 6.962 4.26e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1328 on 5 and 2494 DF, p-value: < 2.2e-16In entrambi i casi non si arriva ad una conclusione uguale al modello 4, ma osservando gli R2 delle due procedure stepwise AIC e BIC, rispettivamente 72.8% e 72.7%, e confrontandolo con l’R2 del modello 4 (72.6%) si preferisce quest’ultimo poichè con un R2 simile ma meno complesso.
Si procede ora col valutare le interazione tra variabili ed eventuali effetti quadratici:
mod5 <- update(mod4,.~+ (Cranio + Lunghezza + Sesso + Gestazione)^2)
summary(mod5)##
## Call:
## lm(formula = Peso ~ Cranio + Lunghezza + Sesso + Gestazione +
## Cranio:Lunghezza + Cranio:Sesso + Cranio:Gestazione + Lunghezza:Sesso +
## Lunghezza:Gestazione + Sesso:Gestazione, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1135.66 -182.34 -12.52 164.15 2540.69
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.160e+02 1.123e+03 -0.281 0.77852
## Cranio -6.736e+00 7.004e+00 -0.962 0.33628
## Lunghezza 1.171e+01 4.991e+00 2.346 0.01908 *
## SessoM -3.258e+01 2.795e+02 -0.117 0.90722
## Gestazione -1.750e+02 6.351e+01 -2.756 0.00589 **
## Cranio:Lunghezza -9.316e-03 1.402e-02 -0.665 0.50634
## Cranio:SessoM -6.507e-01 8.709e-01 -0.747 0.45502
## Cranio:Gestazione 5.822e-01 2.019e-01 2.884 0.00396 **
## Lunghezza:SessoM 9.677e-01 6.068e-01 1.595 0.11090
## Lunghezza:Gestazione 3.781e-02 1.072e-01 0.353 0.72428
## SessoM:Gestazione -3.898e+00 7.691e+00 -0.507 0.61233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.3 on 2489 degrees of freedom
## Multiple R-squared: 0.7301, Adjusted R-squared: 0.729
## F-statistic: 673.3 on 10 and 2489 DF, p-value: < 2.2e-16rmse <- sqrt(mean(residuals(mod5)^2))
rmse## [1] 272.7078mod6 <- lm(Peso ~ (Cranio + Lunghezza + Gestazione + Sesso)^2 +
I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2),
data=dati)
summary(mod6)##
## Call:
## lm(formula = Peso ~ (Cranio + Lunghezza + Gestazione + Sesso)^2 +
## I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2), data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1153.69 -180.99 -10.79 166.96 1315.06
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.402e+02 1.203e+03 -0.200 0.84176
## Cranio -2.284e+01 1.033e+01 -2.212 0.02709 *
## Lunghezza -5.322e+00 5.980e+00 -0.890 0.37355
## Gestazione 1.803e+02 7.639e+01 2.361 0.01832 *
## SessoM -3.857e+01 2.740e+02 -0.141 0.88809
## I(Cranio^2) 5.042e-02 1.879e-02 2.683 0.00733 **
## I(Lunghezza^2) 5.749e-02 6.712e-03 8.566 < 2e-16 ***
## I(Gestazione^2) -4.659e+00 1.631e+00 -2.858 0.00430 **
## Cranio:Lunghezza -8.438e-02 1.687e-02 -5.001 6.10e-07 ***
## Cranio:Gestazione 1.045e+00 2.091e-01 4.996 6.25e-07 ***
## Cranio:SessoM 2.046e-01 8.577e-01 0.239 0.81145
## Lunghezza:Gestazione -2.925e-01 1.977e-01 -1.480 0.13903
## Lunghezza:SessoM -7.240e-01 6.145e-01 -1.178 0.23882
## Gestazione:SessoM 1.031e+01 7.663e+00 1.346 0.17857
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.1 on 2486 degrees of freedom
## Multiple R-squared: 0.7426, Adjusted R-squared: 0.7413
## F-statistic: 551.8 on 13 and 2486 DF, p-value: < 2.2e-16rmse <- sqrt(mean(residuals(mod6)^2))
rmse## [1] 266.3072mod7 <- update(mod6,.~.-Cranio -Lunghezza - Gestazione - Sesso -Cranio:Sesso - Lunghezza:Gestazione - Lunghezza:Sesso - Gestazione:Sesso)
summary(mod7)##
## Call:
## lm(formula = Peso ~ I(Cranio^2) + I(Lunghezza^2) + I(Gestazione^2) +
## Cranio:Lunghezza + Cranio:Gestazione, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1118.18 -176.14 -7.27 162.56 1255.40
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.995e+03 7.021e+01 -28.412 < 2e-16 ***
## I(Cranio^2) 2.777e-02 1.072e-02 2.591 0.00962 **
## I(Lunghezza^2) 4.780e-02 4.434e-03 10.780 < 2e-16 ***
## I(Gestazione^2) -4.498e+00 7.824e-01 -5.749 1.01e-08 ***
## Cranio:Lunghezza -1.059e-01 1.271e-02 -8.328 < 2e-16 ***
## Cranio:Gestazione 1.132e+00 1.774e-01 6.380 2.10e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.9 on 2494 degrees of freedom
## Multiple R-squared: 0.7363, Adjusted R-squared: 0.7357
## F-statistic: 1392 on 5 and 2494 DF, p-value: < 2.2e-16rmse <- sqrt(mean(residuals(mod7)^2))
rmse## [1] 269.583summary(mod4)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso,
## data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1138.2 -184.3 -17.6 163.3 2627.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6651.1188 135.5172 -49.080 < 2e-16 ***
## Gestazione 31.2737 3.7856 8.261 2.31e-16 ***
## Lunghezza 10.2054 0.3007 33.939 < 2e-16 ***
## Cranio 10.6704 0.4245 25.139 < 2e-16 ***
## SessoM 79.1049 11.2117 7.056 2.22e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared: 0.7261, Adjusted R-squared: 0.7257
## F-statistic: 1654 on 4 and 2495 DF, p-value: < 2.2e-16rmse <- sqrt(mean(residuals(mod4)^2))
rmse## [1] 274.728Si può notare che il beneficio su R2 non è rilevante perciò si opterà per tenere il modello 4 con un R2 del 72.6% (RMSE=274.73 g) avendo una performance simile agli altri modelli ottenuti ma una semplicità maggiore. Si sottolinea comunque il fatto che un RMSE di 275 grammi circa significa un errore potenziale di più o meno questa cifra che comparate con il peso medio (3284 g) significa un errore di più o meno 8.4% che non è poco. Tuttavia l’RMSE degli altri modelli testati è molto simile perciò si continua a preferire la semplicità del modello 4.
mod10 <- update(mod4,.~.-Gestazione -Sesso)
summary(mod4)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso,
## data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1138.2 -184.3 -17.6 163.3 2627.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6651.1188 135.5172 -49.080 < 2e-16 ***
## Gestazione 31.2737 3.7856 8.261 2.31e-16 ***
## Lunghezza 10.2054 0.3007 33.939 < 2e-16 ***
## Cranio 10.6704 0.4245 25.139 < 2e-16 ***
## SessoM 79.1049 11.2117 7.056 2.22e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared: 0.7261, Adjusted R-squared: 0.7257
## F-statistic: 1654 on 4 and 2495 DF, p-value: < 2.2e-16summary(mod10)##
## Call:
## lm(formula = Peso ~ Lunghezza + Cranio, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1295.09 -184.36 -11.95 162.85 2792.61
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6306.9134 124.8791 -50.50 <2e-16 ***
## Lunghezza 11.6312 0.2682 43.37 <2e-16 ***
## Cranio 11.2847 0.4298 26.26 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 281.4 on 2497 degrees of freedom
## Multiple R-squared: 0.7129, Adjusted R-squared: 0.7127
## F-statistic: 3101 on 2 and 2497 DF, p-value: < 2.2e-16Provando a rimuovere invece le variabili rimaste si ottengono i seguenti R2:
- Rimuovendo Cranio: 65.5%
- Rimuovendo Lunghezza: 60%
- Rimuovendo Gestazione: 71.8%
- Rimuovendo Sesso: 72%
Di conseguenza “Cranio” e “LUnghezza” sono le features più impattanti ed importanti mentre “Gestazione” e “Sesso” si potrebbero togliere ma si possono considerare variabili di controllo importanti.
Si procede con un analisi dei residui del modello finale (modello 4) e valutazione della presenza di leverage o outliers (come visto inizialmente probabilmente ce ne saranno):
par(mfrow=c(2,2))
plot(mod4)
shapiro.test(residuals(mod4))##
## Shapiro-Wilk normality test
##
## data: residuals(mod4)
## W = 0.9742, p-value < 2.2e-16lmtest::bptest(mod4)##
## studentized Breusch-Pagan test
##
## data: mod4
## BP = 89.148, df = 4, p-value < 2.2e-16lmtest::dwtest(mod4)##
## Durbin-Watson test
##
## data: mod4
## DW = 1.9557, p-value = 0.1337
## alternative hypothesis: true autocorrelation is greater than 0n=length(dati$Peso)
n## [1] 2500# Leverage
lev <- hatvalues(mod4)
plot(lev)
p <- sum(lev)
soglia = 2*p/n
abline(h=soglia,col=2)
lev[lev>soglia]## 15 34 42 61 67 96
## 0.006144457 0.006237758 0.004252678 0.004587185 0.005345322 0.004801899
## 101 106 117 131 151 155
## 0.007185969 0.012812305 0.004746840 0.006964953 0.010847336 0.006670119
## 190 205 206 220 249 295
## 0.005288500 0.005297016 0.009402903 0.007376867 0.004655679 0.004008437
## 304 305 310 312 315 348
## 0.004419003 0.005382162 0.028757757 0.013126063 0.005342858 0.004187962
## 378 383 445 471 486 492
## 0.015366923 0.004305772 0.007094099 0.004289364 0.004740595 0.008175223
## 565 587 592 615 629 638
## 0.004635950 0.008375235 0.006313937 0.004575867 0.004041054 0.006216787
## 656 666 684 697 702 726
## 0.004735898 0.004345904 0.008750225 0.005809934 0.004719868 0.004038656
## 748 750 765 805 821 895
## 0.008236046 0.006712718 0.006049647 0.014303716 0.004039952 0.005203974
## 928 947 956 964 968 991
## 0.022095348 0.007803378 0.007670034 0.004618460 0.004075295 0.004259145
## 1014 1067 1091 1096 1130 1157
## 0.008221844 0.007868942 0.008927908 0.004268008 0.006561457 0.004080725
## 1166 1181 1188 1200 1238 1248
## 0.004020427 0.005586402 0.006404596 0.005445729 0.005374725 0.014573080
## 1273 1283 1293 1294 1356 1357
## 0.007080183 0.004055883 0.005555616 0.004735838 0.005285646 0.006526609
## 1361 1385 1395 1400 1402 1420
## 0.004080786 0.012017154 0.004646054 0.005559147 0.004788164 0.005130876
## 1428 1429 1551 1553 1556 1560
## 0.008191659 0.021037040 0.048521821 0.006810651 0.005868384 0.004588261
## 1593 1606 1610 1619 1628 1634
## 0.004855489 0.004746852 0.007761683 0.014504316 0.004663751 0.004527422
## 1686 1692 1693 1701 1712 1735
## 0.008715871 0.004260736 0.005041490 0.010197488 0.006945920 0.004370535
## 1780 1802 1806 1809 1827 1858
## 0.025528862 0.004005956 0.004090967 0.008388773 0.006008547 0.004328741
## 1868 1893 1911 1920 1977 2037
## 0.004746360 0.004245302 0.004481835 0.004131544 0.005384125 0.004234559
## 2040 2066 2089 2114 2115 2120
## 0.011429685 0.004013637 0.006252340 0.012850646 0.011763363 0.018002540
## 2140 2149 2175 2200 2215 2216
## 0.006090376 0.013033782 0.021167345 0.011030260 0.004833388 0.007548562
## 2220 2224 2225 2257 2307 2318
## 0.004023907 0.005780067 0.005408611 0.005767952 0.013898144 0.004770755
## 2337 2359 2391 2408 2437 2452
## 0.004700477 0.004750437 0.004386621 0.009632048 0.023757012 0.023227398
## 2458 2476 2478
## 0.008002531 0.004070348 0.005745434# Outliers
plot(rstudent(mod4))
abline(h=c(-2,2))
car::outlierTest(mod4)## rstudent unadjusted p-value Bonferroni p
## 1551 9.986149 4.7193e-23 1.1798e-19
## 155 4.951276 7.8654e-07 1.9663e-03
## 1306 4.781188 1.8440e-06 4.6100e-03# Distanza di Cook
cook <- cooks.distance(mod4)
plot(cook)
In seguito allo Shapiro Ttest si rifiuta l’ipotesi di normalità perciò i
residui non sono distribuiti normalmente ed il risultato del
Breusch-Pagan test ci fa rifiutare l’ipotesi di omoschedasticità. Il
Durbin-Watson test ci conferma l’incorrelazione e questo è positivo. Si
nota nel grafico “Residual vs Leverage” un punto che supera la distanza
di Cook (il punto 1551). Si è dunque approfondita la presenza di punti
di leverage e/o outliers individuandone rispettivamente 135 e 3, i quali
potrebbero essere la causa, o parte di essa, dei risultati negativi
ottenuti durante l’analisi dei residui. In ogni caso, con una quantità
grande di dati, lo Shapiro test non è abbastanza per rifiutare il
modello e graficamente (grafico Residuals vs Fitted) si può notare che i
residui sono distribuiti intorno allo 0 (positivo) e non si notano
pattern inoltre (grafico Q-Q Residuals) seguono la diagonale (positivo).
Il grafico Scale-Location (anche se il Breusch-Pagan test è fallito) non
mostra trend e i residui sono distribuiti in una nuvola (casuali) il che
è positivo; si notano però (anche nel grafico Residuals vs Fitted) dei
punti iniziali fuori dalla nuvola il che è dovuto ai vari valori di
leverage e/o outliers. Si procede con la visualizzazione di questi
valori di leverage e outliers in modo da capire se sono dati credibili o
se contengono errori.
index_lev <- which(lev > soglia)
dati[index_lev,]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 15 33 3 0 34 2400 470 298
## 34 27 0 0 39 3150 480 382
## 42 28 1 0 41 2720 500 310
## 61 34 0 0 39 3620 545 334
## 67 29 0 0 33 2400 450 320
## 96 39 3 0 37 3640 510 376
## 101 31 0 0 34 1370 390 287
## 106 29 4 0 30 1340 400 273
## 117 34 1 0 34 2160 435 303
## 131 30 0 0 34 2290 450 285
## 151 20 0 0 41 2280 450 280
## 155 30 0 0 36 3610 410 330
## 190 26 2 0 39 4050 525 390
## 205 45 2 0 38 3850 505 384
## 206 39 1 0 31 1500 405 295
## 220 23 1 0 40 3520 445 363
## 249 36 1 0 34 2500 440 338
## 295 18 0 0 40 1850 460 305
## 304 36 0 0 37 2060 420 326
## 305 41 2 0 33 2300 450 323
## 310 40 3 0 28 1560 420 379
## 312 26 1 0 32 1280 360 276
## 315 33 2 0 35 2910 450 355
## 348 32 0 0 38 3560 460 360
## 378 27 0 0 28 1285 400 274
## 383 22 1 0 39 3600 550 346
## 445 27 0 0 32 1550 410 289
## 471 29 0 0 40 3680 560 346
## 486 22 0 0 33 2250 445 310
## 492 34 2 0 33 1410 380 295
## 565 24 1 0 40 4250 520 386
## 587 16 1 0 31 1900 440 300
## 592 30 1 0 32 2260 440 322
## 615 27 1 0 35 2100 440 345
## 629 23 0 0 34 2450 475 329
## 638 25 0 0 33 1720 420 300
## 656 38 3 0 41 2320 450 330
## 666 32 1 0 40 4240 555 345
## 684 30 1 0 39 3000 475 390
## 697 30 0 0 39 2820 510 300
## 702 25 0 0 33 2220 450 320
## 726 30 0 0 35 2350 446 344
## 748 35 0 0 33 1390 390 277
## 750 24 0 0 35 1450 405 280
## 765 26 1 0 33 1970 445 300
## 805 30 2 0 29 1190 360 272
## 821 22 0 0 41 3050 495 310
## 895 30 1 0 34 2580 470 347
## 928 25 0 0 28 830 310 254
## 947 34 3 0 32 1615 390 297
## 956 25 0 0 41 2210 430 310
## 964 26 0 0 40 4010 550 335
## 968 23 1 0 39 2900 470 366
## 991 35 2 0 40 4050 550 385
## 1014 17 0 0 37 2050 390 295
## 1067 26 3 0 31 1960 420 300
## 1091 30 1 0 33 1770 410 275
## 1096 37 0 0 34 1750 420 306
## 1130 33 11 0 43 3400 475 360
## 1157 26 0 0 40 4060 505 380
## 1166 36 3 0 39 2950 505 307
## 1181 30 1 0 36 4070 500 373
## 1188 21 0 0 40 4140 550 320
## 1200 21 0 0 40 3200 525 310
## 1238 19 0 0 33 2270 440 315
## 1248 26 1 0 30 1170 370 266
## 1273 32 1 0 33 2040 480 307
## 1283 29 0 0 40 4250 550 376
## 1293 30 3 0 38 4600 485 380
## 1294 24 1 0 40 3250 460 360
## 1356 32 1 0 40 4090 525 390
## 1357 22 0 0 32 2340 445 304
## 1361 25 0 0 39 3570 520 315
## 1385 33 0 0 29 1620 410 292
## 1395 30 2 0 39 3790 505 304
## 1400 22 0 0 34 2590 485 314
## 1402 35 1 0 39 2660 460 364
## 1420 32 1 0 38 3770 530 380
## 1428 30 1 0 36 1280 385 292
## 1429 24 4 0 29 1280 390 355
## 1551 35 1 0 38 4370 315 374
## 1553 30 4 0 35 4520 520 360
## 1556 37 0 0 41 2420 490 300
## 1560 17 0 0 41 2800 455 328
## 1593 41 3 0 35 1500 420 304
## 1606 28 0 0 41 3050 460 352
## 1610 37 3 0 33 2000 470 293
## 1619 31 0 0 31 990 340 278
## 1628 31 0 0 35 2120 410 312
## 1634 32 1 0 38 2500 430 328
## 1686 27 0 0 31 1800 430 308
## 1692 15 0 0 35 2700 470 300
## 1693 25 0 0 41 3100 500 305
## 1701 22 0 0 32 1430 380 301
## 1712 28 0 0 39 3800 520 300
## 1735 15 0 0 37 2750 440 345
## 1780 25 2 0 25 900 325 253
## 1802 27 2 0 39 3600 540 330
## 1806 42 2 0 37 3290 455 355
## 1809 35 0 0 32 1780 420 277
## 1827 32 1 0 33 2100 420 310
## 1858 32 1 0 37 3860 520 371
## 1868 29 0 0 40 3470 525 390
## 1893 22 1 0 34 3030 470 312
## 1911 25 1 0 39 3480 520 312
## 1920 26 0 1 39 4930 550 350
## 1977 39 4 0 34 2970 480 350
## 2037 42 3 0 37 4020 525 368
## 2040 27 1 0 38 3240 410 359
## 2066 30 2 0 41 3540 470 355
## 2089 32 1 1 33 1780 400 305
## 2114 36 0 0 31 1180 355 270
## 2115 35 1 0 32 1890 500 309
## 2120 32 0 0 27 1140 370 267
## 2140 30 2 0 33 1600 410 290
## 2149 39 3 0 30 1300 380 276
## 2175 37 8 0 28 930 355 235
## 2200 33 0 0 30 1750 410 294
## 2215 29 1 0 40 2440 465 298
## 2216 22 0 0 32 2580 470 330
## 2220 23 3 1 41 2620 475 313
## 2224 41 1 0 33 2000 425 312
## 2225 27 0 0 35 3140 465 290
## 2257 40 0 0 34 1580 400 300
## 2307 26 1 0 30 1170 370 273
## 2318 17 0 0 41 3100 500 307
## 2337 31 3 1 37 3440 460 362
## 2359 25 6 0 33 2230 430 313
## 2391 36 1 0 41 4400 565 366
## 2408 37 2 0 31 1690 405 290
## 2437 28 1 0 27 980 320 265
## 2452 28 0 0 26 930 345 245
## 2458 31 0 0 31 1730 430 300
## 2476 19 0 0 40 2910 520 315
## 2478 32 1 0 33 2740 475 324
## Tipo.parto Ospedale Sesso
## 15 Ces osp3 M
## 34 Nat osp1 F
## 42 Ces osp2 M
## 61 Nat osp1 F
## 67 Nat osp2 M
## 96 Ces osp2 M
## 101 Nat osp2 F
## 106 Ces osp1 M
## 117 Ces osp3 M
## 131 Nat osp2 M
## 151 Ces osp3 M
## 155 Nat osp1 M
## 190 Nat osp1 M
## 205 Nat osp3 M
## 206 Nat osp3 M
## 220 Nat osp1 F
## 249 Nat osp3 F
## 295 Nat osp3 F
## 304 Nat osp2 F
## 305 Nat osp1 M
## 310 Nat osp3 F
## 312 Nat osp2 M
## 315 Nat osp3 F
## 348 Nat osp3 M
## 378 Nat osp1 F
## 383 Nat osp2 F
## 445 Nat osp1 F
## 471 Nat osp2 M
## 486 Nat osp3 F
## 492 Nat osp2 F
## 565 Nat osp2 F
## 587 Nat osp2 F
## 592 Ces osp3 F
## 615 Ces osp2 F
## 629 Nat osp1 F
## 638 Nat osp1 M
## 656 Nat osp2 F
## 666 Nat osp2 F
## 684 Ces osp2 F
## 697 Ces osp3 F
## 702 Nat osp1 F
## 726 Nat osp1 F
## 748 Nat osp1 F
## 750 Nat osp1 F
## 765 Nat osp3 M
## 805 Nat osp2 F
## 821 Nat osp1 F
## 895 Nat osp2 M
## 928 Nat osp1 F
## 947 Nat osp3 F
## 956 Nat osp3 F
## 964 Nat osp2 F
## 968 Nat osp1 F
## 991 Nat osp1 M
## 1014 Nat osp2 F
## 1067 Nat osp2 F
## 1091 Nat osp3 M
## 1096 Ces osp3 F
## 1130 Nat osp1 M
## 1157 Ces osp1 F
## 1166 Nat osp1 F
## 1181 Nat osp2 M
## 1188 Ces osp1 M
## 1200 Nat osp1 F
## 1238 Nat osp1 M
## 1248 Nat osp2 M
## 1273 Ces osp1 F
## 1283 Ces osp3 F
## 1293 Nat osp1 M
## 1294 Nat osp2 F
## 1356 Ces osp3 F
## 1357 Nat osp1 F
## 1361 Nat osp1 F
## 1385 Nat osp3 F
## 1395 Ces osp3 M
## 1400 Nat osp2 M
## 1402 Ces osp1 F
## 1420 Nat osp3 F
## 1428 Nat osp2 F
## 1429 Nat osp1 F
## 1551 Nat osp3 F
## 1553 Nat osp2 F
## 1556 Ces osp1 M
## 1560 Nat osp3 M
## 1593 Nat osp1 M
## 1606 Nat osp3 F
## 1610 Ces osp1 F
## 1619 Ces osp2 F
## 1628 Nat osp1 F
## 1634 Nat osp1 M
## 1686 Ces osp3 M
## 1692 Nat osp1 F
## 1693 Nat osp2 F
## 1701 Nat osp1 M
## 1712 Nat osp3 F
## 1735 Nat osp2 M
## 1780 Nat osp3 F
## 1802 Nat osp1 M
## 1806 Nat osp1 M
## 1809 Ces osp1 F
## 1827 Nat osp1 M
## 1858 Nat osp1 M
## 1868 Nat osp1 M
## 1893 Nat osp2 F
## 1911 Nat osp1 M
## 1920 Ces osp2 F
## 1977 Ces osp2 F
## 2037 Ces osp1 M
## 2040 Ces osp1 F
## 2066 Nat osp1 M
## 2089 Ces osp1 F
## 2114 Nat osp3 F
## 2115 Nat osp2 F
## 2120 Nat osp3 F
## 2140 Ces osp1 F
## 2149 Nat osp1 M
## 2175 Nat osp1 F
## 2200 Nat osp2 M
## 2215 Nat osp2 F
## 2216 Nat osp1 F
## 2220 Nat osp3 M
## 2224 Ces osp3 M
## 2225 Nat osp2 F
## 2257 Nat osp2 F
## 2307 Nat osp3 M
## 2318 Ces osp1 M
## 2337 Ces osp1 M
## 2359 Nat osp3 F
## 2391 Nat osp1 F
## 2408 Nat osp2 M
## 2437 Nat osp1 M
## 2452 Ces osp3 F
## 2458 Nat osp3 F
## 2476 Nat osp1 F
## 2478 Ces osp2 Fdati[1551,]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1551 35 1 0 38 4370 315 374
## Tipo.parto Ospedale Sesso
## 1551 Nat osp3 Fdati[155,]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 155 30 0 0 36 3610 410 330
## Tipo.parto Ospedale Sesso
## 155 Nat osp1 Mdati[1306,]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1306 23 0 0 41 4900 510 352
## Tipo.parto Ospedale Sesso
## 1306 Nat osp2 FIn seguito alla visualizzazione delle righe dei valori di leverage e outliers non sono stati trovati dati dei quali si potesse affermare con certezza l’erroneità dei valori presenti perciò il modello finale sarà il modello 4 con i limiti risultati dell’analisi dei residui e dalla presenza di parecchi valori tra leverage e outliers di cui uno sulla soglia d’allarme della distanza di Cook. Sicuramente questi valori contribuiscono anche all’RMSE che porta ad un errore di 8.4% rispetto al peso medio e all’R2 di 72.6% che non è male ma neanche ottimo. In linea di massima il modello si potrà utilizzare se l’errore sul peso predetto è accettabile, nel caso serva più precisione si può utilizzare il modello per avere una stima ma va approfondito a parte.
Di seguito alcune previsioni (si utilizzano le 4 variabili con cui è costruito il modello perciò essendo il numero di gravidanze non significativo per il modello non fa parte delle variabili significative da dare in pasto al modello predittivo; inoltre la sola gestazione non è sufficiente).
predict(mod4, newdata = data.frame(Lunghezza=496, Gestazione=39, Cranio=350, Sesso="F"))## 1
## 3365.072predict(mod4, newdata = data.frame(Lunghezza=496, Gestazione=39, Cranio=350, Sesso="M"))## 1
## 3444.177predict(mod4, newdata = data.frame(Lunghezza=496, Gestazione=39, Cranio=380, Sesso="F"))## 1
## 3685.183Si conclude con alcune visualizzazioni grafiche di come varia la nostra variabile target (Peso) in funzione delle variabili ritenute significative per il modello:
ggplot(data=dati)+
geom_point(aes(x=Lunghezza,
y=Peso,
col=Sesso))+
geom_smooth(aes(x=Lunghezza,
y=Peso,
col=Sesso), se=F, method="lm")## `geom_smooth()` using formula = 'y ~ x'
ggplot(data=dati)+
geom_point(aes(x=Cranio,
y=Peso,
col=Sesso))+
geom_smooth(aes(x=Cranio,
y=Peso,
col=Sesso), se=F, method="lm")## `geom_smooth()` using formula = 'y ~ x'
Sia nella rappresentazione grafica del peso rispetto alla lughezza che rispetto al diametro del cranio si nota come le rette (sesso maschile e sesso femminile) ottenute siano molto vicine, ulteriore conferma del fatto che la variabile “Sesso” potrebbe essere rimossa al fine di semplificare il modello (è stato scelto infatti di tenerla come variabile di controllo).
ggplot(data=dati)+
geom_point(aes(x=Gestazione,
y=Peso,
col=Sesso))+
geom_smooth(aes(x=Gestazione,
y=Peso,
col=Sesso), se=F, method="lm")## `geom_smooth()` using formula = 'y ~ x'
X_num <- dati[, c("Lunghezza", "Cranio", "Gestazione")]
cor(X_num, use = "complete.obs")## Lunghezza Cranio Gestazione
## Lunghezza 1.0000000 0.6033410 0.6189204
## Cranio 0.6033410 1.0000000 0.4608289
## Gestazione 0.6189204 0.4608289 1.0000000Nel grafico del peso rispetto alle settimane di gestazione si nota anche qui la poca influenza del sesso del bambino viste le rette molto vicine; c’è da sottolineare però che, nel caso del diametro del cranio e delle settimane di gestazione, le rette sono si vicine ma parallele con la retta del sesso maschile sopra in entrambi i casi a quella del sesso femminile, dunque è una maggiore conferma della scelta di tenere la variabile “Sesso” come variabile di controllo.
Un ulteriore ricontrollo è stato fatto, tramite la matrice di correlazione, sulla possibile dipendenza tra le variabili (dal momento che logicamente lunghezza, diametro del cranio e settimane di gestazione sicuramente hanno una dipendenza anche se piccola): il risultato, come ci si aspettava, è un fattore di correlazione basso tra diametro del cranio e settimane di gestazione ma un fattore di correlazione del 60% tra lunghezza-diametro del cranio e lunghezza-settimane di gestazione. Si mantiene comunque il modello 4 come modello finale considerando “Lunghezza” e “Cranio” vere e proprio variabili del modello mentre “Sesso” e “Gestazione” variabili di controllo del modello.