Modello Statistico per la Previsione del Peso Neonatale
Zoe Cataldi
2025-01-09
Azienda: Neonatal Health Solutions
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.
ANALISI DESCRITTIVA DELLE VARIABILI
I dati raccolti riguardano 2500 neonati.
Iniziamo questo studio con un’analisi descrittiva delle variabili del nostro campione, principalmente per capirne la distribuzione ed evidenziare eventuali anomalie.
data = read.csv("neonati.csv",
sep=",",
stringsAsFactors = T)
attach(data)
table1=summary(data)
n=nrow(data)
kable(table1, digits = 2)| 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 |
Da questa tabella possiamo notare quali sono le variabili prese in esame e per ognuna di esse vengono riportati gli indici di posizione, quali valore minimo e valore massimo, primo e terzo quartile, mediana e media.
Successivamente vedremo meglio ogni variabile per comprenderne l’andamento.
Può essere già evidenziata un’anomalia nella variabile Anni.madre, per la presenza di un valore minimo pari a zero. Iniziamo quindi dallo studio di questa variabile.
STUDIO DELLA VARIABILE ANNI.MADRE
#studio della variabile Anni.madre
boxplot(Anni.madre,
main="Distribuzione età materna")
Anni.madre_CL=cut(Anni.madre, breaks=c(0,11,22,33,44,54), include.lowest=TRUE)
table2=table(Anni.madre_CL)
kable(table2, col.names=c("Classi_Anni.madre", "Freq.Assoluta"))| Classi_Anni.madre | Freq.Assoluta |
|---|---|
| [0,11] | 2 |
| (11,22] | 349 |
| (22,33] | 1769 |
| (33,44] | 378 |
| (44,54] | 2 |
Come mostrato dal boxplot e confermato dalla suddivisione in classi, si rilevano due valori impossibili per l’età materna pari a 0 e 1. Si decide quindi di eliminare tali valori dal campione per evitare che influenzino in maniera errata l’analisi. La suddivisione in classi mostra anche che la fascia di età più frequente è quella tra i 22 e i 33 anni, con un netto distacco dalle fasce 33-44 e 11-22. Solo 2 casi sono stati registrati superiori ai 44 anni.
La tabella seguente contiene i valori puliti.
data_clean=data[data$Anni.madre >= 11, ]
table3=summary(data_clean)
kable(table3, digits=2)| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso | Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso | |
|---|---|---|---|---|---|---|---|---|---|---|
| Min. :13.00 | Min. : 0.0000 | Min. :0.00000 | Min. :25.00 | Min. : 830 | Min. :310.0 | Min. :235 | Ces: 728 | osp1:816 | F:1255 | |
| 1st Qu.:25.00 | 1st Qu.: 0.0000 | 1st Qu.:0.00000 | 1st Qu.:38.00 | 1st Qu.:2990 | 1st Qu.:480.0 | 1st Qu.:330 | Nat:1770 | osp2:848 | M:1243 | |
| Median :28.00 | Median : 1.0000 | Median :0.00000 | Median :39.00 | Median :3300 | Median :500.0 | Median :340 | NA | osp3:834 | NA | |
| Mean :28.19 | Mean : 0.9816 | Mean :0.04163 | Mean :38.98 | Mean :3284 | Mean :494.7 | Mean :340 | NA | NA | NA | |
| 3rd Qu.:32.00 | 3rd Qu.: 1.0000 | 3rd Qu.:0.00000 | 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.00000 | Max. :43.00 | Max. :4930 | Max. :565.0 | Max. :390 | NA | NA | NA |
STUDIO DELLA VARIABILE N.GRAVIDANZE
Per quanto riguarda la variabile N.gravidanze, ci sono dei valori che possono risultare anomali, come mostrato dal boxplot sottostante. Infatti, un numero di gravidanze pregresse superiore a 6 potrebbe risultare anomalo, soprattutto nell’epoca attuale. Nonostante ciò, non possiamo escludere tali valori a priori, ma possiamo comunque tenerne conto nel corso dell’analisi.
#studio della variabile N.gravidanze
boxplot(N.gravidanze,
main="Distribuzione del numero di gravidanze")
STUDIO DELLA VARIABILE FUMATRICI
Per quanto riguarda la variabile Fumatrici, il grafico sottostante mostra che il numero di madri fumatrici nel nostro campione è nettamente inferiore rispetto al numero di madri non fumatrici.
#studio della variabile Fumatrici
si_no=table(Fumatrici)
barplot(si_no,
main="0=Non Fumatrici / 1=Fumatrici",
ylab="Freq.Assolute",
ylim=c(0,2500))
STUDIO DELLA VARIABILE GESTAZIONE
Il boxplot e la suddivisione in classi della variabile, mostrano che la maggior parte delle gravidanze registrate sono state gravidanze a termine, cioè il parto è avvenuto tra le 37 e le 41 settimane di gestazione. Invece, 162 gravidanze sono state pre-termine (parto prematuro) e 58 sono state gravidanze protratte, quando si sono superate la 41 settimane di gestazione.
#studio della variabile Gestazione
boxplot(Gestazione,
main="Distribuzione settimane di gestazione")
Gestazione_CL=cut(Gestazione, breaks=c(25,36,41,43), include.lowest=TRUE)
table4=table(Gestazione_CL)
kable(table4, col.names=c("Classi_Gestazione", "Freq.Assoluta"))| Classi_Gestazione | Freq.Assoluta |
|---|---|
| [25,36] | 162 |
| (36,41] | 2280 |
| (41,43] | 58 |
STUDIO DELLA VARIABILE PESO
Il boxplot riportato di seguito, rappresenta la distribuzione del peso dei neonati divisi per sesso. Evidentemente, le misurazioni registrate dimostrano che la maggior parte dei neonati ha un peso alla nascita tra i 3000 e i 3500 grammi. Questo è l’intervallo di un peso medio in seguito ad una gravidanza a termine. Sono stati registrati anche neonati con un basso o alto peso alla nascita, rispettivamente con un peso inferiore a 2500 grammi o superiore a 4500 grammi.
#studio della variabile Peso
boxplot(Peso~Sesso,
main="Distribuzione del peso dei neonati")
Valutiamo i pesi registrati nel nostro campione in relazione al resto della popolazione. La media del peso neonatale nella popolazione è pari a 3300 grammi. Il test t di Student sottostante ci restituisce un p-value di circa 0.13 e ciò vuol dire che la media del peso del nostro campione di neonati non è significativamente diversa da quella della popolazione.
#test t per la variabile peso
#H0: mu(peso_camp) = mu(peso_pop)
#H1: mu(peso_camp) != mu(peso_pop)
test=t.test(Peso,
mu = 3300,
conf.level = 0.95,
alternative = "two.sided")
library(broom)
results=tidy(test)
print(results)## # A tibble: 1 × 8
## estimate statistic p.value parameter conf.low conf.high method alternative
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 3284. -1.52 0.130 2499 3263. 3305. One Sampl… two.sidedSTUDIO DELLE VARIABILI LUNGHEZZA E CRANIO
La variabile Lunghezza si riferisce alla lunghezza del corpo dei neonati, mentre la variabile Cranio si riferisce al diametro del cranio dei neonati.
In entrambi i grafici, possiamo notare che la maggior parte delle misurazioni registrate nel nostro campione, corrispondono a quello che viene considerato il valore medio di un neonato nato a termine. Stiamo parlando di circa 50 cm per la lunghezza del corpo e 33-36 cm per il diametro del cranio.
#studio della variabile Lunghezza
boxplot(Lunghezza~Sesso,
main="Distribuzione della lunghezza (mm) del corpo dei neonati")
#studio della variabile Cranio
boxplot(Cranio~Sesso,
main="Distribuzione del diametro (mm) del cranio dei neonati")
Possiamo anche confrontare ad esempio la variabile Lunghezza del nostro campione con il resto della popolazione, ricordando che, come detto poco prima, la lunghezza media del cranio neonatale è fissata a circa 50 cm.
Il seguente test t di Student ci restituisce un p-value estremamente piccolo, indicando che la media del nostro campione è significativamente diversa dalla media della popolazione. Ciò può essere dovuto alla presenza di outliers che influenzano la media del nostro campione.
#test t per la variabile lunghezza
#H0: mu(lungh_camp) = mu(lungh_pop)
#H1: mu(lungh_camp) != mu(lungh_pop)
test=t.test(Lunghezza,
mu = 500,
conf.level = 0.95,
alternative = "two.sided")
results=tidy(test)
print(results)## # A tibble: 1 × 8
## estimate statistic p.value parameter conf.low conf.high method alternative
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 495. -10.1 1.81e-23 2499 494. 496. One Samp… two.sidedUn’altra valutazione che possiamo fare sulle misure antropometriche è quella di saggiare l’ipotesi di uguaglianza tra i due sessi, per capire se ci sono delle differenze significative.
I boxplot precedenti ci hanno mostrato che sembra non esserci una grande differenza tra maschi e femmine in termini di peso, lunghezza del corpo e diametro del cranio.
Per avere una risposta più completa è opportuno però eseguire dei test specifici, riportati di seguito. Nelle tabelle risultanti, si intende per estimate 1 la media calcolata per le femmine, mentre per estimate 2 la media calcolata per i maschi. I p-value molto piccoli dimostrano in realtà che le medie delle misure antropometriche sono significativamente diverse tra i due sessi.
#saggiamo l'ipotesi di uguaglianza delle misure antropometriche tra i due sessi
attach(data_clean)
test1=t.test(Peso~Sesso)
results=tidy(test1)
print(results)## # A tibble: 1 × 10
## estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -247. 3161. 3408. -12.1 7.28e-33 2489. -287. -207.
## # ℹ 2 more variables: method <chr>, alternative <chr>test2=t.test(Lunghezza~Sesso)
results=tidy(test2)
print(results)## # A tibble: 1 × 10
## estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -9.91 490. 500. -9.58 2.23e-21 2457. -11.9 -7.88
## # ℹ 2 more variables: method <chr>, alternative <chr>test3=t.test(Cranio~Sesso)
results=tidy(test3)
print(results)## # A tibble: 1 × 10
## estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -4.84 338. 342. -7.44 1.41e-13 2489. -6.11 -3.56
## # ℹ 2 more variables: method <chr>, alternative <chr>STUDIO DELLE VARIABILI TIPO.PARTO E OSPEDALE
Per prima cosa, dalle seguenti tabelle emerge che nel nostro campione sono stati registrati un numero di parti naturali più del doppio del numero di parti cesarei, mentre per quanto riguarda l’ospedale di provenienza notiamo esserci una distribuzione piuttosto equa tra i tre ospedali oggetto di studio.
table5=table(Tipo.parto)
kable(table5)| Tipo.parto | Freq |
|---|---|
| Ces | 728 |
| Nat | 1770 |
table6=table(Ospedale)
kable(table6)| Ospedale | Freq |
|---|---|
| osp1 | 816 |
| osp2 | 848 |
| osp3 | 834 |
Tramite un test del chi-quadrato possiamo studiare una eventuale associazione tra ospedale e tipo di parto, ma il test restituisce un p-value alto, indicando quindi che c’è indipendenza tra le due variabili.
#saggiamo l'ipotesi "in alcuni ospedali si fanno più parti cesarei"
#H0=indipendenza tra ospedale e tipo di parto
#H1=associazione tra ospedale e tipo di parto
table7=table(Ospedale, Tipo.parto)
kable(table7, caption="Frequenze di parti cesarei e naturali nei tre ospedali")| Ces | Nat | |
|---|---|---|
| osp1 | 242 | 574 |
| osp2 | 254 | 594 |
| osp3 | 232 | 602 |
test=chisq.test(table7)
results=tidy(test)
print(results)## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <int> <chr>
## 1 1.08 0.582 2 Pearson's Chi-squared testSTUDIO DELLA CORRELAZIONE E DELLE INTERAZIONI TRA VARIABILI
A questo punto, possiamo valutare l’impatto di ciascuna variabile indipendente sul peso del neonato ed eventuali interazioni.
Infatti, dalla matrice di correlazione tra le variabili possiamo notare già graficamente, che le variabili più rilevanti rispetto al peso dei neonati sono principalmente la lunghezza del corpo, il diametro del cranio e, in misura minore, i mesi di gestazione della madre. In particolare, gli scatterplot mostrano una correlazione positiva, cioè all’aumentare di una di queste variabili, cresce anche la variabile Peso.
Per quanto riguarda le variabili qualitative, tale matrice non è indicativa, quindi saranno eseguiti altri test successivamente.
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(data_clean, upper.panel = panel.smooth, lower.panel = panel.cor)
In particolare, per la variabile N.gravidanze viene eseguito il test del chi-quadrato che dimostra una associazione tra il numero di gravidanze avute dalla madre con il peso del neonato, ma è bene sottolineare che questa associazione non implica una causalità.
#saggiamo l'associazione tra le variabili peso e numero di gravidanze
boxplot(Peso~N.gravidanze)
table8=table(Peso, N.gravidanze)
test=chisq.test(table8)
results=tidy(test)
print(results)## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <int> <chr>
## 1 4025. 5.45e-12 3432 Pearson's Chi-squared testInvece, a proposito della variabile Fumatrici, sulla base del t.test, non c’è una significativa differenza tra le medie del peso dei neonati di madri fumatrici e non fumatrici. Quindi la variabile “Fumatrici” può essere considerata non rilevante nello studio del peso dei neonati.
boxplot(Peso~Fumatrici)
test=t.test(Peso~Fumatrici)
results=tidy(test)
print(results)## # A tibble: 1 × 10
## estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 49.9 3286. 3236. 1.04 0.302 114. -45.5 145.
## # ℹ 2 more variables: method <chr>, alternative <chr>STUDIO DELLA VARIABILE RISPOSTA
Dopo aver analizzato le variabili esplicative, passiamo ora allo studio della nostra variabile risposta, cioè il peso. Mediante test specifici, come il test di Shapiro-Wilk, risulta avere una distribuzione non Normale (p-value molto piccolo), essere leggermente asimmetrica negativa e leptocurtica, quindi con valori alti più frequenti e più allungata rispetto ad una Normale.
#studio della variabile risposta "Peso"
attach(data_clean)
test=shapiro.test(Peso)
results=tidy(test)
print(results)## # A tibble: 1 × 3
## statistic p.value method
## <dbl> <dbl> <chr>
## 1 0.971 3.38e-22 Shapiro-Wilk normality testskewness_value=moments::skewness(Peso)
kurtosis_value=moments::kurtosis(Peso)-3
skewness_value=round(skewness_value, 2)
kurtosis_value=round(kurtosis_value, 2)
results.ske.kur=data.frame(
Index=c("Skewness", "Kurtosis"),
Values=c(skewness_value, kurtosis_value)
)
print(results.ske.kur)## Index Values
## 1 Skewness -0.65
## 2 Kurtosis 2.03plot(density(Peso))
MODELLO 1 DI REGRESSIONE LINEARE MULTIPLA
Costruiamo un Modello 1 di regressione lineare multipla, escludendo le variabili “Fumatrici”, “Tipo.parto” e “Ospedale” in quanto non risultano rilevanti per lo studio del peso dei neonati.
Il riquadro sottostante mostra i risultati del modello ed è emerso che le variabili più significative sono:
“Gestazione”, con un coefficiente di regressione di 32.69 (correlazione positiva, in particolare per ogni incremento unitario della variabile “Gestazione”, si registra un aumento medio del peso dei neonati di 32.69 grammi, tenendo fissate le altre variabili);
“Lunghezza”, con un coefficiente di regressione di 10.25 (da intrpretare come “Gestazione”);
“Cranio”, con un coefficiente di regressione di 10.52 (da intrpretare come “Gestazione”);
“Sesso”, essendo una variabile qualitativa deve essere interpretata in questo modo: tenendo fissate le altre variabili, nei neonati maschi si registra un peso medio di 77.90 grammi in più rispetto alle femmine.
Il modello mostra anche un R2 aggiustato pari a 0.73 circa, che corrisponde ad una variabilità spiegata del 73%, che è un buon valore.
Vediamo, inoltre, che l’età della madre non è significativa e il numero di gravidanze avute dalla madre è molto poco significativo. Quindi, vedremo di seguito come trattare queste due variabili per semplificare ulteriormente il modello.
#costruiamo il primo modello di regressione lineare multipla
mod1 = lm(Peso ~ Anni.madre+N.gravidanze+Gestazione+Lunghezza+Cranio+
Sesso, data=data_clean)
summary(mod1)##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione +
## Lunghezza + Cranio + Sesso, data = data_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1159.84 -181.98 -15.04 164.16 2634.01
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6712.0654 141.3399 -47.489 < 2e-16 ***
## Anni.madre 0.8909 1.1491 0.775 0.4382
## N.gravidanze 11.1203 4.6710 2.381 0.0174 *
## Gestazione 32.6914 3.8219 8.554 < 2e-16 ***
## Lunghezza 10.2461 0.3008 34.058 < 2e-16 ***
## Cranio 10.5239 0.4271 24.642 < 2e-16 ***
## SessoM 77.8998 11.2125 6.948 4.72e-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.7264
## F-statistic: 1106 on 6 and 2491 DF, p-value: < 2.2e-16MODELLO 2 DI REGRESSIONE LINEARE MULTIPLA
Costruiamo un secondo modello.
Il Modello 2 vede l’eliminazione della variabile “N.gravidanze” e le stime significative nel Modello 1 sono rimaste tali. Inoltre, R2 aggiustato è diminuito di un valore estremamente piccolo, quindi questo secondo modello è accettabile e preferibile al primo.
Di seguito sono eseguiti ulteriori test per validare questa operazione.
#eliminazione della variabile "N.gravidanze"
mod2 = update(mod1, ~.- N.gravidanze)
summary(mod2)##
## Call:
## lm(formula = Peso ~ Anni.madre + Gestazione + Lunghezza + Cranio +
## Sesso, data = data_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1163.03 -184.20 -14.07 163.24 2618.69
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6723.2290 141.3943 -47.550 < 2e-16 ***
## Anni.madre 1.8995 1.0692 1.777 0.0758 .
## Gestazione 32.2256 3.8205 8.435 < 2e-16 ***
## Lunghezza 10.2137 0.3008 33.954 < 2e-16 ***
## Cranio 10.6047 0.4261 24.887 < 2e-16 ***
## SessoM 78.6738 11.2182 7.013 2.99e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 275 on 2492 degrees of freedom
## Multiple R-squared: 0.7265, Adjusted R-squared: 0.7259
## F-statistic: 1324 on 5 and 2492 DF, p-value: < 2.2e-16SELEZIONE TRA MODELLO 1 E MODELLO 2
I test riportati di seguito vengono utilizzati come criteri di selezione tra più modelli.
#valutazione dei due modelli
Anova=anova(mod1, mod2)
Anova_values=round(Anova, 2)
kable(Anova_values, caption="Test ANOVA")| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 2491 | 187996688 | NA | NA | NA | NA |
| 2492 | 188424446 | -1 | -427758.5 | 5.67 | 0.02 |
BIC=BIC(mod1, mod2)
kable(BIC, caption="Test BIC")| df | BIC | |
|---|---|---|
| mod1 | 8 | 35200.87 |
| mod2 | 7 | 35198.72 |
#install.packages("car")
library(car)
Vif=vif(mod2)
Vif_values=round(Vif, 2)
kable(Vif_values, caption="Test di multicollinearità" )| x | |
|---|---|
| Anni.madre | 1.03 |
| Gestazione | 1.68 |
| Lunghezza | 2.07 |
| Cranio | 1.62 |
| Sesso | 1.04 |
Nonostante l’analisi della varianza (test ANOVA) suggerisca che l’inclusione della variabile “N.gravidanze” nel modello di regressione sia rilevante in termini di informazione aggiunta, il criterio di selezione BIC suggerisce di preferire il Modello 2. Si è anche verificato che il Modello 2 non presenta problemi di multicollinearità, in quanto tutti i valori sono di molto inferiori a 5.
Per quanto riguarda l’età della madre, che è precedentemente risultata una variabile non significativa, si decide comunque di mantenerla nel modello in quanto l’età può essere una variabile di controllo molto importante in studi medici di questo tipo.
STUDIO DI UN POSSIBILE EFFETTO NON LINEARE
Un’ultima analisi è stata fatta sulla variabile “Gestazione” in quanto dalla matrice di correlazione e dal grafico sottostante emerge un possibile effetto non lineare, la retta infatti tende a diventare una leggera curva. Per questo si costruiscono altri due modelli che verranno interpretati di seguito.
#analisi della variabile gestazione
attach(data_clean)
plot(Gestazione, Peso, pch=20)
mod3 = update(mod2, ~.+ I(Gestazione^2))
summary(mod3)##
## Call:
## lm(formula = Peso ~ Anni.madre + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2), data = data_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1158.58 -182.63 -12.34 164.36 2641.40
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4660.7331 899.5830 -5.181 2.38e-07 ***
## Anni.madre 1.9815 1.0688 1.854 0.0639 .
## Gestazione -83.0515 49.8022 -1.668 0.0955 .
## Lunghezza 10.3202 0.3040 33.945 < 2e-16 ***
## Cranio 10.7010 0.4278 25.017 < 2e-16 ***
## SessoM 76.4062 11.2509 6.791 1.39e-11 ***
## I(Gestazione^2) 1.5394 0.6631 2.322 0.0203 *
## ---
## 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.7264
## F-statistic: 1106 on 6 and 2491 DF, p-value: < 2.2e-16plot(Gestazione, Peso, pch=20)
mod4=update(mod3, ~.+ Gestazione*Anni.madre)
summary(mod4)##
## Call:
## lm(formula = Peso ~ Anni.madre + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2) + Anni.madre:Gestazione, data = data_clean)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1142.39 -183.17 -13.62 166.60 2643.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2578.6145 1138.9010 -2.264 0.02365 *
## Anni.madre -60.5207 21.0483 -2.875 0.00407 **
## Gestazione -143.9749 53.7803 -2.677 0.00748 **
## Lunghezza 10.2958 0.3037 33.905 < 2e-16 ***
## Cranio 10.7121 0.4271 25.081 < 2e-16 ***
## SessoM 76.4631 11.2332 6.807 1.25e-11 ***
## I(Gestazione^2) 1.7328 0.6652 2.605 0.00925 **
## Anni.madre:Gestazione 1.6112 0.5419 2.973 0.00297 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2490 degrees of freedom
## Multiple R-squared: 0.728, Adjusted R-squared: 0.7273
## F-statistic: 952.1 on 7 and 2490 DF, p-value: < 2.2e-16Nel Modello 3 si è aggiunto l’effetto quadratico della variabile, che è risultato essere significativo, dimostrando un cambiamento di correlazione tra la “Gestazione” e il “Peso”. Questo è facilmente intuibile poichè il periodo di gestazione è limitato, pùò variare da caso a caso ma comunque c’è un limite, non aumenterà all’infinito.
Per approfondire l’aspetto della gestazione, è stato creato un Modello 4 in cui si aggiunge in particolare l’interazione tra l’età materna e i mesi di gestazione. L’interazione risulta significativa ed implica che l’impatto dell’età della madre sul peso del neonato varia a seconda di quanto è durata la gravidanza.
Tutto ciò sta ad indicare che la relazione tra età materna, periodo di gestazione e peso del neonato è più complessa di quello che può sembrare e può dipendere da diversi fattori biologici.
A questo punto, valutiamo i tre modelli per cercare di capire quale possa essere il migliore.
#valutazione dei tre modelli
Anova=anova(mod2, mod3, mod4)
Anova_values=round(Anova, 2)
kable(Anova_values, caption="Test ANOVA")| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 2492 | 188424446 | NA | NA | NA | NA |
| 2491 | 188017653 | 1 | 406792.8 | 5.41 | 0.02 |
| 2490 | 187352482 | 1 | 665171.8 | 8.84 | 0.00 |
BIC=BIC(mod2, mod3, mod4)
kable(BIC, caption="Test BIC")| df | BIC | |
|---|---|---|
| mod2 | 7 | 35198.72 |
| mod3 | 8 | 35201.15 |
| mod4 | 9 | 35200.12 |
Si decide di scegliere come modello definitivo il Modello 2, in quanto, nonostante il Modello 4 sembri essere quello che spieghi una maggior variabilità del peso del neonato, è allo stesso tempo il modello più complesso. Infatti, sulla base del BIC il modello migliore risulta essere ancora il 2, in quanto offre un miglior equilibrio tra accuratezza e parsimonia.
ANALISI DEI RESIDUI
Infine, è importante fare un’analisi dei residui sia graficamente che mediante test specifici. Per residui si intendono gli errori del modello, che possono essere calcolati come distanza verticale di ogni punto dalla retta di regressione.
Vediamo di seguito dei grafici che ci permettono di osservare già la distribuzione dei residui e notiamo che ci sono 3 osservazioni, la 1551, la 155 e la 1306, che si discostano molto dal resto delle osservazioni. Successivamente vengono eseguiti test specifici, quali Shapiro-Wilk, Breusch-Pagan e Durbin-Watson.
#analisi dei residui
plot(mod2)
#normalità dei residui
Shapiro=shapiro.test(residuals(mod2))
results_Shapiro=tidy(Shapiro)
print(results)## # A tibble: 1 × 3
## statistic p.value method
## <dbl> <dbl> <chr>
## 1 0.971 3.38e-22 Shapiro-Wilk normality test#install.packages("lmtest")
library(lmtest)
Breusch.Pagan=bptest(mod2)
results_bp=tidy(Breusch.Pagan)
print(results_bp)## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 91.5 3.27e-18 5 studentized Breusch-Pagan testDurbin.Watson=dwtest(mod2)
results_dw=tidy(Durbin.Watson)
print(results_dw)## # A tibble: 1 × 4
## statistic p.value method alternative
## <dbl> <dbl> <chr> <chr>
## 1 1.95 0.114 Durbin-Watson test true autocorrelation is greater than 0Dai test emerge quanto segue:
si rifiuta l’ipotesi di Normalità:
si rifiuta l’ipotesi di omoschedasticità;
non si rifiuta l’ipotesi di incorrelazione.
Questi risultati possono essere interpretati mediante ulteriori test e l’osservazione grafica dei residui.In particolare, la non Normalità e l’eteroschedasticità potrebbero essere dovute ai valori di leva (anche se nessuna osservazione supera una distanza di Cook pari a 1, infatti la massima registrata è pari a 0.82) e a tre osservazioni risultate come outliers significativi, che corrispondono proprio alle tre osservazioni evidenziate precedentemente con l’analisi grafica. Analizzando queste osservazioni non risultano anomalie, ma probabilmente alcuni valori registrati differiscono dall’andamento generale delle osservazioni.
#leverage
leverage=hatvalues(mod2)
plot(leverage)
n=nrow(data_clean)
p=sum(leverage)
soglia=2*p/n
abline(h=soglia, col=2)
leverage[leverage>soglia]## 15 34 61 67 70 96
## 0.006353704 0.006346038 0.005215193 0.005392782 0.004872400 0.005921776
## 101 106 117 131 151 155
## 0.007235067 0.012870847 0.005026404 0.006987052 0.011314488 0.006671764
## 190 205 206 220 230 249
## 0.005512028 0.008596365 0.010364317 0.007874209 0.006087205 0.005082777
## 260 295 304 305 310 312
## 0.005549406 0.005124088 0.005169390 0.006956766 0.029128409 0.013394575
## 315 335 378 383 408 445
## 0.005420796 0.004891525 0.015686242 0.004840026 0.005965561 0.007241700
## 446 486 492 565 587 592
## 0.005482041 0.005747801 0.008414809 0.005022829 0.011684923 0.006332366
## 629 638 656 684 697 702
## 0.004818185 0.006627089 0.006560342 0.008757628 0.005987155 0.005173849
## 724 748 750 765 805 895
## 0.004910972 0.008721493 0.007007600 0.006301965 0.014329109 0.005222084
## 928 947 951 956 1014 1067
## 0.022713537 0.007988137 0.005210745 0.007709457 0.010078143 0.008234789
## 1072 1075 1083 1091 1096 1106
## 0.005005686 0.005952630 0.004941848 0.008940286 0.005075106 0.006619970
## 1130 1151 1166 1181 1188 1194
## 0.007075810 0.005097453 0.005284643 0.005596979 0.006916013 0.005529846
## 1200 1238 1248 1273 1293 1294
## 0.005882067 0.007408052 0.014918330 0.007153963 0.005556821 0.005047183
## 1311 1323 1356 1357 1385 1400
## 0.005549364 0.006013916 0.005412594 0.007601657 0.012049556 0.006488419
## 1402 1420 1428 1429 1551 1553
## 0.005319573 0.005210991 0.008229676 0.022531542 0.048824756 0.006815842
## 1556 1560 1593 1610 1619 1639
## 0.007707020 0.006152994 0.006900067 0.008634450 0.014507340 0.005207855
## 1686 1692 1693 1701 1712 1718
## 0.009008178 0.007106397 0.005054131 0.011405592 0.006979971 0.004807329
## 1735 1743 1780 1806 1809 1827
## 0.007653394 0.004845554 0.026427246 0.006271517 0.008826787 0.006043764
## 1874 1893 1977 2037 2040 2089
## 0.004954263 0.005121269 0.006385291 0.006402534 0.011549071 0.006298709
## 2098 2114 2115 2120 2140 2149
## 0.004859409 0.013351427 0.012109050 0.018008839 0.006092754 0.014006033
## 2175 2200 2215 2216 2224 2225
## 0.021822275 0.011073897 0.004944838 0.008824245 0.007426804 0.005440094
## 2257 2307 2318 2333 2359 2391
## 0.007402116 0.014277773 0.006043710 0.004825852 0.005183660 0.005453547
## 2408 2437 2452 2458 2476 2478
## 0.010205468 0.024057599 0.023455443 0.008005461 0.004937555 0.005782921#outliers
plot(rstudent(mod2))
car::outlierTest(mod2)## rstudent unadjusted p-value Bonferroni p
## 1551 9.955057 6.3845e-23 1.5948e-19
## 155 4.950167 7.9105e-07 1.9760e-03
## 1306 4.812976 1.5759e-06 3.9366e-03#distanza di Cook
cook=cooks.distance(mod2)
plot(cook)
max.cook=round(max(cook), 2)
kable(max.cook, col.names="Distanza di Cook")| Distanza di Cook |
|---|
| 0.82 |
BONTA’ DEL MODELLO E POSSIBILI PREVISIONI
Considerando quindi il Modello 2 un buon modello in base a R2 (circa 0.73), possiamo utilizzarlo per fare previsioni. Ad esempio, sono state fatte di seguito tre previsioni per vedere come cambia il peso del neonato alla 39esima settimana di gestazione, al variare dell’età materna.
- 39 settimane di gestazione, età materna=30, sesso del neonato=femmina
#previsioni
attach(data_clean)
predict(mod2, newdata=data.frame(Gestazione=39,
Anni.madre=30,
Lunghezza=mean(Lunghezza),
Cranio=mean(Cranio),
Sesso="F"))## 1
## 3249.14- 39 settimane di gestazione, età materna=22, sesso del neonato=femmina
predict(mod2, newdata=data.frame(Gestazione=39,
Anni.madre=22,
Lunghezza=mean(Lunghezza),
Cranio=mean(Cranio),
Sesso="F"))## 1
## 3233.943- 39 settimane di gestazione, età materna=44, sesso del neonato=femmina
predict(mod2, newdata=data.frame(Gestazione=39,
Anni.madre=44,
Lunghezza=mean(Lunghezza),
Cranio=mean(Cranio),
Sesso="F"))## 1
## 3275.733RAPPRESENTAZIONI GRAFICHE E CONCLUSIONI
Vediamo infine alcune rappresentazioni grafiche che mostrano alcune relazioni specifiche del nostro modello.
#rappresentazioni grafiche
library(ggplot2)
ggplot(data=data_clean)+
geom_point(aes(x=Gestazione,
y=Peso,
col=Sesso), position="jitter")+
geom_smooth(aes(x=Gestazione,
y=Peso,
col=Sesso), se=F, method="lm")
ggplot(data=data_clean)+
geom_point(aes(x=Anni.madre,
y=Peso,
col=Sesso), position="jitter")+
geom_smooth(aes(x=Anni.madre,
y=Peso,
col=Sesso), se=F, method="lm")
ggplot(data=data_clean)+
geom_boxplot(aes(x=factor(Fumatrici),
y=Peso,
fill=Sesso))+
labs(x="Fumatrici(0=Non fumatrice, 1=Fumatrice)", y="Peso", col="Sesso", fill="Sesso")
Ad esempio, è evidente la correlazione positiva tra le settimane di gestazione e il peso del neonato. L’età materna invece presa singolarmente sembra non influire sul peso del neonato, il grafico mostra infatti una retta orizzontale, leggermente decrescente per le femmine. Anche il fumo sembra non avere un effetto diretto sul peso del neonato.
Quindi questa analisi statistica dimostra che gli unici aspetti che direttamente sembrano influire sul peso dei neonati sono le settimane di gestazione e ovviamente le misure antropometriche. E’ importante tenere presente però che i fenomeni biologici possono essere influenzati da vari fattori, quindi non si può escludere categoricamente che l’età materna o il fumo non possano avere un’incidenza in determinati casi.