Previsione del Peso Neonatale
Eleonora Cagnani
2025-03-19
dati <- read.csv("neonati.csv",stringsAsFactors = TRUE)
str(dati)## 'data.frame': 2500 obs. of 10 variables:
## $ Anni.madre : int 26 21 34 28 20 32 26 25 22 23 ...
## $ N.gravidanze: int 0 2 3 1 0 0 1 0 1 0 ...
## $ Fumatrici : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Gestazione : int 42 39 38 41 38 40 39 40 40 41 ...
## $ Peso : int 3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
## $ Lunghezza : int 490 490 500 515 480 495 480 510 500 510 ...
## $ Cranio : int 325 345 375 365 335 340 345 349 335 362 ...
## $ Tipo.parto : Factor w/ 2 levels "Ces","Nat": 2 2 2 2 2 2 2 2 1 1 ...
## $ Ospedale : Factor w/ 3 levels "osp1","osp2",..: 3 1 2 2 3 2 3 1 2 2 ...
## $ Sesso : Factor w/ 2 levels "F","M": 2 1 2 2 1 1 1 2 1 1 ...attach(dati)
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 Ospedale Sesso
## Min. : 830 Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1256
## 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Nat:1772 osp2:849 M:1244
## Median :3300 Median :500.0 Median :340 osp3:835
## Mean :3284 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :565.0 Max. :390n<-nrow(dati)
n## [1] 2500ANALISI PRELIMINARE
Dopo aver caricato i dati ( n = 2500), si nota che alcune variabili presentano valori anomali:
L’età della madre ha un minimo pari a 0, che è un valore inverosimile.
La lunghezza minima di 310 mm e il peso minimo di 830 g sono compatibili con nati prematuri gravi, ma vanno verificati nel contesto clinico.
Si estraggono i record contenenti i valori anomali: a partire dal dato sull’età materna, con lo scopo di sostituirla con un’età plausibile.
#estraggo caso sospetto eta madre
dati[Anni.madre == 0, ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1380 0 0 0 39 3060 490 330
## Tipo.parto Ospedale Sesso
## 1380 Nat osp3 Mknitr::kable(dati[Anni.madre == 0, ])| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1380 | 0 | 0 | 0 | 39 | 3060 | 490 | 330 | Nat | osp3 | M |
median(Anni.madre) #28## [1] 28dati$Anni.madre[dati$Anni.madre == 0] <- median(dati$Anni.madre[dati$Anni.madre > 0])Il valore anomalo dell’età materna pari a 0 è stato sostituito con la mediana dei valori validi di età. La mediana è stata scelta in quanto misura più robusta rispetto alla media poichè non risente dei valori estremi.
Successivamente si osserva il record che contiene la lunghezza del neonato pari a 310 mm:
#estraggo caso sospetto Lunghezza e Peso
#dati[Lunghezza == 310, ]#| Peso == 830, ]
knitr::kable(dati[Lunghezza == 310, ])| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| 928 | 25 | 0 | 0 | 28 | 830 | 310 | 254 | Nat | osp1 | F |
Nonostante i valori estremi, peso e lunghezza risultano coerenti con una gestazione di 28 settimane e cranio molto piccolo (254 mm). Pertanto, il record è stato mantenuto inalterato nel dataset in quanto rappresentativo di un caso clinico plausibile.
Si prosegue osservando la distribuzione della variabile target peso:
plot(density(Peso))
La curva di densità del peso neonatale mostra una distribuzione
leggermente asimmetrica a destra, con una coda più lunga verso valori
alti.
Per indagare le relazioni tra le variabili quantitative, si costruisce un grafico di correlazione. Questo consente di valutare possibili problemi di multicollinearità tra i regressori, oltre ad identificare le variabili maggiormente associate alla variabile target.
#x11()
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)
}
pairs(dati, upper.panel = panel.smooth, lower.panel = panel.cor)
Le variabili maggiormente correlate al peso neonatale sono la lunghezza (r = 0.80), il cranio (r = 0.70) e la gestazione (r = 0.59), suggerendo una forte associazione lineare positiva con il peso. Questo implica, ad esempio, che all’aumentare della lunghezza il peso tende ad aumentare, così come neonati con cranio più largo o nati più tardi tendono ad avere un peso maggiore. La lunghezza e il diametro cranico risultano fortemente correlate tra loro (r = 0.60).
Variabili come l’età della madre e il numero di gravidanze, invece, mostrano una correlazione debole con il peso. Anche il fumo materno non mostra una chiara relazione lineare, ma il suo effetto potrebbe manifestarsi in modo non lineare o attraverso interazioni. Tuttavia, queste variabili non verrano escluse a priori ma verranno comunque considerate nelle prime fasi dell’analisi per esplorare eventuali effetti non lineari o interazioni.
Per le variabili categoriali, si procederà con boxplot e test t per confrontare i gruppi e identificare eventuali differenze significative nel peso.
I boxplot permettono di evidenziare differenze mediane, dispersione e la presenza di outlier tra i gruppi, mentre i test-t consentono di verificare se tali differenze siano statisticamente significative.
par(mfrow=c(2,2))
boxplot(Peso)
boxplot(Peso~Tipo.parto)
boxplot(Peso~Sesso)
boxplot(Peso~Ospedale)
median(Peso[Sesso=="M"]) #[1] 3430## [1] 3430median(Peso[Sesso=="F"]) #[1] 3160## [1] 3160I boxplot mostrano che, condizionando il peso al sesso, al tipo di parto e all’ospedale, si evidenziano alcune differenze interessanti. I neonati maschi sembrano avere un peso leggermente maggiore rispetto alle femmine. Le differenze tra i tipi di parto e tra i tre ospedali risultano contenute, ma si notano outlier evidenti, in particolare per neonati con peso molto basso.
A confermare la tendenza emersa dal boxplot, la mediana del peso nei neonati maschi (3430 g) è sensibilmente più alta rispetto a quella delle femmine (3160 g).
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.215t.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.916t.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.346par(mfrow=c(1,2))
boxplot(Peso)
boxplot(Peso~Fumatrici)
boxplot(Peso~Tipo.parto)
t.test(Lunghezza~Sesso)##
## 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.6672t.test(Cranio~Sesso)##
## 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.4486boxplot(Lunghezza~Sesso)
boxplot(Cranio~Sesso)
I test t confermano che la differenza tra i sessi è statisticamente significativa (p < 2.2e-16), mentre non emergono differenze significative tra parto cesareo e naturale (p ≈ 0.90) né tra figli di madri fumatrici e non fumatrici.
Anche le misure antropometriche alla nascita (lunghezza e circonferenza cranica) risultano significativamente diverse tra i sessi. I maschi, in media, hanno una lunghezza di circa 499.7 mm contro 489.8 mm nelle femmine, e una circonferenza cranica di 342.4 mm rispetto a 337.6 mm. I p-value molto bassi confermano la significatività statistica di tali differenze.
table(Tipo.parto,Ospedale) ## Ospedale
## Tipo.parto osp1 osp2 osp3
## Ces 242 254 232
## Nat 574 595 603chisq.test(table(Tipo.parto,Ospedale))##
## Pearson's Chi-squared test
##
## data: table(Tipo.parto, Ospedale)
## X-squared = 1.0972, df = 2, p-value = 0.5778Attraverso un test del chi quadrato su tabella di contingenza tra tipo di parto e ospedale, è stato verificato se la frequenza di parti cesarei differisce tra le tre strutture. Il test restituisce un p-value = 0.5778, molto superiore alla soglia di 0.05. I tipi di parto risultano equamente distribuiti.
Inoltre, da letteratura clinica ed in particolare secondo le indicazioni fornite dall’Ospedale Pediatrico Bambino Gesù, è noto che il peso medio alla nascita si aggira intorno ai 3300 grammi, mentre la lunghezza media è comunemente considerata pari a 500 mm. Per verificare se il nostro campione si discosta da questi valori attesi, sono stati eseguiti due test t per un campione.
mean(Peso)## [1] 3284.081mean(Lunghezza)## [1] 494.692t.test(Peso, mu=3300)##
## One Sample t-test
##
## data: Peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081t.test(Lunghezza, mu=500)##
## One Sample t-test
##
## data: Lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 493.6598 495.7242
## sample estimates:
## mean of x
## 494.692L’ipotesi nulla (H₀) assume che non vi sia differenza tra la media campionaria e quella teorica. Per il peso neonatale, non si rifiuta H₀ (p-value = 0.13), quindi la media del campione non differisce significativamente da quella attesa (3300 g). Al contrario, per la lunghezza, si rifiuta H₀ (p-value < 2.2e-16), indicando una differenza statisticamente significativa rispetto al valore ipotizzato (500 mm). Gli intervalli di confidenza al 95% confermano queste conclusioni: il valore teorico di 500 non è incluso per la lunghezza, mentre 3300 sì nel caso del peso.
CREAZIONE DEL MODELLO
Si è scelto di procedere inizialmente con una selezione manuale delle variabili, seguita successivamente da una selezione automatica tramite procedura stepwise, per poter avere un confronto sulla bontà delle variabili selezionate.
Nella fase manuale, sono state incluse tutte le variabili disponibili, al fine di valutarne l’effettivo contributo esplicativo. Alcune di queste, come tipo tipo di ospedale, pur potendo essere escluse a priori per motivi logici (in quanto non legate direttamente alle caratteristiche biologiche del neonato o della gravidanza), sono state comunque mantenute in un primo momento per completezza, e rimosse solo in seguito se non significative o non utili in ottica previsionale.
#selezione manuale del modello
mod1<-lm(Peso~., data=dati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.80 -181.38 -14.49 160.84 2612.20
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6736.8188 141.4196 -47.637 < 2e-16 ***
## Anni.madre 0.8400 1.1393 0.737 0.4610
## N.gravidanze 11.3505 4.6619 2.435 0.0150 *
## Fumatrici -30.1627 27.5391 -1.095 0.2735
## Gestazione 32.5501 3.8189 8.523 < 2e-16 ***
## Lunghezza 10.2945 0.3007 34.236 < 2e-16 ***
## Cranio 10.4722 0.4260 24.582 < 2e-16 ***
## Tipo.partoNat 29.5134 12.0846 2.442 0.0147 *
## Ospedaleosp2 -11.2019 13.4381 -0.834 0.4046
## Ospedaleosp3 28.0836 13.4970 2.081 0.0376 *
## SessoM 77.5246 11.1781 6.935 5.14e-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-16Dal summary del modello iniziale completo (mod1), si osserva che le variabili Anni.madre, Fumatrici e Ospedaleosp2 hanno p-value rispettivamente pari a 0.43, 0.27 e 0.40: valori molto superiori alla soglia di significatività (0.05), indicando che non forniscono un contributo statisticamente rilevante alla spiegazione del peso neonatale e per questo vengono tolte.
# Rimozione delle variabili non significative: età madre, fumo materno e ospedale
mod2<-update(mod1, ~.-Anni.madre - Fumatrici -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-16La variabile numero di gravidanze, pur risultando significativa nel modello precedente (p = 0.003), presenta un effetto molto contenuto (+12.8 g per ogni gravidanza in più) e non sembra contribuire in modo sostanziale alla qualità predittiva del modello.
#rimozione di N.gravidanze
mod3<-update(mod2, ~.-N.gravidanze)
summary(mod3)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1118.43 -185.56 -16.07 161.53 2626.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6675.8084 135.7769 -49.167 < 2e-16 ***
## Gestazione 31.1917 3.7821 8.247 2.59e-16 ***
## Lunghezza 10.2412 0.3008 34.049 < 2e-16 ***
## Cranio 10.6397 0.4242 25.080 < 2e-16 ***
## Tipo.partoNat 29.1062 12.1113 2.403 0.0163 *
## SessoM 79.0670 11.2010 7.059 2.17e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2494 degrees of freedom
## Multiple R-squared: 0.7267, Adjusted R-squared: 0.7262
## F-statistic: 1326 on 5 and 2494 DF, p-value: < 2.2e-16La sua rimozione nel modello mod3 comporta una riduzione minima dell’R² aggiustato e un lieve miglioramento del BIC (da 35221.75 a 35222.59), suggerendo che il modello diventa più parsimonioso senza perdere informazione rilevante.
#Rimozione tipo di parto
mod4<-update(mod3, ~.-Tipo.parto)
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-16#Verifica BIC per i modelli costruiti manualmente
knitr::kable(BIC(mod1,mod2,mod3,mod4))| df | BIC | |
|---|---|---|
| mod1 | 12 | 35241.91 |
| mod2 | 8 | 35221.75 |
| mod3 | 7 | 35222.59 |
| mod4 | 6 | 35220.54 |
Analogamente alla variabile “N.gravidanze”, la variabile tipo di parto presenta un effetto significativo ma contenuto (+29 g per i parti naturali) e la sua esclusione non compromette la bontà del modello: nel modello 4 il BIC migliora ulteriormente (da 35222.59 a 35220.54), mentre l’R² rimane pressoché invariato (R² aggiustato modello 4 = 0.7257).
Il modello 4 spiega circa il 72% della variabilità del peso neonatale.
Dopo aver verificato la presenza di multicollinearità si procede con l’analisi della varianza:
#Verifica multicollinearità
car::vif(mod4)## Gestazione Lunghezza Cranio Sesso
## 1.653502 2.069517 1.606131 1.038813anova(mod4, mod2)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2495 188688687
## 2 2493 187601677 2 1087011 7.2225 0.0007453 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La verifica della multicollinearità tramite il VIF (Variance Inflation Factor) sul modello mod4 mostra valori tutti inferiori alla soglia critica di 5, indicando che non vi è collinearità preoccupante tra le variabili predittive incluse.
Il confronto tramite ANOVA tra mod2 e mod4 permette di testare se le variabili escluse (numero di gravidanze e tipo di parto) apportano un miglioramento significativo al modello.
L’analisi ANOVA restituisce un p-value = 0.00075, indicando che le due variabili rimosse apportano effettivamente un contributo statisticamente significativo alla spiegazione della variabilità del peso neonatale. Tuttavia, come osservato in precedenza, sia l’effetto stimato che l’impatto informativo (in termini di R² aggiustato e BIC) risultano contenuti. In ottica predittiva, si può quindi preferire un modello più parsimonioso, ovvero il modello 4.
A scopo comparativo con la selezione manuale delle variabili, è stata applicata una procedura automatica di selezione stepwise basata sul criterio BIC.
#SCELTA DEL MODELLO CON PROCEDURA STEPWISE
stepwise.mod<-MASS::stepAIC(mod1,direction = "both",
k=log(n))## Start: AIC=28139.4
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 40789 186809099 28132
## - Fumatrici 1 90016 186858326 28133
## - Ospedale 2 685217 187453527 28133
## - N.gravidanze 1 444811 187213121 28138
## - Tipo.parto 1 447562 187215872 28138
## <none> 186768310 28139
## - Sesso 1 3609311 190377621 28179
## - Gestazione 1 5451454 192219764 28204
## - Cranio 1 45344286 232112596 28675
## - Lunghezza 1 87951071 274719381 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 40789 186768310 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 41671 186858326 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 49391 187552286 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 50245 188015301 28125
## - Sesso 1 3655292 191720838 28158
## - Gestazione 1 5464853 193530399 28181
## - Cranio 1 46108583 234174130 28658
## - Lunghezza 1 87632762 275698308 29066summary(stepwise.mod)##
## 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-16La procedura ha selezionato progressivamente un sottoinsieme di variabili che include: N.gravidanze, Gestazione, Lunghezza, Cranio e Sesso. Il modello finale ottenuto tramite stepwise ha un R² aggiustato di 0.7265: migliore rispetto al modello 4 nell’ordine dei centesimi (ricordo R² aggiustato mod4 = 0.7257) .
Dunque anche il modello ottenuto con il metodo stepwise spiega circa il 72% della varianza del peso neonatale.
#Confronto con i modelli selezionati manualmente
knitr::kable(BIC(stepwise.mod, mod4, mod2))| df | BIC | |
|---|---|---|
| stepwise.mod | 7 | 35220.10 |
| mod4 | 6 | 35220.54 |
| mod2 | 8 | 35221.75 |
#anova(mod4, stepwise.mod)
anova(stepwise.mod, mod4)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2494 188065546
## 2 2495 188688687 -1 -623141 8.2637 0.004079 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Il confronto tra i modelli conferma che il modello selezionato tramite stepwise è sostanzialmente equivalente al modello mod4, ottenuto manualmente, con valori di BIC praticamente sovrapponibili (35220.10 contro 35220.54).
Si preferisce il modello mod4 per la sua maggiore semplicità a parità di capacità esplicativa.
Dal grafico di correlazione, alcune relazioni tra le variabili e il peso neonatale mostrano andamenti curvilinei (ad es. Gestazione, Lunghezza e Cranio). Questo comportamento è un indicatore preliminare della possibile presenza di effetti non lineari, che verranno studiati di seguito.
#Introduzione di EFFETTI NON LINEARI
#Aggiunta di termini quadratici per esplorare effetti non lineari
mod_nl <- lm(Peso ~ Gestazione + I(Gestazione^2) +
Lunghezza + I(Lunghezza^2) + Cranio +I(Cranio^2) + Sesso,
data = dati)
summary(mod_nl)##
## Call:
## lm(formula = Peso ~ Gestazione + I(Gestazione^2) + Lunghezza +
## I(Lunghezza^2) + Cranio + I(Cranio^2) + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1173.39 -183.29 -12.16 163.59 1473.80
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.214e+03 1.163e+03 -1.044 0.296
## Gestazione 3.709e+02 6.751e+01 5.494 4.33e-08 ***
## I(Gestazione^2) -4.332e+00 8.852e-01 -4.894 1.05e-06 ***
## Lunghezza -2.865e+01 4.442e+00 -6.448 1.35e-10 ***
## I(Lunghezza^2) 4.008e-02 4.550e-03 8.810 < 2e-16 ***
## Cranio -5.181e+00 9.825e+00 -0.527 0.598
## I(Cranio^2) 2.328e-02 1.448e-02 1.608 0.108
## SessoM 7.366e+01 1.101e+01 6.690 2.75e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269 on 2492 degrees of freedom
## Multiple R-squared: 0.7383, Adjusted R-squared: 0.7376
## F-statistic: 1004 on 7 and 2492 DF, p-value: < 2.2e-16Si toglie il termine quadratico Cranio che risulta non significativo.
#var I(Cranio^2) non sign
mod_nl2<-update(mod_nl, ~.-I(Cranio^2))
summary(mod_nl2)##
## Call:
## lm(formula = Peso ~ Gestazione + I(Gestazione^2) + Lunghezza +
## I(Lunghezza^2) + Cranio + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1177.46 -182.44 -9.66 164.95 1406.21
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.384e+03 9.072e+02 -2.628 0.00865 **
## Gestazione 3.310e+02 6.281e+01 5.270 1.48e-07 ***
## I(Gestazione^2) -3.819e+00 8.261e-01 -4.624 3.97e-06 ***
## Lunghezza -3.161e+01 4.043e+00 -7.820 7.76e-15 ***
## I(Lunghezza^2) 4.310e-02 4.145e-03 10.398 < 2e-16 ***
## Cranio 1.060e+01 4.177e-01 25.376 < 2e-16 ***
## SessoM 7.397e+01 1.101e+01 6.716 2.30e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.1 on 2493 degrees of freedom
## Multiple R-squared: 0.738, Adjusted R-squared: 0.7374
## F-statistic: 1170 on 6 and 2493 DF, p-value: < 2.2e-16Il modello mod_nl2 include effetti quadratici per le variabili Gestazione e Lunghezza, con un R² aggiustato pari a 0.738, leggermente superiore a quello del modello mod4 (che si ricorda ancora essere pari a 0.7257). Questo significa che mod_nl2 spiega circa un punto percentuale in più della variabilità del peso neonatale ma risultando molto più complesso.
knitr::kable(BIC(mod_nl2, mod4))| df | BIC | |
|---|---|---|
| mod_nl2 | 8 | 35124.96 |
| mod4 | 6 | 35220.54 |
Il confronto del BIC mostra che mod_nl2 è preferibile secondo questo criterio (35121.14 vs 35220.54). Tuttavia, si osserva nuovamente che il guadagno di varianza spiegata è minimo. Per questo motivo, si preferisce ancora mantenere come base il modello mod4, più semplice e interpretabile.
Si passa ad osservare la bontà del modello inserendo anche le interazioni tra le variabili:
#Introduzione di INTERAZIONI TRA VARIABILI
mod_inter <- lm(Peso ~ N.gravidanze + Gestazione + I(Gestazione^2) +
Lunghezza + I(Lunghezza^2) + Cranio + Sesso +
Gestazione*Lunghezza + Gestazione*Cranio +
Lunghezza*Cranio + Lunghezza*Sesso,
data = dati)
summary(mod_inter)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + I(Gestazione^2) +
## Lunghezza + I(Lunghezza^2) + Cranio + Sesso + Gestazione *
## Lunghezza + Gestazione * Cranio + Lunghezza * Cranio + Lunghezza *
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1191.72 -181.65 -10.29 166.82 1313.92
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.315e+03 1.105e+03 -1.190 0.234086
## N.gravidanze 1.491e+01 4.222e+00 3.532 0.000419 ***
## Gestazione 1.221e+02 7.337e+01 1.664 0.096201 .
## I(Gestazione^2) -3.677e+00 1.585e+00 -2.320 0.020441 *
## Lunghezza -1.067e+01 5.550e+00 -1.923 0.054634 .
## I(Lunghezza^2) 5.904e-02 6.554e-03 9.008 < 2e-16 ***
## Cranio -2.505e+00 6.739e+00 -0.372 0.710114
## SessoM 1.410e+02 2.133e+02 0.661 0.508823
## Gestazione:Lunghezza -3.817e-01 1.938e-01 -1.969 0.049058 *
## Gestazione:Cranio 1.140e+00 2.064e-01 5.520 3.73e-08 ***
## Lunghezza:Cranio -6.362e-02 1.471e-02 -4.325 1.59e-05 ***
## Lunghezza:SessoM -1.366e-01 4.303e-01 -0.317 0.750922
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.8 on 2488 degrees of freedom
## Multiple R-squared: 0.743, Adjusted R-squared: 0.7418
## F-statistic: 653.8 on 11 and 2488 DF, p-value: < 2.2e-16A partire dal modello completo con interazioni (mod_inter), si è costruito un modello semplificato (mod_inter_simpl) mantenendo solo i termini quadratici e le interazioni risultate significative. Il nuovo modello include: N.gravidanze, I(Gestazione^2), I(Lunghezza^2), Gestazione:Cranio, Gestazione:Lunghezza e Cranio:Lunghezza.
mod_inter_simpl <- lm(Peso ~ N.gravidanze +
I(Gestazione^2) +
I(Lunghezza^2) +
Gestazione:Cranio +
Gestazione:Lunghezza +
Lunghezza:Cranio,
data = dati)
summary(mod_inter_simpl)##
## Call:
## lm(formula = Peso ~ N.gravidanze + I(Gestazione^2) + I(Lunghezza^2) +
## Gestazione:Cranio + Gestazione:Lunghezza + Lunghezza:Cranio,
## data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1157.2 -176.6 -8.2 165.2 1282.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.017e+03 7.015e+01 -28.755 < 2e-16 ***
## N.gravidanze 1.572e+01 4.256e+00 3.693 0.000226 ***
## I(Gestazione^2) -2.572e+00 1.412e+00 -1.821 0.068718 .
## I(Lunghezza^2) 5.594e-02 6.450e-03 8.674 < 2e-16 ***
## Gestazione:Cranio 1.303e+00 1.487e-01 8.763 < 2e-16 ***
## Gestazione:Lunghezza -4.195e-01 1.650e-01 -2.542 0.011087 *
## Cranio:Lunghezza -8.137e-02 1.174e-02 -6.931 5.3e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.2 on 2493 degrees of freedom
## Multiple R-squared: 0.7377, Adjusted R-squared: 0.737
## F-statistic: 1168 on 6 and 2493 DF, p-value: < 2.2e-16L’R² aggiustato (0.737) resta praticamente invariato rispetto al modello completo (0.7418), mentre il numero di parametri stimati si riduce da 11 a 6. Questo suggerisce che il modello semplificato mantiene una buona capacità esplicativa e predittiva, risultando preferibile in ottica di parsimonia.
Per testare un ulteriore modello si prova ad aggiungere le interazioni al mod4:
mod4_inter <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso +
Gestazione:Lunghezza + Gestazione:Cranio +
Lunghezza:Cranio + Lunghezza:Sesso)
summary(mod4_inter)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso +
## Gestazione:Lunghezza + Gestazione:Cranio + Lunghezza:Cranio +
## Lunghezza:Sesso)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1129.26 -182.51 -12.59 163.74 2575.96
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.098e+02 1.118e+03 -0.367 0.71399
## Gestazione -1.747e+02 6.318e+01 -2.764 0.00575 **
## Lunghezza 1.174e+01 4.990e+00 2.352 0.01873 *
## Cranio -6.022e+00 6.890e+00 -0.874 0.38217
## SessoM -2.042e+02 2.151e+02 -0.949 0.34258
## Gestazione:Lunghezza 4.660e-02 1.047e-01 0.445 0.65623
## Gestazione:Cranio 5.631e-01 2.008e-01 2.804 0.00509 **
## Lunghezza:Cranio -9.895e-03 1.379e-02 -0.717 0.47324
## Lunghezza:SessoM 5.603e-01 4.339e-01 1.291 0.19675
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.3 on 2491 degrees of freedom
## Multiple R-squared: 0.73, Adjusted R-squared: 0.7291
## F-statistic: 841.9 on 8 and 2491 DF, p-value: < 2.2e-16Ne risulta che anche il modello mod4_inter, che introduce solo interazioni tra le variabili principali, incrementa l’R² aggiustato in modo minimo (0.7291).
Si esegue nuovamente l’iter per l’analisi sulla bontà del modello:
knitr::kable(BIC(mod4, mod_inter, mod4_inter))| df | BIC | |
|---|---|---|
| mod4 | 6 | 35220.54 |
| mod_inter | 13 | 35116.50 |
| mod4_inter | 10 | 35215.89 |
Dal confronto appena effettuato non risulta che il modello 4 abbia il valore di BIC, ma come ripetutamente sostenuto in precendenza la differenza con gli altri modelli non è tale da giustificare la scelta di un modello piu complesso.
anova(mod4,mod4_inter)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ Gestazione + Lunghezza + Cranio + Sesso + Gestazione:Lunghezza +
## Gestazione:Cranio + Lunghezza:Cranio + Lunghezza:Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2495 188688687
## 2 2491 185994918 4 2693769 9.0193 3.132e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod4, mod_inter_simpl)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + I(Gestazione^2) + I(Lunghezza^2) + Gestazione:Cranio +
## Gestazione:Lunghezza + Lunghezza:Cranio
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2495 188688687
## 2 2493 180712610 2 7976077 55.017 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod_nl, mod_nl2)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + I(Gestazione^2) + Lunghezza + I(Lunghezza^2) +
## Cranio + I(Cranio^2) + Sesso
## Model 2: Peso ~ Gestazione + I(Gestazione^2) + Lunghezza + I(Lunghezza^2) +
## Cranio + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 180290170
## 2 2493 180477124 -1 -186953 2.5841 0.1081Per proseguire con il confronto sono dunque stati eseguiti test ANOVA tra modelli nidificati.
Il test ANOVA risulta significativo nel confronto tra modello semplice (mod4) e modello con interazioni: il Residual Sum of Squares (RSS) passa da 188.688.687 a 185.994.918, e il test F restituisce un valore pari a 9.02 (p-value = 3.13e-07), indicando un miglioramento statisticamente rilevante.
Lo stesso ragionamento si applica al confronto tra mod4 e mod_inter_simpl, dove l’introduzione di interazioni comporta una riduzione dell’RSS significativa. Al contrario, l’output di anova(mod_nl, mod_nl2) mostra che il termine quadratico I(Cranio²) non fornisce un contributo informativo rilevante e può essere rimosso.
Tuttavia, bilanciando l’aumento di variabilità spiegata con l’aumento della complessità, a conferma di quanto sostenuto finora si sceglie come modello finale il modello 4, in quanto più parsimonioso ma comunque efficace dal punto di vista predittivo.
INTERPRETAZIONE DEI COEFFICIENTI DEL MODELLO SELEZIONATO
Concludendo, il modelo selezionato stima il peso neonatale in funzione della durata della gestazione, della lunghezza e del diametro cranico del neonato, e del sesso.
A parità di altre condizioni, ogni settimana aggiuntiva di gestazione è associata in media a un aumento di 31.3 grammi nel peso. Ogni millimetro in più di lunghezza corrisponde a +10.2 grammi, mentre ogni millimetro in più del diametro cranico si associa ad un aumento di 10.7 grammi nel peso del neoato.
A parità di altre variabili, i maschi pesano in media 79 grammi in più rispetto alle femmine.
ANALISI DEI RESIDUI
#test grafici
#x11()
par(mfrow=c(2,2))
plot(mod4)
Per valutare l’adeguatezza del modello, sono stati analizzati i residui attraverso test statistici e diagnostica grafica.
Di seguito il commento della diagnostica grafica.
QQ PLOT: i residui seguono approssimativamente la distribuzione normale, con lievi deviazioni nelle code.
RESIDUALS VS FITTED: non si osservano pattern sistematici, quindi si suppone che l’assunzione di linearità sia ragionevolmente rispettata.
SCALE-LOCATION: la varianza dei residui appare costante lungo i fitted values, non si notano problemi evidenti di eteroschedasticità.
RESIDUALS VD LEVERAGE: sono presenti pochi punti ad alta leva (es. osservazione 1551 che si colloca tra le soglie di Cook di 0.5 e 1 e che potrebbe avere un’influenza moderata sul mmodello), ma nessuno supera la soglia critica di Cook, quindi non emergono influenze eccessive da singole osservazioni.
Nel complesso, la diagnostica grafica non evidenzia violazioni gravi delle ipotesi del modello lineare. Il modello risulta quindi adeguato sotto questo profilo.
#test statistici
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 0#Root Mean Squared Error (RMSE).
sqrt(mean(residuals(mod4)^2)) #[1] 266.1389## [1] 274.728#leverage
lev<-hatvalues(mod4)
plot(lev)
#valore soglia:
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#outlier
plot(rstudent(mod4))
abline(h=c(-2,2),col=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-03Per quanto concerne i test statistici invece di seguito il commento:
Il test di Shapiro-Wilk evidenzia una deviazione dalla normalità dei residui (p < 2.2e-16), ma data la numerosità campionaria elevata, tale violazione non compromette la validità delle stime.
Il test di Breusch-Pagan indica una possibile eteroschedasticità (p < 2.2e-16), confermata da una lieve variazione della varianza nei grafici. Tuttavia, l’impatto pratico è contenuto e non influisce in modo sostanziale sull’interpretazione del modello.
Il test di Durbin-Watson (p = 0.13) non segnala autocorrelazione significativa dei residui.
Infine, l’analisi dei valori di leverage e degli outlier studentizzati mostra pochi punti potenzialmente influenti.
Per valutare l’influenza delle singole osservazioni sul modello finale, è stato esaminato il grafico della distanza di Cook.
cook<-cooks.distance(mod4)
plot(cook)
max(cook)## [1] 0.9783883Questo indice permette di identificare punti che, se rimossi, causerebbero un cambiamento sostanziale nei coefficienti stimati. Nel grafico si osserva un’unica osservazione con valore di Cook elevato (0.978), ma comunque inferiore al valore soglia comunemente utilizzato (1). Questo suggerisce che, sebbene il punto sia influente, non lo è in maniera tale da invalidare il modello.
L’osservazione specifica è la sequnte:
which.max(cook)## 1551
## 1551knitr::kable(dati[which.max(cook), ])| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1551 | 35 | 1 | 0 | 38 | 4370 | 315 | 374 | Nat | osp3 | F |
L’osservazione con il valore di Cook elevato è la 1551: come già specificato risulta sì influente ma non anomala.
PREVISIONE E RISULTATI
Per testare il modello selezionato su un caso pratico, è stato creato un profilo tipo di neonata: madre alla terza gravidanza, parto alla 39ª settimana, con valori medi di lunghezza e cranio (rispettivamente di 494.692 e 340.0292 arrotondati a 495 e 340).
mean(Lunghezza) #494.692## [1] 494.692mean(Cranio) #340.0292## [1] 340.0292#Profilo tipo: neonata femmina nata alla 39a settimana, terza gravidanza
nuova_neonata <- data.frame(
Gestazione = 39,
Lunghezza = 495,
Cranio = 340,
Sesso = factor("F", levels = levels(dati$Sesso))
)
predict(mod4, newdata = nuova_neonata)## 1
## 3248.163# 3248.163 La variabile “numero di gravidanze” era stata esclusa dal modello finale per ragioni di rilevanza, dunque la sua esclusione non compromette la validità della stima.
Il peso previsto per la neonata risulta pari a 3248.163 grammi.
VISUALIZZAZIONI GRAFICHE
Per interpretare visivamente il comportamento del modello, vengono riportate di seguito delle visualizzazioni che mostrano l’effetto della durata della gestazione sul peso previsto.
È stato costruito un dataset di riferimento, fissando le altre variabili (lunghezza e cranio) su valori medi osservati nel campione e simulando un profilo tipo di neonata femmina. In questo modo, si osserva come cambia il peso previsto al variare delle settimane di gestazione, mantenendo costanti gli altri fattori.
#VISUALIZZAZIONI
#effetto della durata della gestazione sul peso
settimane <- seq(min(dati$Gestazione), max(dati$Gestazione), 1)
dati_grafico <- data.frame(
Gestazione = settimane,
Lunghezza = 495, # valore medio approssimato
Cranio = 340, # valore medio
Sesso = factor("F", levels = levels(dati$Sesso))
)
# Aggiunta della previsione
dati_grafico$Peso_previsto <- predict(mod4, newdata = dati_grafico)# Grafico
# Plot base: effetto della gestazione sul peso atteso per neonata tipo (femmina)
plot(dati_grafico$Gestazione, dati_grafico$Peso_previsto,
type = "b", pch = 19,
xlab = "Settimane di gestazione",
ylab = "Peso previsto (g)",
main = "Peso previsto in base alla gestazione (neonata tipo)")
Dal grafico si osserva una relazione crescente e lineare tra le settimane di gestazione e il peso previsto alla nascita. L’andamento è coerente con la matrice di correlazione, in cui l’indice di correlazione tra peso e gestazione è pari a 0.59.
Nel modello, infatti, la variabile Gestazione risulta significativa, con un coefficiente positivo che indica un aumento medio di peso di circa 31 grammi per ogni settimana aggiuntiva.
Nel secondo grafico viene mostrato l’effetto congiunto del sesso e della durata della gestazione sul peso previsto alla nascita, sempre mantenendo costanti le altre variabili (lunghezza e cranio su valori medi). Le due curve presentano un andamento parallelo e crescente, confermando il contributo positivo della gestazione.
library(ggplot2)
# Dataset per entrambi i sessi (con variabili compatibili col mod4)
dati_plot <- do.call(rbind, lapply(c("M", "F"), function(sesso) {
data.frame(
Gestazione = settimane,
Lunghezza = 495,
Cranio = 340,
Sesso = factor(sesso, levels = levels(dati$Sesso)),
SessoLabel = ifelse(sesso == "M", "Maschio", "Femmina")
)
}))
# Previsione peso dal modello 4
dati_plot$Peso_previsto <- predict(mod4, newdata = dati_plot)
# Grafico con ggplot
ggplot(dati_plot, aes(x = Gestazione, y = Peso_previsto, color = SessoLabel)) +
geom_line(size = 1.2) +
geom_point(size = 2) +
scale_color_manual(values = c("Maschio" = "tan2", "Femmina" = "olivedrab4")) +
labs(title = "Peso previsto in base alla gestazione per sesso",
x = "Settimane di gestazione",
y = "Peso previsto (g)",
color = "Sesso") +
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.
La curva relativa ai neonati maschi si colloca sistematicamente al di sopra di quella delle femmine, evidenziando come il sesso maschile incida mediamente per circa 79 grammi in più rispetto al sesso femminile, in linea con quanto indicato dal coefficiente associato alla variabile Sesso (SessoM).
CONCLUSIONI
Il modello finale risulta essere parsimoniamente efficace: spiega una quota consistente della variabilità del peso neonatale, poco più del 72%, e mantiene una struttura semplice con sole variabili principali, tutte significative. Le diagnosi grafiche e statistiche dei residui mostrano una buona adeguatezza del modello.
L’introduzione di interazioni o termini non lineari porta a miglioramenti marginali ed il guadagno in capacità predittiva non ne giustifica la maggiore complessità.