Modello Statistico per la Previsione del Peso Neonatale
Domenico Toscano
2025-06-21
Obiettivo: Creare un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali. Il progetto mira a migliorare la gestione delle gravidanze ad alto rischio, ottimizzare le risorse ospedaliere e garantire migliori risultati per la salute neonatale. Il progetto si inserisce all’interno di un contesto di crescente attenzione verso la prevenzione delle complicazioni neonatali. La possibilità di prevedere il peso alla nascita dei neonati rappresenta un’opportunità fondamentale per migliorare la pianificazione clinica e ridurre i rischi associati a nascite problematiche, come parti prematuri o neonati con basso peso. Di seguito, i principali benefici che questo progetto porterà all’azienda e al settore sanitario.
Raccolta dei Dati
dati <- read.csv("neonati.csv", stringsAsFactors = TRUE, sep=",")
attach(dati)
kable(summary(dati))| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| Min. : 0.00 | Min. : 0.0000 | Min. :0.0000 | Min. :25.00 | Min. : 830 | Min. :310.0 | Min. :235 | Ces: 728 | osp1:816 | F:1256 | |
| 1st Qu.:25.00 | 1st Qu.: 0.0000 | 1st Qu.:0.0000 | 1st Qu.:38.00 | 1st Qu.:2990 | 1st Qu.:480.0 | 1st Qu.:330 | Nat:1772 | osp2:849 | M:1244 | |
| Median :28.00 | Median : 1.0000 | Median :0.0000 | Median :39.00 | Median :3300 | Median :500.0 | Median :340 | NA | osp3:835 | NA | |
| Mean :28.16 | Mean : 0.9812 | Mean :0.0416 | Mean :38.98 | Mean :3284 | Mean :494.7 | Mean :340 | NA | NA | NA | |
| 3rd Qu.:32.00 | 3rd Qu.: 1.0000 | 3rd Qu.:0.0000 | 3rd Qu.:40.00 | 3rd Qu.:3620 | 3rd Qu.:510.0 | 3rd Qu.:350 | NA | NA | NA | |
| Max. :46.00 | Max. :12.0000 | Max. :1.0000 | Max. :43.00 | Max. :4930 | Max. :565.0 | Max. :390 | NA | NA | NA |
Definizione Statistica delle Variabili
Il dataset contiene informazioni su 2500 neonati e include le seguenti variabili:
Variabili Quantitative:
- Anni.madre: Età della madre in anni (variabile quantitativa continua)
- N.gravidanze: Numero di gravidanze precedenti
(variabile quantitativa discreta)
- Gestazione: Durata della gravidanza in settimane (variabile quantitativa continua)
- Peso: Peso del neonato alla nascita in grammi (variabile dipendente - quantitativa continua)
- Lunghezza: Lunghezza del neonato in centimetri (variabile quantitativa continua)
- Cranio: Diametro del cranio del neonato in centimetri (variabile quantitativa continua)
Variabili Qualitative:
- Fumatrici: Indicatore binario del fumo materno (0=non fumatrice, 1=fumatrice) (variabile qualitativa nominale dicotomica)
- Tipo.parto: Modalità del parto (naturale/cesareo) (variabile qualitativa nominale dicotomica)
- Ospedale: Ospedale di nascita (1, 2, 3) (variabile qualitativa nominale politomica)
- Sesso: Sesso del neonato (M/F) (variabile qualitativa nominale dicotomica)
Calcolo Indici di Posizione
calcola_moda <- function(x) {
freq <- table(x)
max_freq <- max(freq)
moda <- as.numeric(names(freq[freq == max_freq]))
return(moda)
}
metrics <- list(
"Anni Madre" = Anni.madre,
"N gravidanze" = N.gravidanze,
"Gestazione [settimane]" = Gestazione,
"Peso [gr]" = Peso,
"Lunghezza [cm]" = Lunghezza,
"Cranio [cm]" = Cranio
)
indici.posizione <- data.frame(
Misure = c("Min", "Moda", "Mediana", "Media","Max")
)
for (name in names(metrics)) {
data <- metrics[[name]]
indici.posizione[[name]] <- c(
formatC(min(data), format = "f", digits = 2),
formatC(calcola_moda(data), format = "f", digits = 2),
formatC(median(data), format = "f", digits = 2),
formatC(mean(data), format = "f", digits = 2),
formatC(max(data), format = "f", digits =2)
)
}
kable(indici.posizione, caption = "Indici di Posizione")%>% kable_styling(position = "center")| Misure | Anni Madre | N gravidanze | Gestazione [settimane] | Peso [gr] | Lunghezza [cm] | Cranio [cm] |
|---|---|---|---|---|---|---|
| Min | 0.00 | 0.00 | 25.00 | 830.00 | 310.00 | 235.00 |
| Moda | 30.00 | 0.00 | 40.00 | 3300.00 | 500.00 | 340.00 |
| Mediana | 28.00 | 1.00 | 39.00 | 3300.00 | 500.00 | 340.00 |
| Media | 28.16 | 0.98 | 38.98 | 3284.08 | 494.69 | 340.03 |
| Max | 46.00 | 12.00 | 43.00 | 4930.00 | 565.00 | 390.00 |
Si può notare come il minimo per la variabile Anni Madre sia 0, questo è evidentemente un errore. Vengono quindi identificate quelle osservazioni per cui l’età della madre è inferiore a 13 anni e vengono eliminate dal dataset. Viene di seguito riportato il summary aggiornato degli indici di posizione.
dati<- dati[Anni.madre>=13,]
attach(dati)metrics <- list(
"Anni Madre" = Anni.madre,
"N gravidanze" = N.gravidanze,
"Gestazione [settimane]" = Gestazione,
"Peso [gr]" = Peso,
"Lunghezza [cm]" = Lunghezza,
"Cranio [cm]" = Cranio
)
indici.posizione <- data.frame(
Misure = c("Min", "Moda", "Mediana", "Media","Max")
)
for (name in names(metrics)) {
data <- metrics[[name]]
indici.posizione[[name]] <- c(
formatC(min(data), format = "f", digits = 2),
formatC(calcola_moda(data), format = "f", digits = 2),
formatC(median(data), format = "f", digits = 2),
formatC(mean(data), format = "f", digits = 2),
formatC(max(data), format = "f", digits =2)
)
}
kable(indici.posizione, caption = "Indici di Posizione")%>% kable_styling(position = "center")| Misure | Anni Madre | N gravidanze | Gestazione [settimane] | Peso [gr] | Lunghezza [cm] | Cranio [cm] |
|---|---|---|---|---|---|---|
| Min | 13.00 | 0.00 | 25.00 | 830.00 | 310.00 | 235.00 |
| Moda | 30.00 | 0.00 | 40.00 | 3300.00 | 500.00 | 340.00 |
| Mediana | 28.00 | 1.00 | 39.00 | 3300.00 | 500.00 | 340.00 |
| Media | 28.19 | 0.98 | 38.98 | 3284.18 | 494.70 | 340.03 |
| Max | 46.00 | 12.00 | 43.00 | 4930.00 | 565.00 | 390.00 |
variabili_numeriche <- c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")
dati_numerici <- dati[, variabili_numeriche]
cat("Summary delle variabili numeriche:\n")## Summary delle variabili numeriche:kable(summary(dati_numerici))| Anni.madre | N.gravidanze | Gestazione | Peso | Lunghezza | Cranio | |
|---|---|---|---|---|---|---|
| Min. :13.00 | Min. : 0.0000 | Min. :25.00 | Min. : 830 | Min. :310.0 | Min. :235 | |
| 1st Qu.:25.00 | 1st Qu.: 0.0000 | 1st Qu.:38.00 | 1st Qu.:2990 | 1st Qu.:480.0 | 1st Qu.:330 | |
| Median :28.00 | Median : 1.0000 | Median :39.00 | Median :3300 | Median :500.0 | Median :340 | |
| Mean :28.19 | Mean : 0.9816 | Mean :38.98 | Mean :3284 | Mean :494.7 | Mean :340 | |
| 3rd Qu.:32.00 | 3rd Qu.: 1.0000 | 3rd Qu.:40.00 | 3rd Qu.:3620 | 3rd Qu.:510.0 | 3rd Qu.:350 | |
| Max. :46.00 | Max. :12.0000 | Max. :43.00 | Max. :4930 | Max. :565.0 | Max. :390 |
CV <- function(x) sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE) * 100
stats_df <- data.frame(
Variabile = variabili_numeriche,
Media = round(sapply(dati_numerici, mean, na.rm = TRUE), 2),
Mediana = round(sapply(dati_numerici, median, na.rm = TRUE), 2),
Varianza = round(sapply(dati_numerici, var, na.rm = TRUE), 2),
Dev_Standard = round(sapply(dati_numerici, sd, na.rm = TRUE), 2),
Coeff_Variazione = round(sapply(dati_numerici, CV), 2),
Asimmetria = round(sapply(dati_numerici, skewness, na.rm = TRUE), 2)
)
kable(stats_df)| Variabile | Media | Mediana | Varianza | Dev_Standard | Coeff_Variazione | Asimmetria | |
|---|---|---|---|---|---|---|---|
| Anni.madre | Anni.madre | 28.19 | 28 | 27.22 | 5.22 | 18.51 | 0.15 |
| N.gravidanze | N.gravidanze | 0.98 | 1 | 1.64 | 1.28 | 130.50 | 2.51 |
| Gestazione | Gestazione | 38.98 | 39 | 3.49 | 1.87 | 4.79 | -2.07 |
| Peso | Peso | 3284.18 | 3300 | 275865.90 | 525.23 | 15.99 | -0.65 |
| Lunghezza | Lunghezza | 494.70 | 500 | 693.21 | 26.33 | 5.32 | -1.51 |
| Cranio | Cranio | 340.03 | 340 | 269.93 | 16.43 | 4.83 | -0.79 |
Interpretazione delle Statistiche Descrittive
Età della madre: - Media di 28.2 anni con distribuzione quasi simmetrica (asimmetria = 0.15) - Coefficiente di variazione del 18.5% indica moderata variabilità
Numero di gravidanze: - Media di 1.0 gravidanze con forte asimmetria positiva (2.51) - Alto coefficiente di variazione (130%) indica distribuzione molto concentrata sui valori bassi
Durata gestazione: - Media di 39.0 settimane con asimmetria negativa (-2.07) - Indica presenza di parti prematuri che creano una coda verso i valori bassi
Peso del neonato (variabile dipendente): - Media di 3284g con leggera asimmetria negativa (-0.65) - Distribuzione relativamente simmetrica con coefficiente di variazione di circa 16%
Lunghezza del neonato e diametro cranio: - Entrambe mostrano distribuzione quasi normale con bassa variabilità (CV < 6%)
Dai boxplot si può notare la presenza di diversi outlier, particolarmente evidenti per il numero di gravidanze e il peso del neonato, la distribuzione asimmetrica del numero di gravidanze e la relativa simmetria dell’età della madre, che confermano quanto valutato attraverso le statistiche descrittive.
Distribuzioni di Frequenza Variabili Qualitative
variabili_categoriche <- c("Fumatrici", "Tipo.parto", "Ospedale", "Sesso")
for (var in variabili_categoriche) {
freq_ass <- table(dati[[var]])
freq_rel <- prop.table(freq_ass)
distr_freq <- cbind(
Frequenza_Assoluta = freq_ass,
Frequenza_Relativa = round(freq_rel, 3),
Percentuale = paste0(round(freq_rel * 100, 1), "%")
)
cat("\n")
print(kable(distr_freq, caption = paste("Distribuzione di Frequenza -", var)))
}##
##
##
## Table: Distribuzione di Frequenza - Fumatrici
##
## | |Frequenza_Assoluta |Frequenza_Relativa |Percentuale |
## |:--|:------------------|:------------------|:-----------|
## |0 |2394 |0.958 |95.8% |
## |1 |104 |0.042 |4.2% |
##
##
##
## Table: Distribuzione di Frequenza - Tipo.parto
##
## | |Frequenza_Assoluta |Frequenza_Relativa |Percentuale |
## |:---|:------------------|:------------------|:-----------|
## |Ces |728 |0.291 |29.1% |
## |Nat |1770 |0.709 |70.9% |
##
##
##
## Table: Distribuzione di Frequenza - Ospedale
##
## | |Frequenza_Assoluta |Frequenza_Relativa |Percentuale |
## |:----|:------------------|:------------------|:-----------|
## |osp1 |816 |0.327 |32.7% |
## |osp2 |848 |0.339 |33.9% |
## |osp3 |834 |0.334 |33.4% |
##
##
##
## Table: Distribuzione di Frequenza - Sesso
##
## | |Frequenza_Assoluta |Frequenza_Relativa |Percentuale |
## |:--|:------------------|:------------------|:-----------|
## |F |1255 |0.502 |50.2% |
## |M |1243 |0.498 |49.8% |Da queste analisi emerge che: solo il 4% delle madri sono fumatrici (2498), per quanto riguarda il tipo di parto, Il 71% dei parti sono naturali e Il 29% sono cesarei. La distribuzione è quasi uniforme tra i tre ospedali (33-34% ciascuno) I grafici a barre confermano visivamente queste distribuzioni di frequenza.
Analisi variabile Peso
La variabile dipendente è il peso, in prima analisi facciamo un test di Shapiro-Wilk per saggiare l’ipotesi di normalità.
shapiro_result <- shapiro.test(Peso)
cat("Test di Shapiro-Wilk per normalità del Peso:\n")## Test di Shapiro-Wilk per normalità del Peso:cat("W =", round(shapiro_result$statistic, 4), ", p-value =", format(shapiro_result$p.value, scientific = TRUE), "\n")## W = 0.9707 , p-value = 3.378772e-22Il risultato è di 0.97, con un p-value prossimo a 0, quindi si rifiuta l’ipotesi nulla di normalità della variabile Peso.
Per avere un riscontro, visualizziamo la distribuzione di densità della variabile peso.
par(mfrow=c(1,2))
plot(density(Peso), main="Distribuzione della Densità del Peso",
xlab="Peso (grammi)", ylab="Densità")
hist(Peso, freq=FALSE, main="Istogramma del Peso",
xlab="Peso (grammi)", ylab="Densità", col="lightblue")
lines(density(Peso), col="red", lwd=2)
La distribuzione del peso mostrata nel grafico appare infatti approssimativamente normale (a campana), anche se si può notare una leggera asimmetria negativa (coda più lunga a sinistra), che è coerente con il coefficiente di asimmetria -0.65 riportato nelle statistiche descrittive. Inoltre c’è da considerare che la dimensione del campione è grande (N circa 2500). Con campioni così grandi, il test di Shapiro-Wilk diventa estremamente sensibile anche a piccole deviazioni dalla normalità. Questo può portare al rifiuto dell’ipotesi nulla anche quando la distribuzione è “abbastanza” normale per scopi pratici.
Matrice di Correlazione
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- (cor(x, y))
txt <- format(c(r, 1), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = 1.5)
}
colonne_interesse <- c("Peso", "Anni.madre", "N.gravidanze", "Gestazione", "Lunghezza", "Cranio")
dati_select <- dati[, colonne_interesse]
pairs(dati_select,lower.panel=panel.cor, upper.panel=panel.smooth)
Dall’analisi della matrice si mostra come il peso abbia una correlazione lineare positiva con le variabili Gestazioni, quindi il numero di settimane della gestazione, con un coefficiente di 0.59, con la lunghezza del neonato (R2 = 0.8) e con il diametro del cranio con un R2 di 0.7. Una leggera correlazione positiva con un R2 di 0.24 è valutata anche con il sesso. Per le restanti variabili i coefficienti di correlazioni sono molto prossimi a zero.
Per le variabili di tipo qualitativo per avere un informazione più di dettaglio costruiamo dei boxplot.
par(mfrow = c(2, 2))
boxplot(Peso, main = "Distribuzione generale del Peso")
boxplot(Peso ~ Sesso, main = "Peso per Sesso")
boxplot(Peso ~ Fumatrici, main = "Peso per Fumatrici")
boxplot(Peso ~ Tipo.parto, main = "Peso per Tipo di Parto")
Dai boxplot osserviamo che: 1. Il peso dei neonati maschi è mediamente più alto rispetto alle femmine. 2. Il fumo materno non sembra avere un impatto significativo sulla distribuzione del peso. 3. Non emergono differenze rilevanti tra i neonati nati con parto cesareo e quelli nati con parto naturale.
Per confermare queste osservazioni, applichiamo un t-test per confrontare le medie.
t.test(Peso ~ Sesso)##
## 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.496t.test(Peso ~ Fumatrici)##
## Welch Two Sample t-test
##
## data: Peso by Fumatrici
## t = 1.0362, df = 114.12, p-value = 0.3023
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -45.5076 145.3399
## sample estimates:
## mean in group 0 mean in group 1
## 3286.262 3236.346t.test(Peso ~ Tipo.parto)##
## Welch Two Sample t-test
##
## data: Peso by Tipo.parto
## t = -0.13626, df = 1494.4, p-value = 0.8916
## alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
## 95 percent confidence interval:
## -46.44246 40.40931
## sample estimates:
## mean in group Ces mean in group Nat
## 3282.047 3285.063chisq.test(Ospedale, Peso)## Warning in chisq.test(Ospedale, Peso): L'approssimazione al Chi-quadrato
## potrebbe essere inesatta##
## Pearson's Chi-squared test
##
## data: Ospedale and Peso
## X-squared = 560.03, df = 572, p-value = 0.6318- Sesso: I maschi pesano significativamente di più delle femmine (+247g, p < 0.001)
- Fumo: Nessuna differenza significativa (p = 0.30) - risultato inaspettato
- Tipo di parto: Nessuna differenza significativa (p = 0.89) - come atteso
- Ospedale: Nessuna associazione (p = 0.63)
Costruzione del Modello di Regressione
Prima della modellazione statistica, escludiamo a priori variabili che non hanno senso predittivo:
dati_modello <- dati[, !names(dati) %in% c("Tipo.parto", "Ospedale")]
cat("\nVariabili incluse nel modello:\n")##
## Variabili incluse nel modello:cat(paste(names(dati_modello), collapse = ", "), "\n")## Anni.madre, N.gravidanze, Fumatrici, Gestazione, Peso, Lunghezza, Cranio, Sessomod1 <- lm(Peso~.-Tipo.parto-Ospedale, data= dati )
summary(mod1)##
## Call:
## lm(formula = Peso ~ . - Tipo.parto - Ospedale, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1160.6 -181.3 -15.7 163.6 2630.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6712.2405 141.3339 -47.492 < 2e-16 ***
## Anni.madre 0.8803 1.1491 0.766 0.444
## N.gravidanze 11.3789 4.6767 2.433 0.015 *
## Fumatrici -30.3958 27.6080 -1.101 0.271
## Gestazione 32.9472 3.8288 8.605 < 2e-16 ***
## Lunghezza 10.2316 0.3011 33.979 < 2e-16 ***
## Cranio 10.5198 0.4271 24.633 < 2e-16 ***
## SessoM 78.0787 11.2132 6.963 4.24e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2490 degrees of freedom
## Multiple R-squared: 0.7272, Adjusted R-squared: 0.7264
## F-statistic: 948.3 on 7 and 2490 DF, p-value: < 2.2e-16Da questo modello sembra che le variabili Anni madre e Fumatrici non abbiano impatto sulla variabile output. La variabile del fumo in gravidanze è una delle variabili di interesse nello studio, quindi per il momento rimuovo solo la variabile anni madre e calcolo un nuovo modello.
mod2 <- update(mod1, ~.-Anni.madre)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1150.24 -181.32 -15.73 162.95 2635.69
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6682.2637 135.7983 -49.207 < 2e-16 ***
## N.gravidanze 12.6996 4.3470 2.921 0.00352 **
## Fumatrici -30.5728 27.6048 -1.108 0.26818
## Gestazione 32.6437 3.8079 8.573 < 2e-16 ***
## Lunghezza 10.2309 0.3011 33.979 < 2e-16 ***
## Cranio 10.5366 0.4265 24.707 < 2e-16 ***
## SessoM 78.1596 11.2117 6.971 4.01e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2491 degrees of freedom
## Multiple R-squared: 0.7271, Adjusted R-squared: 0.7265
## F-statistic: 1106 on 6 and 2491 DF, p-value: < 2.2e-16La rimozione della variabile non influisce sul R2. La variabile Fumatrici continua a non avere un effetto influente, per fare un ulteriore prova prima della rimozione della variabile dal modello provo a considerare l’effetto simultaneo di Fumatrici e Gestazione.
mod3 <- update(mod2, ~.+Fumatrici:Gestazione)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso + Fumatrici:Gestazione, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.83 -181.58 -16.95 163.76 2635.68
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6699.6947 136.7290 -49.000 < 2e-16 ***
## N.gravidanze 12.7452 4.3470 2.932 0.0034 **
## Fumatrici 795.7016 757.5315 1.050 0.2936
## Gestazione 33.2007 3.8418 8.642 < 2e-16 ***
## Lunghezza 10.2252 0.3011 33.957 < 2e-16 ***
## Cranio 10.5313 0.4265 24.693 < 2e-16 ***
## SessoM 78.7436 11.2241 7.016 2.94e-12 ***
## Fumatrici:Gestazione -21.0469 19.2830 -1.091 0.2752
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2490 degrees of freedom
## Multiple R-squared: 0.7273, Adjusted R-squared: 0.7265
## F-statistic: 948.6 on 7 and 2490 DF, p-value: < 2.2e-16L’R2 non subisce variazioni, ma entrambe le variabili che considerano l’effetto del fumo non risultano significative, questo potrebbe dipendere dal numero di osservazioni relative a madri fumatrici molto inferiore rispetto a quello delle madri non fumatrici. Si rimuovono quindi entrambe le variabili.
mod4 <- update(mod3, ~.-Fumatrici-Fumatrici:Gestazione)
summary(mod4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## 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-16L’R2 non subisce variazioni e in questo modello tutti regressori sono significativi. Per la variabile gestazione potrebbe aver senso considerare il suo effetto quadratico, come evidenziato dalla matrice di correlazione.
mod5 <- update(mod4, ~.+I(Gestazione^2))
summary(mod5)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2), data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1144.0 -181.5 -12.9 165.8 2661.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4646.7158 898.6322 -5.171 2.52e-07 ***
## N.gravidanze 12.5489 4.3381 2.893 0.00385 **
## Gestazione -81.2309 49.7402 -1.633 0.10257
## Lunghezza 10.3502 0.3040 34.045 < 2e-16 ***
## Cranio 10.6376 0.4282 24.843 < 2e-16 ***
## SessoM 75.7563 11.2435 6.738 1.99e-11 ***
## I(Gestazione^2) 1.5168 0.6621 2.291 0.02206 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.5 on 2491 degrees of freedom
## Multiple R-squared: 0.7276, Adjusted R-squared: 0.7269
## F-statistic: 1109 on 6 and 2491 DF, p-value: < 2.2e-16Si vede un aumento di R2 trascurabile e la variabile Gestazione perde la sua influenza sul modello, l’effetto quadratico di Gestazione non viene quindi considerato. Anche la variabile Cranio potrebbe avere un effetto quadratico, provo quindi a considerarla nel modello.
mod6 <- update(mod5, ~.-I(Gestazione^2)+I(Cranio^2))
summary(mod6)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso + I(Cranio^2), data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1138.6 -179.4 -14.8 163.4 2622.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.10118 1151.77280 0.073 0.94180
## N.gravidanze 12.76356 4.31259 2.960 0.00311 **
## Gestazione 38.90540 3.93291 9.892 < 2e-16 ***
## Lunghezza 10.48745 0.30157 34.776 < 2e-16 ***
## Cranio -31.79371 7.16973 -4.434 9.63e-06 ***
## SessoM 73.10236 11.16590 6.547 7.11e-11 ***
## I(Cranio^2) 0.06262 0.01059 5.915 3.77e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 272.8 on 2491 degrees of freedom
## Multiple R-squared: 0.7308, Adjusted R-squared: 0.7301
## F-statistic: 1127 on 6 and 2491 DF, p-value: < 2.2e-16Selezione del Modello Finale
Procedo ora alla selezione del modello migliore tra i sei modelli creati. Utilizzo i criteri BIC e AIC per scegliere il modello migliore.
AIC(mod1,mod2,mod3,mod4,mod5,mod6)## df AIC
## mod1 9 35155.07
## mod2 8 35153.66
## mod3 9 35154.46
## mod4 7 35152.89
## mod5 8 35149.63
## mod6 8 35120.04BIC(mod1,mod2,mod3,mod4,mod5,mod6)## df BIC
## mod1 9 35207.48
## mod2 8 35200.24
## mod3 9 35206.87
## mod4 7 35193.65
## mod5 8 35196.21
## mod6 8 35166.63Per entrambi i criteri il modello 6 è quello migliore (AIC e BIC minore), seguito dal modello 4. Gli R2 dei due modelli differivano di un meno di un punto percentuale, utilizzo il test anova per verificare se c’è stato un aumento significativo della varianza spiegata con l’aggiunta di una variabile nel modello 6.
anova(mod6,mod4)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso +
## I(Cranio^2)
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2491 185437522
## 2 2492 188042054 -1 -2604532 34.987 3.775e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1C’è un aumento significativo di varianza spiegata. Verifico ora la multicollinearità tra le variabili.
vif(mod6)## N.gravidanze Gestazione Lunghezza Cranio Sesso I(Cranio^2)
## 1.023611 1.812253 2.114666 465.421862 1.045890 452.902551vif(mod4)## N.gravidanze Gestazione Lunghezza Cranio Sesso
## 1.023462 1.669779 2.075747 1.624568 1.040184In presenza del termine quadratico i VIF aumentano leggermente per le variabili comuni a entrambi i modelli, ma hanno un valore >> 5 per la variabile cranio e il suo termine quadratico. Il modello 6 ha un aumento di meno di un punto percentuale su R2, quindi viene preferito il più semplice (modello 4).
Diagnostica del Modello
Analisi dei Residui
mod_finale = mod4
par(mfrow = c(2,2))
plot(mod_finale, main = "Diagnostica Modello")
cat("\n=== TEST DIAGNOSTICI ===\n")##
## === TEST DIAGNOSTICI ===shapiro_res <- shapiro.test(residuals(mod_finale))
cat("Test Shapiro-Wilk (normalità residui):\n")## Test Shapiro-Wilk (normalità residui):cat("W =", round(shapiro_res$statistic, 4), ", p-value =", format(shapiro_res$p.value, scientific = TRUE), "\n")## W = 0.9741 , p-value = 7.318386e-21bp_test <- bptest(mod_finale)
cat("\nTest Breusch-Pagan (omoschedasticità):\n") ##
## Test Breusch-Pagan (omoschedasticità):cat("BP =", round(bp_test$statistic, 4), ", p-value =", format(bp_test$p.value, scientific = TRUE), "\n")## BP = 90.2966 , p-value = 5.821337e-18dw_test <- dwtest(mod_finale)
cat("\nTest Durbin-Watson (autocorrelazione):\n")##
## Test Durbin-Watson (autocorrelazione):cat("DW =", round(dw_test$statistic, 4), ", p-value =", format(dw_test$p.value, scientific = TRUE), "\n")## DW = 1.9532 , p-value = 1.209109e-01par(mfrow = c(1,1))
plot(density(residuals(mod_finale)), main="Distribuzione dei Residui",
xlab="Residui", ylab="Densità")
Il test di Breusch-Pagan rifiuta l’ipotesi nulla quindi si rifiuta l’ipotesi nulla di omoschedasticità, questo potrebbe influenzare l’efficienza delle stime dei coefficienti. Non viene invece rifiutata l’ipotesi nulla di Durbin-Watson, quindi i residui non sono autocorrelati. Anche il test di Shapiro-Wilk rifiuta l’ipotesi nulla di normalità, ma dalla visualizzazione grafica la distribuzione dei residui si avvicina a quella normale.
Identificazione punti di leva e outliers Influenti
lev <- hatvalues(mod_finale)
cook_dist <- cooks.distance(mod_finale)
stud_res <- rstudent(mod_finale)
soglia_lev <- 2 * mean(lev)
soglia_cook <- 4 / nrow(dati_modello)
par(mfrow = c(2,2))
plot(lev, main="Leverage", ylab="Hat Values")
abline(h = soglia_lev, col="red", lty=2)
high_lev <- which(lev > soglia_lev)
plot(stud_res, main="Residui Studentizzati", ylab="Residui Studentizzati")
abline(h = c(-2, 2), col="blue", lty=2)
outliers <- which(abs(stud_res) > 2)
plot(cook_dist, main="Distanza di Cook", ylab="Cook's Distance")
abline(h = soglia_cook, col="red", lty=2)
influential <- which(cook_dist > soglia_cook)
Previsioni con il Modello
Esempio di Previsione
Facciamo una previsione considerando questi dati:
Sesso = F Numero gravidanze = 2 (terza gravidanza) Settimane Gestazione = 39 Fumatrice = NO Non essendoci forniti dati per la lunghezza o il diametro del cranio si usano come dati le medie di lunghezza e diametro per le neonate.
dati_lunghezza <- dati %>%
group_by(Sesso) %>%
summarise(LunghezzaMedia = mean(Lunghezza, na.rm = TRUE))
dati_diametro <- dati %>%
group_by(Sesso) %>%
summarise(DiametroMedio = mean(Cranio, na.rm = TRUE))
new_data <- data.frame(N.gravidanze=2,Gestazione=39,Lunghezza=489.7641, Cranio=339.6231,Sesso="F")
peso_new <- predict(mod_finale,new_data)Il modello prevede un peso della neonata di 3203.92 g.
Visualizzazioni Finali
Visualizziamo graficamente i risultati del modello per mostrare le relazioni tra le variabili. Viene di seguito mostrata la relazione tra settimane di gestazione e peso del neonato, tenendo in considerazione anche il sesso.
ggplot(data = dati) +
geom_point(aes(x = Gestazione, y = Peso, col = Sesso), position = "jitter") +
geom_smooth(aes(x = Gestazione, y = Peso, col = Sesso), se = FALSE, method = "lm") +
labs(
title = "Relazione tra Settimane di Gestazione e Peso del Neonato",
x = "Settimane di Gestazione",
y = "Peso del Neonato (g)",
color = "Sesso"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title.x = element_text(size = 12, face = "bold"),
axis.title.y = element_text(size = 12, face = "bold"),
legend.title = element_text(size = 12, face = "bold"),
legend.position = "right"
)
Si può notare come la retta di regressione per maschi e femmine abbia la stessa tendenza crescente positiva, ma che i pesi delle neonate siano inferiori a quelli dei neonati. Le osservazioni per una durata minore della gestazione si trovano al di sotto delle retta di regressione quindi il modello potrebbe sovrastimare il peso del neonato quando questo nasce prematuramente. La nuvola di punti sembra invece più concentrata sulla retta di regressioni per gravidanze più lunghe. Indaghiamo ora la relazione tra il numero di gravidanze e il peso del neonato.
ggplot(data = dati) +
geom_point(aes(x = N.gravidanze, y = Peso, col = Sesso), position = "jitter") +
geom_smooth(aes(x = N.gravidanze, y = Peso, col = Sesso), se = FALSE, method = "lm") +
labs(
title = "Relazione tra Numero di Gravidanze e Peso del Neonato",
x = "Numero di gravidanze",
y = "Peso del Neonato (g)",
color = "Sesso"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title.x = element_text(size = 12, face = "bold"),
axis.title.y = element_text(size = 12, face = "bold"),
legend.title = element_text(size = 12, face = "bold"),
legend.position = "right"
)## `geom_smooth()` using formula = 'y ~ x'
Il numero di gravidanze non sembra avere influenza sul peso del neonato. Ancora una volta il peso delle neonate risulta essere inferiore a quello dei neonati. Visualizziamo la relazione tra lunghezza e peso.
ggplot(data = dati) +
geom_point(aes(x = Lunghezza, y = Peso, col = Sesso), position = "jitter") +
geom_smooth(aes(x = Lunghezza, y = Peso, col = Sesso), se = FALSE, method = "lm") +
labs(
title = "Relazione tra Lunghezza e Peso del Neonato",
x = "Lunghezza (cm)",
y = "Peso del Neonato (g)",
color = "Sesso"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title.x = element_text(size = 12, face = "bold"),
axis.title.y = element_text(size = 12, face = "bold"),
legend.title = element_text(size = 12, face = "bold"),
legend.position = "right"
)## `geom_smooth()` using formula = 'y ~ x'
Anche in questo caso, come per le settimane di gestazione, si ha una dipendenza lineare positiva della variabile output per ambo i sessi. Si nota come sia meno visibile la differenza di peso tra maschi e femmine al variare della lunghezza. Analizziamo anche la relazione tra il diametro del cranio e il peso del neonato.
ggplot(data = dati) +
geom_point(aes(x = Cranio, y = Peso, col = Sesso), position = "jitter") +
geom_smooth(aes(x = Cranio, y = Peso, col = Sesso), se = FALSE, method = "lm") +
labs(
title = "Relazione tra Diametro del Cranio e Peso del Neonato",
x = "Diametro del cranio",
y = "Peso del Neonato (g)",
color = "Sesso"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
axis.title.x = element_text(size = 12, face = "bold"),
axis.title.y = element_text(size = 12, face = "bold"),
legend.title = element_text(size = 12, face = "bold"),
legend.position = "right"
)## `geom_smooth()` using formula = 'y ~ x'
Si può osservare in questo caso la stessa tendenza evidenziata per le settimane di gestazione. Si nota una crescita lineare, ma il modello sembra sovrastimare i pesi per diametri del cranio minori. Nuovamente si nota la differenza in peso tra maschi e femmine.