Modello previsione peso neonati v2
Riccardo Siega
2024-10-31
1. Raccolta dati e struttura dataset
In questo documento verrà creato un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali. Il dataframe ‘neonati’ contiene informazioni su 2,500 neonati provenienti da tre ospedali. Le variabili raccolte includono:
Età della madre (anni).
Numero di gravidanze (antecedenti l’attuale).
Fumo materno (0=non fumatrice, 1=fumatrice).
Durata della gravidanza (settimane).
Peso neonato alla nascita (grammi).
Lunghezza e cranio (cm), misurabili anche durante la gravidanza (ecografie).
Tipo di parto (Nat, Ces).
Ospedale di nascita (osp1, osp2, osp3).
Sesso del neonato (M, F)
Per prima cosa, dobbiamo importare il dataframe in R. Il working path dovrà essere modificato dall’utente finale per personalizzarlo in base alle proprie esigenze.
path = "/Users/Riccardo/My Drive/Master Data Science/Mod3_Statistica_inferenziale/Cap9_Project"
setwd(path)
neonati <- read.csv("https://drive.google.com/uc?export=download&id=1ChfwftuOSH-WLIto_1AvV-_sQIksGeTq", sep = ",")
attach(neonati)Successivamente, importiamo le librerie e definiamo le funzioni che ci serviranno.
# Librerie
install.packages("ggplot2")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpUF0O3i/downloaded_packagesinstall.packages("moments")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpUF0O3i/downloaded_packagesinstall.packages("dplyr")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpUF0O3i/downloaded_packagesinstall.packages("tidyr")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpUF0O3i/downloaded_packagesinstall.packages("clipr")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpUF0O3i/downloaded_packagesinstall.packages("lubridate")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpUF0O3i/downloaded_packagesinstall.packages("RColorBrewer")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpUF0O3i/downloaded_packagesinstall.packages("knitr")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpUF0O3i/downloaded_packagesinstall.packages("kableExtra")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpUF0O3i/downloaded_packagesinstall.packages("car")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpUF0O3i/downloaded_packagesinstall.packages("lmtest")##
## The downloaded binary packages are in
## /var/folders/sj/npkvv0kx4vq6sn0_l_5fzpfh0000gn/T//RtmpUF0O3i/downloaded_packageslibrary(ggplot2)
library(moments)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(clipr)## Welcome to clipr. See ?write_clip for advisories on writing to the clipboard in R.library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(RColorBrewer)
library(knitr)
library(kableExtra)##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric# Funzioni
CV <- function(x){
return(sd(x)/mean(x)*100)
}2. Analisi descrittiva
Verichiamo per prima cosa la presenza di eventuali valori mancanti nel dataframe importato.
sum(is.na(neonati))## [1] 02.1 Variabili quantitative
var_cont = data.frame(
Variabile = c("Peso","Età Madre", "Gravidanze", "Gestazione","Lunghezza","Cranio"),
Min = c(min(Peso), min(Anni.madre), min(N.gravidanze), min(Gestazione),
min(Lunghezza), min(Cranio)),
Mass = c(max(Peso), max(Anni.madre), max(N.gravidanze), max(Gestazione),
max(Lunghezza), max(Cranio)),
Media = c(round(mean(Peso),0), round(mean(Anni.madre),0), round(mean(N.gravidanze),0),
round(mean(Gestazione),0), round(mean(Lunghezza),0), round(mean(Cranio),0)),
Mediana = c(median(Peso), median(Anni.madre), median(N.gravidanze), median(Gestazione),
median(Lunghezza), median(Cranio)),
CVar = c(round(CV(Peso),0), round(CV(Anni.madre),0), round(CV(N.gravidanze),0),
round(CV(Gestazione),0), round(CV(Lunghezza),0), round(CV(Cranio),0)),
Q1 = c(quantile(Peso, 0.25), quantile(Anni.madre, 0.25), quantile(N.gravidanze, 0.25),
quantile(Gestazione, 0.25), quantile(Lunghezza, 0.25), quantile(Cranio, 0.25)),
Q3 = c(quantile(Peso, 0.75), quantile(Anni.madre, 0.75), quantile(N.gravidanze, 0.75),
quantile(Gestazione, 0.75), quantile(Lunghezza, 0.75), quantile(Cranio, 0.75)),
IQR = c(IQR(Peso), IQR(Anni.madre), IQR(N.gravidanze), IQR(Gestazione),
IQR(Lunghezza), IQR(Cranio)),
Simmetria = c(round(skewness(Peso),1), round(skewness(Anni.madre),1),
round(skewness(N.gravidanze),1), round(skewness(Gestazione),1),
round(skewness(Lunghezza),1), round(skewness(Cranio),1)),
Curtosi = c(round(kurtosis(Peso)-3,1), round(kurtosis(Anni.madre)-3,1),
round(kurtosis(N.gravidanze)-3,1), round(kurtosis(Gestazione)-3,1),
round(kurtosis(Lunghezza)-3,1), round(kurtosis(Cranio)-3,1))
)
kable(var_cont, caption = "Statistiche delle variabili continue") %>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "center")| Variabile | Min | Mass | Media | Mediana | CVar | Q1 | Q3 | IQR | Simmetria | Curtosi |
|---|---|---|---|---|---|---|---|---|---|---|
| Peso | 830 | 4930 | 3284 | 3300 | 16 | 2990 | 3620 | 630 | -0.6 | 2.0 |
| Età Madre | 0 | 46 | 28 | 28 | 19 | 25 | 32 | 7 | 0.0 | 0.4 |
| Gravidanze | 0 | 12 | 1 | 1 | 131 | 0 | 1 | 1 | 2.5 | 11.0 |
| Gestazione | 25 | 43 | 39 | 39 | 5 | 38 | 40 | 2 | -2.1 | 8.3 |
| Lunghezza | 310 | 565 | 495 | 500 | 5 | 480 | 510 | 30 | -1.5 | 6.5 |
| Cranio | 235 | 390 | 340 | 340 | 5 | 330 | 350 | 20 | -0.8 | 2.9 |
2.1.1 Peso
Trattandosi della variabile risposta del modello previsionale, procediamo a verificarne la normalità . L’indice di simmetria è molto basso mentre è evidente la leptocurtosi; procediamo in ogni caso ad effettuare il test di Shapiro-Wilk.
shapiro.test(Peso)##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16Il p-value praticamente nullo ci fa rifiutare l’ipotesi di normalità .
Rappresentando graficamente la curva di densità della variabile rispetto alla normale notiamo effettivamente l’asimmetria negativa e la leptocurtosi.
ggplot(data = neonati, aes(x = Peso)) +
geom_density(fill = "lightblue", color = "black", alpha = 0.7) +
labs(title = "Peso neonato", x = "Peso (g)",
y = "Densità ") +
stat_function(fun = dnorm, args = list(mean = mean(Peso), sd = sd(Peso)),
color = "red", size = 1) +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
2.1.2 Età della madre
Il valore minimo della variabile, zero, è evidentemente errato. In particolare, analizzando il dataframe risultano errate le osservazioni O1152 e O1380.
neonati[Anni.madre < 13, ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1152 1 1 0 41 3250 490 350
## 1380 0 0 0 39 3060 490 330
## Tipo.parto Ospedale Sesso
## 1152 Nat osp2 F
## 1380 Nat osp3 MLe correggiamo entrambe facendo uso del valore mediano della variabile.
neonati$Anni.madre[neonati$Anni.madre < 13] = median(Anni.madre)2.2 Variabili qualitative
var_cat = data.frame(
Variabile = c("Fumatrici", "Parto", "Ospedale", "Sesso"),
Livelli = c(paste(names(table(Fumatrici)), collapse = ", "),
paste(names(table(Tipo.parto)), collapse = ", "),
paste(names(table(Ospedale)), collapse = ", "),
paste(names(table(Sesso)), collapse = ", ")),
Freq_abs = c(paste(table(Fumatrici), collapse = ", "),
paste(table(Tipo.parto), collapse = ", "),
paste(table(Ospedale), collapse = ", "),
paste(table(Sesso), collapse = ", ")),
Freq_rel = c(paste(round(prop.table(table(Fumatrici)) * 100, 1), collapse = ", "),
paste(round(prop.table(table(Tipo.parto)) * 100, 1), collapse = ", "),
paste(round(prop.table(table(Ospedale)) * 100, 1), collapse = ", "),
paste(round(prop.table(table(Sesso)) * 100, 1), collapse = ", "))
)
kable(var_cat, caption = "Statistiche principali variabili categoriche") %>%
kable_styling(bootstrap_options = "striped", full_width = F, position = "center")| Variabile | Livelli | Freq_abs | Freq_rel |
|---|---|---|---|
| Fumatrici | 0, 1 | 2396, 104 | 95.8, 4.2 |
| Parto | Ces, Nat | 728, 1772 | 29.1, 70.9 |
| Ospedale | osp1, osp2, osp3 | 816, 849, 835 | 32.6, 34, 33.4 |
| Sesso | F, M | 1256, 1244 | 50.2, 49.8 |
3. Modello di regressione
3.1 Normalità della variabile risposta
Già analizzata in precedenza, la variabile ‘Peso’ non è normale. L’asimmetria è minima ma è presente leptocurtosi: procederemo nell’analisi nonostante queste deviazioni.
3.2 Matrice di correlazione
neonati_numeric = neonati[sapply(neonati, is.numeric)]
matrix_corr = round(cor(neonati_numeric, use = "complete.obs"),2)
print(matrix_corr)## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## Anni.madre 1.00 0.38 0.01 -0.13 -0.02 -0.06
## N.gravidanze 0.38 1.00 0.05 -0.10 0.00 -0.06
## Fumatrici 0.01 0.05 1.00 0.03 -0.02 -0.02
## Gestazione -0.13 -0.10 0.03 1.00 0.59 0.62
## Peso -0.02 0.00 -0.02 0.59 1.00 0.80
## Lunghezza -0.06 -0.06 -0.02 0.62 0.80 1.00
## Cranio 0.02 0.04 -0.01 0.46 0.70 0.60
## Cranio
## Anni.madre 0.02
## N.gravidanze 0.04
## Fumatrici -0.01
## Gestazione 0.46
## Peso 0.70
## Lunghezza 0.60
## Cranio 1.00Le uniche variabili quantitative che presentano significativi valori di correlazione con ‘Peso’ sono:
Gestazione; correlazione positiva +0.59
Lunghezza: correlazione positiva +0.80
Cranio: correlazione positiva +0.70
Per le variabili qualitative bisogna fare uso di strumenti diversi, come i t-test che consentono di verificare se la media del peso tra differenti modalità della variabile è statisticamente significativa.
Partiamo alla variabile ‘Fumatrici’: con un p-value di 0.30 non rifiutiamo l’iposi nulla di uguaglianza delle medie del peso neonato tra fumatrici e non.
t.test(Peso~Fumatrici)##
## 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.346Passiamo quindi alla variabile ’Tipo.parto’: con un p-value di 0.89 non rifiutiamo l’iposi nulla di uguaglianza delle medie del peso neonato tra parto naturale e cesareo.
t.test(Peso~Tipo.parto)##
## 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.916Per la variabile ‘Ospedale’, che presenta tre modalità , dobbiamo fare uso del pairwise t.test: con un p-value della matrice pari a 1/0.33 non rifiutiamo l’iposi nulla di uguaglianza delle medie del peso neonato tra i vari ospedali.
pairwise.t.test(Peso, Ospedale, paired = F,
pool.sd = T,
p.adjust.method = "bonferroni")##
## Pairwise comparisons using t tests with pooled SD
##
## data: Peso and Ospedale
##
## osp1 osp2
## osp2 1.00 -
## osp3 0.33 0.33
##
## P value adjustment method: bonferroniPassiamo infine alla variabile ‘Sesso’: con un p-value praticamente nullo rifiutiamo l’iposi nulla di uguaglianza delle medie del peso neonato maschio e femmina.
t.test(Peso~Sesso)##
## 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.2153.3 Modello 1: tutte le variabili
mod1 = lm(Peso ~ . , data = neonati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = neonati)
##
## 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 *
## Fumatrici -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-16Variabili maggiormente significative:
Gestazione, Lunghezza, Cranio - coerenti con le considerazioni emerse dall’analisi della matrice di correlazione.
Sesso - coerente con il risultato del t-test sulla differenza delle medie nei due gruppi.
Variabili abbastanza significative:
- Gravidanze, tipo parto e ospedale - non coerenti con l’analisi della matrice di correlazione e i risultati dei t-test precedenti.
3.4 Modello 2: prima semplificazione
Rimuoviamo dal modello:
Anni madre (correlazione molto bassa con peso)
Fumatrici (t-test non rifiuta l’ipotesi di medie del peso uguali nei due gruppi)
mod2 = update(mod1, ~. -Anni.madre -Fumatrici)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = neonati)
##
## 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-163.5 Modello 3: seconda semplificazione
Considerate le evidenze precedenti dei t-test in merito alla differenza di peso tra neonati soggetti a differente tipo parto od ospediale, eliminiamo anche queste due variabili.
mod3 = update(mod2, ~. -Tipo.parto -Ospedale)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = neonati)
##
## 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-16La perdita di capacità predittiva del modello è minima (72.65% vs 72.78%) e tutte le variabili rimaste presentano elevata significatività . In particolare, considerati costanti tutti gli altri fattori:
Ogni gravidanza addizionale aggiunge 12.47g di peso al neonato.
Ogni settimana di gestazione addizionale aggiunge 32.33g di peso al neonato.
Ogni centimetro di lunghezza addizionale aggiunge 10.24g di peso al neonato.
Ogni centimetro di cranio addizionale aggiunge 10.54g di peso al neonato.
Il sesso maschile aggiunge 77.99g di peso al neonato rispetto al femminile.
vif_mod3 = vif(mod3)
vif_df_mod3 = data.frame(Variabile = names(vif_mod3), VIF = vif_mod3)
ggplot(vif_df_mod3, aes(x = reorder(Variabile, VIF), y = VIF, fill = VIF > 5)) +
geom_bar(stat = "identity", color = "black") +
scale_fill_manual(values = c("FALSE" = "blue", "TRUE" = "red"),
labels = c("Accettabile", "Elevato")) +
labs(title = "Valori VIF Mod3",
x = "Variabile",
y = "VIF") +
theme_minimal() +
coord_flip() 
Tutti i VIF sono <5, pertanto non si evidenziano fenomeni di multicollinearità .
3.5 Interazioni & effetti non lineari
Partiamo da un’analisi delle interazioni tra variabili. Solo quella Gestazione-Cranio pare significativa.
mod_interactions = lm(Peso ~ N.gravidanze*Gestazione +N.gravidanze*Lunghezza
+N.gravidanze*Cranio +N.gravidanze*Sesso +Gestazione*Lunghezza
+Gestazione*Cranio +Gestazione*Sesso +Lunghezza*Cranio +Lunghezza*Sesso
+Cranio*Sesso,
data = neonati)
summary(mod_interactions)##
## Call:
## lm(formula = Peso ~ N.gravidanze * Gestazione + N.gravidanze *
## Lunghezza + N.gravidanze * Cranio + N.gravidanze * Sesso +
## Gestazione * Lunghezza + Gestazione * Cranio + Gestazione *
## Sesso + Lunghezza * Cranio + Lunghezza * Sesso + Cranio *
## Sesso, data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1165.82 -180.78 -12.87 163.43 2535.30
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.120e+02 1.185e+03 0.601 0.548092
## N.gravidanze -1.955e+02 8.555e+01 -2.285 0.022377 *
## Gestazione -2.058e+02 6.435e+01 -3.198 0.001402 **
## Lunghezza 1.380e+01 5.158e+00 2.675 0.007517 **
## Cranio -1.197e+01 7.676e+00 -1.560 0.118923
## SessoM -2.164e+01 2.810e+02 -0.077 0.938604
## N.gravidanze:Gestazione -2.579e+00 2.743e+00 -0.940 0.347153
## N.gravidanze:Lunghezza 1.983e-01 2.343e-01 0.847 0.397343
## N.gravidanze:Cranio 6.115e-01 3.568e-01 1.714 0.086721 .
## N.gravidanze:SessoM 7.344e+00 9.101e+00 0.807 0.419770
## Gestazione:Lunghezza 2.884e-03 1.123e-01 0.026 0.979519
## Gestazione:Cranio 7.363e-01 2.149e-01 3.427 0.000621 ***
## Gestazione:SessoM -4.129e+00 7.724e+00 -0.535 0.592996
## Lunghezza:Cranio -1.207e-02 1.413e-02 -0.855 0.392871
## Lunghezza:SessoM 1.169e+00 6.078e-01 1.923 0.054646 .
## Cranio:SessoM -9.735e-01 8.773e-01 -1.110 0.267235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 272.5 on 2484 degrees of freedom
## Multiple R-squared: 0.7322, Adjusted R-squared: 0.7306
## F-statistic: 452.9 on 15 and 2484 DF, p-value: < 2.2e-16Passiamo poi ad una valutazione di eventuali effetti non lineari. Solo gli effetti quadratici di Gestazione e Lunghezza appaiono significativi.
mod_non_linear = lm(Peso ~ N.gravidanze + I(N.gravidanze^2) +
Gestazione + I(Gestazione^2) +
Lunghezza + I(Lunghezza^2) +
Cranio + I(Cranio^2) +
Sesso,
data = neonati)
summary(mod_non_linear)##
## Call:
## lm(formula = Peso ~ N.gravidanze + I(N.gravidanze^2) + Gestazione +
## I(Gestazione^2) + Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^2) +
## Sesso, data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1195.58 -181.99 -11.69 162.33 1469.96
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.075e+03 1.162e+03 -0.925 0.354814
## N.gravidanze 2.822e+01 7.897e+00 3.574 0.000359 ***
## I(N.gravidanze^2) -2.646e+00 1.276e+00 -2.074 0.038228 *
## Gestazione 3.692e+02 6.740e+01 5.478 4.74e-08 ***
## I(Gestazione^2) -4.282e+00 8.839e-01 -4.845 1.35e-06 ***
## Lunghezza -2.912e+01 4.434e+00 -6.567 6.21e-11 ***
## I(Lunghezza^2) 4.061e-02 4.542e-03 8.942 < 2e-16 ***
## Cranio -5.376e+00 9.800e+00 -0.549 0.583360
## I(Cranio^2) 2.328e-02 1.444e-02 1.612 0.107101
## SessoM 7.240e+01 1.099e+01 6.588 5.41e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.2 on 2490 degrees of freedom
## Multiple R-squared: 0.7399, Adjusted R-squared: 0.739
## F-statistic: 787.2 on 9 and 2490 DF, p-value: < 2.2e-16Proviamo a questo punto a raffinare mod3 per tenere conto di quanto scoperto. Purtroppo l’analisi VIF fa emergere evidenti effetti di multicollinearità che rendono il modello instabile.
mod4 = update(mod3, ~. + Gestazione*Cranio +I(Gestazione^2) +I(Lunghezza^2))
summary(mod4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2) + I(Lunghezza^2) + Gestazione:Cranio,
## data = neonati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1184.88 -182.72 -12.83 164.10 1543.27
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.200e+02 1.091e+03 -0.752 0.452317
## N.gravidanze 1.463e+01 4.243e+00 3.448 0.000574 ***
## Gestazione 3.282e+02 6.272e+01 5.233 1.80e-07 ***
## Lunghezza -2.769e+01 4.402e+00 -6.292 3.70e-10 ***
## Cranio -4.173e+00 5.807e+00 -0.719 0.472400
## SessoM 7.179e+01 1.099e+01 6.533 7.81e-11 ***
## I(Gestazione^2) -5.426e+00 1.028e+00 -5.279 1.41e-07 ***
## I(Lunghezza^2) 3.914e-02 4.514e-03 8.670 < 2e-16 ***
## Gestazione:Cranio 3.797e-01 1.504e-01 2.524 0.011656 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.2 on 2491 degrees of freedom
## Multiple R-squared: 0.7399, Adjusted R-squared: 0.7391
## F-statistic: 885.8 on 8 and 2491 DF, p-value: < 2.2e-16vif_mod4 = vif(mod4)## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vifvif_df_mod4 = data.frame(Variabile = names(vif_mod4), VIF = vif_mod4)
ggplot(vif_df_mod4, aes(x = reorder(Variabile, VIF), y = VIF, fill = VIF > 5)) +
geom_bar(stat = "identity", color = "black") +
scale_fill_manual(values = c("FALSE" = "blue", "TRUE" = "red"),
labels = c("Accettabile", "Elevato")) +
labs(title = "Valori VIF Mod4",
x = "Variabile",
y = "VIF") +
theme_minimal() +
coord_flip() 
3.6 Scelta del modello
Non resta che analizzare la bontà dei modelli ottenuti tramite i criteri di informazione Akaike e Bayes. Il primo propende per mod2, il secondo per mod3. Seguiamo la logica del rasoio di Occam e scegliamo mod3 per la sua semplicità .
AIC(mod1,mod2,mod3)## df AIC
## mod1 12 35172.09
## mod2 10 35169.79
## mod3 7 35179.33BIC(mod1,mod2,mod3)## df BIC
## mod1 12 35241.97
## mod2 10 35228.03
## mod3 7 35220.103.6 Qualità del modello
Adjusted R-squared è pario a 72.65%, onestamente non molto alto.
Per la valutazione del parametro MSE meglio riferirsi alla sua radice quadrata, RMSE, da confrontare con la varianza della variabile risposta. RMSE (274)<<sd(Peso) (525), segno che il modello ha una buona capacità predittiva.
Peso_est = predict(mod3, neonati)
RMSE = sqrt(mean((Peso - Peso_est)^2))
RMSE## [1] 274.274sd(Peso)## [1] 525.03873.7 Residui, leverage, outlier e revisione del modello
Passiamo ad un’analisi dei residui, partendo da una valutazione grafica:
I residui si dispongono per la gran parte attorno ad una media nulla.
I residui si allineano per la gran parte al QQplot della normale.
L’andamento della varianza dei residui è abbastanza costante.
par(mfrow=c(2,2))
plot(mod3)
I test statistici di verifica ci dicono che:
Normalità (Shapiro-Wilk): il p-value praticamente nullo ci fa rifiutare l’ipotesi di normalità dei residui. Non essendo normale la stessa variabile ‘Peso’, si tratta di un risultato atteso.
Omoschedasticità (Breusch-Pagan): il p-value praticamente nullo ci fa rifiutare l’ipotesi di omoschedasticità dei residui.
Autocorrelazione (Durbin-Watson): il p-value non ci fa rifiutare l’ipotesi di non correlazione dei residui.
shapiro.test(residuals(mod3))##
## Shapiro-Wilk normality test
##
## data: residuals(mod3)
## W = 0.97408, p-value < 2.2e-16bptest(mod3)##
## studentized Breusch-Pagan test
##
## data: mod3
## BP = 90.253, df = 5, p-value < 2.2e-16dwtest(mod3)##
## Durbin-Watson test
##
## data: mod3
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0L’ultimo grafico della rappresentazione precedente identificava eventuali osservazioni critiche in termini di distanza di Cook, che misura l’influenza combinata di leverage (distanza dei valori indipendenti) e residuo (distanza della predizione dal valore reale). Il superamento della soglia 0.5 avviene solo per O-1551.
cook = cooks.distance(mod3)
plot(cook)
neonati[which(cook > 0.5), ] ## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1551 35 1 0 38 4370 315 374
## Tipo.parto Ospedale Sesso
## 1551 Nat osp3 FProviamo a rimuovere questo leverage dal modello: Adjusted R-squared sale a 73.67%.
neonati_cln = neonati[-1551, ]
mod3_cln = update(mod3, data = neonati_cln)
summary(mod3_cln)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = neonati_cln)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1165.74 -179.59 -12.74 162.89 1410.88
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6683.4142 133.0802 -50.221 < 2e-16 ***
## N.gravidanze 13.1652 4.2557 3.094 0.002 **
## Gestazione 29.5891 3.7340 7.924 3.43e-15 ***
## Lunghezza 10.8927 0.3017 36.109 < 2e-16 ***
## Cranio 9.9187 0.4225 23.476 < 2e-16 ***
## SessoM 78.1348 10.9840 7.114 1.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.3 on 2493 degrees of freedom
## Multiple R-squared: 0.7372, Adjusted R-squared: 0.7367
## F-statistic: 1399 on 5 and 2493 DF, p-value: < 2.2e-16Rivediamo la considerazioni precedenti alla luce della modifica fatta. RMSE (268)<<sd(Peso) (524), segno che il modello ha ancora una buona capacità predittiva.
Peso_est = predict(mod3_cln, neonati_cln)
RMSE = sqrt(mean((neonati_cln$Peso - Peso_est)^2))
RMSE## [1] 268.9331sd(neonati_cln$Peso)## [1] 524.694I test di normalità e non autocorrelazione danno gli stessi risultati del caso precedente, ma quello di omoschedasticità stavolta è al limite del 5% di p-value, un grosso miglioramento.
shapiro.test(residuals(mod3_cln))##
## Shapiro-Wilk normality test
##
## data: residuals(mod3_cln)
## W = 0.98886, p-value = 4.764e-13bptest(mod3_cln)##
## studentized Breusch-Pagan test
##
## data: mod3_cln
## BP = 11.393, df = 5, p-value = 0.04411dwtest(mod3_cln)##
## Durbin-Watson test
##
## data: mod3_cln
## DW = 1.954, p-value = 0.1251
## alternative hypothesis: true autocorrelation is greater than 0Passiamo infine all’analisi degli outlier. Una prima valutazione grafica mostra che ce ne sono in numero significativo.
plot(rstudent(mod3_cln))
abline(h=c(-2,2), col=2)
Il numero totale di outlier è 109, quelli statisticamente significativi sono O-155, O-1306, O-1399.
res_stud = rstudent(mod3_cln)
obs_outlier = which(res_stud < -2 | res_stud > 2)
length(obs_outlier)## [1] 109outlierTest(mod3_cln)## rstudent unadjusted p-value Bonferroni p
## 155 5.287902 1.3448e-07 0.00033607
## 1306 4.942809 8.2119e-07 0.00205210
## 1399 -4.350348 1.4139e-05 0.03533400neonati_cln[c(155,1306,1399), ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 155 30 0 0 36 3610 410 330
## 1306 23 0 0 41 4900 510 352
## 1399 42 2 0 38 2560 525 349
## Tipo.parto Ospedale Sesso
## 155 Nat osp1 M
## 1306 Nat osp2 F
## 1399 Ces osp2 M4. Previsioni e risultati
Il modello finale scelto è quindi mod3_cln, che fa riferimento al dataframe neonati_cln.
Vogliamo ora stimare il peso di una neonata considerando una madre alla terza gravidanza, non fumatrice, che partorirà alla 39ma settimana. Per le variabili regressori non specificate, faremo uso di valori mediani o modali a seconda del caso.
Il modello ci fornisce un valore di peso atteso pari a 3316g.
neonata = data.frame(N.gravidanze = 2, Fumatrici = 0, Gestazione = 39,
Anni.madre = median(neonati_cln$Anni.madre),
Lunghezza = median(neonati_cln$Lunghezza),
Cranio = median(neonati_cln$Cranio), Sesso = "F")
predict(mod3_cln, neonata)## 1
## 3315.6125. Visualizzazioni
Mostriamo infine alcune visualizzazioni relative ai risultati del modello selezionato.
5.1 Gestazione vs Peso (scatterplot)
Visualizza le osservazioni del peso vs settimane di gestazione e sovrappone la linea di regressione del modello, per mostrare come si adatta ai dati.
neonati_cln$Peso_est = predict(mod3_cln, neonati_cln)
ggplot(neonati_cln, aes(x = Gestazione)) +
geom_point(aes(y = Peso), alpha = 0.6, color = "blue") +
geom_smooth(aes(y = Peso_est), method = "lm", col = "red", se = FALSE) +
labs(title = "Gestazione vs Peso",
x = "Settimane di Gestazione", y = "Peso (g)") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
5.2 Valori osservati vs predetti (scatterplot)
Mostra i valori predetti dal modello vs quelli osservati. Se il modello è buono, i punti dovrebbero allinearsi lungo la linea diagonale (y = x).
ggplot(neonati_cln, aes(x = Peso, y = Peso_est)) +
geom_point(color = "blue", alpha = 0.6) +
geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
labs(title = "Valori Predetti vs Osservati",
x = "Peso Osservato", y = "Peso Predetto") +
theme_minimal()
5.2 Distribuzione residui (istogramma)
Un modello buono dovrebbe avere dei residui distribuiti simmetricamente attorno allo zero.
neonati_cln$res = neonati_cln$Peso - neonati_cln$Peso_est
ggplot(neonati_cln, aes(x = res)) +
geom_histogram(bins = 30, fill = "blue", alpha = 0.7) +
labs(title = "Distribuzione dei residui",
x = "Residuo", y = "Frequenza") +
theme_minimal()
5.3 Sesso vs Peso (boxplot)
Confrontiamo i pesi predetti dal modello in base a delle variabili qualitative, come ad esempio il sesso del neonato.
ggplot(neonati_cln, aes(x = Sesso, y = Peso_est, fill = Sesso)) +
geom_boxplot() +
stat_summary(fun = median, geom = "text", aes(label = round(..y.., 0)),
vjust = -0.5, color = "black", size = 3) +
stat_summary(fun = function(x) quantile(x, 0.25),
geom = "text", aes(label = round(..y.., 0)),
vjust = 1.5, color = "blue", size = 3) +
stat_summary(fun = function(x) quantile(x, 0.75),
geom = "text", aes(label = round(..y.., 0)),
vjust = -1.5, color = "blue", size = 3) +
labs(title = "Previsioni Peso per Sesso", x = "Sesso", y = "Peso Predetto (g)") +
theme_minimal()## Warning: The dot-dot notation (`..y..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(y)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.