Modello Statistico per la Previsione del Peso Neonatale
Giulia Calvino
2025-05-20
Analisi del Peso dei Neonati
Questo documento presenta l’analisi del peso dei neonati utilizzando R. L’obiettivo è costruire un modello di regressione lineare multipla per stimare il peso del neonato in base a diverse variabili materne e neonatali.
1. Raccolta Dati e Struttura del Dataset
Il dataset “neonati.csv” è stato caricato in R. Le prime righe del dataset e la sua struttura sono state esaminate per comprendere le variabili e i loro tipi.
setwd("C:/Users/giuli/OneDrive/Desktop/MasterAI") # Imposta la directory
dati <- read.csv("neonati.csv", stringsAsFactors = TRUE, sep = ",") # Carica i dati
kable(summary(dati)) # Riepilogo statistico del dataset| 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 |
str(dati) # Struttura del dataset## '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) # Rende le variabili direttamente accessibili
# Verifica asimmetria e curtosi della variabile Peso
moments::skewness(Peso)## [1] -0.6470308moments::kurtosis(Peso) - 3## [1] 2.031532n <- nrow(dati) # Numero di osservazioni
print(n)## [1] 25002. Analisi e Modellizzazione
2.1 Analisi Preliminare
Variabili Numeriche
Si definisce una funzione analisi_descrittiva() per calcolare statistiche descrittive (media, mediana, deviazione standard, quartili, asimmetria, curtosi) per le variabili numeriche.
analisi_descrittiva <- function(x, nome_var = "Variabile") {
# Calcolo delle statistiche
media <- mean(x, na.rm = TRUE)
mediana <- median(x, na.rm = TRUE)
deviazione <- sd(x, na.rm = TRUE)
quartili <- quantile(x, probs = c(0.25, 0.75), na.rm = TRUE)
asimmetria <- moments::skewness(x, na.rm = TRUE)
curtosi <- moments::kurtosis(x, na.rm = TRUE) - 3
# Creazione di un data frame con i risultati
risultati <- data.frame(
Variabile = nome_var,
Media = media,
Mediana = mediana,
DeviazioneStandard = deviazione,
Quartile_25 = unname(quartili[1]), # Usa unname per rimuovere i nomi
Quartile_75 = unname(quartili[2]), # Usa unname per rimuovere i nomi
Asimmetria = asimmetria,
Curtosi = curtosi
)
return(risultati)
}
# Calcolo delle statistiche per ogni variabile e combinazione dei risultati
risultati_peso <- analisi_descrittiva(dati$Peso, "Peso")
risultati_anni_madre <- analisi_descrittiva(dati$Anni.madre, "Anni.madre")
risultati_n_gravidanze <- analisi_descrittiva(dati$N.gravidanze, "Numero Gravidanze")
risultati_gestazione <- analisi_descrittiva(dati$Gestazione, "Gestazione")
risultati_lunghezza <- analisi_descrittiva(dati$Lunghezza, "Lunghezza")
risultati_cranio <- analisi_descrittiva(dati$Cranio, "Diametro Cranio")
# Combinazione di tutti i risultati in un unico data frame
analisi_completa <- rbind(risultati_peso, risultati_anni_madre, risultati_n_gravidanze,
risultati_gestazione, risultati_lunghezza, risultati_cranio)
# Kable per creare la tabella
kable(analisi_completa)| Variabile | Media | Mediana | DeviazioneStandard | Quartile_25 | Quartile_75 | Asimmetria | Curtosi |
|---|---|---|---|---|---|---|---|
| Peso | 3284.0808 | 3300 | 525.038744 | 2990 | 3620 | -0.6470308 | 2.0315318 |
| Anni.madre | 28.1640 | 28 | 5.273578 | 25 | 32 | 0.0428115 | 0.3804165 |
| Numero Gravidanze | 0.9812 | 1 | 1.280587 | 0 | 1 | 2.5142541 | 10.9894064 |
| Gestazione | 38.9804 | 39 | 1.868639 | 38 | 40 | -2.0653133 | 8.2581496 |
| Lunghezza | 494.6920 | 500 | 26.318644 | 480 | 510 | -1.5146991 | 6.4871739 |
| Diametro Cranio | 340.0292 | 340 | 16.425330 | 330 | 350 | -0.7850527 | 2.9462063 |
- Le variabili mostrano diverse caratteristiche di distribuzione, con alcune che presentano skewness e curtosi significative, suggerendo la presenza di outliers o distribuzioni non normali.
- La variabilità all’interno dei dati è ben rappresentata dalle deviazioni standard e dai quartili.
- Le analisi di skewness e curtosi possono essere utili per ulteriori indagini statistiche e per comprendere la forma delle distribuzioni delle variabili.
Questi risultati forniscono un quadro chiaro e dettagliato delle caratteristiche statistiche delle variabili analizzate, utile per interpretazioni più approfondite e per guidare ulteriori analisi dei dati.
Variabili Categoriche
analisi_categorica <- function(x, nome_var = "Variabile") {
# Frequenze assolute e relative
freq_ass <- table(x)
freq_rel <- prop.table(freq_ass)
# Creazione del dataframe per la tabella
freq_df <- data.frame(
Categoria = names(freq_ass),
Fr_Assoluta = as.integer(freq_ass),
Fr_Relativa = round(as.numeric(freq_rel), 3)
)
# Formattazione dei nomi delle colonne
colnames(freq_df) <- c(nome_var, "Fr. Assoluta", "Fr. Relativa")
# Stampare la tabella
kable(freq_df, caption = paste("Distribuzione della variabile", nome_var), align = "lcc", escape = FALSE)
}
analisi_categorica(dati$Sesso, "Sesso")| Sesso | Fr. Assoluta | Fr. Relativa |
|---|---|---|
| F | 1256 | 0.502 |
| M | 1244 | 0.498 |
analisi_categorica(dati$Fumatrici, "Fumatrici")| Fumatrici | Fr. Assoluta | Fr. Relativa |
|---|---|---|
| 0 | 2396 | 0.958 |
| 1 | 104 | 0.042 |
analisi_categorica(dati$Tipo.parto, "Tipo_Parto")| Tipo_Parto | Fr. Assoluta | Fr. Relativa |
|---|---|---|
| Ces | 728 | 0.291 |
| Nat | 1772 | 0.709 |
analisi_categorica(dati$Ospedale, "Ospedale")| Ospedale | Fr. Assoluta | Fr. Relativa |
|---|---|---|
| osp1 | 816 | 0.326 |
| osp2 | 849 | 0.340 |
| osp3 | 835 | 0.334 |
Test di Ipotesi
# Calcolo del test
tab_parto_ospedale <- table(dati$Tipo.parto, dati$Ospedale)
chi_test <- chisq.test(tab_parto_ospedale)
# Estrazione valori
x2_value <- round(chi_test$statistic, 3)
df_value <- chi_test$parameter
p_value <- format.pval(chi_test$p.value)
# Output in formato Markdown
cat(paste0("**Risultato Test del Chi-quadrato**\n\n",
"- Statistica Chi-quadrato: ", x2_value, "\n",
"- Gradi di libertà: ", df_value, "\n",
"- Valore p: ", p_value, "\n\n"))Risultato Test del Chi-quadrato
- Statistica Chi-quadrato: 1.097
- Gradi di libertà: 2
- Valore p: 0.57776
Il test del Chi-quadrato non ha evidenziato un’associazione significativa tra il tipo di parto e l’ospedale (p = 0.578). Pertanto, non ci sono differenze statisticamente significative nella distribuzione del tipo di parto tra i diversi ospedali.
Tramite ricerca sul web, è stato possibile trovare una statistica sulla media del peso e lunghezza neonatale che useremo come confronto con i nostri dati:
Peso medio della popolazione alla nascita: 3.300g con una dev.st. di 500g
Lunghezza media della popolazione alla nascita: 500mm con una dev.st. di 20mm
# Test t per Peso
pop_media_peso <- 3300
test_peso <- t.test(dati$Peso, mu = pop_media_peso, alternative = "two.sided")
# Test t per Lunghezza
pop_media_lunghezza <- 500
test_lunghezza <- t.test(dati$Lunghezza, mu = pop_media_lunghezza, alternative = "two.sided")
# Pulizia output
peso_mean <- round(test_peso$estimate, 1)
peso_ci <- round(test_peso$conf.int, 1)
peso_t <- round(test_peso$statistic, 3)
peso_p <- format.pval(test_peso$p.value, digits = 4, eps = 0.0001)
lung_mean <- round(test_lunghezza$estimate, 1)
lung_ci <- round(test_lunghezza$conf.int, 1)
lung_t <- round(test_lunghezza$statistic, 3)
lung_p <- format.pval(test_lunghezza$p.value, digits = 4, eps = 0.0001)
# Stampa in Markdown
cat("**Test t per confronto con la popolazione**\n\n")Test t per confronto con la popolazione
cat("Peso alla nascita:\n")Peso alla nascita:
cat(paste0("- Media campionaria: ", peso_mean, "\n",
"- Valore t: ", peso_t, "n\n",
"- Intervallo di confidenza al 95%: [", peso_ci[1], ", ", peso_ci[2], "]\n",
"- Valore p: ", peso_p, "\n\n"))- Media campionaria: 3284.1
- Valore t: -1.516n
- Intervallo di confidenza al 95%: [3263.5, 3304.7]
- Valore p: 0.1296
cat("Lunghezza alla nascita:\n")Lunghezza alla nascita:
cat(paste0("- Media campionaria: ", lung_mean, " mm\n",
"- Valore t: ", lung_t, "\n",
"- Intervallo di confidenza al 95%: [", lung_ci[1], ", ", lung_ci[2], "]\n",
"- Valore p: ", lung_p, "\n"))- Media campionaria: 494.7 mm
- Valore t: -10.084
- Intervallo di confidenza al 95%: [493.7, 495.7]
- Valore p: < 1e-04
La media del peso dei neonati (3284.1 g) non differisce significativamente dalla media teorica della popolazione, con p = 0.13. Non si rifiuta l’ipotesi nulla.
La lunghezza media dei neonati (494.7 mm) è significativamente inferiore rispetto alla media ipotizzata di 500 mm (p < 0.001). Si rifiuta l’ipotesi nulla.
# Test t tra sessi
test_peso_sesso <- t.test(Peso ~ Sesso, data = dati)
test_lung_sesso <- t.test(Lunghezza ~ Sesso, data = dati)
test_cranio_sesso <- t.test(Cranio ~ Sesso, data = dati)
# Estrazione e pulizia risultati
estrai_info <- function(test) {
list(
media_F = round(test$estimate["mean in group F"], 1),
media_M = round(test$estimate["mean in group M"], 1),
t = round(test$statistic, 3),
p = format.pval(test$p.value, digits = 4, eps = 0.0001),
ci = round(test$conf.int, 2)
)
}
peso <- estrai_info(test_peso_sesso)
lung <- estrai_info(test_lung_sesso)
cranio <- estrai_info(test_cranio_sesso)
# Stampa Markdown
cat("**Test t per confronto tra sessi**\n\n")Test t per confronto tra sessi
cat("Peso alla nascita:\n")Peso alla nascita:
cat(paste0("- Media femmine: ", peso$media_F, " g\n",
"- Media maschi: ", peso$media_M, " g\n",
"- Valore t: ", peso$t, "\n",
"- Intervallo di confidenza al 95%: [", peso$ci[1], ", ", peso$ci[2], "]\n",
"- Valore p: ", peso$p, "\n\n"))- Media femmine: 3161.1 g
- Media maschi: 3408.2 g
- Valore t: -12.106
- Intervallo di confidenza al 95%: [-287.11, -207.06]
- Valore p: < 1e-04
cat("Lunghezza alla nascita:\n")Lunghezza alla nascita:
cat(paste0("- Media femmine: ", lung$media_F, " mm\n",
"- Media maschi: ", lung$media_M, " mm\n",
"- Valore t: ", lung$t, "\n",
"- Intervallo di confidenza al 95%: [", lung$ci[1], ", ", lung$ci[2], "]\n",
"- Valore p: ", lung$p, "\n\n"))- Media femmine: 489.8 mm
- Media maschi: 499.7 mm
- Valore t: -9.582
- Intervallo di confidenza al 95%: [-11.93, -7.88]
- Valore p: < 1e-04
cat("Circonferenza cranica:\n")Circonferenza cranica:
cat(paste0("- Media femmine: ", cranio$media_F, " mm\n",
"- Media maschi: ", cranio$media_M, " mm\n",
"- Valore t: ", cranio$t, "\n",
"- Intervallo di confidenza al 95%: [", cranio$ci[1], ", ", cranio$ci[2], "]\n",
"- Valore p: ", cranio$p, "\n"))- Media femmine: 337.6 mm
- Media maschi: 342.4 mm
- Valore t: -7.41
- Intervallo di confidenza al 95%: [-6.09, -3.54]
- Valore p: < 1e-04
I maschi hanno un peso medio significativamente maggiore rispetto alle femmine (p < 0.001). Differenza media: circa 247 g.
Anche la lunghezza media è maggiore nei maschi rispetto alle femmine, con una differenza significativa (p < 0.001).
La circonferenza cranica media è maggiore nei maschi, con differenza significativa (p < 0.001).
Gestione dei Dati Mancanti
È stato verificato il numero di dati mancanti per ciascuna variabile. In questo caso, non è stato riscontrato nessun dato mancante.
kable(colSums(is.na(dati))) # Numero di dati mancanti per colonna| x | |
|---|---|
| Anni.madre | 0 |
| N.gravidanze | 0 |
| Fumatrici | 0 |
| Gestazione | 0 |
| Peso | 0 |
| Lunghezza | 0 |
| Cranio | 0 |
| Tipo.parto | 0 |
| Ospedale | 0 |
| Sesso | 0 |
2.2 Creazione del Modello di Regressione
È stato costruito un modello di regressione lineare multipla per stimare il peso del neonato.
Correlazione: È stata visualizzata la matrice di correlazione tra le variabili per identificare eventuali forti correlazioni.
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) {
par(usr = c(0, 1, 0, 1))
r <- 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, col = "red") # Aggiunto colore rosso
}
pairs(dati, lower.panel = panel.cor, upper.panel = panel.smooth)
Verifica di Normalità: È stato effettuato uno Shapiro-Wilk test per verificare la normalità della variabile Peso.
shapiro.test(Peso)##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16Il test ha evidenziato che la variabile non è normalmente distribuita.
Costruzione e Selezione del Modello: Sono stati creati e confrontati diversi modelli (mod1, mod2, mod3) utilizzando ANOVA e BIC per selezionare il modello più parsimonioso e con buona capacità predittiva. P.S. Qui ho deciso di lasciare l’output della console per semplicità di codice
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-16car::vif(mod1) # Verifica multicollinearità## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.187454 1 1.089704
## N.gravidanze 1.186428 1 1.089233
## Fumatrici 1.007392 1 1.003689
## Gestazione 1.695810 1 1.302233
## Lunghezza 2.085755 1 1.444214
## Cranio 1.630796 1 1.277026
## Tipo.parto 1.004242 1 1.002119
## Ospedale 1.004071 2 1.001016
## Sesso 1.040643 1 1.020119Il modello iniziale (mod1) è stato costruito per esaminare l’effetto di tutte le variabili (Anni.madre, N.gravidanze, Fumatrici, Gestazione, Lunghezza, Cranio, Tipo.parto, Ospedale, e Sesso) sul Peso del neonato. L’output di summary(mod1) fornisce informazioni dettagliate sulla significatività di ciascuna variabile e sulla bontà di adattamento del modello.
Variabili Significative: N.gravidanze, Gestazione, Lunghezza, Cranio, Tipo.partoNat, Ospedaleosp3, e SessoM mostrano valori p significativi (p < 0.05), indicando che queste variabili hanno un impatto statisticamente significativo sul Peso del neonato. In particolare, si osserva che: 1) Un aumento nel numero di gravidanze è associato a un aumento del peso del neonato (11.2665 unità). 2) La durata della gestazione e le dimensioni (Lunghezza e Cranio) hanno un forte effetto positivo sul peso. 3) I neonati nati con parto naturale hanno un peso diverso rispetto ad altri tipi di parto. 4) I neonati nati nell’ospedale 3 hanno un peso diverso rispetto all’ospedale di riferimento. 5) I neonati maschi hanno un peso diverso rispetto alle femmine.
Anni.madre, Fumatrici, e Ospedaleosp2 non risultano significativi, suggerendo che, in questo modello, non hanno un impatto statisticamente rilevante sul peso del neonato.
Bontà di Adattamento del Modello: R-squared: Il valore di R-squared (0.7289) indica che il modello spiega circa il 72.89% della variabilità nel Peso del neonato. L’Adjusted R-squared (0.7278) è simile, suggerendo che il modello non è eccessivamente complesso e che le variabili incluse contribuiscono effettivamente alla spiegazione della variabilità.
F-statistic: Il valore elevato dell’F-statistic e il p-value molto basso (p < 2.2e-16) confermano la significatività complessiva del modello.
Verifica della Multicollinearità: L’output di car::vif(mod1) mostra i valori del Variance Inflation Factor (VIF). Valori di VIF superiori a 5 o 10 sono generalmente considerati indicativi di multicollinearità, che può gonfiare gli errori standard dei coefficienti e rendere difficile l’interpretazione dei singoli effetti delle variabili. In questo caso, tutti i valori di GVIF sono inferiori a 2.1, suggerendo che la multicollinearità non è un problema grave nel modello.
2.3 Selezione del Modello Ottimale
mod2 <- update(mod1, ~ . -Anni.madre -Fumatrici -Ospedale)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1129.31 -181.70 -16.31 161.07 2638.85
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.2971 135.9911 -49.322 < 2e-16 ***
## N.gravidanze 12.7558 4.3366 2.941 0.0033 **
## Gestazione 32.2713 3.7941 8.506 < 2e-16 ***
## Lunghezza 10.2864 0.3007 34.207 < 2e-16 ***
## Cranio 10.5057 0.4260 24.659 < 2e-16 ***
## Tipo.partoNat 30.0342 12.0969 2.483 0.0131 *
## SessoM 77.9285 11.1905 6.964 4.22e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2493 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.727
## F-statistic: 1110 on 6 and 2493 DF, p-value: < 2.2e-16car::vif(mod2)## N.gravidanze Gestazione Lunghezza Cranio Tipo.parto Sesso
## 1.024171 1.669258 2.080039 1.626199 1.003438 1.040060#confronto tra mod1 e mod2
anova(mod1, mod2) ## Analysis of Variance Table
##
## Model 1: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2489 186762521
## 2 2493 187601677 -4 -839155 2.7959 0.02479 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1BIC(mod1, mod2)## df BIC
## mod1 12 35241.84
## mod2 8 35221.75Il modello mod2 è stato ottenuto da mod1 rimuovendo le variabili non significative identificate nell’analisi precedente: Anni.madre, Fumatrici e Ospedale. L’obiettivo era semplificare il modello mantenendo una buona capacità predittiva.
Variabili Significative: Tutte le variabili rimaste nel modello (N.gravidanze, Gestazione, Lunghezza, Cranio, Tipo.parto, e Sesso) sono altamente significative (p < 0.01 per la maggior parte). Questo conferma la loro importanza nel predire il peso del neonato. I segni e la grandezza degli effetti delle variabili significative rimangono sostanzialmente simili a quelli osservati in mod1, suggerendo che la rimozione delle variabili non significative non ha alterato drasticamente l’impatto delle variabili importanti.
Bontà di Adattamento del Modello: R-squared: L’R-squared è leggermente diminuito (da 0.7289 a 0.7277), ma l’Adjusted R-squared è solo leggermente inferiore (da 0.7278 a 0.727). Questo indica che la perdita di potere esplicativo dovuta alla rimozione delle variabili è minima, e il modello semplificato (mod2) è ancora in grado di spiegare una grande parte della variabilità nel peso del neonato.
F-statistic: L’F-statistic è aumentato (da 669.2 a 1110) e il p-value rimane estremamente basso, confermando la significatività complessiva del modello semplificato.
Verifica della Multicollinearità in mod2: L’output di car::vif(mod2) mostra che tutti i valori di VIF sono molto bassi (inferiori a 2.1). Questo conferma che la multicollinearità non è un problema nel modello semplificato. Anzi, la rimozione delle variabili ha potenzialmente ridotto ulteriormente qualsiasi rischio di multicollinearità.
Confronto tra mod1 e mod2 usando anova(mod1, mod2): L’analisi ANOVA confronta formalmente i due modelli. Il risultato mostra un p-value di 0.02479, che è inferiore a 0.05. Questo indica che c’è una differenza statisticamente significativa tra i due modelli. Tuttavia, la significatività è marginale, e la differenza nella somma dei quadrati degli errori (Sum of Sq) è relativamente piccola rispetto alla variabilità totale dei dati. Ciò suggerisce che, sebbene mod1 sia leggermente migliore in termini di adattamento puro, mod2 non è significativamente peggiore e offre una maggiore parsimonia. Confronto tra mod1 e mod2 usando BIC(mod1, mod2):
Il Bayesian Information Criterion (BIC) penalizza la complessità del modello. I risultati di BIC(mod1, mod2) mostrano un valore di BIC inferiore per mod2 (35221.75) rispetto a mod1 (35241.84). Questo è un forte indicatore che mod2 è preferibile, poiché il BIC favorisce i modelli più semplici che spiegano i dati in modo adeguato. In Sintesi:
Il confronto tra mod1 e mod2 suggerisce che la semplificazione del modello, rimuovendo le variabili non significative, è una scelta appropriata. Sebbene mod1 possa avere un leggermente migliore adattamento ai dati (come indicato dall’ANOVA), mod2 offre una parsimonia significativa con una perdita minima di potere predittivo. Il BIC conferma questa conclusione, supportando la selezione di mod2 come modello preferito.
mod3 <- update(mod2, ~ . -Tipo.parto)
summary(mod3)##
## 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-16anova(mod2, mod3) #confronto tra mod2 e mod3## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2493 187601677
## 2 2494 188065546 -1 -463870 6.1643 0.0131 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1BIC(mod1, mod2, mod3)## df BIC
## mod1 12 35241.84
## mod2 8 35221.75
## mod3 7 35220.10Il modello finale scelto, il mod3, è una versione semplificata del mod2, da cui è stata rimossa la variabile Tipo.parto. Questa scelta è stata fatta perché l’obiettivo è prevedere il peso del neonato prima della nascita, momento in cui Tipo.parto non è ancora noto. Anche se risultava significativa nel mod2, il suo contributo era minimo, e si è preferito ottenere un modello più semplice ed essenziale.
Bontà di Adattamento del Modello: I valori di Multiple R-squared (0.727) e Adjusted R-squared (0.7265) del modello mod3 sono molto simili a quelli di mod2 (0.7277 e 0.727). Questo mostra che, anche se la variabile Tipo.parto era statisticamente significativa in mod2, la sua rimozione ha avuto un effetto trascurabile sulla capacità del modello di spiegare il peso del neonato. Inoltre, l’F-statistic è aumentato (1328) e il p-value è rimasto molto basso, indicando che mod3 è un modello nel complesso altamente significativo.
Confronto tra mod2 e mod3 usando anova(mod2, mod3): Il confronto tra mod2 e mod3 tramite ANOVA mostra un p-value di 0.0131, quindi una differenza statisticamente significativa tra i due modelli. Tuttavia, la differenza nella somma dei quadrati degli errori è piccola, segnalando che la perdita di informazione dovuta all’eliminazione di Tipo.parto è minima in termini di varianza spiegata.
Confronto tra mod1, mod2 e mod3 usando BIC(mod1, mod2, mod3): I valori di BIC mostrano una riduzione costante da mod1 (35241.84) a mod2 (35221.75), fino a mod3 (35220.10). Questo conferma che mod3 è il modello più adatto, poiché il BIC premia i modelli più semplici che riescono comunque a spiegare bene i dati.
Conclusione: La scelta di mod3 come modello finale è giustificata dalla sua semplicità (meno predittori) e dal valore di BIC più basso, che segnala un miglior equilibrio tra accuratezza e complessità. Anche se Tipo.parto era significativa in mod2, la sua esclusione non ha compromesso la capacità predittiva del modello, che resta adatto all’obiettivo di prevedere il peso alla nascita prima del parto. L’ANOVA segnala una differenza significativa rispetto a mod2, ma il lieve calo nell’R-squared e il BIC più favorevole supportano mod3 come scelta più efficiente e pratica per la predizione prenatale.
Analisi della Qualità del Modello
Sono stati analizzati i residui del modello finale per verificare le assunzioni della regressione lineare (linearità, omoschedasticità, normalità). Sono stati utilizzati grafici dei residui, test di Shapiro-Wilk, Breusch-Pagan e Durbin-Watson.
# Analisi dei residui
par(mfrow=c(2,2))
plot(mod3)
par(mfrow=c(2,2))
shapiro.test(residuals(mod3))##
## Shapiro-Wilk normality test
##
## data: residuals(mod3)
## W = 0.97408, p-value < 2.2e-16lmtest::bptest(mod3)##
## studentized Breusch-Pagan test
##
## data: mod3
## BP = 90.253, df = 5, p-value < 2.2e-16lmtest::dwtest(mod3)##
## Durbin-Watson test
##
## data: mod3
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0#leverage
lev<-hatvalues(mod3)
plot(lev)
p<-sum(lev)
n<-length(lev)
soglia=2*p/n
abline(h=soglia,col=2)
lev[lev>soglia]## 13 15 34 67 89 96
## 0.005630918 0.007050422 0.006744143 0.005891590 0.012816470 0.005351586
## 101 106 131 134 151 155
## 0.007526751 0.014479871 0.007229339 0.007552819 0.010883937 0.007207682
## 161 189 190 204 205 206
## 0.020335483 0.004893297 0.005366557 0.014490604 0.005351634 0.009476652
## 220 294 305 310 312 315
## 0.007393997 0.005912765 0.005442061 0.028812123 0.013169272 0.005385800
## 378 440 442 445 486 492
## 0.015934061 0.005404736 0.007723662 0.007509382 0.005164446 0.008274018
## 497 516 582 587 592 614
## 0.005166306 0.013079851 0.011665555 0.008412325 0.006384116 0.005299262
## 638 656 657 684 697 702
## 0.006688287 0.005927777 0.005322685 0.008818987 0.005863826 0.005202259
## 729 748 750 757 765 805
## 0.005023115 0.008565543 0.006942097 0.008145491 0.006070298 0.014356657
## 828 893 895 913 928 946
## 0.007179817 0.005075205 0.005295896 0.005571144 0.022742332 0.006909044
## 947 956 985 1008 1014 1049
## 0.008409465 0.007784123 0.007039416 0.005343037 0.008470133 0.004956169
## 1067 1091 1106 1130 1166 1181
## 0.008465430 0.008933360 0.005967317 0.031728597 0.005513559 0.005677676
## 1188 1200 1219 1238 1248 1273
## 0.006477203 0.005492370 0.030694311 0.005908078 0.014622914 0.007085831
## 1291 1293 1311 1321 1325 1356
## 0.006117497 0.006073639 0.009625908 0.009293111 0.004857169 0.005303442
## 1357 1385 1395 1400 1402 1411
## 0.006965051 0.012636943 0.005126697 0.005925069 0.004811441 0.008048184
## 1420 1428 1429 1450 1505 1551
## 0.005155654 0.008192811 0.021757172 0.015104831 0.013330439 0.048769569
## 1553 1556 1573 1593 1606 1610
## 0.008504889 0.005919673 0.005047204 0.005623758 0.005001812 0.008722184
## 1617 1619 1628 1686 1693 1701
## 0.004866796 0.015067498 0.005069731 0.009349313 0.005077858 0.010842957
## 1712 1718 1727 1735 1780 1781
## 0.006992084 0.006958857 0.013300523 0.004884846 0.025538678 0.016832361
## 1809 1827 1868 1892 1962 1967
## 0.008707504 0.006065698 0.005205637 0.005332812 0.005540442 0.005337356
## 1977 2037 2040 2046 2086 2089
## 0.006927281 0.004889127 0.011494872 0.005471670 0.013193090 0.006293550
## 2098 2114 2115 2120 2140 2146
## 0.005094455 0.013316875 0.011772090 0.018659995 0.006244232 0.005802168
## 2148 2149 2157 2175 2200 2215
## 0.007926839 0.013583436 0.005907225 0.032527273 0.011670024 0.004892265
## 2216 2220 2221 2224 2225 2244
## 0.008117864 0.005414040 0.021628717 0.005838076 0.005591261 0.006929217
## 2257 2307 2317 2318 2337 2359
## 0.006170254 0.013965608 0.007673614 0.004831118 0.005230450 0.010067364
## 2408 2422 2436 2437 2452 2458
## 0.009696691 0.021532808 0.004986522 0.023943328 0.023838489 0.008506087
## 2471 2478
## 0.020903740 0.005775173#outliers
plot(rstudent(mod3))
abline(h=c(-2,2))
car::outlierTest(mod3)## rstudent unadjusted p-value Bonferroni p
## 1551 10.051908 2.4906e-23 6.2265e-20
## 155 5.027798 5.3138e-07 1.3285e-03
## 1306 4.827238 1.4681e-06 3.6702e-03#distanza di cook
cook<-cooks.distance(mod3)
plot(cook,ylim = c(0,1)) 
L’analisi dei residui del modello mod3 è stata condotta per verificarne le assunzioni. I grafici plot(mod3) (Immagine 1) mostrano:
Residuals vs Fitted: Non si osservano pattern evidenti, il che suggerisce che la relazione tra predittori e risposta è abbastanza lineare e che l’eteroschedasticità è contenuta, anche se ci sono alcuni punti estremi.
Normal Q-Q: I residui si discostano leggermente dalla linea ideale, segnalando una possibile violazione dell’assunzione di normalità. Questo è confermato dal test di Shapiro-Wilk (W = 0.974, p < 2.2e-16), che rigetta l’ipotesi di normalità dei residui.
Scale-Location: La dispersione dei residui standardizzati rispetto ai valori fitted non è del tutto omogenea, suggerendo una lieve eteroschedasticità.
Residuals vs Leverage: Alcune osservazioni mostrano valori di leverage elevati, in particolare la numero 1551, che potrebbe influenzare in modo rilevante il modello.
Il test di Breusch-Pagan (bptest, BP = 90.253, df = 5, p < 2.2e-16) indica una significativa eteroschedasticità. Il test di Durbin-Watson (dwtest, DW = 1.9535, p = 0.1224) non fornisce evidenza significativa di autocorrelazione dei residui.
Punti di Leva (Leverage): Il grafico dei leverage (plot(lev)) mostra alcuni punti con valori di leverage superiori alla soglia (2*p/n ≈ 0.004). Le osservazioni con leverage elevato sono state identificate e includono, tra le altre, 15, 155, e 161 (Immagine 2 - grafico in alto a sinistra).
Outlier: Il grafico degli outliers mostra che l’osservazione 1551 supera ampiamente la soglia di ±2, configurandosi come outlier significativo, come confermato anche dal test di Bonferroni (t = 10.05, p < 2e-16). Anche le osservazioni 155 e 1306 risultano potenzialmente influenti, con p-value corretti significativi.
Distanza di Cook: Il grafico della distanza di Cook (plot(cook)) mostra che alcune osservazioni, in particolare la 1551, esercitano un’influenza considerevole sul modello (Immagine 2 - grafico in basso).
Conclusioni: L’analisi dei residui suggerisce alcune violazioni delle assunzioni del modello lineare, in particolare la normalità e l’omoschedasticità. Sono stati identificati diversi punti con elevato leverage e outlier influenti (in primis l’osservazione 1551), che potrebbero influenzare le stime dei coefficienti e la validità delle inferenze del modello.”
Trasformazione della Variabile Risposta
È stata applicata una trasformazione di Box-Cox per migliorare la normalità e l’omoschedasticità dei residui. La trasformazione di radice quadrata è stata ritenuta appropriata.
library(MASS)
boxcox(mod3)
boxcox_risultati <- boxcox(mod3)
lambda_ottimale <- boxcox_risultati$x[which.max(boxcox_risultati$y)]
print(paste("Il valore ottimale di lambda è:", lambda_ottimale))## [1] "Il valore ottimale di lambda è: 0.545454545454546"PesoNeonatale_trasformato <- sqrt(dati$Peso)
mod3_tr <- lm(PesoNeonatale_trasformato ~ . -Peso -Anni.madre -Fumatrici -Ospedale -Tipo.parto, data = dati)
summary(mod3_tr)##
## Call:
## lm(formula = PesoNeonatale_trasformato ~ . - Peso - Anni.madre -
## Fumatrici - Ospedale - Tipo.parto, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.2551 -1.5609 -0.0605 1.4752 23.3685
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -37.548308 1.176930 -31.904 < 2e-16 ***
## N.gravidanze 0.110906 0.037631 2.947 0.00324 **
## Gestazione 0.389886 0.032934 11.838 < 2e-16 ***
## Lunghezza 0.094180 0.002607 36.126 < 2e-16 ***
## Cranio 0.095434 0.003696 25.819 < 2e-16 ***
## SessoM 0.611264 0.097140 6.293 3.68e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.381 on 2494 degrees of freedom
## Multiple R-squared: 0.7574, Adjusted R-squared: 0.7569
## F-statistic: 1557 on 5 and 2494 DF, p-value: < 2.2e-16par(mfrow = c(2, 2))
plot(mod3_tr)
BIC(mod3, mod3_tr)## df BIC
## mod3 7 35220.10
## mod3_tr 7 11481.57Sebbene un confronto diretto tramite ANOVA tra mod3 e mod3_tr non sia appropriato a causa della diversa variabile risposta (Peso in mod3 e la sua radice quadrata in mod3_tr), il BIC fornisce un criterio valido per confrontare i due modelli. Il valore BIC notevolmente inferiore per mod3_tr (11481.57) rispetto a mod3 (35220.10) suggerisce che mod3_tr è il modello preferibile per descrivere la relazione tra le variabili predittive e il peso del neonato (trasformato).
Infine, per valutare l’influenza dell’osservazione 1551, identificata come outlier influente, è stato costruito un nuovo modello (mod_no_1551) escludendo tale osservazione. Inoltre, come per mod3_tr, la variabile dipendente Peso è stata trasformata utilizzando la radice quadrata (sqrt(dati_no_1551$Peso))
#creazione nuovo dataframe senza l'osservazione 1551
dati_no_1551 <- dati[-1551, ]
PesoNeonatale_trasformato <- sqrt(dati_no_1551$Peso)
#creazione del modello
mod_no_1551 <- lm(PesoNeonatale_trasformato ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso,
data = dati_no_1551)
summary(mod_no_1551)##
## Call:
## lm(formula = PesoNeonatale_trasformato ~ N.gravidanze + Gestazione +
## Lunghezza + Cranio + Sesso, data = dati_no_1551)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.399 -1.531 -0.052 1.468 13.248
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -37.568401 1.153027 -32.582 < 2e-16 ***
## N.gravidanze 0.117017 0.036872 3.174 0.00152 **
## Gestazione 0.365604 0.032352 11.301 < 2e-16 ***
## Lunghezza 0.099882 0.002614 38.216 < 2e-16 ***
## Cranio 0.089933 0.003661 24.568 < 2e-16 ***
## SessoM 0.612522 0.095167 6.436 1.46e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.333 on 2493 degrees of freedom
## Multiple R-squared: 0.7669, Adjusted R-squared: 0.7664
## F-statistic: 1640 on 5 and 2493 DF, p-value: < 2.2e-16BIC(mod3_tr, mod_no_1551)## df BIC
## mod3_tr 7 11481.57
## mod_no_1551 7 11374.43L’output di summary(mod_no_1551) rivela che tutte le variabili predittive (N.gravidanze, Gestazione, Lunghezza, Cranio, e Sesso) rimangono altamente significative (p < 0.01). L’adjusted R-squared, invece, è migliorato (0.7664 vs 0.7569).
Il confronto dei valori BIC tra il modello mod3_tr (stimato sull’intero dataset e con con la variabile risposta trasformata) e il modello mod_no_1551 (stimato escludendo l’osservazione 1551 e con la variabile risposta trasformata) rivela che mod_no_1551 presenta un valore BIC inferiore (11374.43) rispetto a mod3_tr (11481.57). Questa riduzione nel BIC suggerisce che il modello mod_no_1551, stimato senza l’osservazione influente e con la trasformazione della variabile risposta, è leggermente preferibile a mod3_tr.
Procediamo con l’analisi dei residui:
par(mfrow = c(2, 2))
plot(mod_no_1551)
par(mfrow=c(2,2))
shapiro.test(residuals(mod_no_1551))##
## Shapiro-Wilk normality test
##
## data: residuals(mod_no_1551)
## W = 0.9914, p-value = 4.771e-11lmtest::bptest(mod_no_1551)##
## studentized Breusch-Pagan test
##
## data: mod_no_1551
## BP = 15.888, df = 5, p-value = 0.00717lmtest::dwtest(mod_no_1551)##
## Durbin-Watson test
##
## data: mod_no_1551
## DW = 1.9485, p-value = 0.099
## alternative hypothesis: true autocorrelation is greater than 0#leverage
lev<-hatvalues(mod_no_1551)
plot(lev)
p<-sum(lev)
n<-length(lev)
soglia=2*p/n
abline(h=soglia,col=2)
lev[lev>soglia]## 13 15 34 61 67 89
## 0.005650394 0.007130572 0.006935698 0.004869002 0.005891963 0.012816487
## 96 101 106 131 134 151
## 0.005357728 0.007596954 0.014485358 0.007270462 0.007556342 0.010884155
## 155 161 189 190 204 205
## 0.007446908 0.020391756 0.004964398 0.005385960 0.014525113 0.005397518
## 206 220 294 304 305 310
## 0.009487528 0.007738686 0.005922509 0.004963924 0.005442684 0.029052940
## 312 315 348 378 440 442
## 0.013282873 0.005490179 0.004859888 0.015943467 0.005417654 0.007731162
## 445 486 492 497 516 582
## 0.007513589 0.005164643 0.008396610 0.005168514 0.013093538 0.011688847
## 587 592 614 638 656 657
## 0.008427004 0.006389106 0.005313060 0.006698260 0.006033902 0.005406699
## 684 697 702 729 748 750
## 0.009094955 0.006010925 0.005202699 0.005144702 0.008591238 0.006959863
## 757 765 805 828 893 895
## 0.008148112 0.006079582 0.014398868 0.007181297 0.005083104 0.005303785
## 913 928 946 947 949 956
## 0.005573692 0.022931675 0.006921754 0.008471606 0.004821777 0.007910332
## 964 985 1008 1014 1049 1067
## 0.004869054 0.007045075 0.005343464 0.008645855 0.004958181 0.008465605
## 1091 1106 1130 1166 1181 1188
## 0.008933569 0.005981876 0.031853467 0.005615493 0.005691116 0.006736247
## 1200 1219 1238 1248 1273 1291
## 0.005645850 0.030694454 0.005911473 0.014638558 0.007180250 0.006155082
## 1293 1294 1311 1321 1325 1356
## 0.006180971 0.004948504 0.009628045 0.009319849 0.004862079 0.005336228
## 1357 1385 1395 1400 1402 1411
## 0.006973085 0.012637033 0.005238687 0.005990413 0.005013423 0.008100110
## 1420 1428 1429 1450 1505 1553
## 0.005155699 0.008348777 0.022077783 0.015113725 0.013331272 0.008538588
## 1556 1560 1573 1593 1606 1610
## 0.005948770 0.004863437 0.005064179 0.005652255 0.005190559 0.008844546
## 1617 1619 1628 1686 1693 1701
## 0.004869744 0.015297674 0.005177326 0.009349658 0.005121069 0.010984425
## 1712 1718 1727 1735 1780 1781
## 0.007202578 0.006968691 0.013306762 0.005048297 0.025579769 0.016846099
## 1809 1827 1868 1892 1962 1967
## 0.008717017 0.006092729 0.005240785 0.005349151 0.005635253 0.005338037
## 1977 2037 2040 2046 2086 2089
## 0.006927445 0.004897704 0.012068488 0.005483403 0.013196335 0.006372111
## 2098 2114 2115 2120 2140 2146
## 0.005114483 0.013412224 0.011998021 0.018662104 0.006251984 0.005806194
## 2148 2149 2157 2175 2200 2215
## 0.007969158 0.013594736 0.005910374 0.032538514 0.011671278 0.004893518
## 2216 2220 2221 2224 2225 2244
## 0.008125783 0.005414196 0.021638156 0.005858043 0.005651613 0.006929239
## 2257 2307 2317 2318 2337 2359
## 0.006250963 0.013996551 0.007714734 0.004869147 0.005355071 0.010072587
## 2408 2422 2436 2437 2452 2458
## 0.009700281 0.021552185 0.004990057 0.024101713 0.023840276 0.008507518
## 2471 2478
## 0.020920615 0.005796548#outliers
plot(rstudent(mod_no_1551))
abline(h=c(-2,2))
car::outlierTest(mod_no_1551)## rstudent unadjusted p-value Bonferroni p
## 155 5.736393 1.0841e-08 2.7092e-05
## 1399 -4.480241 7.7924e-06 1.9473e-02
## 1694 4.459686 8.5719e-06 2.1421e-02
## 1306 4.297473 1.7939e-05 4.4829e-02#distanza di cook
cook<-cooks.distance(mod_no_1551)
plot(cook,ylim = c(0,1)) 
L’analisi visiva dei residui di mod_no_1551 (dai grafici plot(mod_no_1551)) mostra alcuni miglioramenti rispetto a mod3:
Residuals vs Fitted: La dispersione dei residui appare più omogenea. Normal Q-Q: I punti si avvicinano maggiormente alla linea retta, suggerendo un miglioramento nella normalità dei residui. Questo è supportato dal test di Shapiro-Wilk (shapiro.test, W = 0.9914, p = 4.771e-11), il cui p-value è aumentato rispetto a mod3, sebbene rimanga significativo, indicando che la normalità perfetta non è ancora raggiunta. Scale-Location: La dispersione dei residui standardizzati sembra più uniforme. Il test di Breusch-Pagan (bptest, BP = 15.888, df = 5, p = 0.00717) indica ancora una presenza significativa di eteroschedasticità, sebbene il valore del test sia diminuito rispetto a mod3. Il test di Durbin-Watson (dwtest, DW = 1.9485, p = 0.099) continua a non fornire evidenza significativa di autocorrelazione.
Punti di Leva (Leverage) e Outlier nel Modello Trasformato: Il grafico dei leverage (plot(lev)) mostra pattern simili a mod3 per i punti con leverage elevato. Il grafico degli studentized residuals (plot(rstudent(mod_no_1551))) e il test di outlier di Bonferroni (outlierTest) identificano ancora alcuni outlier significativi (osservazioni 155, 1399, 1694, 1306), sebbene l’estrema outlierness dell’osservazione 1551 sia stata eliminata per costruzione.
Il grafico della distanza di Cook (plot(cook)) mostra che l’influenza delle altre osservazioni sul modello è generalmente bassa dopo la rimozione di 1551.
Conclusioni sull’Influenza dell’Osservazione 1551: L’esclusione dell’osservazione 1551 e la trasformazione della variabile dipendente hanno portato a un miglioramento della bontà di adattamento del modello (R-squared più elevato) e a un parziale miglioramento della normalità e dell’omoschedasticità dei residui. Tuttavia, il test di Shapiro-Wilk e di Breusch-Pagan rimangono significativi, suggerendo che le assunzioni del modello lineare non sono ancora pienamente soddisfatte.
#creazione dataframe senza outliers
dati_no_outliers <- dati[-c(1551, 155, 1306), ]
PesoNeonatale_trasformato_no_out <- sqrt(dati_no_outliers$Peso)
mod_no_outliers <- lm(PesoNeonatale_trasformato_no_out ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso,
data = dati_no_outliers)
summary(mod_no_outliers)##
## Call:
## lm(formula = PesoNeonatale_trasformato_no_out ~ N.gravidanze +
## Gestazione + Lunghezza + Cranio + Sesso, data = dati_no_outliers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.4214 -1.5189 -0.0542 1.4775 10.3839
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -37.709131 1.142754 -32.998 < 2e-16 ***
## N.gravidanze 0.124054 0.036522 3.397 0.000693 ***
## Gestazione 0.362632 0.032038 11.319 < 2e-16 ***
## Lunghezza 0.100985 0.002594 38.923 < 2e-16 ***
## Cranio 0.089047 0.003627 24.552 < 2e-16 ***
## SessoM 0.604314 0.094303 6.408 1.75e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.31 on 2491 degrees of freedom
## Multiple R-squared: 0.7709, Adjusted R-squared: 0.7705
## F-statistic: 1677 on 5 and 2491 DF, p-value: < 2.2e-16par(mfrow = c(2, 2))
plot(mod_no_outliers)
# Confronto BIC
BIC(mod_no_1551, mod_no_outliers)## df BIC
## mod_no_1551 7 11374.43
## mod_no_outliers 7 11315.95#Analisi residui
shapiro.test(residuals(mod_no_outliers))##
## Shapiro-Wilk normality test
##
## data: residuals(mod_no_outliers)
## W = 0.99462, p-value = 6.803e-08lmtest::bptest(mod_no_outliers)##
## studentized Breusch-Pagan test
##
## data: mod_no_outliers
## BP = 12.426, df = 5, p-value = 0.02939lmtest::dwtest(mod_no_outliers)##
## Durbin-Watson test
##
## data: mod_no_outliers
## DW = 1.9487, p-value = 0.09977
## alternative hypothesis: true autocorrelation is greater than 0#leverage
lev<-hatvalues(mod_no_outliers)
plot(lev)
p<-sum(lev)
n<-length(lev)
soglia=2*p/n
abline(h=soglia,col=2)
lev[lev>soglia]## 13 15 34 61 67 89
## 0.005650754 0.007135223 0.006951881 0.004887490 0.005906436 0.012820261
## 96 101 106 131 134 151
## 0.005358916 0.007624834 0.014507340 0.007275889 0.007559014 0.010887329
## 161 189 190 204 205 206
## 0.020430712 0.004970145 0.005387164 0.014547073 0.005402605 0.009520238
## 220 294 304 305 310 312
## 0.007766771 0.005929618 0.004991429 0.005454507 0.029096990 0.013349168
## 315 348 378 440 442 445
## 0.005502863 0.004887643 0.015959180 0.005422279 0.007731498 0.007527907
## 486 492 497 516 582 587
## 0.005169314 0.008430002 0.005170852 0.013096421 0.011695556 0.008430863
## 592 614 638 656 657 684
## 0.006397171 0.005313321 0.006724033 0.006039275 0.005414839 0.009114526
## 697 702 729 748 750 757
## 0.006022135 0.005207930 0.005158386 0.008613529 0.006973918 0.008149543
## 765 805 828 893 895 913
## 0.006089025 0.014436903 0.007183933 0.005087575 0.005314691 0.005575157
## 928 946 947 949 956 964
## 0.023025539 0.006925025 0.008494534 0.004822973 0.007926018 0.004890468
## 985 1008 1014 1049 1067 1091
## 0.007048319 0.005343996 0.008680452 0.004958544 0.008473680 0.008952845
## 1106 1130 1166 1181 1188 1200
## 0.005986907 0.031854299 0.005628823 0.005696511 0.006752567 0.005661675
## 1219 1238 1248 1273 1291 1293
## 0.030705628 0.005930019 0.014686990 0.007182620 0.006159357 0.006191574
## 1294 1311 1321 1325 1356 1357
## 0.004964365 0.009631204 0.009319961 0.004863849 0.005339153 0.006976712
## 1385 1395 1400 1402 1411 1420
## 0.012652950 0.005243948 0.005992427 0.005030067 0.008102884 0.005157693
## 1428 1429 1450 1505 1553 1556
## 0.008381003 0.022132668 0.015115212 0.013338679 0.008545037 0.005948863
## 1560 1573 1593 1606 1610 1617
## 0.004880111 0.005076781 0.005671387 0.005207356 0.008850640 0.004871165
## 1619 1628 1686 1693 1701 1712
## 0.015372557 0.005203434 0.009372427 0.005126615 0.011051528 0.007219694
## 1718 1727 1735 1780 1781 1809
## 0.006976554 0.013308609 0.005083897 0.025643216 0.016854859 0.008723177
## 1827 1868 1892 1962 1967 1977
## 0.006120500 0.005245142 0.005356059 0.005653250 0.005342543 0.006927927
## 2037 2040 2046 2086 2089 2098
## 0.004898554 0.012126418 0.005488012 0.013198725 0.006397889 0.005114497
## 2114 2115 2120 2140 2146 2148
## 0.013461997 0.012007149 0.018695854 0.006262405 0.005813267 0.007969789
## 2149 2157 2175 2200 2215 2216
## 0.013633239 0.005918830 0.032558581 0.011702790 0.004893892 0.008127433
## 2220 2221 2224 2225 2244 2257
## 0.005414568 0.021654583 0.005882968 0.005652065 0.006930660 0.006277938
## 2307 2317 2318 2337 2359 2408
## 0.014049577 0.007715169 0.004869664 0.005372914 0.010076902 0.009727853
## 2422 2436 2437 2452 2458 2471
## 0.021554960 0.004992074 0.024215698 0.023884475 0.008515571 0.020924633
## 2478
## 0.005796981#outliers
plot(rstudent(mod_no_outliers))
abline(h=c(-2,2))
car::outlierTest(mod_no_outliers)## rstudent unadjusted p-value Bonferroni p
## 1399 -4.534830 6.038e-06 0.015077
## 1694 4.517251 6.557e-06 0.016373#distanza di cook
cook<-cooks.distance(mod_no_outliers)
plot(cook,ylim = c(0,1))
Dopo aver rimosso anche le osservazioni 155 e 1306, il modello è stato ricalcolato come mod_no_outliers. L’analisi del summary() rivela un leggero aumento dell’R-squared (0.7709) rispetto a mod_no_1551 (0.7669), suggerendo un miglioramento marginale nella varianza spiegata. I coefficienti stimati presentano piccole variazioni, ma la loro significatività rimane generalmente elevata.
L’analisi dei residui mostra che il test di Shapiro-Wilk (W = 0.99462, p = 6.803e-08) indica ancora una deviazione significativa dalla normalità, sebbene il p-value sia aumentato rispetto a mod_no_1551. Il test di Breusch-Pagan (BP = 12.426, df = 5, p = 0.02939) continua a suggerire la presenza di eteroschedasticità, con un p-value leggermente superiore rispetto a prima. Il test di Durbin-Watson (DW = 1.9487, p = 0.09977) non evidenzia autocorrelazione significativa.
Il confronto del BIC mostra un valore inferiore per mod_no_outliers (11315.95) rispetto a mod_no_1551 (11374.43). Questo suggerisce che, secondo il criterio BIC, il modello stimato senza i tre outlier è leggermente preferibile.
Infine, l’esecuzione di car::outlierTest() su mod_no_outliers identifica ancora le osservazioni 1399 e 1694 come outlier significativi dopo la correzione di Bonferroni. Questo indica che, nonostante la rimozione di tre outlier, potrebbero persistere ulteriori osservazioni con residui insolitamente grandi.
In conclusione, la rimozione degli outlier ha portato a un modesto miglioramento nella bontà di adattamento del modello (R-squared) e a una leggera preferenza secondo il criterio BIC. Tuttavia, le assunzioni di normalità ed omoschedasticità dei residui non sono pienamente soddisfatte e persistono ulteriori outlier significativi secondo il test di Bonferroni. La decisione finale sul modello da utilizzare dovrebbe bilanciare il miglioramento nell’adattamento con la potenziale perdita di generalizzabilità dovuta alla rimozione di osservazioni.
3. Previsioni e Risultati
Per effettuare una previsione pratica, si è considerato lo scenario di una neonata (Sesso Femminile) alla terza gravidanza della madre, con una gestazione di 39 settimane. Poiché i valori di Lunghezza e Diametro Cranio non erano forniti nella richiesta, sono state utilizzate le medie di queste variabili, calcolate precedentemente tramite la funzione analisi_descrittiva, per ottenere una stima del peso neonatale.
# Estrazione delle medie di lunghezza e cranio
media_lunghezza_pratica <- analisi_completa[analisi_completa$Variabile == "Lunghezza", "Media"]
media_cranio_pratica <- analisi_completa[analisi_completa$Variabile == "Diametro Cranio", "Media"]
# Creazione di un nuovo dataframe con le caratteristiche specificate (usando le medie calcolate)
nuovi_dati_previsione_pratica_medie <- data.frame(
N.gravidanze = 3,
Gestazione = 39,
Lunghezza = media_lunghezza_pratica,
Cranio = media_cranio_pratica,
Sesso = factor("F", levels = c("F", "M"))
)
nuovi_dati_previsione_pratica_medie$SessoM <- ifelse(nuovi_dati_previsione_pratica_medie$Sesso == "M", 1, 0)
# Previsione con il modello mod_no_outliers
previsione_pratica_medie_sqrt <- predict(mod_no_outliers, newdata = nuovi_dati_previsione_pratica_medie)
# Calcolo del peso in grammi
previsione_pratica_peso_grammi_medie <- (previsione_pratica_medie_sqrt)^2
cat(paste("Il peso previsto (usando le medie calcolate) è di circa:", round(previsione_pratica_peso_grammi_medie, 2), "grammi.\n"))## Il peso previsto (usando le medie calcolate) è di circa: 3253.64 grammi.Dunque, utilizzando il modello validato mod_no_outliers, il peso previsto per la neonata è di circa 3253.64 grammi
4. Visualizzazioni
Per illustrare visivamente i risultati del modello e le relazioni chiave tra le variabili, sono stati generati diversi grafici. Questi aiuteranno a comprendere come i fattori predittivi influenzano il peso alla nascita. N.B. La variabile “fumatrici” non è inclusa nel modello finale, dunque, come variabile di stratificazione è stata scelta la variabile “Sesso”.
Grafico 1: Relazione tra Durata della Gestazione e Peso alla Nascita, stratificata per Sesso
ggplot(data = dati_no_outliers) +
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 Durata della Gestazione e Peso alla Nascita",
x = "Durata della Gestazione (settimane)",
y = "Peso alla Nascita (grammi)",
color = "Sesso") +
theme_minimal()
Si osserva chiaramente una relazione positiva tra gestazione e peso per entrambi i sessi. Le linee di tendenza mostrano che, a parità di gestazione, i neonati maschi (azzurro) tendono a un peso leggermente superiore rispetto alle femmine (rosso). Questo riflette la significatività di Sesso nel modello.
Grafico 2: Relazione tra Numero di Gravidanze e Peso alla Nascita, stratificata per Sesso
ggplot(data = dati_no_outliers) +
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 N.gravidanze e Peso alla Nascita",
x = "N.gravidanze",
y = "Peso alla Nascita (grammi)",
color = "Sesso") +
theme_minimal()
La visualizzazione della relazione tra il numero di gravidanze e il peso, differenziata per sesso, mostra una tendenza molto debole. In questo caso, si osserva un leggero divario nel peso medio tra maschi e femmine, consistente con i risultati del modello, ma l’influenza del numero di gravidanze sul peso appare marginale.
Grafico 3: Relazione tra Lunghezza e Peso alla Nascita, stratificata per Sesso
ggplot(data = dati_no_outliers) +
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 del Neonato e Peso alla Nascita",
x = "Lunghezza del Neonato (cm)",
y = "Peso alla Nascita (grammi)",
color = "Sesso") +
theme_minimal()
Questo grafico illustra la forte relazione positiva tra la lunghezza del neonato e il suo peso alla nascita, differenziata per sesso. Si osserva che a parità di lunghezza, i neonati maschi tendono a pesare leggermente di più delle femmine, una tendenza coerente con la significatività di Sesso nel modello finale.
Grafico 4: Relazione tra Diametro Cranico e Peso alla Nascita, stratificata per Sesso
ggplot(data = dati_no_outliers) +
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 Cranico e Peso alla Nascita",
x = "Diametro Cranico (cm)",
y = "Peso alla Nascita (grammi)",
color = "Sesso") +
theme_minimal()
La visualizzazione mostra una chiara correlazione positiva tra il diametro cranico e il peso alla nascita. Analogamente alla lunghezza, anche qui si rileva una tendenza dei maschi a presentare un peso leggermente superiore rispetto alle femmine per lo stesso diametro cranico, confermando il ruolo di Sesso come predittore nel modello.
5. Considerazioni Finali
Questo progetto ha mirato alla costruzione e selezione di un modello predittivo robusto per il peso neonatale, basandosi su un’ampia gamma di variabili demografiche e cliniche. Il percorso di modellazione ha incluso diverse fasi di raffinamento, partendo da un modello completo e procedendo verso una versione più parsimoniosa. Il modello finale selezionato, mod_no_outliers, che prevede la radice quadrata del peso neonatale, rappresenta il miglior compromesso tra bontà di adattamento e semplicità. In conclusione, questo progetto ha fornito un modello predittivo valido e interpretabile per il peso neonatale, identificando i fattori chiave che ne influenzano la variabilità. Le intuizioni derivanti da questo modello possono essere un punto di partenza prezioso per studi clinici e interventi volti a monitorare e migliorare gli esiti della nascita.