Modello Statistico per la Previsione del Peso Neonato
Andrea Travascio
2025-12-21
Dettagli del Progetto
1) Raccolta dei Dati e Struttura del Dataset
DATASET: Per costruire il modello predittivo, abbiamo raccolto dati su 2500 neonati provenienti da tre ospedali.
- Età della madre: Misura dell’età in anni.
- Numero di gravidanze: Quante gravidanze ha avuto la madre.
- Fumo materno: Un indicatore binario (0=non fumatrice, 1=fumatrice).
- Durata della gravidanza: Numero di settimane di gestazione.
- Peso del neonato: Peso alla nascita in grammi.
- Lunghezza e diametro del cranio: Lunghezza del neonato e diametro craniale, misurabili anche durante la gravidanza tramite ecografie.
- Tipo di parto: Naturale o cesareo.
- Ospedale di nascita: Ospedale 1, 2 o 3.
- Sesso del neonato: Maschio (M) o femmina (F).
OBIETTIVO: L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature
2) Analisi e Modellizzazione
2.1) Analisi Preliminare
Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.
#------------ Esploro Variabile Risposta -------------
# funzione analisi variabili "continue"
DescribeVar <- function(var, LabelVar = "Variabile", alpha = 0.05) {
#----- Analisi distribuzione variabile ---------
print(LabelVar)
print(summary(var))
# Shapiro–Wilk test
st <- shapiro.test(var)
cat("Shapiro–Wilk test\n")
st
if (st$p.value <= alpha) {
cat("H0 rifiutata: i dati NON sono compatibili con una distribuzione normale.\n\n")
} else {
cat("H0 NON rifiutata: i dati sono compatibili con una distribuzione normale.\n\n")
}
# Skewness
sk <- moments::skewness(var, na.rm = TRUE)
cat("Skewness =", round(sk, 3), "\n")
if (abs(sk) <= 0.1) {
cat("Distribuzione approssimativamente simmetrica.\n")
} else if (sk < -0.1) {
cat("Distribuzione asimmetrica negativa (coda a sinistra).\n\n")
} else {
cat("Distribuzione asimmetrica positiva (coda a destra).\n\n")
}
# Kurtosis
kt <- moments::kurtosis(var)-3
cat("Kurtosis =", round(kt, 3), "\n")
if (abs(kt) == 0) {
cat("Mesocurtica: Simile alla distribuzione normale\n")
} else if (kt > 0) {
cat("Leptocurtica: Picco più alto, code più spesse, maggiore probabilità di outliers \n\n")
} else if (kt < 0){
cat("Platicurtica: Picco più basso, code più sottili, distribuzione più piatta e valori estremi meno probabili \n\n")
}
#Plot visualizzazione distribuzione valori
par(mfrow=c(1,2))
plot(density(var),
main = paste("Distribuzione di", LabelVar),
xlab = LabelVar)
abline(v=mean(var),col=2)
abline(v=median(var),col=3)
abline(v=quantile(var, probs = c(0.16, 0.84)),col = 4, lty = 2)
legend("topright",
legend = c("Media", "Mediana", "16perc", "84perc"),col = c(2, 3, 4, 4),
lty = c(1, 1, 2, 2),lwd = c(2, 2, 1, 1), bty = "n")
#----- Identificare eventuali outlier o anomalie ---------
Q1 <- quantile(var, 0.25, na.rm = TRUE)
Q3 <- quantile(var, 0.75, na.rm = TRUE)
IQR_val <- Q3 - Q1
outlier_idx <- which(
var < Q1 - 1.5*IQR_val |
var > Q3 + 1.5*IQR_val)
boxplot(var,
main = "Boxplot Check valori anomali",
ylab = LabelVar)
cat("Numero di valori anomali:", length(outlier_idx)," su ",length(var),"\n\n")
}
# funzione analisi variabili "binarie"
DescribeVar2 <- function(var, LabelVar = "Variabile"){
cat(LabelVar, "\n")
cat(rep("-", nchar(LabelVar)), "\n", sep = "")
tab <- prop.table(table(var))
print(tab)
barplot(table(var),
ylab="Frequenza",
main = LabelVar)
}
#------------ Esploro Variabile -------------
# inizio da Anni.madre perchè ho notato un problema
DescribeVar(Anni.madre,LabelVar="Anni madre (anni)")## [1] "Anni madre (anni)"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 25.00 28.00 28.16 32.00 46.00
## Shapiro–Wilk test
## H0 rifiutata: i dati NON sono compatibili con una distribuzione normale.
##
## Skewness = 0.043
## Distribuzione approssimativamente simmetrica.
## Kurtosis = 0.38
## Leptocurtica: Picco più alto, code più spesse, maggiore probabilità di outliers
## Numero di valori anomali: 13 su 2500# infatti Anni.madre ha un minimo di zero, non è plausibile.
dati[dati$Anni.madre <= 10, ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1152 1 1 0 41 3250 490 350
## 1380 0 0 0 39 3060 490 330
## Tipo.parto Ospedale Sesso
## 1152 Nat osp2 F
## 1380 Nat osp3 M# sotto i 10 anni ne ho solo due: anni=1,0 non reali, elimino queste righe per tutti i casi
dati_clean <- dati[Anni.madre > 12, ]
n <- nrow(dati_clean)
summary(dati_clean$Anni.madre)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13.00 25.00 28.00 28.19 32.00 46.00detach(dati)
attach(dati_clean)
DescribeVar(Peso,LabelVar="Peso (g)")## [1] "Peso (g)"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 2990 3300 3284 3620 4930
## Shapiro–Wilk test
## H0 rifiutata: i dati NON sono compatibili con una distribuzione normale.
##
## Skewness = -0.647
## Distribuzione asimmetrica negativa (coda a sinistra).
##
## Kurtosis = 2.029
## Leptocurtica: Picco più alto, code più spesse, maggiore probabilità di outliers
## Numero di valori anomali: 69 su 2498DescribeVar(Lunghezza,LabelVar="Lunghezza (mm)")## [1] "Lunghezza (mm)"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 310.0 480.0 500.0 494.7 510.0 565.0
## Shapiro–Wilk test
## H0 rifiutata: i dati NON sono compatibili con una distribuzione normale.
##
## Skewness = -1.515
## Distribuzione asimmetrica negativa (coda a sinistra).
##
## Kurtosis = 6.481
## Leptocurtica: Picco più alto, code più spesse, maggiore probabilità di outliers
## Numero di valori anomali: 59 su 2498DescribeVar(Cranio,LabelVar="Diamtero Cranio (mm)")## [1] "Diamtero Cranio (mm)"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 235 330 340 340 350 390
## Shapiro–Wilk test
## H0 rifiutata: i dati NON sono compatibili con una distribuzione normale.
##
## Skewness = -0.785
## Distribuzione asimmetrica negativa (coda a sinistra).
##
## Kurtosis = 2.945
## Leptocurtica: Picco più alto, code più spesse, maggiore probabilità di outliers
## Numero di valori anomali: 48 su 2498DescribeVar(Anni.madre,LabelVar="Anni madre (anni)")## [1] "Anni madre (anni)"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13.00 25.00 28.00 28.19 32.00 46.00
## Shapiro–Wilk test
## H0 rifiutata: i dati NON sono compatibili con una distribuzione normale.
##
## Skewness = 0.151
## Distribuzione asimmetrica positiva (coda a destra).
##
## Kurtosis = -0.106
## Platicurtica: Picco più basso, code più sottili, distribuzione più piatta e valori estremi meno probabili
## Numero di valori anomali: 11 su 2498#DescribeVar(Gestazione,LabelVar="Gestazione [settimane]")
# per vairibili binarie
Fumatrici_2 <- factor(Fumatrici,levels = c(0, 1),labels = c("No", "Sì"))
Tipo.parto_2 <- factor(Tipo.parto,levels = c("Nat", "Ces"),labels = c("Naturale", "Cesario"))
par(mfrow=c(1,3))
DescribeVar2(Fumatrici_2,LabelVar="Fumatrici")## Fumatrici
## ---------
## var
## No Sì
## 0.95836669 0.04163331DescribeVar2(Tipo.parto_2,LabelVar="Tipo di parto")## Tipo di parto
## -------------
## var
## Naturale Cesario
## 0.7085669 0.2914331DescribeVar2(Ospedale,LabelVar="Ospedale")## Ospedale
## --------
## var
## osp1 osp2 osp3
## 0.3266613 0.3394716 0.3338671
par(mfrow=c(1,3))
DescribeVar2(Sesso,LabelVar="Sesso")## Sesso
## -----
## var
## F M
## 0.5024019 0.4975981DescribeVar2(N.gravidanze,LabelVar="Numero gravitanze")## Numero gravitanze
## -----------------
## var
## 0 1 2 3 4 5
## 0.4383506805 0.3270616493 0.1361088871 0.0600480384 0.0192153723 0.0084067254
## 6 7 8 9 10 11
## 0.0044035228 0.0004003203 0.0032025620 0.0008006405 0.0012009608 0.0004003203
## 12
## 0.0004003203DescribeVar2(Gestazione,LabelVar="Gestazione")## Gestazione
## ----------
## var
## 25 26 27 28 29 30
## 0.0004003203 0.0004003203 0.0008006405 0.0016012810 0.0012009608 0.0020016013
## 31 32 33 34 35 36
## 0.0032025620 0.0036028823 0.0072057646 0.0064051241 0.0132105685 0.0248198559
## 37 38 39 40 41 42
## 0.0768614892 0.1749399520 0.2321857486 0.2966373098 0.1313050440 0.0224179343
## 43
## 0.0008006405
Inoltre si saggeranno le seguenti ipotesi con i test adatti: 1. in alcuni ospedali si fanno più parti cesarei 2. La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione 3. Le misure antropometriche sono significativamente diverse tra i due sessi
#----------- test -> in alcuni ospedali si fanno più parti cesarei?
# creare una tabella incrociata tipo di parto vs Ospedale e creare una matrice dei valori attesi
tab <- table(Tipo.parto_2, Ospedale)
test.indipendenza <- chisq.test(tab)
test.indipendenza##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 1.083, df = 2, p-value = 0.5819Il test Chi2 Pearson ci dice che non si può statisticamente rifiutare l’ipotesi di indipendenza (pvalue>0.05) suggerendo che la distribuzione di frequenze di parti naturali/cesarei non differisce in modo significativo tra i tre ospedali.
———– test -> media peso/lunghezza significativamente = quelli pop mondiale Ho trovato questo studio come riferimento https://www.scielo.br/j/spmj/a/5HBjGPShrJ4nYpPgQRDKt7m/?format=pdf&lang=en (Table 1) peso medio pop mondiale: mean= 3215.4 g sd= 488.5 g lughezza media pop mondiale: mean= 493 mm sd= 22 mm
——-Test z/t, H0: Peso = Peso mondiale ——–
BW_mondiale <- rnorm(10000000,mean=3215.4,sd=488.5)
plot(density(Peso),col = 1,lwd = 2, xlab = "Peso alla nascita (g)",
main = "Confronto densità peso alla nascita")
lines(density(BW_mondiale), col = 4,lwd = 2)
legend("topright",legend = c("Pop. campione","Pop. mondiale"),
col = c(1, 4),lwd = 2,bty = "n")
# H0: media = 3215.4 g
z.test(Peso,mu=3215.4,stdev=488.5,alternative="two.sided",conf.level=0.95)##
## One Sample z-test
##
## data: Peso
## z = 7.0375, n = 2498.0000, Std. Dev. = 488.5000, Std. Dev. of the
## sample mean = 9.7739, p-value = 1.957e-12
## alternative hypothesis: true mean is not equal to 3215.4
## 95 percent confidence interval:
## 3265.028 3303.341
## sample estimates:
## mean of Peso
## 3284.184# Z~7 e p-value ~ 0, quindi rifiuto l'ipotesi nulla = la media del campione è significativamente diversa da quella mondiale
# più appropriato usare il t-student per casi del mondo reale
t.test(Peso, mu=3215.4, alternative="two.sided", conf.level=0.95)##
## One Sample t-test
##
## data: Peso
## t = 6.5454, df = 2497, p-value = 7.176e-11
## alternative hypothesis: true mean is not equal to 3215.4
## 95 percent confidence interval:
## 3263.577 3304.791
## sample estimates:
## mean of x
## 3284.184# stessa conclusione di prima con t ~ 6.5z-test e t-test ci dicono che la media della distribuzione di valori ‘Peso’ di questo campione è significativamente diversi da quella mondiale
——-Test z/t, H0: Lunghezza = Lunghezza mondiale ——–
BL_mondiale <- rnorm(10000000,mean=493,sd=22)
plot(density(Lunghezza),col = 1,lwd = 2, xlab = "Lunghezza (mm)",
main = "Confronto densità lunghezza alla nascita")
lines(density(BL_mondiale), col = 4,lwd = 2)
legend("topright",legend = c("Pop. campione","Pop. mondiale"),
col = c(1, 4),lwd = 2,bty = "n")
# H0: media = 3215.4 g
z.test(Lunghezza,mu=493,stdev=22,alternative="two.sided",conf.level=0.95)##
## One Sample z-test
##
## data: Lunghezza
## z = 3.8525, n = 2498.00000, Std. Dev. = 22.00000, Std. Dev. of the
## sample mean = 0.44018, p-value = 0.0001169
## alternative hypothesis: true mean is not equal to 493
## 95 percent confidence interval:
## 493.8330 495.5585
## sample estimates:
## mean of Lunghezza
## 494.6958# z~3.9 e p-value=1.2e-4 -> riufiuto l'ipotesi ma entro solo 3 sigma rispetto a prima
t.test(Lunghezza, mu=493, alternative="two.sided", conf.level=0.95)##
## One Sample t-test
##
## data: Lunghezza
## t = 3.2191, df = 2497, p-value = 0.001303
## alternative hypothesis: true mean is not equal to 493
## 95 percent confidence interval:
## 493.6628 495.7287
## sample estimates:
## mean of x
## 494.6958#Anche in questo caso non cambia il risultato: z~3.2 e p-value=1.3e-3z-test e t-test ci dicono che la media della distribuzione di valori ‘Lunghezza’ di questo campione è significativamente diversi da quella mondiale
#------ misure antropometriche diverse tra i due sessi? ---------
# misure= Peso Lunghezza Cranio
# sessi = Sesso
print("-------Test t, H0: Peso(M) = Peso(F) --------")## [1] "-------Test t, H0: Peso(M) = Peso(F) --------"par(mfrow=c(1,3))
boxplot(Peso ~ Sesso)
t.test(Peso ~ Sesso, alternative = "two.sided", conf.level = 0.95)##
## Welch Two Sample t-test
##
## data: Peso by Sesso
## t = -12.115, df = 2488.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.4841 -207.3844
## sample estimates:
## mean in group F mean in group M
## 3161.061 3408.496print("-------Test t, H0: Lunghezza(M) = Lunghezza(F) --------")## [1] "-------Test t, H0: Lunghezza(M) = Lunghezza(F) --------"boxplot(Lunghezza ~ Sesso)
t.test(Lunghezza ~ Sesso, alternative = "two.sided", conf.level = 0.95)##
## Welch Two Sample t-test
##
## data: Lunghezza by Sesso
## t = -9.5823, df = 2457.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.939001 -7.882672
## sample estimates:
## mean in group F mean in group M
## 489.7641 499.6750print("-------Test t, H0: Cranio(M) = Cranio(F) --------")## [1] "-------Test t, H0: Cranio(M) = Cranio(F) --------"boxplot(Cranio ~ Sesso)
t.test(Cranio ~ Sesso, alternative = "two.sided", conf.level = 0.95)##
## Welch Two Sample t-test
##
## data: Cranio by Sesso
## t = -7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -6.110504 -3.560417
## sample estimates:
## mean in group F mean in group M
## 337.6231 342.4586In tutti i casi p-valu~0, quindi rifiutiamo H0, ed esiste una differenza significativa tra M e F. In particolare, vediamo che Peso, Lunghezza e Diametro Cranio dei Maschi è sempre maggiore di quello delle Femmine.
2.2) Creazione del Modello di Regressione
Verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo, potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato ed eventuali interazioni. Ad esempio, ci aspettiamo che una maggiore durata della gestazione aumenterebbe in media il peso del neonato.
2.3) Selezione del Modello Ottimale
Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC), selezioneremo il modello più parsimonioso, eliminando le variabili non significative. Verranno considerati anche modelli con interazioni tra le variabili e possibili effetti non lineari.
# Voglio solo valori numerici
dati_num <- dati_clean[,sapply(dati_clean, is.numeric)]
#covarianza =sum((Peso-mean(Peso))*(Gestazione-mean(Gestazione))/(n-1)
cov(Peso,Gestazione) ## [1] 581.0835#coefficente di correlazione lineare di Pearson
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)}
par(mfrow=c(1,1))
pairs(dati_num,upper.panel = panel.smooth,lower.panel = panel.cor)
round(cor(dati_num),2)## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza
## Anni.madre 1.00 0.38 0.01 -0.13 -0.02 -0.06
## N.gravidanze 0.38 1.00 0.05 -0.10 0.00 -0.06
## Fumatrici 0.01 0.05 1.00 0.03 -0.02 -0.02
## Gestazione -0.13 -0.10 0.03 1.00 0.59 0.62
## Peso -0.02 0.00 -0.02 0.59 1.00 0.80
## Lunghezza -0.06 -0.06 -0.02 0.62 0.80 1.00
## Cranio 0.02 0.04 -0.01 0.46 0.70 0.60
## Cranio
## Anni.madre 0.02
## N.gravidanze 0.04
## Fumatrici -0.01
## Gestazione 0.46
## Peso 0.70
## Lunghezza 0.60
## Cranio 1.00Dalla matrice di correlazione noto delle correlazioni positive tra peso e gestazione, lunghezza e diametro del cranio, ed una lieve correlazione positiva tra anni della madre e numero di gravidanze, che però è abbastanza ovvia come correlazione e non inerente al nostro scopo. Per il resto non abbiamo correlazioni evidenti.
par(mfrow=c(1,3))
plot(Peso, Gestazione, pch=20,
xlab="Peso (g)", ylab="Gestazione (settimane)", main="Peso vs Gestazione")
plot(Peso, Lunghezza, pch=20,
xlab="Peso (g)", ylab="Lunghezza (mm)", main="Peso vs Lunghezza")
plot(Peso, Cranio, pch=20,
xlab="Peso (g)", ylab="Cranio (mm)", main="Peso vs Cranio")
# 2.2) modello che include tutto
mod_all <- lm(Peso ~., data=dati_clean)
summary(mod_all)##
## Call:
## lm(formula = Peso ~ ., data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.26 -181.53 -14.45 161.05 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6735.7960 141.4790 -47.610 < 2e-16 ***
## Anni.madre 0.8018 1.1467 0.699 0.4845
## N.gravidanze 11.3812 4.6686 2.438 0.0148 *
## Fumatrici -30.2741 27.5492 -1.099 0.2719
## Gestazione 32.5773 3.8208 8.526 < 2e-16 ***
## Lunghezza 10.2922 0.3009 34.207 < 2e-16 ***
## Cranio 10.4722 0.4263 24.567 < 2e-16 ***
## Tipo.partoNat 29.6335 12.0905 2.451 0.0143 *
## Ospedaleosp2 -11.0912 13.4471 -0.825 0.4096
## Ospedaleosp3 28.2495 13.5054 2.092 0.0366 *
## SessoM 77.5723 11.1865 6.934 5.18e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274 on 2487 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 668.7 on 10 and 2487 DF, p-value: < 2.2e-16Anni.madre, Fumatrici e Ospedale sembrano non incluira sul parametro peso. Le caratteristiche che più sono legate al peso sono: Gestazione, che varia di ~32 ogni grammo, Lunghezza e Cranio, che variano circa 10 mm ogni grammo, e il sesso.
2.3) —————– Selezione del Modello Ottimale —————–
# Applico il metodo automatico step wise
mod_lin <- lm(Peso ~., data=dati_clean)
stepwise.mod_lin<-MASS::stepAIC(mod_lin,direction="both",k=log(n))## Start: AIC=28118.61
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 36710 186779904 28111
## - Fumatrici 1 90677 186833870 28112
## - Ospedale 2 687555 187430749 28112
## - N.gravidanze 1 446244 187189438 28117
## - Tipo.parto 1 451073 187194266 28117
## <none> 186743194 28119
## - Sesso 1 3610705 190353899 28159
## - Gestazione 1 5458852 192202046 28183
## - Cranio 1 45318506 232061700 28654
## - Lunghezza 1 87861708 274604902 29074
##
## Step: AIC=28111.28
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 91599 186871503 28105
## - Ospedale 2 693914 187473818 28105
## - Tipo.parto 1 452049 187231953 28110
## <none> 186779904 28111
## - N.gravidanze 1 631082 187410986 28112
## + Anni.madre 1 36710 186743194 28119
## - Sesso 1 3617809 190397713 28151
## - Gestazione 1 5424800 192204704 28175
## - Cranio 1 45569477 232349381 28649
## - Lunghezza 1 87852027 274631931 29066
##
## Step: AIC=28104.68
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 702925 187574428 28098
## - Tipo.parto 1 444404 187315907 28103
## <none> 186871503 28105
## - N.gravidanze 1 608136 187479640 28105
## + Fumatrici 1 91599 186779904 28111
## + Anni.madre 1 37633 186833870 28112
## - Sesso 1 3601860 190473363 28144
## - Gestazione 1 5358199 192229702 28168
## - Cranio 1 45613331 232484834 28642
## - Lunghezza 1 88259386 275130889 29063
##
## Step: AIC=28098.41
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 467626 188042054 28097
## <none> 187574428 28098
## - N.gravidanze 1 648873 188223301 28099
## + Ospedale 2 702925 186871503 28105
## + Fumatrici 1 100610 187473818 28105
## + Anni.madre 1 44184 187530244 28106
## - Sesso 1 3644818 191219246 28139
## - Gestazione 1 5457887 193032315 28162
## - Cranio 1 45747094 233321522 28636
## - Lunghezza 1 87955701 275530129 29051
##
## Step: AIC=28096.81
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 188042054 28097
## - N.gravidanze 1 621053 188663107 28097
## + Tipo.parto 1 467626 187574428 28098
## + Ospedale 2 726146 187315907 28103
## + Fumatrici 1 92548 187949505 28103
## + Anni.madre 1 45366 187996688 28104
## - Sesso 1 3650790 191692844 28137
## - Gestazione 1 5477493 193519547 28161
## - Cranio 1 46098547 234140601 28637
## - Lunghezza 1 87532691 275574744 29044print("Metodo STEP WISE \n")## [1] "Metodo STEP WISE \n"summary(stepwise.mod_lin)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.37 -180.98 -15.57 163.69 2639.09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.7251 135.8036 -49.201 < 2e-16 ***
## N.gravidanze 12.4554 4.3416 2.869 0.00415 **
## Gestazione 32.3827 3.8008 8.520 < 2e-16 ***
## Lunghezza 10.2455 0.3008 34.059 < 2e-16 ***
## Cranio 10.5410 0.4265 24.717 < 2e-16 ***
## SessoM 77.9807 11.2111 6.956 4.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1327 on 5 and 2492 DF, p-value: < 2.2e-16#Modello preferito: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
#Adjusted R-squared: 0.7265
# Verifico tra diversi casi manualmente
# Facciamo un check
# iniziamo il modello base
summary(mod_lin)##
## Call:
## lm(formula = Peso ~ ., data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.26 -181.53 -14.45 161.05 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6735.7960 141.4790 -47.610 < 2e-16 ***
## Anni.madre 0.8018 1.1467 0.699 0.4845
## N.gravidanze 11.3812 4.6686 2.438 0.0148 *
## Fumatrici -30.2741 27.5492 -1.099 0.2719
## Gestazione 32.5773 3.8208 8.526 < 2e-16 ***
## Lunghezza 10.2922 0.3009 34.207 < 2e-16 ***
## Cranio 10.4722 0.4263 24.567 < 2e-16 ***
## Tipo.partoNat 29.6335 12.0905 2.451 0.0143 *
## Ospedaleosp2 -11.0912 13.4471 -0.825 0.4096
## Ospedaleosp3 28.2495 13.5054 2.092 0.0366 *
## SessoM 77.5723 11.1865 6.934 5.18e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274 on 2487 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 668.7 on 10 and 2487 DF, p-value: < 2.2e-16# Adjusted R-squared: 0.7278
# Effatticamente Anni.madre,Fumatrici,Ospedale sono ininfluenti Pr(>|t|)>0.05
# N.gravidanze e Tipo.partoNat non sono così significative
# Provo a non rimuovere N.gravidanze e Tipo.partoNat (anche se penso non influisca sul peso della nascita, ma verifichiamolo scientificamente)
mod_lin1 <- update(mod_lin,~.-Anni.madre-Fumatrici-Ospedale)
summary(mod_lin1)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso, data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1129.14 -181.97 -16.26 160.95 2638.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6708.0171 136.0715 -49.298 < 2e-16 ***
## N.gravidanze 12.7356 4.3385 2.935 0.00336 **
## Gestazione 32.3253 3.7969 8.514 < 2e-16 ***
## Lunghezza 10.2833 0.3009 34.177 < 2e-16 ***
## Cranio 10.5063 0.4263 24.648 < 2e-16 ***
## Tipo.partoNat 30.1601 12.1027 2.492 0.01277 *
## SessoM 77.9171 11.1994 6.957 4.42e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.4 on 2491 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.727
## F-statistic: 1109 on 6 and 2491 DF, p-value: < 2.2e-16# Rimuovo Tipo.partoNat che ha un p-value minore.
mod_lin2 <- update(mod_lin1,~.-Tipo.partoNat)
summary(mod_lin2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso, data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1129.14 -181.97 -16.26 160.95 2638.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6708.0171 136.0715 -49.298 < 2e-16 ***
## N.gravidanze 12.7356 4.3385 2.935 0.00336 **
## Gestazione 32.3253 3.7969 8.514 < 2e-16 ***
## Lunghezza 10.2833 0.3009 34.177 < 2e-16 ***
## Cranio 10.5063 0.4263 24.648 < 2e-16 ***
## Tipo.partoNat 30.1601 12.1027 2.492 0.01277 *
## SessoM 77.9171 11.1994 6.957 4.42e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.4 on 2491 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.727
## F-statistic: 1109 on 6 and 2491 DF, p-value: < 2.2e-16#Infine il modello predetto migliore
mod_lin3 <- update(mod_lin2,~.-N.gravidanze)
summary(mod_lin3)##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso, data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1118.29 -185.67 -15.86 161.47 2625.55
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6676.5072 135.8546 -49.145 < 2e-16 ***
## Gestazione 31.2475 3.7849 8.256 2.42e-16 ***
## Lunghezza 10.2381 0.3009 34.019 < 2e-16 ***
## Cranio 10.6399 0.4245 25.067 < 2e-16 ***
## Tipo.partoNat 29.2394 12.1171 2.413 0.0159 *
## SessoM 79.0655 11.2097 7.053 2.25e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.8 on 2492 degrees of freedom
## Multiple R-squared: 0.7268, Adjusted R-squared: 0.7262
## F-statistic: 1326 on 5 and 2492 DF, p-value: < 2.2e-16# Verifico se effetivamente mod3 è il migliore
BIC(mod_lin,mod_lin1,mod_lin2,mod_lin3)## df BIC
## mod_lin 12 35215.45
## mod_lin1 8 35195.25
## mod_lin2 8 35195.25
## mod_lin3 7 35196.05# mod_lin3 è il più semplice e BIC è vicino al minore (di ~1)! Preferisco premiare il modello più semplice
# Nei plot noto che Peso~Lunghezza,Cranio,Gestazione hanno relazioni non proprio lineari
mod2_4 <- update(mod_lin3,~.+I(Lunghezza^2))
mod2_5 <- update(mod_lin3,~.+I(Cranio^2))
mod2_6 <- update(mod_lin3,~.+I(Gestazione^2))
BIC(mod_lin3,mod2_4,mod2_5,mod2_6)## df BIC
## mod_lin3 7 35196.05
## mod2_4 8 35114.89
## mod2_5 8 35170.19
## mod2_6 8 35198.94# mod2_4 migliora per il BIC
# vediamo l'f-test
anova(mod_lin3,mod2_4)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + Tipo.parto + Sesso
## Model 2: Peso ~ Gestazione + Lunghezza + Cranio + Tipo.parto + Sesso +
## I(Lunghezza^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 188223301
## 2 2491 181635999 1 6587302 90.34 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Pr(>F) < 2.2e-16 MIGLIORA! Evidentemente abbiamo bisogno di Lunghezza^2 nel modello
# dal BIC vediamo che Cranio^2 ha un impatto più forte su Gestazione^2
mod2_7 <- update(mod2_4,~.+I(Cranio^2))
anova(mod2_4,mod2_7)## Analysis of Variance Table
##
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + Tipo.parto + Sesso +
## I(Lunghezza^2)
## Model 2: Peso ~ Gestazione + Lunghezza + Cranio + Tipo.parto + Sesso +
## I(Lunghezza^2) + I(Cranio^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2491 181635999
## 2 2490 181633161 1 2838 0.0389 0.8437# Pr(>F) = 0.8429 -> non vale la pena aggiungere Cranio^2relazioni non lineari
# Verifichiamo con step wise anche l'aggiunta di relazioni non lineari
mod_nonlin <- lm(Peso ~.+I(Lunghezza^2)+I(Cranio^2)+I(Gestazione^2), data=dati_clean)
stepwise.mod_nonlin <-MASS::stepAIC(mod_nonlin,direction="both",k=log(n))## Start: AIC=28024.85
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso + I(Lunghezza^2) +
## I(Cranio^2) + I(Gestazione^2)
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 7980 178189981 28017
## - Cranio 1 22334 178204334 28017
## - Fumatrici 1 56274 178238274 28018
## - I(Cranio^2) 1 188204 178370204 28020
## - Ospedale 2 791297 178973297 28020
## - Tipo.parto 1 380030 178562030 28022
## <none> 178182000 28025
## - N.gravidanze 1 659734 178841735 28026
## - I(Gestazione^2) 1 1857057 180039057 28043
## - Gestazione 1 2338034 180520034 28050
## - Sesso 1 3085862 181267862 28060
## - Lunghezza 1 3113588 181295588 28060
## - I(Lunghezza^2) 1 5778914 183960914 28097
##
## Step: AIC=28017.14
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + I(Lunghezza^2) + I(Cranio^2) +
## I(Gestazione^2)
##
## Df Sum of Sq RSS AIC
## - Cranio 1 22818 178212798 28010
## - Fumatrici 1 56604 178246584 28010
## - I(Cranio^2) 1 189810 178379790 28012
## - Ospedale 2 795220 178985201 28013
## - Tipo.parto 1 380418 178570399 28015
## <none> 178189981 28017
## - N.gravidanze 1 827253 179017233 28021
## + Anni.madre 1 7980 178182000 28025
## - I(Gestazione^2) 1 1875592 180065572 28036
## - Gestazione 1 2356644 180546624 28042
## - Sesso 1 3089555 181279536 28052
## - Lunghezza 1 3125072 181315053 28053
## - I(Lunghezza^2) 1 5795552 183985532 28089
##
## Step: AIC=28009.64
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Tipo.parto +
## Ospedale + Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2)
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 57311 178270110 28003
## - Ospedale 2 788143 179000941 28005
## - Tipo.parto 1 383021 178595820 28007
## <none> 178212798 28010
## - N.gravidanze 1 829371 179042169 28013
## + Cranio 1 22818 178189981 28017
## + Anni.madre 1 8464 178204334 28017
## - I(Gestazione^2) 1 1992242 180205040 28030
## - Gestazione 1 2537680 180750478 28037
## - Sesso 1 3098406 181311205 28045
## - Lunghezza 1 4041673 182254471 28058
## - I(Lunghezza^2) 1 7329961 185542760 28102
## - I(Cranio^2) 1 44354850 222567648 28557
##
## Step: AIC=28002.62
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Tipo.parto + Ospedale +
## Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2)
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 795733 179065843 27998
## - Tipo.parto 1 377311 178647421 28000
## <none> 178270110 28003
## - N.gravidanze 1 809494 179079603 28006
## + Fumatrici 1 57311 178212798 28010
## + Cranio 1 23525 178246584 28010
## + Anni.madre 1 8813 178261297 28010
## - I(Gestazione^2) 1 1989309 180259418 28022
## - Gestazione 1 2531989 180802098 28030
## - Sesso 1 3085257 181355366 28038
## - Lunghezza 1 4055525 182325634 28051
## - I(Lunghezza^2) 1 7354099 185624209 28096
## - I(Cranio^2) 1 44396326 222666436 28550
##
## Step: AIC=27998.1
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Tipo.parto + Sesso +
## I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2)
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 399869 179465712 27996
## <none> 179065843 27998
## - N.gravidanze 1 858225 179924067 28002
## + Ospedale 2 795733 178270110 28003
## + Fumatrici 1 64902 179000941 28005
## + Cranio 1 16327 179049516 28006
## + Anni.madre 1 12845 179052998 28006
## - I(Gestazione^2) 1 1889559 180955402 28016
## - Gestazione 1 2423950 181489792 28024
## - Sesso 1 3119397 182185239 28033
## - Lunghezza 1 3977509 183043352 28045
## - I(Lunghezza^2) 1 7243620 186309463 28089
## - I(Cranio^2) 1 44578202 223644045 28546
##
## Step: AIC=27995.84
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Sesso + I(Lunghezza^2) +
## I(Cranio^2) + I(Gestazione^2)
##
## Df Sum of Sq RSS AIC
## <none> 179465712 27996
## + Tipo.parto 1 399869 179065843 27998
## - N.gravidanze 1 829418 180295130 28000
## + Ospedale 2 818291 178647421 28000
## + Fumatrici 1 58757 179406955 28003
## + Cranio 1 18449 179447262 28003
## + Anni.madre 1 13436 179452276 28004
## - I(Gestazione^2) 1 1880915 181346627 28014
## - Gestazione 1 2415573 181881285 28021
## - Sesso 1 3120261 182585973 28031
## - Lunghezza 1 4018162 183483874 28043
## - I(Lunghezza^2) 1 7286352 186752064 28087
## - I(Cranio^2) 1 44919559 224385270 28546summary(stepwise.mod_nonlin)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Sesso +
## I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2), data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1188.23 -181.87 -12.06 163.33 1447.93
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.584e+03 9.066e+02 -1.747 0.080769 .
## N.gravidanze 1.441e+01 4.247e+00 3.392 0.000704 ***
## Gestazione 3.624e+02 6.260e+01 5.789 7.96e-09 ***
## Lunghezza -3.016e+01 4.039e+00 -7.467 1.13e-13 ***
## SessoM 7.239e+01 1.100e+01 6.580 5.73e-11 ***
## I(Lunghezza^2) 4.167e-02 4.145e-03 10.055 < 2e-16 ***
## I(Cranio^2) 1.542e-02 6.179e-04 24.965 < 2e-16 ***
## I(Gestazione^2) -4.206e+00 8.233e-01 -5.109 3.49e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 268.5 on 2490 degrees of freedom
## Multiple R-squared: 0.7395, Adjusted R-squared: 0.7387
## F-statistic: 1010 on 7 and 2490 DF, p-value: < 2.2e-16# Adjusted R-squared: 0.7387
# Peso ~ N.gravidanze + Gestazione + Lunghezza + Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2)
# La nuova formula mostra che Cranio entra solo al quadrato e Gestazione^2 migliora in questo modo
BIC(stepwise.mod_nonlin,mod2_4)## df BIC
## stepwise.mod_nonlin 9 35092.68
## mod2_4 8 35114.89AIC(stepwise.mod_nonlin,mod2_4)## df AIC
## stepwise.mod_nonlin 9 35040.28
## mod2_4 8 35068.30# step wise sia dal BIC e AIC sembra miglioreInterazioni tra variabili?
# ho provato altre combinazioni meno aspettate, anche con combinazioni
mod_nonlin_inter <- lm(Peso ~.+I(Lunghezza^2)+I(Cranio^2)+I(Gestazione^2)+Lunghezza*Cranio + Lunghezza*Gestazione + Cranio*Gestazione, data=dati_clean)
stepwise.mod_nonlin_inter <-MASS::stepAIC(mod_nonlin_inter,direction="both",k=log(n))## Start: AIC=28007.6
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso + I(Lunghezza^2) +
## I(Cranio^2) + I(Gestazione^2) + Lunghezza * Cranio + Lunghezza *
## Gestazione + Cranio * Gestazione
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 10064 175311114 28000
## - Fumatrici 1 50250 175351300 28000
## - Gestazione:Lunghezza 1 124152 175425203 28002
## - Ospedale 2 770639 176071690 28003
## - Tipo.parto 1 399770 175700820 28006
## - I(Cranio^2) 1 509902 175810953 28007
## <none> 175301051 28008
## - I(Gestazione^2) 1 629043 175930093 28009
## - N.gravidanze 1 692836 175993886 28010
## - Lunghezza:Cranio 1 1851489 177152539 28026
## - Gestazione:Cranio 1 1873762 177174813 28026
## - Sesso 1 3155808 178456859 28044
## - I(Lunghezza^2) 1 5215228 180516279 28073
##
## Step: AIC=27999.92
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + I(Lunghezza^2) + I(Cranio^2) +
## I(Gestazione^2) + Lunghezza:Cranio + Gestazione:Lunghezza +
## Gestazione:Cranio
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 50592 175361707 27993
## - Gestazione:Lunghezza 1 125002 175436116 27994
## - Ospedale 2 775142 176086256 27995
## - Tipo.parto 1 400171 175711286 27998
## - I(Cranio^2) 1 512981 175824096 27999
## <none> 175311114 28000
## - I(Gestazione^2) 1 631337 175942451 28001
## - N.gravidanze 1 874864 176185979 28004
## + Anni.madre 1 10064 175301051 28008
## - Lunghezza:Cranio 1 1849800 177160915 28018
## - Gestazione:Cranio 1 1869926 177181040 28019
## - Sesso 1 3160306 178471420 28037
## - I(Lunghezza^2) 1 5233000 180544115 28066
##
## Step: AIC=27992.82
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) +
## Lunghezza:Cranio + Gestazione:Lunghezza + Gestazione:Cranio
##
## Df Sum of Sq RSS AIC
## - Gestazione:Lunghezza 1 124864 175486570 27987
## - Ospedale 2 782529 176144236 27988
## - Tipo.parto 1 394695 175756402 27991
## - I(Cranio^2) 1 518697 175880403 27992
## <none> 175361707 27993
## - I(Gestazione^2) 1 630766 175992472 27994
## - N.gravidanze 1 855818 176217524 27997
## + Fumatrici 1 50592 175311114 28000
## + Anni.madre 1 10406 175351300 28000
## - Lunghezza:Cranio 1 1857054 177218761 28011
## - Gestazione:Cranio 1 1869066 177230773 28012
## - Sesso 1 3148061 178509768 28029
## - I(Lunghezza^2) 1 5246721 180608427 28059
##
## Step: AIC=27986.77
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) +
## Lunghezza:Cranio + Gestazione:Cranio
##
## Df Sum of Sq RSS AIC
## - Ospedale 2 793335 176279906 27982
## - Tipo.parto 1 402656 175889226 27985
## <none> 175486570 27987
## - I(Cranio^2) 1 642784 176129354 27988
## - N.gravidanze 1 858843 176345413 27991
## + Gestazione:Lunghezza 1 124864 175361707 27993
## + Fumatrici 1 50454 175436116 27994
## + Anni.madre 1 11269 175475302 27994
## - Gestazione:Cranio 1 1763042 177249612 28004
## - Lunghezza:Cranio 1 2454257 177940828 28014
## - I(Gestazione^2) 1 2779344 178265914 28018
## - Sesso 1 3110845 178597416 28023
## - I(Lunghezza^2) 1 7122586 182609156 28078
##
## Step: AIC=27982.4
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) +
## Lunghezza:Cranio + Gestazione:Cranio
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 426654 176706559 27981
## <none> 176279906 27982
## - I(Cranio^2) 1 621528 176901433 27983
## + Ospedale 2 793335 175486570 27987
## - N.gravidanze 1 906927 177186833 27987
## + Gestazione:Lunghezza 1 135670 176144236 27988
## + Fumatrici 1 57917 176221989 27989
## + Anni.madre 1 16199 176263707 27990
## - Gestazione:Cranio 1 1750411 178030317 27999
## - Lunghezza:Cranio 1 2473678 178753584 28009
## - I(Gestazione^2) 1 2650017 178929923 28012
## - Sesso 1 3145874 179425780 28019
## - I(Lunghezza^2) 1 7104425 183384331 28073
##
## Step: AIC=27980.61
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) + Lunghezza:Cranio +
## Gestazione:Cranio
##
## Df Sum of Sq RSS AIC
## <none> 176706559 27981
## - I(Cranio^2) 1 629701 177336260 27982
## + Tipo.parto 1 426654 176279906 27982
## + Ospedale 2 817333 175889226 27985
## - N.gravidanze 1 875974 177582533 27985
## + Gestazione:Lunghezza 1 144400 176562160 27986
## + Fumatrici 1 51930 176654629 27988
## + Anni.madre 1 16822 176689737 27988
## - Gestazione:Cranio 1 1730936 178437495 27997
## - Lunghezza:Cranio 1 2448504 179155063 28007
## - I(Gestazione^2) 1 2640382 179346941 28010
## - Sesso 1 3146204 179852764 28017
## - I(Lunghezza^2) 1 7109371 183815931 28071summary(stepwise.mod_nonlin_inter)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) +
## Lunghezza:Cranio + Gestazione:Cranio, data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1167.2 -182.7 -11.5 168.3 1318.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.014e+02 1.191e+03 -0.169 0.865697
## N.gravidanze 1.481e+01 4.218e+00 3.511 0.000454 ***
## Gestazione 1.795e+02 7.436e+01 2.414 0.015847 *
## Lunghezza -7.196e+00 5.659e+00 -1.272 0.203647
## Cranio -2.057e+01 1.006e+01 -2.044 0.041019 *
## SessoM 7.273e+01 1.093e+01 6.654 3.49e-11 ***
## I(Lunghezza^2) 5.045e-02 5.044e-03 10.003 < 2e-16 ***
## I(Cranio^2) 5.477e-02 1.840e-02 2.977 0.002939 **
## I(Gestazione^2) -6.294e+00 1.033e+00 -6.096 1.26e-09 ***
## Lunghezza:Cranio -9.274e-02 1.580e-02 -5.870 4.93e-09 ***
## Gestazione:Cranio 1.015e+00 2.057e-01 4.936 8.51e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.6 on 2487 degrees of freedom
## Multiple R-squared: 0.7435, Adjusted R-squared: 0.7424
## F-statistic: 720.8 on 10 and 2487 DF, p-value: < 2.2e-16# Adjusted R-squared: 0.7424
# Il modello è: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) + Lunghezza:Cranio + Gestazione:Cranio
# il best model include anche Lunghezza:Cranio e Gestazione:Cranio
# con ottime probabilità, a discapito delle relazioni lineari con le stesse caratteristiche.Tuttavia, in questo caso Lunghezza Pr(>|t|) >0.05
stepwise.mod_nonlin_inter2 <- update(stepwise.mod_nonlin_inter,~.-Lunghezza)
summary(stepwise.mod_nonlin_inter2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Cranio + Sesso +
## I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) + Lunghezza:Cranio +
## Gestazione:Cranio, data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1160.25 -181.37 -12.23 168.14 1316.71
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.939e+01 1.177e+03 0.025 0.980079
## N.gravidanze 1.470e+01 4.218e+00 3.484 0.000502 ***
## Gestazione 1.289e+02 6.284e+01 2.052 0.040281 *
## Cranio -2.655e+01 8.895e+00 -2.985 0.002861 **
## SessoM 7.290e+01 1.093e+01 6.670 3.14e-11 ***
## I(Lunghezza^2) 4.731e-02 4.397e-03 10.759 < 2e-16 ***
## I(Cranio^2) 6.521e-02 1.647e-02 3.959 7.73e-05 ***
## I(Gestazione^2) -6.191e+00 1.029e+00 -6.014 2.08e-09 ***
## Cranio:Lunghezza -1.049e-01 1.260e-02 -8.318 < 2e-16 ***
## Gestazione:Cranio 1.140e+00 1.807e-01 6.313 3.23e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.6 on 2488 degrees of freedom
## Multiple R-squared: 0.7433, Adjusted R-squared: 0.7424
## F-statistic: 800.5 on 9 and 2488 DF, p-value: < 2.2e-16# Gestazione ha P(r) ~0.04, che accade se lo rimuovo?
stepwise.mod_nonlin_inter3 <- update(stepwise.mod_nonlin_inter2,~.-Gestazione)
summary(stepwise.mod_nonlin_inter3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Cranio + Sesso + I(Lunghezza^2) +
## I(Cranio^2) + I(Gestazione^2) + Cranio:Lunghezza + Gestazione:Cranio,
## data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1157.80 -178.88 -9.86 168.26 1304.88
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.343e+02 1.152e+03 0.464 0.642752
## N.gravidanze 1.472e+01 4.221e+00 3.487 0.000496 ***
## Cranio -1.522e+01 6.976e+00 -2.181 0.029241 *
## SessoM 7.179e+01 1.092e+01 6.572 6.03e-11 ***
## I(Lunghezza^2) 4.795e-02 4.389e-03 10.924 < 2e-16 ***
## I(Cranio^2) 4.578e-02 1.348e-02 3.395 0.000697 ***
## I(Gestazione^2) -4.812e+00 7.803e-01 -6.167 8.11e-10 ***
## Cranio:Lunghezza -1.065e-01 1.259e-02 -8.465 < 2e-16 ***
## Cranio:Gestazione 1.212e+00 1.774e-01 6.833 1.04e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 266.8 on 2489 degrees of freedom
## Multiple R-squared: 0.7429, Adjusted R-squared: 0.742
## F-statistic: 898.9 on 8 and 2489 DF, p-value: < 2.2e-16# Cranio ha P(r) ~0.03, che accade se lo rimuovo?
stepwise.mod_nonlin_inter4 <- update(stepwise.mod_nonlin_inter3,~.-Cranio)
summary(stepwise.mod_nonlin_inter4) ##
## Call:
## lm(formula = Peso ~ N.gravidanze + Sesso + I(Lunghezza^2) + I(Cranio^2) +
## I(Gestazione^2) + Cranio:Lunghezza + Cranio:Gestazione, data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1161.94 -179.49 -7.93 168.32 1310.33
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.973e+03 7.017e+01 -28.124 < 2e-16 ***
## N.gravidanze 1.458e+01 4.223e+00 3.453 0.000564 ***
## SessoM 7.341e+01 1.091e+01 6.730 2.09e-11 ***
## I(Lunghezza^2) 4.836e-02 4.388e-03 11.021 < 2e-16 ***
## I(Cranio^2) 2.757e-02 1.060e-02 2.601 0.009342 **
## I(Gestazione^2) -4.590e+00 7.742e-01 -5.928 3.49e-09 ***
## Cranio:Lunghezza -1.080e-01 1.258e-02 -8.583 < 2e-16 ***
## Cranio:Gestazione 1.155e+00 1.755e-01 6.577 5.82e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267 on 2490 degrees of freedom
## Multiple R-squared: 0.7424, Adjusted R-squared: 0.7417
## F-statistic: 1025 on 7 and 2490 DF, p-value: < 2.2e-16# Sono rimasti tutti valori estremamente dipendenti con P(r) <<0.01
# Adjusted R-squared 74% (non è cambiato molto)
BIC(stepwise.mod_nonlin_inter,stepwise.mod_nonlin,stepwise.mod_nonlin_inter2,stepwise.mod_nonlin_inter3,stepwise.mod_nonlin_inter4)## df BIC
## stepwise.mod_nonlin_inter 12 35077.45
## stepwise.mod_nonlin 9 35092.68
## stepwise.mod_nonlin_inter2 11 35071.25
## stepwise.mod_nonlin_inter3 10 35067.65
## stepwise.mod_nonlin_inter4 9 35064.60AIC(stepwise.mod_nonlin_inter,stepwise.mod_nonlin,stepwise.mod_nonlin_inter2,stepwise.mod_nonlin_inter3,stepwise.mod_nonlin_inter4)## df AIC
## stepwise.mod_nonlin_inter 12 35007.57
## stepwise.mod_nonlin 9 35040.28
## stepwise.mod_nonlin_inter2 11 35007.20
## stepwise.mod_nonlin_inter3 10 35009.42
## stepwise.mod_nonlin_inter4 9 35012.19# BIC/AIC stepwise.mod_nonlin_inter4 più bassi
anova(stepwise.mod_nonlin_inter,stepwise.mod_nonlin_inter4)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) + Lunghezza:Cranio +
## Gestazione:Cranio
## Model 2: Peso ~ N.gravidanze + Sesso + I(Lunghezza^2) + I(Cranio^2) +
## I(Gestazione^2) + Cranio:Lunghezza + Cranio:Gestazione
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2487 176706559
## 2 2490 177459324 -3 -752765 3.5315 0.01427 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(stepwise.mod_nonlin_inter2,stepwise.mod_nonlin_inter4)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Cranio + Sesso + I(Lunghezza^2) +
## I(Cranio^2) + I(Gestazione^2) + Lunghezza:Cranio + Gestazione:Cranio
## Model 2: Peso ~ N.gravidanze + Sesso + I(Lunghezza^2) + I(Cranio^2) +
## I(Gestazione^2) + Cranio:Lunghezza + Cranio:Gestazione
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2488 176821442
## 2 2490 177459324 -2 -637882 4.4877 0.01134 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(stepwise.mod_nonlin_inter3,stepwise.mod_nonlin_inter4)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Cranio + Sesso + I(Lunghezza^2) + I(Cranio^2) +
## I(Gestazione^2) + Cranio:Lunghezza + Gestazione:Cranio
## Model 2: Peso ~ N.gravidanze + Sesso + I(Lunghezza^2) + I(Cranio^2) +
## I(Gestazione^2) + Cranio:Lunghezza + Cranio:Gestazione
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2489 177120675
## 2 2490 177459324 -1 -338649 4.7589 0.02924 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# anova: in tutti i casi l'aumento dei parametri non significativamentre giustificato
# Domanda logica, e poco scientifica: veramente il numero di gravidanze è importante ai fini del peso del neonato?
# Vediamo se la rimuoviamo che accade
stepwise.mod_nonlin_inter5 <- update(stepwise.mod_nonlin_inter4,~.-N.gravidanze)
summary(stepwise.mod_nonlin_inter5) ##
## Call:
## lm(formula = Peso ~ Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) +
## Cranio:Lunghezza + Cranio:Gestazione, data = dati_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1148.65 -183.62 -10.63 167.47 1298.28
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.948e+03 6.994e+01 -27.855 < 2e-16 ***
## SessoM 7.475e+01 1.092e+01 6.843 9.75e-12 ***
## I(Lunghezza^2) 4.788e-02 4.396e-03 10.892 < 2e-16 ***
## I(Cranio^2) 2.817e-02 1.062e-02 2.652 0.00806 **
## I(Gestazione^2) -4.508e+00 7.755e-01 -5.813 6.93e-09 ***
## Cranio:Lunghezza -1.067e-01 1.260e-02 -8.468 < 2e-16 ***
## Cranio:Gestazione 1.132e+00 1.758e-01 6.441 1.42e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.5 on 2491 degrees of freedom
## Multiple R-squared: 0.7411, Adjusted R-squared: 0.7405
## F-statistic: 1189 on 6 and 2491 DF, p-value: < 2.2e-16# Adjusted R-squared praticamente ~ invariato
BIC(stepwise.mod_nonlin_inter5,stepwise.mod_nonlin_inter4)## df BIC
## stepwise.mod_nonlin_inter5 8 35068.71
## stepwise.mod_nonlin_inter4 9 35064.60anova(stepwise.mod_nonlin_inter5,stepwise.mod_nonlin_inter4)## Analysis of Variance Table
##
## Model 1: Peso ~ Sesso + I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) +
## Cranio:Lunghezza + Cranio:Gestazione
## Model 2: Peso ~ N.gravidanze + Sesso + I(Lunghezza^2) + I(Cranio^2) +
## I(Gestazione^2) + Cranio:Lunghezza + Cranio:Gestazione
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2491 178308981
## 2 2490 177459324 1 849657 11.922 0.000564 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# sia BIC che anova mi dice che è meglio tenere N.gravidanze nel modello. Probabilmente esiste una dipendenza n-Dim che non visualizziamo ad occhio
# giusto per sicurezza lo confronto con i modelli best che non includevano i parametri di non linearità e di interazione ottenuti da step wise
BIC(stepwise.mod_nonlin_inter,stepwise.mod_lin,stepwise.mod_nonlin_inter4)## df BIC
## stepwise.mod_nonlin_inter 12 35077.45
## stepwise.mod_lin 7 35193.65
## stepwise.mod_nonlin_inter4 9 35064.60anova(stepwise.mod_nonlin_inter,stepwise.mod_nonlin_inter4)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## I(Lunghezza^2) + I(Cranio^2) + I(Gestazione^2) + Lunghezza:Cranio +
## Gestazione:Cranio
## Model 2: Peso ~ N.gravidanze + Sesso + I(Lunghezza^2) + I(Cranio^2) +
## I(Gestazione^2) + Cranio:Lunghezza + Cranio:Gestazione
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2487 176706559
## 2 2490 177459324 -3 -752765 3.5315 0.01427 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# continua ad essere questo il miglior modelloConclusione best model:
print(stepwise.mod_nonlin_inter4$call)## lm(formula = Peso ~ N.gravidanze + Sesso + I(Lunghezza^2) + I(Cranio^2) +
## I(Gestazione^2) + Cranio:Lunghezza + Cranio:Gestazione, data = dati_clean)bestmodel <- stepwise.mod_nonlin_inter4
#Adjusted R-squared del 74% va abbastanza bene. 2.4) Analisi della Qualità del Modello
Una volta ottenuto il modello finale, valuteremo la sua capacità predittiva utilizzando metriche come R² e il Root Mean Squared Error (RMSE). Un’attenzione particolare sarà rivolta all’analisi dei residui e alla presenza di valori influenti, che potrebbero distorcere le previsioni, indagando su di essi.
par(mfrow=c(2,2))
plot(bestmodel)
Residuals vs Fitted: i punti sono sparsi grosso modo casualmente tranne
per un pattern che si trova per un numero piccolo di punti per bassi
valori di Fitted values Q-Q Residuals: tra -2 e +2 del theoretical
quantiles abbiamo i punti ben allineati, ma alle ali abbiamo dei valori
dei quantili che discostano da quelli aspettati per una normale, in
particolare i punti 1399 per l’ala inferiore, e 155 e 1306 per l’ala
superiore. Scale-Location: anche qui non vi è alcun pattern, e sembra
esserci una varianza costante, con un fit piatto. Residuals vs Leverage:
in questo plot non ci sono valori outliers o leverage eccetto un valore
(1551) che supera il valore 1 (Allarme!) della distanza di cook. 1551 è
un residuo di osservazione potenzialmente influente sulle stime di
regressione.
#Usiamo un metodo numerico per indagare valori di Leverage e Outliers
#Leverage
par(mfrow=c(1,1))
lev <- hatvalues(bestmodel)
plot(lev)
p=sum(lev)
soglia <- 2*p/n
abline(h=soglia, col=2)
lev[lev>soglia]## 1 3 15 34 42 61
## 0.007795012 0.006931860 0.007492752 0.011768702 0.010421340 0.009144867
## 89 93 96 101 106 113
## 0.012923339 0.006499817 0.012405822 0.007078324 0.013578751 0.007347929
## 131 134 151 155 161 190
## 0.007142078 0.007747457 0.057828728 0.010192614 0.021238035 0.011039910
## 204 205 206 220 295 310
## 0.014746297 0.011732195 0.008909822 0.015470634 0.008667683 0.371814436
## 312 315 328 348 353 378
## 0.012131730 0.009573681 0.007744695 0.006423532 0.006422903 0.013741055
## 383 442 445 471 492 516
## 0.006537369 0.007878217 0.007461676 0.008718369 0.008433961 0.013297839
## 540 566 582 587 592 638
## 0.007108068 0.006465566 0.011867639 0.007871116 0.007918075 0.006810504
## 657 666 684 697 701 748
## 0.009542239 0.007765327 0.021091550 0.013308701 0.009128566 0.008107601
## 750 757 805 821 828 895
## 0.006679316 0.008321821 0.012574130 0.010138390 0.007304080 0.010048629
## 928 946 947 949 956 959
## 0.019960989 0.007636703 0.008137928 0.006499057 0.014263201 0.006767627
## 964 985 991 1014 1067 1091
## 0.009928789 0.008243689 0.006413592 0.007376424 0.008087891 0.008442248
## 1105 1130 1155 1166 1181 1188
## 0.006444593 0.033023289 0.008332296 0.008349269 0.016626241 0.020889242
## 1200 1218 1219 1248 1273 1293
## 0.013583439 0.008064255 0.030923670 0.013476010 0.006723945 0.011037878
## 1294 1311 1318 1321 1356 1357
## 0.007222811 0.009954726 0.007165393 0.010289277 0.007966786 0.006644323
## 1361 1385 1395 1402 1411 1420
## 0.007830674 0.012312717 0.009091072 0.007831738 0.008383580 0.010743346
## 1428 1429 1450 1505 1551 1553
## 0.006972053 0.130404156 0.015372832 0.013493882 0.650854940 0.017158491
## 1556 1586 1610 1619 1686 1693
## 0.019128023 0.007324893 0.009361823 0.015328655 0.009763020 0.014444620
## 1697 1701 1712 1718 1727 1780
## 0.006704751 0.010935192 0.019189977 0.012662286 0.013379251 0.019845192
## 1781 1802 1809 1858 1868 1911
## 0.016997553 0.008905700 0.008300971 0.009557750 0.007893529 0.009377636
## 1977 2037 2040 2086 2114 2115
## 0.013410868 0.008779528 0.039275720 0.013349900 0.012068380 0.011890402
## 2120 2136 2148 2149 2175 2200
## 0.015779328 0.006690562 0.008102874 0.012630036 0.028800271 0.011046432
## 2215 2216 2220 2221 2244 2307
## 0.012622863 0.013380851 0.010331689 0.021786510 0.007033774 0.012824098
## 2317 2318 2337 2359 2408 2422
## 0.007954581 0.012821514 0.007668370 0.009987638 0.009226119 0.021820700
## 2437 2452 2458 2471 2476
## 0.020851158 0.019390703 0.008140225 0.021123045 0.008077471# Vediamo che percentuali di valori sono sopra la soglia
infl <- which(lev > soglia)
nlev <- length(infl)
perclev <- round(length(infl)*100/n,1)
cat("% Leverage =", perclev, "% \n")## % Leverage = 5.7 %cat("Num Leverage =", nlev, "\n")## Num Leverage = 143# 143 valori, che però rappresentano il 5.7% dei valori sono sopra la soglia. In realtà non sono significativi
#Outliers
plot(rstudent(bestmodel))
abline(h=c(-2,2), col=2)
outlierTest(bestmodel)## rstudent unadjusted p-value Bonferroni p
## 1306 4.935596 8.5189e-07 0.002128
## 155 4.595439 4.5339e-06 0.011326
## 1399 -4.373741 1.2716e-05 0.031765
## 1694 4.331685 1.5384e-05 0.038429cat("Num Outliers =", length(outlierTest(bestmodel))-1, "\n")## Num Outliers = 4# 4 outliers: 1306, 155, 1399, 1694
#distanza di cook
cookD <- cooks.distance(bestmodel)
plot(cookD)
cat("Max CookDistance =", max(cookD), "\n")## Max CookDistance = 3.482531# Test sui residui
bptest(bestmodel)##
## studentized Breusch-Pagan test
##
## data: bestmodel
## BP = 47.611, df = 7, p-value = 4.241e-08cat("bptest: rifiutata H0 -> varianza non costante")## bptest: rifiutata H0 -> varianza non costantedwtest(bestmodel)##
## Durbin-Watson test
##
## data: bestmodel
## DW = 1.9582, p-value = 0.1476
## alternative hypothesis: true autocorrelation is greater than 0cat("dwtest: non rifiuta H0 -> residui non sono autocorrelati")## dwtest: non rifiuta H0 -> residui non sono autocorrelatishapiro.test(residuals(bestmodel))##
## Shapiro-Wilk normality test
##
## data: residuals(bestmodel)
## W = 0.99106, p-value = 2.474e-11cat("shapiro.test: rifiuta H0 -> residui non sono normali")## shapiro.test: rifiuta H0 -> residui non sono normalicat("-------check visivo PDF dei residui-------")## -------check visivo PDF dei residui-------plot(density(residuals(bestmodel)))
# alle ali non è normale, però la maggior parte lo sono.
# Deviazione dalla normalità, principalmente nelle code della distribuzione.
# CASO ANOMALO VISUALIZZATO in Peso~4370 g Lunghezza~315 mm
plot(Peso, Lunghezza, pch=20, xlab="Peso (g)", ylab="Lunghezza (mm)", main="Peso vs Lunghezza")
# in basso a destra
x_out <- 4370
y_out <- 315
symbols(x_out, y_out,circles = 100,inches = FALSE,add = TRUE,fg = "red",lwd = 2)
text(x_out, y_out, labels = "Leverage", pos = 3, col = "red")
Il modello e l’interpretazione dei coefficienti sono affidabili,
specialmente alla luce della dimensione campionaria!
3) Previsioni e Risultati
Una volta validato il modello, lo useremo per fare previsioni pratiche. Ad esempio, potremo stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.
# Dall'esempio ora capisco che probabilmente il numero di gravidanze è un parametro importante per il fit.
L39 <- median(Lunghezza[Gestazione == 39])
C39 <- median(Cranio[Gestazione == 39])
#Predizione: N.gravidanze=3,Sesso=F,Gestazione=39,Lunghezza=500 mm,Cranio=340 mm
caso_test <- data.frame(N.gravidanze=3,Sesso="F",Gestazione=39,Lunghezza=L39,Cranio=C39)
#PESO (g):
predict(bestmodel, newdata=caso_test, interval="prediction", level=0.95)## fit lwr upr
## 1 3324.878 2800.847 3848.9094) Visualizzazioni
Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili. Ad esempio, potremmo visualizzare l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the visreg package.
## Please report the issue at <https://github.com/pbreheny/visreg/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggplot2 package.
## Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Conclusioni
Il progetto di previsione del peso neonatale è un’iniziativa fondamentale per Neonatal Health Solutions. Attraverso l’uso di dati clinici dettagliati e strumenti di analisi statistica avanzati, possiamo contribuire a migliorare la qualità della cura prenatale, ridurre i rischi per i neonati e ottimizzare l’efficienza delle risorse ospedaliere. Questo progetto rappresenta un punto di svolta per l’azienda, consentendo non solo un miglioramento della pratica clinica ma anche l’implementazione di politiche sanitarie più informate e proattive.
Conclusioni: la relazione che più permette di predire il peso del bambino è una relazione non lineare in Lunghezza Cranio e Gestazione, come si evince dai rispettivi trend nei grafici, e con interazione tra dimensioni Cranio:Lunghezza e Cranio:Gestazione che aiutano a migliorare il fit N-dimensionale del modello. Il Sesso e il numero di gravidanze sono fondamentali per comprendere il peso effettivo. Rimane poco chiaro come il numero di gravidanze sia legato effettivamente al peso direttamente o indirettamente. Nella matrice di correlazione il fit appare orizzontale/piatto. Si può ipotizzare che un approccio della madre diverso durante la gravitanza possa contribuire ad un diverso nutrimento e acquisizione del peso da parte del nenoato, ma necessita ulteriori studi.
Miglior modello di dipendenza dal Peso: Peso ~ N.gravidanze + Sesso + I(Lunghezza^2) + I(Cranio^2) +I(Gestazione^2) + Cranio:Lunghezza + Cranio:Gestazione
Le distribuzioni di queste variabili non seguono una distribuzione normale, e tanto meno i residui che se ne ricavano dai fit. Tuttavia, questo solo in una piccola percentuale distribuita nelle ali.