Neonatal Health Solutions - Modello Statistico per la Previsione del Peso Neonatale
2025-09-09
Importazione ed analisi dati Neonatal Health Solutions
dati <- read.csv("https://docs.google.com/uc?id=1ChfwftuOSH-WLIto_1AvV-_sQIksGeTq&export=download",
sep = ",")
kable(head(dati, 10), digits = 2, caption = "national health solution dataset")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso |
|---|---|---|---|---|---|---|---|---|---|
| 26 | 0 | 0 | 42 | 3380 | 490 | 325 | Nat | osp3 | M |
| 21 | 2 | 0 | 39 | 3150 | 490 | 345 | Nat | osp1 | F |
| 34 | 3 | 0 | 38 | 3640 | 500 | 375 | Nat | osp2 | M |
| 28 | 1 | 0 | 41 | 3690 | 515 | 365 | Nat | osp2 | M |
| 20 | 0 | 0 | 38 | 3700 | 480 | 335 | Nat | osp3 | F |
| 32 | 0 | 0 | 40 | 3200 | 495 | 340 | Nat | osp2 | F |
| 26 | 1 | 0 | 39 | 3100 | 480 | 345 | Nat | osp3 | F |
| 25 | 0 | 0 | 40 | 3580 | 510 | 349 | Nat | osp1 | M |
| 22 | 1 | 0 | 40 | 3670 | 500 | 335 | Ces | osp2 | F |
| 23 | 0 | 0 | 41 | 3700 | 510 | 362 | Ces | osp2 | F |
Ho alcune variabili qualitative, queste le devo trattare come factor, facendo così R le identifica come “variabili label”, nei vari summari e nei modelli se necessario si crea in automatico variabili dummy scegliendo un riferimento e stimando un coefficiente rispetto a tale riferimento
dati$Fumatrici <- factor(dati$Fumatrici, levels=c(0,1), labels=c("No","Sì"))
dati$Tipo.parto <- as.factor(dati$Tipo.parto)
dati$Ospedale <- as.factor(dati$Ospedale)
dati$Sesso <- as.factor(dati$Sesso)Faccio un summary del dataset e calcolo il numero di osservazioni totali
kable(summary(dati),digits = 2, caption = "summary del dataset")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| Min. : 0.00 | Min. : 0.0000 | No:2396 | Min. :25.00 | Min. : 830 | Min. :310.0 | Min. :235 | Ces: 728 | osp1:816 | F:1256 | |
| 1st Qu.:25.00 | 1st Qu.: 0.0000 | Sì: 104 | 1st Qu.:38.00 | 1st Qu.:2990 | 1st Qu.:480.0 | 1st Qu.:330 | Nat:1772 | osp2:849 | M:1244 | |
| Median :28.00 | Median : 1.0000 | NA | Median :39.00 | Median :3300 | Median :500.0 | Median :340 | NA | osp3:835 | NA | |
| Mean :28.16 | Mean : 0.9812 | NA | Mean :38.98 | Mean :3284 | Mean :494.7 | Mean :340 | NA | NA | NA | |
| 3rd Qu.:32.00 | 3rd Qu.: 1.0000 | NA | 3rd Qu.:40.00 | 3rd Qu.:3620 | 3rd Qu.:510.0 | 3rd Qu.:350 | NA | NA | NA | |
| Max. :46.00 | Max. :12.0000 | NA | Max. :43.00 | Max. :4930 | Max. :565.0 | Max. :390 | NA | NA | NA |
nrow <- nrow(dati)Vale la pena di analizzare i dati di anni.madre, per farlo partendo mettendoli in un boxplot
boxplot(dati$Anni.madre) 
Analizzando il boxplot trovo alcuni valori che non sono plausibili, parto da quella più lontana a vado ad analizzarla per vericare se può esser un errore o un valore mancante
kable(dati[dati$Anni.madre == 0,],digits = 2, caption = "osservazioni anni.madre = 0")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1380 | 0 | 0 | No | 39 | 3060 | 490 | 330 | Nat | osp3 | M |
Siccome ho tutti gli altri valori ipotizzo sia un valore mancante e vado a sostituirlo con la media della colonna anni.madre
dati$Anni.madre[dati$Anni.madre == 0] <- mean(dati$Anni.madre)Vado a rifare il summari per controllare che ora tutti i valori siano plausibili
boxplot(dati$Anni.madre)
kable(summary(dati),digits = 2, caption = "summary del dataset")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| Min. : 1.00 | Min. : 0.0000 | No:2396 | Min. :25.00 | Min. : 830 | Min. :310.0 | Min. :235 | Ces: 728 | osp1:816 | F:1256 | |
| 1st Qu.:25.00 | 1st Qu.: 0.0000 | Sì: 104 | 1st Qu.:38.00 | 1st Qu.:2990 | 1st Qu.:480.0 | 1st Qu.:330 | Nat:1772 | osp2:849 | M:1244 | |
| Median :28.00 | Median : 1.0000 | NA | Median :39.00 | Median :3300 | Median :500.0 | Median :340 | NA | osp3:835 | NA | |
| Mean :28.18 | Mean : 0.9812 | NA | Mean :38.98 | Mean :3284 | Mean :494.7 | Mean :340 | NA | NA | NA | |
| 3rd Qu.:32.00 | 3rd Qu.: 1.0000 | NA | 3rd Qu.:40.00 | 3rd Qu.:3620 | 3rd Qu.:510.0 | 3rd Qu.:350 | NA | NA | NA | |
| Max. :46.00 | Max. :12.0000 | NA | Max. :43.00 | Max. :4930 | Max. :565.0 | Max. :390 | NA | NA | NA |
Ho ancora un min che non ha senso quindi faccio un controllo sulle osservazioni che hanno anni.madre < 10
kable(dati[dati$Anni.madre < 10,],digits = 2, caption = "osservazioni anni.madre < 10")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1152 | 1 | 1 | No | 41 | 3250 | 490 | 350 | Nat | osp2 | F |
Siccome ho tutti gli altri valori ipotizzo sia un valore mancante e vado a sostituirlo con la media della colonna anni.madre
dati$Anni.madre[dati$Anni.madre < 10] <- mean(dati$Anni.madre)Vado a rifare il summari per controllare che ora tutti i valori siano plausibili
boxplot(dati$Anni.madre)
kable(summary(dati),digits = 2, caption = "summary del dataset")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| Min. :13.00 | Min. : 0.0000 | No:2396 | Min. :25.00 | Min. : 830 | Min. :310.0 | Min. :235 | Ces: 728 | osp1:816 | F:1256 | |
| 1st Qu.:25.00 | 1st Qu.: 0.0000 | Sì: 104 | 1st Qu.:38.00 | 1st Qu.:2990 | 1st Qu.:480.0 | 1st Qu.:330 | Nat:1772 | osp2:849 | M:1244 | |
| Median :28.00 | Median : 1.0000 | NA | Median :39.00 | Median :3300 | Median :500.0 | Median :340 | NA | osp3:835 | NA | |
| Mean :28.19 | Mean : 0.9812 | NA | Mean :38.98 | Mean :3284 | Mean :494.7 | Mean :340 | NA | NA | NA | |
| 3rd Qu.:32.00 | 3rd Qu.: 1.0000 | NA | 3rd Qu.:40.00 | 3rd Qu.:3620 | 3rd Qu.:510.0 | 3rd Qu.:350 | NA | NA | NA | |
| Max. :46.00 | Max. :12.0000 | NA | Max. :43.00 | Max. :4930 | Max. :565.0 | Max. :390 | NA | NA | NA |
Ho un minimo di 13 anni, prematuro ma possibile, controllo per srupolo osservazioni inferiori a 15
kable(dati[dati$Anni.madre < 15,],digits = 2, caption = "osservazioni anni.madre < 15")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| 138 | 13 | 0 | No | 38 | 2760 | 470 | 325 | Nat | osp2 | F |
| 1075 | 14 | 1 | No | 39 | 3510 | 490 | 365 | Nat | osp2 | M |
| 1532 | 14 | 0 | No | 39 | 3550 | 500 | 355 | Ces | osp1 | M |
Analisi di Anni.madre ha evidenziato tre valori <15 anni, essi sono rari ma rientrano nella plausibilità biologica e non mostrano incongruenze con le altre variabili, scelgo quindi di mantenerli.
Creo una funzione che mi permette graficamente di analizzare tutto il dataset a colpo d’occhio
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}Effettuo un analisi grafica del dataset per verificare eventuali relazioni presenti fra peso nascita e e altre variabili numeriche
num_vars <- names(Filter(is.numeric, dati))
pairs(dati[ , num_vars],
upper.panel = panel.smooth,
lower.panel = panel.cor) 
Analisi feature presenti
Verifico indipendenza fra tipo parto ed ospedale, per poterlo fare devo condurre un test chi^2 ipotesi nulla è che ci sia indipendenza fra le due variabili
ospparto <- table(dati$Ospedale, dati$Tipo.parto)
res <- tidy(chisq.test(ospparto))
kable(res,digits = 2, caption = "check indipendenza parto ospedale")| statistic | p.value | parameter | method |
|---|---|---|---|
| 1.1 | 0.58 | 2 | Pearson’s Chi-squared test |
Il p-value ha un valore di 0.57 per cui non posso rifiutare H0 che è ipotesi di indipendenza, proviamo a visualizzare graficamente
ggpubr::ggballoonplot(data=as.data.frame(ospparto), fill="blue")
Per dimostrare che la media del campione relativa a peso e lunghezza è significativamente uguale a quella della popolazione posso eseguire dei t.test, ipotesi nulla è che la media del campione sia la media della popolazione.
Non posso utilizzare la media del mio campione ma devo andare a recuperare una media del peso e lunghezza dei neonati da un altra fonte.
Da fonte WHO World Health Organization la stima del peso globale alla nascita è di 3300gr, la lungezza media secondo uno studio a livello mondiale portato avanti da Lancet/Villar è di 49.5cm
mu_peso <- 3300
res <- tidy(t.test(dati$Peso,
mu = mu_peso,
conf.level = 1-0.05,
alternative = "two.sided"))
kable(res,digits = 2, caption = "peso: media campione vs media popolazione")| estimate | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|
| 3284.08 | -1.52 | 0.13 | 2499 | 3263.49 | 3304.67 | One Sample t-test | two.sided |
Il t.test ci da un p.value di 0.12 tale valore ci permette di non rifiutare H0 pertanto non ci sono prove statistiche di differenza dalla media del peso della popolazione
my_lun <- 495
res <- tidy(t.test(dati$Lunghezza,
mu = my_lun,
conf.level = 1-0.05,
alternative = "two.sided"))
kable(res,digits = 2, caption = "lunghezza: media campione vs media popolazione")| estimate | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|
| 494.69 | -0.59 | 0.56 | 2499 | 493.66 | 495.72 | One Sample t-test | two.sided |
Il t.test ci da un p.value di 0.55 tale valore ci permette di non rifiutare H0 pertanto non ci sono prove statistiche di differenza dalla media della lunghezza della popolazione
Occorre verificare dipendenza fra misure antropometriche (peso, lunghezza) rispetto al sesso, di fatto siamo abbastanza convinti di questa dipendenza già dalla prima rappresentazione grafica.
par(mfrow=c(1,2))
boxplot(dati$Peso~dati$Sesso)
boxplot(dati$Lunghezza~dati$Sesso) 
Vado inoltre ad eseguire un t.test per verificare che i due gruppi M/F siano significativamente diversi rispetto a pesa e lunghezza, H0 di questo test è che non vi siano differenze fra i due gruppi
res <- tidy(t.test(dati$Peso ~ dati$Sesso, data = dati))
kable(res,digits = 2, caption = "t.test peso in funzione sesso")| estimate | estimate1 | estimate2 | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|
| -247.08 | 3161.13 | 3408.22 | -12.11 | 0 | 2490.72 | -287.11 | -207.06 | Welch Two Sample t-test | two.sided |
res <- tidy(t.test(dati$Lunghezza ~ dati$Sesso, data = dati))
kable(res,digits = 2, caption = "t.test lunghezza in funzione sesso")| estimate | estimate1 | estimate2 | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|
| -9.9 | 489.76 | 499.67 | -9.58 | 0 | 2459.3 | -11.93 | -7.88 | Welch Two Sample t-test | two.sided |
Entrambi i pvalue sono prossimi allo zero, si rifiuta in entrambi i casi l’ipotesi nulla cioè che non vi siano differenze fra i due sessi.
Modelli regressione
Modalità Manuale
Andiamo a creare un modello di regressione utilizzando come primo step tutte le variabili
mod <- lm(dati$Peso~ ., data=dati)
res <- tidy(summary(mod))
fit_stats <- glance(mod)
kable(res,digits = 2, caption = "model features")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -6735.14 | 141.40 | -47.63 | 0.00 |
| Anni.madre | 0.80 | 1.15 | 0.70 | 0.49 |
| N.gravidanze | 11.41 | 4.67 | 2.45 | 0.01 |
| FumatriciSì | -30.16 | 27.54 | -1.10 | 0.27 |
| Gestazione | 32.53 | 3.82 | 8.52 | 0.00 |
| Lunghezza | 10.30 | 0.30 | 34.24 | 0.00 |
| Cranio | 10.47 | 0.43 | 24.58 | 0.00 |
| Tipo.partoNat | 29.50 | 12.08 | 2.44 | 0.01 |
| Ospedaleosp2 | -11.22 | 13.44 | -0.84 | 0.40 |
| Ospedaleosp3 | 28.10 | 13.50 | 2.08 | 0.04 |
| SessoM | 77.55 | 11.18 | 6.94 | 0.00 |
kable(fit_stats,digits = 3, caption = "model summary")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.729 | 0.728 | 273.933 | 669.137 | 0 | 10 | -17574.04 | 35172.09 | 35241.97 | 186772772 | 2489 | 2500 |
Con questo modello si evidenzia che la feature che ha più influenza sul peso è il Sesso, la variabile che non ha pressochè nessuna influenza invece è Anni.madre che andrò a togliere. L’R^2 è di 0.7289
mod2 <- update(mod, ~ .-Anni.madre)
res <- tidy(summary(mod2))
fit_stats <- glance(mod2)
kable(res,digits = 2, caption = "model features, all variables except Anni.madre")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -6708.11 | 135.94 | -49.35 | 0.00 |
| N.gravidanze | 12.61 | 4.34 | 2.91 | 0.00 |
| FumatriciSì | -30.31 | 27.54 | -1.10 | 0.27 |
| Gestazione | 32.25 | 3.80 | 8.49 | 0.00 |
| Lunghezza | 10.29 | 0.30 | 34.24 | 0.00 |
| Cranio | 10.49 | 0.43 | 24.65 | 0.00 |
| Tipo.partoNat | 29.54 | 12.08 | 2.44 | 0.01 |
| Ospedaleosp2 | -11.08 | 13.44 | -0.82 | 0.41 |
| Ospedaleosp3 | 28.37 | 13.49 | 2.10 | 0.04 |
| SessoM | 77.62 | 11.18 | 6.95 | 0.00 |
kable(fit_stats,digits = 3, caption = "model summary, all variables except Anni.madre")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.729 | 0.728 | 273.905 | 743.586 | 0 | 9 | -17574.29 | 35170.57 | 35234.64 | 186809099 | 2490 | 2500 |
R^2 rimane pressochè identico per cui è stato corretto togliere la variabile dato che ho semplificato il modello, proseguo togliendo ospedale e tipo parto che hanno un p-value di significatività troppo alto e che anche logicamente non dovrebbero aver nessuna influenza sul peso.
mod3 <- update(mod2, ~ .-Ospedale)
res <- tidy(summary(mod3))
fit_stats <- glance(mod3)
kable(res,digits = 2, caption = "model features, all variables except Anni.madre and Ospedale")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -6708.07 | 135.98 | -49.33 | 0.00 |
| N.gravidanze | 13.01 | 4.34 | 3.00 | 0.00 |
| FumatriciSì | -31.76 | 27.57 | -1.15 | 0.25 |
| Gestazione | 32.54 | 3.80 | 8.56 | 0.00 |
| Lunghezza | 10.27 | 0.30 | 34.13 | 0.00 |
| Cranio | 10.50 | 0.43 | 24.65 | 0.00 |
| Tipo.partoNat | 30.30 | 12.10 | 2.50 | 0.01 |
| SessoM | 78.11 | 11.19 | 6.98 | 0.00 |
kable(fit_stats,digits = 3, caption = "model summary, all variables except Anni.madre and Ospedale")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.728 | 0.727 | 274.302 | 951.957 | 0 | 7 | -17578.91 | 35175.83 | 35228.24 | 187501837 | 2492 | 2500 |
mod4 <- update(mod3, ~ .-Tipo.parto)
res <- tidy(summary(mod4))
fit_stats <- glance(mod4)
kable(res,digits = 2, caption = "model features, all variables except Anni.madre, Tipo.Parto and Ospedale")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -6681.67 | 135.72 | -49.23 | 0.00 |
| N.gravidanze | 12.72 | 4.35 | 2.93 | 0.00 |
| FumatriciSì | -30.46 | 27.59 | -1.10 | 0.27 |
| Gestazione | 32.59 | 3.81 | 8.57 | 0.00 |
| Lunghezza | 10.23 | 0.30 | 34.01 | 0.00 |
| Cranio | 10.54 | 0.43 | 24.72 | 0.00 |
| SessoM | 78.17 | 11.20 | 6.98 | 0.00 |
kable(fit_stats,digits = 3, caption = "model summary, all variables except Anni.madre, Tipo.Parto and Ospedale")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.727 | 0.726 | 274.592 | 1107.23 | 0 | 6 | -17582.05 | 35180.11 | 35226.7 | 187973654 | 2493 | 2500 |
R^2 scende leggermente 0.7265 per cui ho semplificato il modello a discapito di un minimo di spiegabilità dello stesso.
Una variabile che poc meno influenza sul peso è N.Gravidanze, dato che anche logicamente non dovrebbe aver impatto, provo a toglierla.
mod5 <- update(mod4, ~ .-N.gravidanze)
res <- tidy(summary(mod5))
fit_stats <- glance(mod5)
kable(res,digits = 2, caption = "model features, all variables except Anni.madre, Tipo.Parto, N.gravidanze and Ospedale")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -6651.07 | 135.52 | -49.08 | 0.00 |
| FumatriciSì | -26.36 | 27.60 | -0.96 | 0.34 |
| Gestazione | 31.48 | 3.79 | 8.30 | 0.00 |
| Lunghezza | 10.19 | 0.30 | 33.86 | 0.00 |
| Cranio | 10.67 | 0.42 | 25.13 | 0.00 |
| SessoM | 79.28 | 11.21 | 7.07 | 0.00 |
kable(fit_stats,digits = 3, caption = "model summary, all variables except Anni.madre, Tipo.Parto, N.gravidanze and Ospedale")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.726 | 0.726 | 275.008 | 1322.948 | 0 | 5 | -17586.34 | 35186.69 | 35227.45 | 188619694 | 2494 | 2500 |
R^2 passa a 0.7262 di fatto non perdo spiegabilità per cui non scarto il modello trovato.
Modalità Automatica
modAic <- MASS::stepAIC(
mod,
direction="both",
k=2) ## Start: AIC=28075.39
## dati$Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 36327 186809099 28074
## - Fumatrici 1 89979 186862751 28075
## <none> 186772772 28075
## - Tipo.parto 1 447229 187220001 28079
## - N.gravidanze 1 448852 187221624 28079
## - Ospedale 2 686401 187459173 28081
## - Sesso 1 3611594 190384366 28121
## - Gestazione 1 5446480 192219252 28145
## - Cranio 1 45338159 232110931 28617
## - Lunghezza 1 87959834 274732606 29038
##
## Step: AIC=28073.88
## dati$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 36327 186772772 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
## dati$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 37246 186862751 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 29039res <- tidy(summary(modAic))
fit_stats <- glance(modAic)
kable(res,digits = 2, caption = "model summary autogenerato metodo AIC")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -6707.43 | 135.94 | -49.34 | 0.00 |
| N.gravidanze | 12.36 | 4.33 | 2.85 | 0.00 |
| Gestazione | 31.99 | 3.79 | 8.44 | 0.00 |
| Lunghezza | 10.31 | 0.30 | 34.32 | 0.00 |
| Cranio | 10.49 | 0.43 | 24.66 | 0.00 |
| Tipo.partoNat | 29.28 | 12.08 | 2.42 | 0.02 |
| Ospedaleosp2 | -11.02 | 13.44 | -0.82 | 0.41 |
| Ospedaleosp3 | 28.64 | 13.49 | 2.12 | 0.03 |
| SessoM | 77.44 | 11.18 | 6.93 | 0.00 |
kable(fit_stats,digits = 3, caption = "model statistic model autogenerato metodo AIC")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.729 | 0.728 | 273.916 | 836.312 | 0 | 8 | -17574.89 | 35169.79 | 35228.03 | 186899997 | 2491 | 2500 |
modBic <- MASS::stepAIC(
mod,
direction="both",
k=log(nrow)) ## Start: AIC=28139.46
## dati$Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 36327 186809099 28132
## - Fumatrici 1 89979 186862751 28133
## - Ospedale 2 686401 187459173 28133
## - Tipo.parto 1 447229 187220001 28138
## - N.gravidanze 1 448852 187221624 28138
## <none> 186772772 28140
## - Sesso 1 3611594 190384366 28180
## - Gestazione 1 5446480 192219252 28204
## - Cranio 1 45338159 232110931 28675
## - Lunghezza 1 87959834 274732606 29096
##
## Step: AIC=28132.12
## dati$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 36327 186772772 28140
## - Sesso 1 3618736 190427835 28172
## - Gestazione 1 5412879 192221978 28196
## - Cranio 1 45588236 232397335 28670
## - Lunghezza 1 87950050 274759149 29089
##
## Step: AIC=28125.51
## dati$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 37246 186862751 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
## dati$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 43775 187557901 28127
## - Sesso 1 3649259 191250936 28160
## - Gestazione 1 5444109 193045786 28183
## - Cranio 1 45758101 233359778 28657
## - Lunghezza 1 88054432 275656108 29074
##
## Step: AIC=28117.58
## dati$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 44978 188020568 28125
## - Sesso 1 3655292 191720838 28158
## - Gestazione 1 5464853 193530399 28181
## - Cranio 1 46108583 234174130 28658
## - Lunghezza 1 87632762 275698308 29066res <- tidy(summary(modBic))
fit_stats <- glance(modBic)
kable(res,digits = 2, caption = "model summary autogenerato metodo BIC")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -6681.14 | 135.72 | -49.23 | 0 |
| N.gravidanze | 12.47 | 4.34 | 2.87 | 0 |
| Gestazione | 32.33 | 3.80 | 8.51 | 0 |
| Lunghezza | 10.25 | 0.30 | 34.09 | 0 |
| Cranio | 10.54 | 0.43 | 24.73 | 0 |
| SessoM | 77.99 | 11.20 | 6.96 | 0 |
kable(fit_stats,digits = 3, caption = "model statistic model autogenerato metodo BIC")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.727 | 0.726 | 274.604 | 1328.316 | 0 | 5 | -17582.67 | 35179.33 | 35220.1 | 188065546 | 2494 | 2500 |
Scelta del modello
vado ora a confrontare il migliore fra i due trovati in automatico, testo solo ANOVA, non hanno senso BIA ed AIC
kable(anova(modAic, modBic), caption = "valutazione anova modelli autogenerati")| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 2491 | 186899996 | NA | NA | NA | NA |
| 2494 | 188065546 | -3 | -1165550 | 5.178143 | 0.0014428 |
Dal confronto tramite ANOVA tra i due modelli si osserva che:
Il modello BIC presenta un RSS maggiore rispetto al modello AIC, quindi fornisce un adattamento peggiore ai dati.
Il test ANOVA restituisce un valore di F = 5.18 con p-value = 0.0014, indicando che la differenza fra i modelli non è attribuibile al caso.
Dato che il peggioramento introdotto dal modello BIC è significativo, si rifiuta l’ipotesi nulla di equivalenza tra i due modelli.
IL modello AIC è da preferire in quanto fornisce una migliore descrizione dei dati.
Vado ora a valutare i modelli lineari trovati manualmente
kable(anova(mod, mod2, mod3, mod4, mod5), caption = "valutazione anova modelli trovati manualmente")| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 2489 | 186772772 | NA | NA | NA | NA |
| 2490 | 186809099 | -1 | -36327.1 | 0.4841078 | 0.4866325 |
| 2492 | 187501837 | -2 | -692738.1 | 4.6158362 | 0.0099788 |
| 2493 | 187973654 | -1 | -471817.2 | 6.2876032 | 0.0122217 |
| 2494 | 188619694 | -1 | -646039.3 | 8.6093483 | 0.0033750 |
kable(BIC(mod, mod2, mod3, mod4, mod5), caption = "valutazione BIC modelli trovati manualmente")| df | BIC | |
|---|---|---|
| mod | 12 | 35241.97 |
| mod2 | 11 | 35234.64 |
| mod3 | 9 | 35228.24 |
| mod4 | 8 | 35226.70 |
| mod5 | 7 | 35227.45 |
kable(AIC(mod, mod2, mod3, mod4, mod5), caption = "valutazione AIC modelli trovati manualmente")| df | AIC | |
|---|---|---|
| mod | 12 | 35172.09 |
| mod2 | 11 | 35170.57 |
| mod3 | 9 | 35175.83 |
| mod4 | 8 | 35180.11 |
| mod5 | 7 | 35186.69 |
Dall’analisi comparativa, mod2 risulta il modello statisticamente migliore, in quanto rappresenta un buon compromesso tra adattamento e parsimonia (ANOVA p = 0.486, AIC molto vicino a mod). Dato che però include variabili come tipo di parto e ospedale che, non hanno senso in ottica previsionale e riducono la generalizzabilità del modello scelgo di escluderlo. La scelta finale ricade su mod4, che pur mostrando un lieve peggioramento rispetto a mod2, presenta valori più favorevoli di AIC e BIC rispetto a mod5.
manualmod <- mod4Confronto infine i due miglior candidati
kable(anova(modAic, manualmod), caption = "valutazione anova best candidates")| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 2491 | 186899996 | NA | NA | NA | NA |
| 2493 | 187973654 | -2 | -1073658 | 7.154847 | 0.0007972 |
kable(BIC(modAic, manualmod), caption = "valutazione BIC best candidates")| df | BIC | |
|---|---|---|
| modAic | 10 | 35228.03 |
| manualmod | 8 | 35226.70 |
kable(AIC(modAic, manualmod), caption = "valutazione AIC best candidates")| df | AIC | |
|---|---|---|
| modAic | 10 | 35169.79 |
| manualmod | 8 | 35180.11 |
Di fatto confermo che il modello migliore è modAic, sembra che fumatrici non sia così significativo
Effetti non lineari
Faccio un ultima prova di modello andando ad introdurre possibili effetti non lineari
plot(dati$N.gravidanze, dati$Peso, output.pch=20)
dai grafici la variabile di cui non riesco ad identificare la relazione è NGravidanze per cui provo su di essa una quadratica
modquad <- update(modAic,~. + I(N.gravidanze^2))
res <- tidy(summary(modquad))
fit_stats <- glance(modquad)
kable(res,digits = 2, caption = "model summary valutazione effetto quadratico")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -6716.37 | 135.96 | -49.40 | 0.00 |
| N.gravidanze | 25.08 | 8.06 | 3.11 | 0.00 |
| Gestazione | 32.57 | 3.80 | 8.57 | 0.00 |
| Lunghezza | 10.30 | 0.30 | 34.29 | 0.00 |
| Cranio | 10.45 | 0.43 | 24.53 | 0.00 |
| Tipo.partoNat | 30.09 | 12.08 | 2.49 | 0.01 |
| Ospedaleosp2 | -10.44 | 13.43 | -0.78 | 0.44 |
| Ospedaleosp3 | 28.06 | 13.49 | 2.08 | 0.04 |
| SessoM | 77.66 | 11.17 | 6.95 | 0.00 |
| I(N.gravidanze^2) | -2.43 | 1.30 | -1.87 | 0.06 |
kable(fit_stats,digits = 3, caption = "model statistic valutazione effetto quadratico")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.729 | 0.728 | 273.779 | 744.522 | 0 | 9 | -17573.14 | 35168.28 | 35232.35 | 186637964 | 2490 | 2500 |
Confronto con il best candidate
kable(anova(modAic, modquad), caption = "valutazione anova modello migliore vs modello quadratico")| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 2491 | 186899996 | NA | NA | NA | NA |
| 2490 | 186637964 | 1 | 262032.7 | 3.495867 | 0.0616395 |
kable(BIC(modAic, modquad), caption = "valutazione anova BIC migliore vs modello quadratico")| df | BIC | |
|---|---|---|
| modAic | 10 | 35228.03 |
| modquad | 11 | 35232.35 |
kable(AIC(modAic, modquad), caption = "valutazione anova AIC migliore vs modello quadratico")| df | AIC | |
|---|---|---|
| modAic | 10 | 35169.79 |
| modquad | 11 | 35168.28 |
L’inclusione del termine quadratico non produce un miglioramento significativo del modello per cui anche se statisticamente si evidenzia un lieve vantaggio per il modello quadratico, la differenza osservata non è sufficiente a giustificare l’aggiunta di complessità. Pertanto, in base all’insieme delle evidenze, si conferma la scelta del modello più parsimonioso come il candidato preferibile.
Ultimo scrupolo è provare il quadratico su gestazione, dalla mia esperienza l’aumento in gravidanza non è lineare ed aumenta con l’aumento delle settimane
modquad2 <- update(modAic,~. + I(Gestazione^2))
res <- tidy(summary(modquad2))
fit_stats <- glance(modquad2)
kable(res,digits = 2, caption = "model summary valutazione secondo modello effetto quadratico")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -4804.22 | 896.97 | -5.36 | 0.00 |
| N.gravidanze | 12.45 | 4.33 | 2.88 | 0.00 |
| Gestazione | -74.24 | 49.63 | -1.50 | 0.13 |
| Lunghezza | 10.41 | 0.30 | 34.28 | 0.00 |
| Cranio | 10.58 | 0.43 | 24.77 | 0.00 |
| Tipo.partoNat | 28.85 | 12.07 | 2.39 | 0.02 |
| Ospedaleosp2 | -10.28 | 13.43 | -0.77 | 0.44 |
| Ospedaleosp3 | 28.31 | 13.48 | 2.10 | 0.04 |
| SessoM | 75.39 | 11.21 | 6.73 | 0.00 |
| I(Gestazione^2) | 1.42 | 0.66 | 2.15 | 0.03 |
kable(fit_stats,digits = 3, caption = "model statistic valutazione secondo modello effetto quadratico")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.729 | 0.728 | 273.718 | 744.977 | 0 | 9 | -17572.58 | 35167.17 | 35231.23 | 186554769 | 2490 | 2500 |
Confronto con il best candidate
kable(anova(modAic, modquad2), caption = "valutazione anova modello migliore vs secondo modello quadratico")| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 2491 | 186899996 | NA | NA | NA | NA |
| 2490 | 186554769 | 1 | 345227.4 | 4.607849 | 0.0319223 |
kable(BIC(modAic, modquad2), caption = "valutazione anova BIC migliore vs secondo modello quadratico")| df | BIC | |
|---|---|---|
| modAic | 10 | 35228.03 |
| modquad2 | 11 | 35231.23 |
kable(AIC(modAic, modquad2), caption = "valutazione anova AIC migliore vs secondo modello quadratico")| df | AIC | |
|---|---|---|
| modAic | 10 | 35169.79 |
| modquad2 | 11 | 35167.17 |
L’aggiunta del termine quadratico sulla variabile gestazione produce un miglioramento statisticamente significativo del modello tuttavia, il guadagno informativo non è sostanziale. Pertanto, la preferenza rimane per il modello più parsimonioso, in quanto garantisce un equilibrio migliore tra semplicità e capacità esplicativa.
Dato che nel modello sono presenti variabili come ospedale e tipo parto che logicamente non dovrebbero avere influenza vado a toglierle e fare un rapido confronto fra i modelli.
modAicClear <- update(modAic,~. - Tipo.parto -Ospedale )
res <- tidy(summary(modAicClear))
fit_stats <- glance(modAicClear)
kable(res,digits = 2, caption = "model summary modello pulito da variabili ininfluenti")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -6681.14 | 135.72 | -49.23 | 0 |
| N.gravidanze | 12.47 | 4.34 | 2.87 | 0 |
| Gestazione | 32.33 | 3.80 | 8.51 | 0 |
| Lunghezza | 10.25 | 0.30 | 34.09 | 0 |
| Cranio | 10.54 | 0.43 | 24.73 | 0 |
| SessoM | 77.99 | 11.20 | 6.96 | 0 |
kable(fit_stats,digits = 3, caption = "model statistic modello pulito da variabili ininfluenti")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.727 | 0.726 | 274.604 | 1328.316 | 0 | 5 | -17582.67 | 35179.33 | 35220.1 | 188065546 | 2494 | 2500 |
Confronto con il modello originale
kable(anova(modAic, modAicClear), caption = "valutazione anova modello migliore vs modello migliore pulito")| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 2491 | 186899996 | NA | NA | NA | NA |
| 2494 | 188065546 | -3 | -1165550 | 5.178143 | 0.0014428 |
kable(BIC(modAic, modAicClear), caption = "valutazione anova BIC modello migliore vs modello migliore pulito")| df | BIC | |
|---|---|---|
| modAic | 10 | 35228.03 |
| modAicClear | 7 | 35220.10 |
kable(AIC(modAic, modAicClear), caption = "valutazione anova AIC modello migliore vs modello migliore pulito")| df | AIC | |
|---|---|---|
| modAic | 10 | 35169.79 |
| modAicClear | 7 | 35179.33 |
Poiché i modelli presentano un R² simile, e quindi una capacità esplicativa pressoché equivalente, scelgo di mantenere il modello pulito per ragioni di semplicità e maggiore interpretabilità.
Modello finale
mymodel <- modAicClear| Variabile | Coefficiente (β) | Interpretazione |
|---|---|---|
| N.gravidanze | +12.47 g | Ogni gravidanza in più è associata a un aumento medio di circa 12 g nel peso del neonato. |
| Gestazione | +32.33 g | Ogni settimana in più di gestazione corrisponde a un incremento medio di circa 32 g. |
| Lunghezza | +10.25 g | Ogni unità in più di lunghezza è associata a circa 10 g in più di peso. |
| Cranio | +10.54 g | Ogni unità in più di circonferenza cranica è associata a circa 11 g in più di peso. |
| Sesso (M) | +77.99 g | I neonati maschi pesano in media circa 78 g in più rispetto alle femmine. |
Di fatto emerge che tutte le variabili cliniche considerate nel modello hanno un effetto positivo sul peso neonatale, confermando che i fattori biologici come il sesso sono i più influenti.
Valutazione del modello
Calcoliamo R^2 ed RMSE sul modello
mymodelR2 <- summary(mymodel)$r.squared
mymodelMSE <- mean(residuals(mymodel)^2)
mymodelRMSE <- sqrt(mymodelMSE)
mymodeltargetmean <- mean(dati$Peso)
kable(mymodelR2,digits = 3, caption = "Calcolo R^2")| x |
|---|
| 0.727 |
kable(mymodelMSE,digits = 3, caption = "Calcolo MSE")| x |
|---|
| 75226.22 |
kable(mymodelRMSE,digits = 3, caption = "Calcolo RMSE")| x |
|---|
| 274.274 |
kable(mymodeltargetmean,digits = 3, caption = "Calcolo Media Target")| x |
|---|
| 3284.081 |
Il modello riesce a spiegare in modo piuttosto buono la variabilità dei dati: circa il 72% delle differenze osservate è giustificato dalle variabili incluse. L’errore medio di previsione è di circa 274 grammi, che confrontato con il peso medio osservato pari a 3.284 grammi rappresenta solo l’8% circa del valore medio, il che indica che le stime prodotte dal modello sono ragionevolmente precise e che l’affidabilità complessiva può essere considerata soddisfacente.
Parte Erratica
Andiamo ora ad analizzare che la parte erratica rispetti le assunzioni sui residui - omoschedasticità - incorrelazione - normalità
per prima cosa rappresento graficamente i residui e le loro caratteristiche
par(mfrow=c(2,2))
plot(mymodel) 
eseguiamo qualche test in maniera automatica sulle assunzioni
res <- tidy(shapiro.test(residuals(mymodel)))
kable(res,digits = 6, caption = "summary test normalità sui residui")| statistic | p.value | method |
|---|---|---|
| 0.97408 | 0 | Shapiro-Wilk normality test |
dai primi test sulla normalità si evince una distribuzione intorno ad una media di 1
effettuando uno shapiro test però sono costretto a rifiutare H0 di normalità
res <- tidy(lmtest::bptest(mymodel))
kable(res,digits = 6, caption = "summary test omoschedasticità sui residui")| statistic | p.value | parameter | method |
|---|---|---|---|
| 90.25279 | 0 | 5 | studentized Breusch-Pagan test |
di fatto il pvalue è molto piccolo per cui sono costretto a rifiutare H0 cioè che i residui godono di omoschedasticità (varianza costante)
res <- tidy(lmtest::dwtest(mymodel))
kable(res,digits = 2, caption = "summary test autocorrelazione sui residui")| statistic | p.value | method | alternative |
|---|---|---|---|
| 1.95 | 0.12 | Durbin-Watson test | true autocorrelation is greater than 0 |
il test evidenzia un pvalue di 0.12 e quindi non rifiutiamo H0 cioè che l’assenza di autocorrelazione
L’analisi dei residui evidenzia alcune violazioni delle ipotesi classiche del modello lineare:
- i residui non risultano normalmente distribuiti,
- è presente eteroschedasticità,
- non vi sono invece segnali di autocorrelazione.
Questi risultati suggeriscono che, pur avendo un modello discretamente predittivo, la validità delle inferenze statistiche potrebbe essere parzialmente compromessa.
Outlaiers e Laverage
Andiamo infine a verificare alcuni outlaiers o eventuali valori di leva, andiamo ad identificare i leverage e, tramite la regola standard 2*sum(p)/nrow andiamo a calcolare la soglia di attenzione
leva <- hatvalues(mymodel)
par(mfrow=c(1,1))
soglia <- 2*sum(leva)/nrow
kable(leva[leva>soglia],digits = 2, caption = "valori di leva sopra la soglia")| x | |
|---|---|
| 13 | 0.01 |
| 15 | 0.01 |
| 34 | 0.01 |
| 67 | 0.01 |
| 89 | 0.01 |
| 96 | 0.01 |
| 101 | 0.01 |
| 106 | 0.01 |
| 131 | 0.01 |
| 134 | 0.01 |
| 151 | 0.01 |
| 155 | 0.01 |
| 161 | 0.02 |
| 189 | 0.00 |
| 190 | 0.01 |
| 204 | 0.01 |
| 205 | 0.01 |
| 206 | 0.01 |
| 220 | 0.01 |
| 294 | 0.01 |
| 305 | 0.01 |
| 310 | 0.03 |
| 312 | 0.01 |
| 315 | 0.01 |
| 378 | 0.02 |
| 440 | 0.01 |
| 442 | 0.01 |
| 445 | 0.01 |
| 486 | 0.01 |
| 492 | 0.01 |
| 497 | 0.01 |
| 516 | 0.01 |
| 582 | 0.01 |
| 587 | 0.01 |
| 592 | 0.01 |
| 614 | 0.01 |
| 638 | 0.01 |
| 656 | 0.01 |
| 657 | 0.01 |
| 684 | 0.01 |
| 697 | 0.01 |
| 702 | 0.01 |
| 729 | 0.01 |
| 748 | 0.01 |
| 750 | 0.01 |
| 757 | 0.01 |
| 765 | 0.01 |
| 805 | 0.01 |
| 828 | 0.01 |
| 893 | 0.01 |
| 895 | 0.01 |
| 913 | 0.01 |
| 928 | 0.02 |
| 946 | 0.01 |
| 947 | 0.01 |
| 956 | 0.01 |
| 985 | 0.01 |
| 1008 | 0.01 |
| 1014 | 0.01 |
| 1049 | 0.00 |
| 1067 | 0.01 |
| 1091 | 0.01 |
| 1106 | 0.01 |
| 1130 | 0.03 |
| 1166 | 0.01 |
| 1181 | 0.01 |
| 1188 | 0.01 |
| 1200 | 0.01 |
| 1219 | 0.03 |
| 1238 | 0.01 |
| 1248 | 0.01 |
| 1273 | 0.01 |
| 1291 | 0.01 |
| 1293 | 0.01 |
| 1311 | 0.01 |
| 1321 | 0.01 |
| 1325 | 0.00 |
| 1356 | 0.01 |
| 1357 | 0.01 |
| 1385 | 0.01 |
| 1395 | 0.01 |
| 1400 | 0.01 |
| 1402 | 0.00 |
| 1411 | 0.01 |
| 1420 | 0.01 |
| 1428 | 0.01 |
| 1429 | 0.02 |
| 1450 | 0.02 |
| 1505 | 0.01 |
| 1551 | 0.05 |
| 1553 | 0.01 |
| 1556 | 0.01 |
| 1573 | 0.01 |
| 1593 | 0.01 |
| 1606 | 0.01 |
| 1610 | 0.01 |
| 1617 | 0.00 |
| 1619 | 0.02 |
| 1628 | 0.01 |
| 1686 | 0.01 |
| 1693 | 0.01 |
| 1701 | 0.01 |
| 1712 | 0.01 |
| 1718 | 0.01 |
| 1727 | 0.01 |
| 1735 | 0.00 |
| 1780 | 0.03 |
| 1781 | 0.02 |
| 1809 | 0.01 |
| 1827 | 0.01 |
| 1868 | 0.01 |
| 1892 | 0.01 |
| 1962 | 0.01 |
| 1967 | 0.01 |
| 1977 | 0.01 |
| 2037 | 0.00 |
| 2040 | 0.01 |
| 2046 | 0.01 |
| 2086 | 0.01 |
| 2089 | 0.01 |
| 2098 | 0.01 |
| 2114 | 0.01 |
| 2115 | 0.01 |
| 2120 | 0.02 |
| 2140 | 0.01 |
| 2146 | 0.01 |
| 2148 | 0.01 |
| 2149 | 0.01 |
| 2157 | 0.01 |
| 2175 | 0.03 |
| 2200 | 0.01 |
| 2215 | 0.00 |
| 2216 | 0.01 |
| 2220 | 0.01 |
| 2221 | 0.02 |
| 2224 | 0.01 |
| 2225 | 0.01 |
| 2244 | 0.01 |
| 2257 | 0.01 |
| 2307 | 0.01 |
| 2317 | 0.01 |
| 2318 | 0.00 |
| 2337 | 0.01 |
| 2359 | 0.01 |
| 2408 | 0.01 |
| 2422 | 0.02 |
| 2436 | 0.00 |
| 2437 | 0.02 |
| 2452 | 0.02 |
| 2458 | 0.01 |
| 2471 | 0.02 |
| 2478 | 0.01 |
andiamo ora ad evidenziare le osservazioni che presentano leva alta, fondamentale andare a fare ulteriori analisi
plot(rstudent(mymodel))
abline(h=c(-2,2),col=2)
out <- car::outlierTest(mymodel)
out_df <- as.data.frame(out$rstudent)
kable(out_df, caption = "Outlier test (Bonferroni)")| out$rstudent | |
|---|---|
| 1551 | 10.051908 |
| 155 | 5.027798 |
| 1306 | 4.827238 |
a questo punto vengono restituite le tre osservazioni che sono considerate outlaiers all’interno del dataset fornito, verifichiamo tramite distanza di cook se influenza il modello
cook <- cooks.distance(mymodel)
plot(cook)
kable(max(cook),digits = 2, caption = "Cook distance maggiore")| x |
|---|
| 0.83 |
con un valore di 0.83 di fatto l’osservazione influenza il modello, quello che sarebbe corretto fare è andare a verificarla
which.max(cook)## 1551
## 1551kable(dati[which.max(cook),],digits = 2, caption = "osservazione con Cook distance maggiore")| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1551 | 35 | 1 | No | 38 | 4370 | 315 | 374 | Nat | osp3 | F |
Analizzando questa osservazione si evince una particolarità, il sesso è F ma il peso è 4370, questo può rappresentare sia un outlaier influente che un errore di battitura rispetto al sesso.
Non avendo accesso ai dati grezzi provo ad andare ad aggiornare il modello eliminando l’osservazione problematica
mymodel_clean <- update(mymodel, subset = -which.max(cook))
res <- tidy(summary(mymodel_clean))
fit_stats <- glance(mymodel_clean)
kable(res,digits = 2, caption = "model summary modello pulito da osservazione 1551")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -6683.41 | 133.08 | -50.22 | 0 |
| N.gravidanze | 13.17 | 4.26 | 3.09 | 0 |
| Gestazione | 29.59 | 3.73 | 7.92 | 0 |
| Lunghezza | 10.89 | 0.30 | 36.11 | 0 |
| Cranio | 9.92 | 0.42 | 23.48 | 0 |
| SessoM | 78.13 | 10.98 | 7.11 | 0 |
kable(fit_stats,digits = 3, caption = "model statistic modello pulito da osservazione 1551")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.737 | 0.737 | 269.257 | 1398.552 | 0 | 5 | -17526.49 | 35066.98 | 35107.74 | 180740178 | 2493 | 2499 |
andiamo a ripetere le analisi fatte sul modello aggiornato
mymodel <- mymodel_clean
par(mfrow=c(2,2))
plot(mymodel)
leva <- hatvalues(mymodel)
par(mfrow=c(1,1))
soglia <- 2*sum(leva)/nrow
abline(h=soglia, col=2)
plot(rstudent(mymodel))
abline(h=c(-2,2),col=2)
out <- car::outlierTest(mymodel)
out_df <- as.data.frame(out$rstudent)
kable(out_df, caption = "Outlier test (Bonferroni)")| out$rstudent | |
|---|---|
| 155 | 5.287902 |
| 1306 | 4.942809 |
| 1399 | -4.350348 |
cook <- cooks.distance(mymodel)
plot(cook)
max(cook)## [1] 0.06511638nonostante vi siano ancora degli outlaiers on distanza di cook massima di 0.06 non ho nessuna osservazione che possa influenzare pesantemente il modello
res <- tidy(shapiro.test(residuals(mymodel)))
kable(res,digits = 6, caption = "summary test normalità sui residui")| statistic | p.value | method |
|---|---|---|
| 0.98886 | 0 | Shapiro-Wilk normality test |
continuo a dover rifiutare H0 cioè che i residui siano normalmente distribuiti
res <- tidy(lmtest::bptest(mymodel))
kable(res,digits = 6, caption = "summary test omoschedasticità sui residui")| statistic | p.value | parameter | method |
|---|---|---|---|
| 11.3935 | 0.044113 | 5 | studentized Breusch-Pagan test |
il pvalue è di poco al di sotto del limite, per cui dovrei rifiutare H0, il modello soffre ancora di un minimo di eteroschedasticità
res <- tidy(lmtest::dwtest(mymodel))
kable(res,digits = 2, caption = "summary test autocorrelazione sui residui")| statistic | p.value | method | alternative |
|---|---|---|---|
| 1.95 | 0.13 | Durbin-Watson test | true autocorrelation is greater than 0 |
il pvalue di 0.13 è sopra il limite quindi non rifiutiamo H0 cioè che l’assenza di autocorrelazione
In conclusione scelgo di adottare il modello pulito dall’osservazione 1551, essa è stata esclusa dall’analisi in quanto caratterizzata da un valore elevato di Cook’s distance, tale da esercitare un’influenza sproporzionata sulla stima dei coefficienti. La presenza di un peso alla nascita anomalo rispetto al sesso suggerisce la possibilità di un errore di rilevazione o, di un outlier estremo. La sua rimozione consente di ottenere un modello più stabile, privo di osservazioni influenti, pur mantenendo invariata la capacità esplicativa complessiva
Stima e valutazioni risultati
proviamo ad andare in stima considerando una madre alla terza gravidanza che partorirà alla 39esima settimana
record <- data.frame(
N.gravidanze = 3,
Gestazione = 39,
Lunghezza = mean(dati$Lunghezza),
Cranio = mean(dati$Cranio),
Tipo.parto = "Nat",
Ospedale = "osp1",
Sesso = "F"
)punto di stima, intervallo di confidenza e predizione
stima <- predict(mymodel, newdata = record) visualizzo la stima graficamente rispetto alle variabili significative del modello
par(mfrow=c(3,3))
plot(dati$N.gravidanze, dati$Peso,
xlab="N.gravidanze", ylab="Peso (g)",
pch=20, col="grey")
points(record$N.gravidanze, stima, col="red", pch=20, cex=1.5)
plot(dati$Gestazione, dati$Peso,
xlab="Gestazione", ylab="Peso (g)",
pch=20, col="grey")
points(record$Gestazione, stima, col="red", pch=20, cex=1.5)
plot(dati$Lunghezza, dati$Peso,
xlab="Lunghezza", ylab="Peso (g)",
pch=20, col="grey")
points(record$Lunghezza, stima, col="red", pch=20, cex=1.5)
plot(dati$Cranio, dati$Peso,
xlab="Cranio", ylab="Peso (g)",
pch=20, col="grey")
points(record$Cranio, stima, col="red", pch=20, cex=1.5)
boxplot(dati$Peso ~ dati$Tipo.parto, data = dati,
xlab = "Tipo Parto", ylab = "Peso (g)",
col = c("lightblue","pink"))
points(record$Tipo.parto, stima, col="red", pch=20, cex=1.5)
boxplot(dati$Peso ~ dati$Ospedale, data = dati,
xlab = "Ospedale", ylab = "Peso (g)",
col = c("lightblue","pink"))
points(record$Ospedale, stima, col="red", pch=20, cex=1.5)
boxplot(dati$Peso ~ dati$Sesso, data = dati,
xlab = "Sesso", ylab = "Peso (g)",
col = c("lightblue","pink"))
points(record$Sesso, stima, col="red", pch=20, cex=1.5)
boxplot(dati$Peso ~ dati$Fumatrici, data = dati,
xlab = "Fumatrici", ylab = "Peso (g)",
col = c("lightblue","pink"))
points(record$Sesso, stima, col="red", pch=20, cex=1.5)
dopo aver visualizzato il trend di tutte le features rispetto a quella target vado a riassumerle graficamente in un unico grafico
car::crPlots(mymodel)
come ulteriore riassunto vado a graficare in forma tabellare l’influenza delle features rispetto al target
par(mfrow=c(1,1))
corrplot::corrplot(
cor(dati[, c("Peso","N.gravidanze","Gestazione","Lunghezza","Cranio")]),
method="color", addCoef.col="black"
)
Extra Valutazione Fumatrici
Nel mio modello non ho incluso fumatori come features, tutte le strategie suggerivano di toglierlo, provo comunque a reinserirlo dato che è una delle features protagoniste dello studio
mymodel_fumatrici <- update(mymodel,~. +dati$Fumatrici)
res <- tidy(summary(mymodel_fumatrici))
fit_stats <- glance(mymodel_fumatrici)
kable(res,digits = 2, caption = "model summary modello con fumatrici")| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | -6659.81 | 135.50 | -49.15 | 0.00 |
| N.gravidanze | 12.84 | 4.33 | 2.96 | 0.00 |
| Gestazione | 30.90 | 3.82 | 8.08 | 0.00 |
| Lunghezza | 10.20 | 0.30 | 33.95 | 0.00 |
| Cranio | 10.72 | 0.43 | 25.04 | 0.00 |
| SessoM | 77.66 | 11.18 | 6.95 | 0.00 |
| dati$FumatriciSì | -30.37 | 27.52 | -1.10 | 0.27 |
kable(fit_stats,digits = 3, caption = "model statistic modello con fumatrici")| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.727 | 0.727 | 273.893 | 1108.571 | 0 | 6 | -17568.65 | 35153.3 | 35199.89 | 186943065 | 2492 | 2499 |
Confronto modello con o senza fumatrici
kable(anova(mymodel, mymodel_fumatrici), caption = "valutazione anova modello senza e con fumatrici")| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 2493 | 180740178 | NA | NA | NA | NA |
| 2492 | 186943065 | 1 | -6202887 | NA | NA |
kable(BIC(mymodel, mymodel_fumatrici), caption = "valutazione BIC modello senza e con fumatrici")| df | BIC | |
|---|---|---|
| mymodel | 7 | 35107.74 |
| mymodel_fumatrici | 8 | 35199.89 |
kable(AIC(mymodel, mymodel_fumatrici), caption = "valutazione AIC modello senza e con fumatrici")| df | AIC | |
|---|---|---|
| mymodel | 7 | 35066.98 |
| mymodel_fumatrici | 8 | 35153.30 |
anche rieffettuando l’analisi non ho nessuna evidenza che esser fumatrice abbia un impatto sul peso del neonato per cui non aggiorno il modello
Conclusioni
Il modello sviluppato, basato su variabili cliniche raccolte da tre
ospedali, ha mostrato una buona capacità predittiva del peso neonatale
alla nascita (R^2 = 0.737).
L’analisi ha evidenziato che fattori quali:
- Lunghezza
- Dimensioni craniche
- Durata della gestazione
- Numero di gravidanze pregresse
- Sesso del neonato
hanno un impatto significativo e positivo sul peso alla
nascita.
Al contrario, la variabile Fumatrici non è risultata
influente nel campione analizzato: un risultato particolarmente
interessante, poiché in controtendenza rispetto a quanto frequentemente
riportato in letteratura, il risultato potrebbe anche esser legato ad un
problema sul campione, infatti nel campione la percentuale di fumatrici
è circa il 4% mentre la media delle fumatrici femminili >15 anni è
circa il 16.8% (https://globalactiontoendsmoking.org/research/tobacco-around-the-world/italy/?utm_source=chatgpt.com)