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] 2500shapiro.test(Peso)##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16ANALISI PRELIMINARE
Dopo aver verificato la struttura del dataset e il numero totale di osservazioni (n = 2500), è stato eseguito il test di Shapiro-Wilk sulla variabile risposta Peso per valutarne la normalità. Il test restituisce un p-value molto basso pari a 2.2e-16, permettendo di rifiutare l’ipotesi nulla: si conclude quindi che il peso alla nascita non segue una distribuzione normale nel campione analizzato.
plot(density(Peso))
qqnorm(Peso)
qqline(Peso, col="red")
I due grafici, il Q-Q plot e la curva di densità, confermano che la variabile Peso non segue una distribuzione normale. Come si vede nel Q-Q plot i punti si discostano visibilmente dalla linea teorica nelle code. La curva di densità mostra una distribuzione leggermente asimmetrica a destra.
Data la dimensione del campione, la regressione lineare è comunque applicabile, sarà necessario osservare la distribuzione dei residui.
#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)
Dal grafico di correlazione emerge che le variabili più strettamente associate al peso neonatale sono la lunghezza (r = 0.80), il diametro cranico (r = 0.70) e le settimane di gestazione (r = 0.59). Questi risultati indicano una relazione lineare positiva tra peso e caratteristiche antropometriche del neonato. In particolare, la forte correlazione tra lunghezza e cranio (r = 0.60) suggerisce che entrambe potrebbero spiegare parte della variabilità del peso. 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.
Per analizzare come il peso alla nascita varia in relazione a diverse variabili categoriali, sono stati utilizzati boxplot e test t per campioni indipendenti, prima condizionando il peso al tipo di parto, al sesso e all’ospedale, e successivamente all’avere una madre fumatrice.
par(mfrow=c(2,2))
boxplot(Peso)
boxplot(Peso~Tipo.parto)
boxplot(Peso~Sesso)
boxplot(Peso~Ospedale)
median(Peso[Sesso=="M"])## [1] 3430median(Peso[Sesso=="F"])## [1] 3160I boxplot mostrano che, condizionando il peso al sesso, al tipo di parto e all’ospedale, si evidenziano alcune differenze interessanti. I 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)
Per valutare se il peso medio alla nascita differisca tra gruppi, sono stati effettuati test t a due campioni indipendenti. Per quanto riguarda il sesso il t-test mostra una differenza statisticamente significativa tra maschi e femmine (p < 2.2e-16). I neonati maschi hanno un peso medio maggiore rispetto alle femmine. Mentre per il tipo di parto il test non evidenzia una differenza significativa tra neonati nati da parto cesareo e naturale (p ≈ 0.90).
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 è 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, non si rifiuta H₀ (p value = 0.13), mentre per la lunghezza si rifiuta H₀ (p < 2.2e-16), indicando una differenza significativa. Dunque la media del peso del campione non differisce significativamente da quella di popolazione, mentre la lunghezza media del neonato è significativamente diversa. Gli intervalli di confidenza supportano queste conclusioni.
CREAZIONE DEL MODELLO
#selezione manuale del modello
mod1<-lm(Peso~., data=dati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1124.40 -181.66 -14.42 160.91 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6738.4762 141.3087 -47.686 < 2e-16 ***
## Anni.madre 0.8921 1.1323 0.788 0.4308
## N.gravidanze 11.2665 4.6608 2.417 0.0157 *
## Fumatrici -30.1631 27.5386 -1.095 0.2735
## Gestazione 32.5696 3.8187 8.529 < 2e-16 ***
## Lunghezza 10.2945 0.3007 34.236 < 2e-16 ***
## Cranio 10.4707 0.4260 24.578 < 2e-16 ***
## Tipo.partoNat 29.5254 12.0844 2.443 0.0146 *
## Ospedaleosp2 -11.2095 13.4379 -0.834 0.4043
## Ospedaleosp3 28.0958 13.4957 2.082 0.0375 *
## SessoM 77.5409 11.1776 6.937 5.08e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 669.2 on 10 and 2489 DF, p-value: < 2.2e-16# Rimozione delle variabili non significative: età madre e fumo materno
mod2<-update(mod1, ~.-Anni.madre - Fumatrici )
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1113.18 -181.16 -16.58 161.01 2620.19
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.4293 135.9438 -49.340 < 2e-16 ***
## N.gravidanze 12.3619 4.3325 2.853 0.00436 **
## Gestazione 31.9909 3.7896 8.442 < 2e-16 ***
## Lunghezza 10.3086 0.3004 34.316 < 2e-16 ***
## Cranio 10.4922 0.4254 24.661 < 2e-16 ***
## Tipo.partoNat 29.2803 12.0817 2.424 0.01544 *
## Ospedaleosp2 -11.0227 13.4363 -0.820 0.41209
## Ospedaleosp3 28.6408 13.4886 2.123 0.03382 *
## SessoM 77.4412 11.1756 6.930 5.36e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared: 0.7287, Adjusted R-squared: 0.7278
## F-statistic: 836.3 on 8 and 2491 DF, p-value: < 2.2e-16# Esclusione della variabile Ospedale
mod3<-update(mod2, ~.-Ospedale)
summary(mod3)##
## 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-16#rimozione di N.gravidanze
mod4<-update(mod3, ~.-N.gravidanze)
summary(mod4)##
## 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-16#Rimozione anche del tipo di parto
mod5<-update(mod4, ~.-Tipo.parto)
summary(mod5)##
## 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#R-squared: 0.7257 con modello 5 si perde variabilita spiegata
#Verifica BIC per i modelli costruiti manualmente
BIC(mod1,mod2,mod3,mod4,mod5) ## df BIC
## mod1 12 35241.84
## mod2 10 35228.03
## mod3 8 35221.75
## mod4 7 35222.59
## mod5 6 35220.54#Verifica multicollinearità
car::vif(mod5)## Gestazione Lunghezza Cranio Sesso
## 1.653502 2.069517 1.606131 1.038813anova(mod5, mod3)## 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 ' ' 1#SCELTA DEL MODELLO CON PROCEDURA STEPWISE
stepwise.mod<-MASS::stepAIC(mod1,direction = "both",
k=log(n))## Start: AIC=28139.32
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 46578 186809099 28132
## - Fumatrici 1 90019 186852540 28133
## - Ospedale 2 685979 187448501 28133
## - N.gravidanze 1 438452 187200974 28137
## - Tipo.parto 1 447929 187210450 28138
## <none> 186762521 28139
## - Sesso 1 3611021 190373542 28179
## - Gestazione 1 5458403 192220925 28204
## - Cranio 1 45326172 232088693 28675
## - Lunghezza 1 87951062 274713583 29096
##
## Step: AIC=28132.12
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 90897 186899996 28126
## - Ospedale 2 692738 187501837 28126
## - Tipo.parto 1 448222 187257321 28130
## <none> 186809099 28132
## - N.gravidanze 1 633756 187442855 28133
## + Anni.madre 1 46578 186762521 28139
## - Sesso 1 3618736 190427835 28172
## - Gestazione 1 5412879 192221978 28196
## - Cranio 1 45588236 232397335 28670
## - Lunghezza 1 87950050 274759149 29089
##
## Step: AIC=28125.51
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 701680 187601677 28119
## - Tipo.parto 1 440684 187340680 28124
## <none> 186899996 28126
## - N.gravidanze 1 610840 187510837 28126
## + Fumatrici 1 90897 186809099 28132
## + Anni.madre 1 47456 186852540 28133
## - Sesso 1 3602797 190502794 28165
## - Gestazione 1 5346781 192246777 28188
## - Cranio 1 45632149 232532146 28664
## - Lunghezza 1 88355030 275255027 29086
##
## Step: AIC=28119.23
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 463870 188065546 28118
## <none> 187601677 28119
## - N.gravidanze 1 651066 188252743 28120
## + Ospedale 2 701680 186899996 28126
## + Fumatrici 1 99840 187501837 28126
## + Anni.madre 1 54392 187547285 28126
## - Sesso 1 3649259 191250936 28160
## - Gestazione 1 5444109 193045786 28183
## - Cranio 1 45758101 233359778 28657
## - Lunghezza 1 88054432 275656108 29074
##
## Step: AIC=28117.58
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 188065546 28118
## - N.gravidanze 1 623141 188688687 28118
## + Tipo.parto 1 463870 187601677 28119
## + Ospedale 2 724866 187340680 28124
## + Fumatrici 1 91892 187973654 28124
## + Anni.madre 1 54816 188010731 28125
## - Sesso 1 3655292 191720838 28158
## - Gestazione 1 5464853 193530399 28181
## - Cranio 1 46108583 234174130 28658
## - Lunghezza 1 87632762 275698308 29066summary(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-16#Confronto con i modelli selezionati manualmente
BIC(stepwise.mod, mod5, mod3) ## df BIC
## stepwise.mod 7 35220.10
## mod5 6 35220.54
## mod3 8 35221.75anova(mod5, stepwise.mod)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2495 188688687
## 2 2494 188065546 1 623141 8.2637 0.004079 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(stepwise.mod, mod3)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + 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 2494 188065546
## 2 2493 187601677 1 463870 6.1643 0.0131 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1stepwise.mod[["call"]]## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)#Introduzione di EFFETTI NON LINEARI
mod_nl <- lm(Peso ~ N.gravidanze + Gestazione + I(Gestazione^2) +
Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^2) + Sesso,
data = dati)
summary(mod_nl)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + I(Gestazione^2) +
## Lunghezza + I(Lunghezza^2) + Cranio + I(Cranio^2) + Sesso,
## data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1187.03 -182.50 -12.03 162.60 1469.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.217e+03 1.160e+03 -1.049 0.2945
## N.gravidanze 1.441e+01 4.246e+00 3.394 0.0007 ***
## Gestazione 3.754e+02 6.738e+01 5.572 2.79e-08 ***
## I(Gestazione^2) -4.374e+00 8.834e-01 -4.951 7.86e-07 ***
## Lunghezza -2.925e+01 4.436e+00 -6.592 5.28e-11 ***
## I(Lunghezza^2) 4.075e-02 4.544e-03 8.967 < 2e-16 ***
## Cranio -4.979e+00 9.804e+00 -0.508 0.6116
## I(Cranio^2) 2.276e-02 1.445e-02 1.575 0.1154
## SessoM 7.232e+01 1.100e+01 6.577 5.83e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.4 on 2491 degrees of freedom
## Multiple R-squared: 0.7395, Adjusted R-squared: 0.7387
## F-statistic: 883.9 on 8 and 2491 DF, p-value: < 2.2e-16#togliere var I(Cranio^2) non sign
mod_nl2<-update(mod_nl, ~.-I(Cranio^2))
summary(mod_nl2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + I(Gestazione^2) +
## Lunghezza + I(Lunghezza^2) + Cranio + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1191.07 -182.28 -13.74 163.34 1403.53
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.360e+03 9.053e+02 -2.607 0.00919 **
## N.gravidanze 1.448e+01 4.247e+00 3.410 0.00066 ***
## Gestazione 3.365e+02 6.270e+01 5.367 8.74e-08 ***
## I(Gestazione^2) -3.873e+00 8.245e-01 -4.698 2.77e-06 ***
## Lunghezza -3.215e+01 4.037e+00 -7.963 2.53e-15 ***
## I(Lunghezza^2) 4.370e-02 4.140e-03 10.556 < 2e-16 ***
## Cranio 1.045e+01 4.192e-01 24.922 < 2e-16 ***
## SessoM 7.261e+01 1.100e+01 6.602 4.93e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.5 on 2492 degrees of freedom
## Multiple R-squared: 0.7392, Adjusted R-squared: 0.7385
## F-statistic: 1009 on 7 and 2492 DF, p-value: < 2.2e-16BIC(mod_nl2, stepwise.mod)## df BIC
## mod_nl2 9 35121.14
## stepwise.mod 7 35220.10#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-16mod_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)#Modello finale (semplificato)
mod_finale <- lm(Peso ~ N.gravidanze + Gestazione + I(Gestazione^2) +
Lunghezza + I(Lunghezza^2) +
Gestazione*Lunghezza + Gestazione*Cranio +
Lunghezza*Cranio,
data = dati)
summary(mod_finale)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + I(Gestazione^2) +
## Lunghezza + I(Lunghezza^2) + Gestazione * Lunghezza + Gestazione *
## Cranio + Lunghezza * Cranio, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1161.62 -175.53 -6.97 163.81 1271.03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.082e+02 1.111e+03 -0.727 0.467050
## N.gravidanze 1.598e+01 4.255e+00 3.756 0.000177 ***
## Gestazione 1.065e+02 7.395e+01 1.440 0.150114
## I(Gestazione^2) -3.795e+00 1.598e+00 -2.375 0.017605 *
## Lunghezza -1.138e+01 5.565e+00 -2.046 0.040899 *
## I(Lunghezza^2) 5.796e-02 6.511e-03 8.902 < 2e-16 ***
## Cranio -3.007e+00 6.796e+00 -0.442 0.658175
## Gestazione:Lunghezza -3.356e-01 1.952e-01 -1.720 0.085627 .
## Gestazione:Cranio 1.150e+00 2.082e-01 5.525 3.64e-08 ***
## Lunghezza:Cranio -6.316e-02 1.484e-02 -4.257 2.15e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.1 on 2490 degrees of freedom
## Multiple R-squared: 0.7383, Adjusted R-squared: 0.7374
## F-statistic: 780.5 on 9 and 2490 DF, p-value: < 2.2e-16# Confronto finale dei modelli
BIC(mod5, mod_nl2, mod_inter, mod_finale, stepwise.mod)## df BIC
## mod5 6 35220.54
## mod_nl2 9 35121.14
## mod_inter 13 35116.50
## mod_finale 11 35145.69
## stepwise.mod 7 35220.10Sono stati creati e confrontati manualmente diversi modelli di regressione lineare, rimuovendo progressivamente le variabili meno significative. Successivamente, a scopo comparativo, è stata applicata una procedura automatica di selezione stepwise basata sul criterio BIC, al fine di individuare la combinazione più parsimoniosa ed efficace di predittori: stepwise.mod risulta il miglior compromesso tra semplicità e accuratezza. La formula del modello lineare è: lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso, data = dati).
Sono stati inoltre valutati possibili effetti non lineari e interazioni tra le variabili. Rispetto al modello stepwise.mod il modello che include le relazioni quadratiche per Gestazione e Lunghezza risulta migliore. Inserire relazioni quadratiche aiuta a spiegare circa il 74% della variabilità del peso neonatale.
Per migliorare ulteriormente la capacità predittiva del modello, sono state introdotte interazioni tra variabili continue (es. Gestazione × Cranio, Lunghezza × Cranio), che potrebbero influenzare congiuntamente il peso alla nascita. Il modello con interazioni (mod_inter) mostra i migliori risultati in termini di BIC e R² aggiustato, risultando quindi preferibile.
Variabili significative incluse nel modello predittivo: 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)
Si è testato un modello semplificato (mod_finale) eliminando alcune interazioni non significative, ma ha mostrato un BIC più alto e un R² adj inferiore. Il modello completo mod_inter resta quindi il più performante.
ANALISI DEI RESIDUI
#test grafici
#x11()
par(mfrow=c(2,2))
plot(mod_inter)
plot(mod_nl2)
plot(stepwise.mod)
Per valutare l’adeguatezza del modello mod_inter, sono stati analizzati i residui attraverso test statistici e diagnostica grafica.
Di seguito il commento per la diagnostica grafica.
QQ PLOT: I residui seguono approssimativamente la distribuzione normale, con leggere deviazioni nelle code estreme. RESIDUALS VS FITTED: Non si osservano pattern sistematici evidenti, suggerendo che l’ipotesi di linearità sia soddisfatta. SCALE LOCATION: La varianza dei residui appare costante, non si notano gravi problemi di eteroschedasticità. RESIDUALS VS LEVERAGE: Non ci sono punti che superano la distanza di Cook, quindi non sembra esserci un’influenza eccessiva da parte di alcune osservazioni - dalla seguente analisi statistica dei reisui emergerà un solo valore che supera la distanza di Cook.
#test statistici
shapiro.test(residuals(mod_inter)) ##
## Shapiro-Wilk normality test
##
## data: residuals(mod_inter)
## W = 0.99099, p-value = 2.152e-11lmtest::bptest(mod_inter) ##
## studentized Breusch-Pagan test
##
## data: mod_inter
## BP = 50.374, df = 11, p-value = 5.362e-07lmtest::dwtest(mod_inter) ##
## Durbin-Watson test
##
## data: mod_inter
## DW = 1.9577, p-value = 0.1451
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(residuals(stepwise.mod)) ##
## Shapiro-Wilk normality test
##
## data: residuals(stepwise.mod)
## W = 0.97408, p-value < 2.2e-16lmtest::bptest(stepwise.mod) ##
## studentized Breusch-Pagan test
##
## data: stepwise.mod
## BP = 90.253, df = 5, p-value < 2.2e-16lmtest::dwtest(stepwise.mod) ##
## Durbin-Watson test
##
## data: stepwise.mod
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0#Root Mean Squared Error (RMSE).
sqrt(mean(residuals(mod_inter)^2)) #[1] 266.1389## [1] 266.1389sqrt(mean(residuals(stepwise.mod)^2))## [1] 274.274#leverage
lev<-hatvalues(mod_inter)
plot(lev)
#valore soglia:
p=sum(lev)
soglia=2*p/n
abline(h=soglia, col=2)
lev[lev>soglia]## 1 15 36 61 67 89
## 0.009832526 0.018740781 0.011441416 0.012531723 0.011485033 0.013012843
## 96 101 106 119 131 151
## 0.012459167 0.025721741 0.042893496 0.010466988 0.017680418 0.075347042
## 155 161 190 204 206 220
## 0.024908005 0.026235289 0.012058009 0.015412340 0.022786598 0.018793941
## 249 295 304 305 310 312
## 0.012169674 0.013296590 0.010740034 0.011243306 0.382817149 0.073682418
## 315 317 322 328 378 383
## 0.014024299 0.010467450 0.009975805 0.013146308 0.079545791 0.013212511
## 445 471 492 516 582 587
## 0.013378930 0.013308410 0.026783102 0.014159935 0.013743171 0.030120231
## 592 615 638 656 666 684
## 0.016263113 0.012502926 0.013060803 0.019251458 0.013909174 0.012114652
## 701 702 726 729 748 750
## 0.011399807 0.009834915 0.010478666 0.014772317 0.026429544 0.027281518
## 757 765 805 895 928 946
## 0.009659273 0.012987557 0.043224624 0.013510991 0.138754605 0.011106763
## 947 949 956 958 964 988
## 0.019200282 0.011783068 0.060262187 0.010155875 0.013782652 0.011506174
## 991 1014 1067 1091 1105 1130
## 0.016850569 0.053714270 0.018771084 0.026064439 0.012003529 0.057059118
## 1181 1188 1200 1219 1222 1238
## 0.015490907 0.019766053 0.010773682 0.030972922 0.010238392 0.010677255
## 1248 1267 1273 1275 1283 1311
## 0.054127866 0.010503834 0.031814790 0.014983549 0.015386482 0.012157390
## 1321 1323 1356 1357 1370 1385
## 0.011187340 0.013729809 0.010434167 0.018316635 0.009891310 0.048069352
## 1400 1420 1428 1429 1450 1505
## 0.018599550 0.017785620 0.047865263 0.218347092 0.016855007 0.013557936
## 1513 1551 1553 1556 1560 1573
## 0.014033419 0.673499824 0.040402962 0.012369353 0.017715801 0.009634622
## 1593 1596 1606 1610 1619 1628
## 0.014356481 0.010915453 0.012400598 0.034789975 0.076761302 0.012126708
## 1634 1639 1686 1692 1701 1705
## 0.013086963 0.012431200 0.025300530 0.009939538 0.043186576 0.010533104
## 1712 1718 1727 1735 1743 1780
## 0.012796606 0.012484673 0.014785101 0.011513407 0.010889920 0.140597572
## 1781 1802 1809 1827 1858 1893
## 0.017185468 0.009851903 0.023309353 0.013420943 0.014354647 0.010128721
## 1920 1948 1977 2006 2008 2016
## 0.013214532 0.010039018 0.017492987 0.009888862 0.011336911 0.010472014
## 2037 2040 2086 2089 2111 2112
## 0.016264438 0.038601188 0.013517016 0.014684791 0.010861129 0.012167585
## 2114 2115 2120 2136 2140 2148
## 0.052642783 0.100480378 0.079710706 0.012385821 0.011917261 0.009724121
## 2149 2175 2200 2215 2216 2220
## 0.040944299 0.129779918 0.033334607 0.014386110 0.034652252 0.012703720
## 2221 2224 2225 2245 2257 2307
## 0.021861471 0.012196152 0.013119094 0.012000104 0.015198051 0.049368265
## 2353 2359 2391 2408 2422 2437
## 0.010689093 0.013422871 0.018268240 0.023020955 0.022090482 0.159569293
## 2438 2452 2458 2471 2478
## 0.011413306 0.122981538 0.023046838 0.022061248 0.019650719#outlier
plot(rstudent(mod_inter))
abline(h=c(-2,2),col=2)
car::outlierTest(mod_inter)## rstudent unadjusted p-value Bonferroni p
## 1306 4.954608 7.7345e-07 0.0019336
## 155 4.537838 5.9535e-06 0.0148840
## 1399 -4.494155 7.3044e-06 0.0182610
## 1694 4.393681 1.1612e-05 0.0290300#Gli outlier sono pochi ma significativi
cook<-cooks.distance(mod_inter)
plot(cook)
max(cook)## [1] 3.110081which.max(cook)## 1551
## 1551dati[which.max(cook), ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1551 35 1 0 38 4370 315 374
## Tipo.parto Ospedale Sesso
## 1551 Nat osp3 Fdati_senza_outlier <- dati[-which.max(cook), ] # Esclusione il valore influente
mod_inter_clean <- lm(Peso ~ N.gravidanze + Gestazione + I(Gestazione^2) +
Lunghezza + I(Lunghezza^2) + Cranio + Sesso +
Gestazione*Lunghezza + Gestazione*Cranio +
Lunghezza*Cranio + Lunghezza*Sesso,
data = dati_senza_outlier)
summary(mod_inter_clean)##
## 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_senza_outlier)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1174.62 -181.04 -11.06 165.32 1311.90
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -475.23780 1118.88344 -0.425 0.671061
## N.gravidanze 14.65433 4.20793 3.483 0.000505 ***
## Gestazione 145.69193 73.33007 1.987 0.047055 *
## I(Gestazione^2) -4.52483 1.59227 -2.842 0.004523 **
## Lunghezza -9.18038 5.54229 -1.656 0.097762 .
## I(Lunghezza^2) 0.02544 0.01023 2.487 0.012941 *
## Cranio -12.23884 7.09221 -1.726 0.084530 .
## SessoM 154.22147 212.61613 0.725 0.468305
## Gestazione:Lunghezza -0.01569 0.21134 -0.074 0.940837
## Gestazione:Cranio 0.73014 0.22699 3.217 0.001314 **
## Lunghezza:Cranio -0.01190 0.01902 -0.626 0.531601
## Lunghezza:SessoM -0.16547 0.42890 -0.386 0.699685
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 265.9 on 2487 degrees of freedom
## Multiple R-squared: 0.7444, Adjusted R-squared: 0.7433
## F-statistic: 658.4 on 11 and 2487 DF, p-value: < 2.2e-16Il test di Shapiro-Wilk ha restituito un p-value molto basso, portando al rifiuto dell’ipotesi nulla di normalità dei residui. Il test di Breusch-Pagan ha evidenziato eteroschedasticità significativa, mentre il test di Durbin-Watson non ha mostrato autocorrelazione.
È stata poi calcolata la Root Mean Squared Error (RMSE), pari a circa 266.1, inferiore rispetto a quella del modello stepwise.mod (274.3) - utilizzato nuovamente come confronto per essere sicuri della bontà del modello scelto. E’ confermata una maggiore precisione predittiva di mod_inter.
L’analisi del leverage ha mostrato poche osservazioni con valori elevati, suggerendo una bassa influenza individuale della maggior parte dei punti. Il test sugli outlier ha invece identificato alcune osservazioni potenzialmente problematiche, e l’indice di Cook ha evidenziato un’unica osservazione con influenza marcata.
Dopo aver escluso l’osservazione più influente e ricreato il modello (mod_inter_clean), si è osservato un miglioramento trascurabile. Inoltre, alcune variabili precedentemente significative hanno perso rilevanza, suggerendo che la rimozione non apporta reali benefici. Per questo motivo, l’outlier è stato mantenuto nel modello finale, che si riconferma essere mod_inter.
ANALISI DEI RESIDUI
#ANALISI RESIDUI
library(lmtest)## Warning: il pacchetto 'lmtest' è stato creato con R versione 4.2.3## Caricamento del pacchetto richiesto: zoo## Warning: il pacchetto 'zoo' è stato creato con R versione 4.2.3##
## Caricamento pacchetto: 'zoo'## I seguenti oggetti sono mascherati da 'package:base':
##
## as.Date, as.Date.numericbptest(mod_inter) #eteroschedasticità ##
## studentized Breusch-Pagan test
##
## data: mod_inter
## BP = 50.374, df = 11, p-value = 5.362e-07dwtest(mod_inter) ##
## Durbin-Watson test
##
## data: mod_inter
## DW = 1.9577, p-value = 0.1451
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(residuals(mod_inter)) ##
## Shapiro-Wilk normality test
##
## data: residuals(mod_inter)
## W = 0.99099, p-value = 2.152e-11#rendere il modello robusto per rimediare all'eteroschedasticita
#install.packages("sandwich")
library(sandwich)
library(lmtest)
coeftest(mod_inter, vcov = vcovHC(mod_inter, type = "HC1"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3154e+03 1.1432e+03 -1.1507 0.2499699
## N.gravidanze 1.4914e+01 4.4935e+00 3.3190 0.0009164 ***
## Gestazione 1.2211e+02 8.7728e+01 1.3919 0.1640777
## I(Gestazione^2) -3.6773e+00 2.1416e+00 -1.7171 0.0860926 .
## Lunghezza -1.0672e+01 7.0634e+00 -1.5108 0.1309583
## I(Lunghezza^2) 5.9038e-02 1.2932e-02 4.5652 5.232e-06 ***
## Cranio -2.5051e+00 8.4708e+00 -0.2957 0.7674590
## SessoM 1.4096e+02 2.2400e+02 0.6293 0.5292228
## Gestazione:Lunghezza -3.8167e-01 3.0535e-01 -1.2499 0.2114416
## Gestazione:Cranio 1.1396e+00 2.6673e-01 4.2725 2.005e-05 ***
## Lunghezza:Cranio -6.3619e-02 2.2518e-02 -2.8252 0.0047632 **
## Lunghezza:SessoM -1.3661e-01 4.5232e-01 -0.3020 0.7626591
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(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-16#modello semplificato rispetto a mod_inter
#con quanto di rilevante dal coeftest()
mod_inter_simpl <- lm(Peso ~ N.gravidanze +
I(Gestazione^2) +
I(Lunghezza^2) +
Cranio +
Gestazione:Cranio +
Lunghezza:Cranio,
data = dati)
summary(mod_inter_simpl)##
## Call:
## lm(formula = Peso ~ N.gravidanze + I(Gestazione^2) + I(Lunghezza^2) +
## Cranio + Gestazione:Cranio + Lunghezza:Cranio, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1132.35 -174.55 -8.47 163.28 1293.38
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.547e+03 9.219e+02 -1.678 0.093539 .
## N.gravidanze 1.581e+01 4.262e+00 3.710 0.000212 ***
## I(Gestazione^2) -5.983e+00 7.236e-01 -8.268 < 2e-16 ***
## I(Lunghezza^2) 4.208e-02 4.139e-03 10.167 < 2e-16 ***
## Cranio -2.657e+00 5.529e+00 -0.481 0.630822
## Cranio:Gestazione 1.474e+00 1.654e-01 8.911 < 2e-16 ***
## Cranio:Lunghezza -8.932e-02 1.190e-02 -7.509 8.26e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269.6 on 2493 degrees of freedom
## Multiple R-squared: 0.737, Adjusted R-squared: 0.7364
## F-statistic: 1164 on 6 and 2493 DF, p-value: < 2.2e-16BIC(mod_inter_simpl, mod_inter)## df BIC
## mod_inter_simpl 8 35134.46
## mod_inter 13 35116.50sqrt(mean(residuals(mod_inter_simpl)^2)) ## [1] 269.1945# Analisi dei residui
#x11()
par(mfrow=c(2,2))
plot(mod_inter_simpl)
bptest(mod_inter_simpl)##
## studentized Breusch-Pagan test
##
## data: mod_inter_simpl
## BP = 46.76, df = 6, p-value = 2.089e-08shapiro.test(residuals(mod_inter_simpl))##
## Shapiro-Wilk normality test
##
## data: residuals(mod_inter_simpl)
## W = 0.99073, p-value = 1.3e-11dwtest(mod_inter_simpl)##
## Durbin-Watson test
##
## data: mod_inter_simpl
## DW = 1.9678, p-value = 0.2099
## alternative hypothesis: true autocorrelation is greater than 0# Confronto con modello completo
anova(mod_inter_simpl, mod_inter)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + I(Gestazione^2) + I(Lunghezza^2) + Cranio +
## Gestazione:Cranio + Lunghezza:Cranio
## Model 2: Peso ~ N.gravidanze + Gestazione + I(Gestazione^2) + Lunghezza +
## I(Lunghezza^2) + Cranio + Sesso + Gestazione * Lunghezza +
## Gestazione * Cranio + Lunghezza * Cranio + Lunghezza * Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2493 181164171
## 2 2488 177074776 5 4089395 11.492 5.458e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1I test statistici sull’analisi dei residui hanno restituito i seguenti risultati: Shapiro-Wilk: residui non perfettamente normali (p < 2.2e-16). Breusch-Pagan: eteroschedasticità significativa (p < 0.001), che viola l’ipotesi di varianza costante dei residui. Durbin-Watson: nessuna autocorrelazione rilevata (p > 0.1).
Per ovviare all’eteroschedasticità, è stata utilizzata la funzione coeftest() del pacchetto sandwich con stime robuste degli errori standard (vcovHC), confermando che le variabili principali rimangono significative anche con questa correzione.
È stato inoltre testato un modello semplificato (mod_inter_simpl), che conserva solo i termini quadratici e le interazioni più significative emerse dalla funzione coeftest(). Questo modello presenta:
BIC più alto (35134.46) rispetto a mod_inter (35116.50) R² aggiustato più basso (0.7364 vs 0.7418) RMSE più elevato (indicando una minore precisione nelle previsioni) ANOVA tra modelli: differenza statisticamente significativa (p < 0.001), a favore del modello completo
Concludendo, dopo tutti gli accertamenti statistici, confronti tra modelli e valutazioni di bontà dell’adattamento, mod_inter è risultato il modello migliore per prevedere il peso neonatale. Include termini quadratici e interazioni significative tra Gestazione, Lunghezza e Cranio. La diagnostica dei residui mostra leggera eteroschedasticità, ma non autocorrelazione.
PREVISIONE E RISULTATI
mean(Lunghezza) #494.692## [1] 494.692mean(Cranio) #340.0292## [1] 340.0292nuova_neonata <- data.frame(
N.gravidanze = 3,
Gestazione = 39,
Lunghezza = 495,
Cranio = 340,
Sesso = factor("F", levels = levels(dati$Sesso))
)
nuova_neonata$`I(Gestazione^2)` <- 39^2
nuova_neonata$`I(Lunghezza^2)` <- 500^2
nuova_neonata$`Gestazione:Lunghezza` <- 39 * 500
nuova_neonata$`Gestazione:Cranio` <- 39 * 340
nuova_neonata$`Lunghezza:Cranio` <- 500 * 340
nuova_neonata$`Lunghezza:SessoM` <- ifelse(nuova_neonata$Sesso == "M", 500, 0)
predict(mod_inter, newdata = nuova_neonata)## 1
## 3265.883# 3265.88 gPer 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. Il peso previsto risulta pari a circa 3266 grammi, coerente con quanto ci si aspettava.
#VISUALIZZAZIONI
#effetto della durata della gestazione sul peso
settimane <- seq(min(dati$Gestazione), max(dati$Gestazione), 1)
#nuovo dataset con le altre variabili fissate
dati_grafico <- data.frame(
N.gravidanze = 1,
Gestazione = settimane,
Lunghezza = 500,
Cranio = 340,
Sesso = factor("F", levels = levels(dati$Sesso)),
`Gestazione:Lunghezza` = settimane * 500,
`Gestazione:Cranio` = settimane * 340,
`Lunghezza:Cranio` = 500 * 340,
`Lunghezza:SessoM` = 0
)
# Peso previsto
dati_grafico$Peso_previsto <- predict(mod_inter, newdata = dati_grafico)
# Grafico
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)")
library(ggplot2)
# Costruzione del dataset per entrambi i sessi
dati_plot <- do.call(rbind, lapply(c("M", "F"), function(sesso) {
data.frame(
N.gravidanze = 1,
Gestazione = settimane,
Lunghezza = 500,
Cranio = 340,
Sesso = factor(sesso, levels = levels(dati$Sesso)),
`Gestazione:Lunghezza` = settimane * 500,
`Gestazione:Cranio` = settimane * 340,
`Lunghezza:Cranio` = 500 * 340,
`Lunghezza:SessoM` = ifelse(sesso == "M", 500, 0),
SessoLabel = ifelse(sesso == "M", "Maschio", "Femmina")
)
}))
# Calcolo del peso previsto
dati_plot$Peso_previsto <- predict(mod_inter, newdata = dati_plot)
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") +
scale_x_continuous(limits = c(min(dati$Gestazione), max(dati$Gestazione))) +
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.
Infine, è stato visualizzato un grafico che mostra le previsioni del peso atteso tra la 25ª e la 43ª settimana di gestazione, distinguendo tra maschi e femmine. Si evidenzia una crescita non lineare del peso con l’aumentare della gestazione e una differenza tra i sessi: a parità di settimane, i neonati maschi tendono ad avere un peso leggermente superiore. Questo risultato conferma l’efficacia predittiva del modello e l’importanza delle interazioni incluse.