Progetto di Statistica Inferenziale
Lavinia Solimani
23/04/2026
Questo progetto ha come obiettivo la creazione di un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali. Lo scopo è di migliorare la gestione delle gravidanze ad alto rischio, ottimizzare le risorse ospedaliere e garantire migliori risultati per la salute neonatale.
1. Raccolta dei Dati e Struttura del Dataset
newborn_data <- read.csv("neonati.csv",
sep=",",
stringsAsFactors = T)
newborn_data$Fumatrici <- factor(newborn_data$Fumatrici,
levels = c(0, 1))
glimpse(newborn_data)## Rows: 2,500
## Columns: 10
## $ Anni.madre <int> 26, 21, 34, 28, 20, 32, 26, 25, 22, 23, 29, 21, 36, 24, 3…
## $ N.gravidanze <int> 0, 2, 3, 1, 0, 0, 1, 0, 1, 0, 2, 2, 5, 0, 3, 2, 2, 3, 0, …
## $ Fumatrici <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ Gestazione <int> 42, 39, 38, 41, 38, 40, 39, 40, 40, 41, 38, 40, 38, 40, 3…
## $ Peso <int> 3380, 3150, 3640, 3690, 3700, 3200, 3100, 3580, 3670, 370…
## $ Lunghezza <int> 490, 490, 500, 515, 480, 495, 480, 510, 500, 510, 480, 51…
## $ Cranio <int> 325, 345, 375, 365, 335, 340, 345, 349, 335, 362, 330, 34…
## $ Tipo.parto <fct> Nat, Nat, Nat, Nat, Nat, Nat, Nat, Nat, Ces, Ces, Ces, Na…
## $ Ospedale <fct> osp3, osp1, osp2, osp2, osp3, osp2, osp3, osp1, osp2, osp…
## $ Sesso <fct> M, F, M, M, F, F, F, M, F, F, M, F, F, F, M, M, M, M, F, …Il dataset è costituito da 2500 osservazioni delle seguenti 10 variabili.
- Anni.madre: misura dell’età in anni (variabile quantitativa discreta);
- N.gravidanze: quante gravidanze ha avuto la madre (variabile quantitativa discreta);
- Fumatrici: indicatore se la madre è non fumatrice (0) o fumatrice (1) (variabile qualitativa codificata in dummy);
- Gestazione: numero di settimane di gestazione (variabile quantitativa continua);
- Peso: peso alla nascita in grammi (variabile quantitativa continua);
- Lunghezza: lunghezza alla nascita in millimetri (variabile quantitativa continua);
- Cranio: diametro craniale alla nascita in millimetri (variabile quantitativa continua);
- Tipo.parto: tipo di parto, se naturale o cesareo (variabile qualitativa su scala nominale);
- Ospedale: ospedale di nascita tra 1, 2 e 3 (variabile qualitativa su scala nominale);
- Sesso: sesso del neonato (variabile qualitativa su scala nominale).
L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature. Pertanto, la variabile risposta è Peso, mentre le altre sono le variabili esplicative.
Definiamo la variabile N, che descriverà la numerosità del campione, e visualizziamo le statistiche di sintesi delle variabili.
attach(newborn_data)
N <- nrow(newborn_data)
pander(summary(newborn_data))| Anni.madre | N.gravidanze | Fumatrici | Gestazione | Peso |
|---|---|---|---|---|
| Min. : 0.00 | Min. : 0.0000 | 0:2396 | Min. :25.00 | Min. : 830 |
| 1st Qu.:25.00 | 1st Qu.: 0.0000 | 1: 104 | 1st Qu.:38.00 | 1st Qu.:2990 |
| Median :28.00 | Median : 1.0000 | NA | Median :39.00 | Median :3300 |
| Mean :28.16 | Mean : 0.9812 | NA | Mean :38.98 | Mean :3284 |
| 3rd Qu.:32.00 | 3rd Qu.: 1.0000 | NA | 3rd Qu.:40.00 | 3rd Qu.:3620 |
| Max. :46.00 | Max. :12.0000 | NA | Max. :43.00 | Max. :4930 |
| Lunghezza | Cranio | Tipo.parto | Ospedale | Sesso |
|---|---|---|---|---|
| Min. :310.0 | Min. :235 | Ces: 728 | osp1:816 | F:1256 |
| 1st Qu.:480.0 | 1st Qu.:330 | Nat:1772 | osp2:849 | M:1244 |
| Median :500.0 | Median :340 | NA | osp3:835 | NA |
| Mean :494.7 | Mean :340 | NA | NA | NA |
| 3rd Qu.:510.0 | 3rd Qu.:350 | NA | NA | NA |
| Max. :565.0 | Max. :390 | NA | NA | NA |
2. Analisi e Modellizzazione
Analisi Preliminare
Innanzitutto, verifichiamo che la distribuzione della variabile risposta sia approssimativamente normale, ricavandone gli indici di forma ed effettuando un test di Shapiro-Wilk.
round(moments::skewness(Peso), digits = 2)## [1] -0.65round(moments::kurtosis(Peso)-3, digits = 2)## [1] 2.03shapiro.test(Peso)##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16Con un indice di asimmetria di Fisher di -0.65 e un coefficiente di curtosi di 2.03, la distribuzione è, rispettivamente, leggermente asimmetrica negativa e leptocurtica. Pertanto, la distribuzione è un po’ più allungata rispetto alla normale, e sono più frequenti modalità alte.
Per quanto riguarda il test di Shapiro-Wilk, la statistica W riporta un valore di 0.97. Essendo il valore molto vicino a 1, la distribuzione risulta molto simile a una normale. Tuttavia, con un p-value molto vicino allo zero, l’ipotesi nulla di normalità della distribuzione andrebbe rifiutata.
Considerato, però, che il test di Shapiro-Wilk diventa più sensibile all’aumentare della numerosità del campione, e il dataset in questione è costituito da 2500 osservazioni, è opportuno appurare graficamente se la distribuzione della variabile risposta è sufficientemente normale da poter comunque condurre l’analisi con un modello di regressione lineare.
Rappresentiamo graficamente la funzione di densità della variabile Peso.
ggplot(newborn_data, aes(x = Peso))+
geom_density(fill = "lightblue")+
geom_vline(xintercept = mean(Peso),
color = "red",
linetype = "dashed")+
labs(title = "Funzione di densità della variabile Peso",
x = "Peso (g)",
y = "Densità")+
theme_minimal()
Coerentemente con il valore ottenuto per l’indice di asimmetria di Fisher e il coefficiente di curtosi, la distribuzione risulta avere una coda allungata a sinistra ed essere un po’ più allungata rispetto a una normale. Tuttavia, la forma in generale è molto simile alla classica “campana” della gaussiana.
Ora visualizziamo un Q-Q plot, che confronta i quantili dei dati della variabile Peso con i quantili di una distribuzione normale.
ggplot(newborn_data, aes(sample = Peso))+
geom_qq()+
geom_qq_line(color = "red")+
labs(title = "Q-Q plot della variabile Peso",
x = "Quantili di riferimento",
y = "Quantili del campione")+
theme_minimal()
Da questo grafico si evince chiaramente il problema di questa distribuzione, che si concentra nelle osservazioni più basse. Nella parte sinistra del grafico, infatti, i punti giacciono al di sotto della retta, mentre si dispongo correttamente sopra di essa nel resto del grafico.
Alla luce di questa breve analisi statistica e grafica, possiamo ritenere che la distribuzione della variabile Peso sia sufficientemente normale da poter utilizzare un modello di regressione lineare.
Ora visualizziamo la matrice di correlazione, con lo scopo di identificare eventuali variabili molto correlate alla variabile risposta, che sono quindi utili ai fini dello studio. Allo stesso tempo, la matrice ci permette di verificare la presenza di multicollinearità, ovvero se ci sono regressori molto correlati tra loro che, quindi, possono causare problemi nel modello.
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr = usr))
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 <- 1.2
text(0.5, 0.5, txt, cex = cex.cor)
}
pairs(newborn_data, upper.panel = panel.smooth, lower.panel = panel.cor)
Prendiamo intanto in esame la variabile risposta Peso e vediamo con quali variabili esplicative è maggiormente correlata:
- Le variabili che presentano una forte correlazione con la variabile Peso sono Lunghezza, con un coefficiente di correlazione lineare di 0.80, Cranio (0.70) e Gestazione (0.59). Il coefficiente è positivo in tutti questi casi, indicando quindi una concordanza tra queste variabili e la variabile risposta. Anche nel caso della variabile Sesso c’è una correlazione positiva, ma è molto più bassa (0.24);
- La correlazione delle variabili Anni.madre, N.gravidanze, Fumatrici, Tipo.parto e Ospedale con la variabile risposta, invece, è praticamente nulla.
Si può notare una leggera non linearità nell’interazione delle variabili Gestazione, Lunghezza e Cranio con la variabile risposta, che potremo esplorare in fase di creazione del modello di regressione.
Ora verifichiamo la presenza o meno di multicollinearità tra i regressori:
- Le variabili Anni.madre e N.gravidanze presentano un coefficiente di correlazione di 0.38, quindi non eccessivamente alto;
- Le variabili Gestazione, Lunghezza e Cranio risultano abbastanza correlate tra di loro, con coefficienti che vanno da 0.46 a 0.62.
Per analizzare graficamente le variabili qualitative, invece, utilizziamo i boxplot.
Variabile Fumatrici
ggplot(newborn_data, aes(x = Fumatrici,
y = Peso))+
geom_boxplot(fill = "steelblue",
alpha = 0.7)+
labs(title = "Peso del neonato in assenza e presenza di fumo materno",
x = "Fumatrici",
y = "Peso (g)")+
theme_minimal()+
theme(legend.position = "none")
Dai boxplot si nota una leggera differenza nelle distribuzioni della variabile risposta Peso condizionate per le modalità della variabile esplicativa Fumatrici. In particolare, il peso dei neonati risulta essere leggermente più basso in media, se la madre è fumatrice. Tuttavia, è necessario verificarlo tramite un test statistico t, per testare l’ipotesi di uguaglianza tra medie per gruppi indipendenti. Un altra osservazione che si può fare da questo grafico è la presenza di molti outliers, soprattutto di neonati con un peso più basso in assenza di fumo materno. Questa concentrazione di outliers nella parte inferiore della distribuzione è coerente con quanto rilevato nel Q-Q plot.
t.test(Peso ~ Fumatrici)##
## Welch Two Sample t-test
##
## data: Peso by Fumatrici
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -45.61354 145.22674
## sample estimates:
## mean in group 0 mean in group 1
## 3286.153 3236.346Con un p-value di circa 0.30, quindi maggiore di 0.05, non si rifiuta l’ipotesi nulla di uguaglianza tra medie. Pertanto, nonostante si noti una leggera differenza tra le due distribuzioni a livello grafico, non si tratta di una differenza significativa. Potenzialmente, quindi, il fatto che la madre sia o meno fumatrice non incide sul peso del neonato.
Variabile Tipo.parto
ggplot(newborn_data, aes(x = Tipo.parto,
y = Peso))+
geom_boxplot(fill = "steelblue",
alpha = 0.7)+
labs(title = "Peso del neonato in caso di parto cesareo e naturale",
y = "Peso (g)")+
theme_minimal()+
theme(legend.position = "none")
Visualizzando le distribuzioni della variabile Peso condizionate dalle modalità della variabile Tipo.parto, notiamo che esse risultano quasi identiche. A prima vista, quindi, il tipo di parto non sembra avere influenza sul peso del neonato, ma eseguiamo anche un t-test per verificare l’ipotesi di uguaglianza tra medie per gruppi indipendenti. Ad ogni modo, anche in questo grafico constatiamo la presenza di diversi outliers nella parte inferiore della distribuzione, soprattutto nel caso di parto naturale.
t.test(Peso ~ Tipo.parto)##
## Welch Two Sample t-test
##
## data: Peso by Tipo.parto
## t = -0.12968, df = 1493, p-value = 0.8968
## alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
## 95 percent confidence interval:
## -46.27992 40.54037
## sample estimates:
## mean in group Ces mean in group Nat
## 3282.047 3284.916Con un p-value di circa 0.90, non si rifiuta l’ipotesi nulla e, pertanto, possiamo confermare l’ipotesi di uguaglianza tra le medie. Probabilmente, il coefficiente di regressione per questa variabile non sarà significativo e, quindi, non occorrerà inserire la variabile nel modello.
Variabile Ospedale
ggplot(newborn_data, aes(x = Ospedale,
y = Peso))+
geom_boxplot(fill = "steelblue",
alpha = 0.7)+
labs(title = "Peso del neonato in base all'ospedale di nascita",
x = "Ospedale",
y = "Peso (g)")+
theme_minimal()+
theme(legend.position = "none")
Anche nel caso della variabile Ospedale, le distribuzioni della variabile risposta condizionate dalle modalità della variabile esplicativa risultano quasi identiche. Anche in questo grafico notiamo un notevole numero di outliers, sempre nella parte inferiore della distribuzione, per tutti e tre gli ospedali.
Eseguiamo anche il t-test per verificare l’ipotesi di uguaglianza tra medie per gruppi indipendenti.
pairwise.t.test(Peso, Ospedale,
paired = F,
pool.sd = T,
p.adjust.method = "bonferroni")##
## Pairwise comparisons using t tests with pooled SD
##
## data: Peso and Ospedale
##
## osp1 osp2
## osp2 1.00 -
## osp3 0.33 0.33
##
## P value adjustment method: bonferroniIn tutti i t-test tra le varie combinazioni di modalità, il p-value risulta decisamente maggiore di 0.05. Pertanto, anche in questo caso non si rifiuta l’ipotesi nulla e non constatiamo alcuna differenza significativa tra gli ospedali sul peso dei neonati. Probabilmente, quindi, neanche questa variabile sarà necessaria all’interno del modello.
Variabile Sesso
ggplot(newborn_data, aes(x = Sesso,
y = Peso))+
geom_boxplot(fill = "steelblue",
alpha = 0.7)+
labs(title = "Peso del neonato in base al sesso",
x = "Sesso",
y = "Peso (g)")+
theme_minimal()+
theme(legend.position = "none")
In questo caso, notiamo che il boxplot della distribuzione condizionata per il sesso maschile è leggermente più alta nel grafico, e ciò significa che il peso dei neonati è mediamente maggiore se sono maschi. Verifichiamo l’ipotesi nulla di uguaglianza tra medie di gruppi indipendenti con il t-test. Da notare la costante presenza di outliers, maggiormente concentrati nella parte inferiore della distribuzione.
t.test(Peso ~ Sesso)##
## Welch Two Sample t-test
##
## data: Peso by Sesso
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -287.1051 -207.0615
## sample estimates:
## mean in group F mean in group M
## 3161.132 3408.215Con un p-value quasi pari a zero, si rifiuta l’ipotesi nulla di uguaglianza tra medie e si conclude, quindi, che le due medie sono significativamente diverse. Pertanto, sarà opportuno inserire questa variabile nel modello di regressione lineare.
Verifichiamo ora le seguenti ipotesi relative al dataset:
- In alcuni ospedali si fanno più parti cesarei.
- La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione.
- Le misure antropometriche sono significativamente diverse tra i due sessi.
Prima ipotesi
Testiamo l’ipotesi che in alcuni ospedali si fanno più parti cesarei, quindi consideriamo le variabili Tipo.parto e Ospedale. Eseguiremo un test chi-quadrato di Pearson, per testare l’ipotesi nulla di indipendenza tra le due variabili.
\[ \begin{cases} H_0: \text{Indipendenza} \\ H_1: \text{Associazione} \end{cases} \]
Innanzitutto, creiamo una tabella di contingenza con le osservazioni delle due variabili, e la rappresentiamo graficamente per identificare eventuali pattern e, quindi, possibili correlazioni.
table_birth_hospital <- table(Tipo.parto, Ospedale)
ggballoonplot(data = as.data.frame(table_birth_hospital),
fill = "steelblue")+
labs(title = "Frequenza del tipo di parto per ospedale")
Dal grafico si nota una leggera maggior frequenza di parti cesarei nell’ospedale 2, ma non sembra esserci un pattern evidente. Ad ogni modo, ora eseguiremo il test chi-quadrato per verificare se ci sia effettivamente un’associazione tra il tipo di parto e l’ospedale di nascita.
chisq.test(table_birth_hospital)##
## Pearson's Chi-squared test
##
## data: table_birth_hospital
## X-squared = 1.0972, df = 2, p-value = 0.5778Con un p-value di 0.58, non si rifiuta l’ipotesi nulla e, pertanto, non si può dire che in alcuni ospedali si facciano più parti cesarei.
Seconda ipotesi
Testiamo l’ipotesi che le medie stimate per il peso e la lunghezza in questo campione di neonati sono significativamente uguali a quelle della popolazione.
Considerato che la numerosità del campione è grande, ma non conosciamo la varianza della popolazione, per testare quest’ipotesi procederemo con un t-test, che è più prudente nel rifiuto dell’ipotesi nulla rispetto a un test Z. Nello specifico, dovremo eseguire due test, uno per la variabile Peso e uno per la variabile Lunghezza.
Mediamente, i neonati pesano circa 3300 grammi e misurano 50 centimetri di lunghezza, pertanto useremo questi valori per la media del peso e della lunghezza dei neonati nella popolazione, che indicheremo con \(\mu\). Le medie del nostro campione, invece, che sono pari a 3284 grammi e 494.7 millimetri (quindi circa 49.5 centimetri), saranno indicate con \(\hat{\mu}\).
\[ \left\{ \begin{aligned} H_0&: \hat{\mu} = \mu \\ H_1&: \hat{\mu} \neq \mu \end{aligned} \right. \]
Eseguiamo il t-test sulla variabile Peso, impostando un livello di confidenza al 95%.
t.test(Peso,
mu = 3300,
conf.level = 0.95,
alternative = "two.sided")##
## One Sample t-test
##
## data: Peso
## t = -1.516, df = 2499, p-value = 0.1296
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
## 3263.490 3304.672
## sample estimates:
## mean of x
## 3284.081La statistica test t ha un valore di -1.516, in una distribuzione t di Student con n-1 = 2499 gradi di libertà, con un p-value di circa 0.13 che è, quindi, maggiore del valore soglia di 0.05. Quindi, non si rifiuta l’ipotesi nulla e si può confermare che la media del peso di questo campione di neonati è significativamente uguale a quella della popolazione.
Ora eseguiamo il t-test sulla variabile Lunghezza, impostando anche in questo caso un livello di confidenza al 95%.
t.test(Lunghezza,
mu = 500,
conf.level = 0.95,
alternative = "two.sided")##
## One Sample t-test
##
## data: Lunghezza
## t = -10.084, df = 2499, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
## 493.6598 495.7242
## sample estimates:
## mean of x
## 494.692In questo caso, con un p-value estremamente vicino allo zero, si rifiuta l’ipotesi nulla di uguaglianza significativa della media del campione alla media della popolazione. In altre parole, la lunghezza media dei neonati in questo campione è significativamente diversa dalla lunghezza media della popolazione di riferimento. Si può notare, infatti, che l’intervallo di confidenza è tra 493.7 e 495.7 millimetri, quindi il valore di 500 della media della popolazione non rientra nell’intervallo.
Terza ipotesi
Testiamo, infine, l’ipotesi che le misure antropometriche sono significativamente diverse tra i due sessi. Le variabili da prendere in esame, quindi, sono Peso, Lunghezza e Cranio, condizionate per la variabile Sesso.
\[ \left\{ \begin{aligned} H_0&: \mu_M = \mu_F \\ H_1&: \mu_M \neq \mu_F \end{aligned} \right. \]
Abbiamo già verificato con un t-test che il peso del neonato è significativamente diverso tra i due sessi. Ora, quindi, procederemo con la verifica dell’ipotesi per la variabile Lunghezza, rappresentando graficamente la distribuzione di questa variabile rispetto al sesso maschile e al sesso femminile, tramite un boxplot condizionato.
ggplot(newborn_data, aes(x = Sesso,
y = Lunghezza))+
geom_boxplot(fill = "steelblue",
alpha = 0.7)+
labs(title = "Lunghezza del neonato in base al sesso",
x = "Sesso",
y = "Lunghezza (mm)")+
theme_minimal()+
theme(legend.position = "none")
Notiamo subito come entrambe le scatole nel grafico siano molto sottili, indicando quindi che il 50% dei dati che si trova tra il primo e il terzo quartile, ovvero il corpo centrale della distribuzione, sia molto concentrato su pochi valori attorno ai 50 centimetri. Tuttavia, si notano anche molti outliers, quindi osservazioni che si discostano molto dal corpo centrale della distribuzione. Questo è coerente con quanto mostrato nei precedenti boxplot, dal momento che sono state rilevate molte osservazioni che riportavano un peso più basso nella distribuzione, e considerato che il peso e la lunghezza sono variabili positivamente correlate.
Dal boxplot risulta esserci una leggera differenza tra i due sessi; nello specifico, i neonati presentano una lunghezza maggiore rispetto alle neonate. Verifichiamo l’ipotesi con un test statistico t per confronti tra medie di gruppi indipendenti.
t.test(Lunghezza ~ Sesso)##
## Welch Two Sample t-test
##
## data: Lunghezza by Sesso
## t = -9.582, df = 2459.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -11.929470 -7.876273
## sample estimates:
## mean in group F mean in group M
## 489.7643 499.6672Con un p-value estremamente vicino allo zero, si rifiuta l’ipotesi nulla di uguaglianza tra medie della variabile Lunghezza condizionate per Sesso. Pertanto, la media della lunghezza per le neonate, pari a circa 49.98 centimetri, è significativamente inferiore alla media della lunghezza per i neonati, pari a circa 49.97 centimetri.
Concludiamo con l’ultima misura antropometrica, ovvero il diametro del cranio. Rappresentiamo, quindi, la variabile Cranio condizionata per Sesso.
ggplot(newborn_data, aes(x = Sesso,
y = Cranio))+
geom_boxplot(fill = "steelblue",
alpha = 0.7)+
labs(title = "Diametro craniale del neonato in base al sesso",
x = "Sesso",
y = "Diametro del cranio (mm)")+
theme_minimal()+
theme(legend.position = "none")
La distribuzione condizionata per sesso maschile risulta leggermente più alta rispetto alla distribuzione condizionata per sesso femminile. Anche in questo caso, coerentemente col fatto che Cranio è positivamente correlata sia alla variabile risposta Peso che alla variabile esplicativa Lunghezza, risultano esserci molti outliers nella parte inferiore della distribuzione. Verifichiamo se le medie delle due distribuzioni sono significativamente diverse con un test statistico t.
t.test(Cranio ~ Sesso)##
## Welch Two Sample t-test
##
## data: Cranio by Sesso
## t = -7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -6.089912 -3.541270
## sample estimates:
## mean in group F mean in group M
## 337.6330 342.4486Con un p-value molto vicino allo zero, si rifiuta l’ipotesi nulla di uguaglianza tra medie per gruppi indipendenti e si può confermare, quindi, che la media del diametro craniale delle neonate stimata per questo campione (33.76 centimetri) è significativamente inferiore alla media del diametro craniale dei neonati, sempre stimata per questo campione (34.24 centimetri).
In conclusione, le misure antropometriche sono significativamente diverse tra i due sessi. Nello specifico, nelle neonate si rilevano misure più basse rispetto ai neonati.
Creazione del Modello di Regressione
Per il primo modello di regressione lineare multipla, includiamo tutte le variabili del dataset.
mod1 <- lm(Peso ~ ., data = newborn_data)
summary(mod1)##
## Call:
## lm(formula = Peso ~ ., data = newborn_data)
##
## 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 *
## Fumatrici1 -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-16Analizziamo le variabili e i coefficienti: - Il numero di gravidanze è positivo e significativo. Per ogni gravidanza aggiuntiva della madre, i neonati pesano circa 11.3 grammi in più; - La durata della gestazione è positiva ed estremamente significativa, a favore dell’ipotesi oggetto di studio. Per ogni settimana di gestazione in più, i neonati pesano circa 32.6 grammi in più; - La lunghezza e diametro del cranio sono positivi ed estremamente significativi, con un coefficiente di 10.29 e 10.47, rispettivamente; - L’età della madre non è una variabile significativa; - La variabile dummy sul fumo materno presenta un coefficiente negativo, il che significa che, se la madre è fumatrice, il neonato pesa in media 30.1 grammi in meno. Tuttavia, contrariamente all’ipotesi oggetto di studio, la variabile non risulta essere significativa; - La variabile dummy sulla tipologia di parto è positiva e significativa. I neonati che nascono con parto naturale pesano mediamente 29.5 grammi in più rispetto ai neonati che nascono con parto cesareo; - La variabile dummy sull’ospedale di nascita è positiva e significativa solamente nel caso in cui si confronta l’ospedale 3 con l’ospedale 1. In media, i neonati nati nell’ospedale 3 pesano circa 28.1 grammi in più; - La variabile dummy sul sesso del neonato è positiva ed estremamente significativa. Un neonato pesa, mediamente, 77.5 grammi in più rispetto a una neonata.
Questo modello presenta un \(R^2\) aggiustato di 0.7278, quindi una variabilità spiegata del 72.8%. Proviamo ora a creare dei modelli alternativi, per selezionare il modello migliore.
Selezione del Modello Ottimale
Proviamo a rimuovere dal modello le variabili che non sono risultate significative.
mod2 <- update(mod1, ~ . - (Anni.madre + Fumatrici))
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso, data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1113.18 -181.16 -16.58 161.01 2620.19
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6707.4293 135.9438 -49.340 < 2e-16 ***
## N.gravidanze 12.3619 4.3325 2.853 0.00436 **
## Gestazione 31.9909 3.7896 8.442 < 2e-16 ***
## Lunghezza 10.3086 0.3004 34.316 < 2e-16 ***
## Cranio 10.4922 0.4254 24.661 < 2e-16 ***
## Tipo.partoNat 29.2803 12.0817 2.424 0.01544 *
## Ospedaleosp2 -11.0227 13.4363 -0.820 0.41209
## Ospedaleosp3 28.6408 13.4886 2.123 0.03382 *
## SessoM 77.4412 11.1756 6.930 5.36e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared: 0.7287, Adjusted R-squared: 0.7278
## F-statistic: 836.3 on 8 and 2491 DF, p-value: < 2.2e-16Notiamo che la variabile relativa al numero di gravidanze ora è molto significativa, mentre la significatività delle altre variabili è rimasta la stessa. Inoltre, dato che l’\(R^2\) aggiustato è rimasto invariato, questo modello è a prima vista quello da preferire, avendo meno parametri ed essendo, quindi, più semplice.
Considerato il coefficiente di correlazione lineare tra le due variabili relative alle misure antropometriche, pari a 0.60, possiamo provare a togliere una delle due variabili dal modello. Iniziamo rimuovendo la variabile Cranio.
mod3 <- update(mod2, ~ . - Cranio)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Tipo.parto +
## Ospedale + Sesso, data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1114.3 -199.3 -21.1 188.7 3593.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5356.0469 138.7415 -38.605 < 2e-16 ***
## N.gravidanze 23.7416 4.8041 4.942 8.25e-07 ***
## Gestazione 45.7346 4.1802 10.941 < 2e-16 ***
## Lunghezza 13.6620 0.2987 45.735 < 2e-16 ***
## Tipo.partoNat 38.9531 13.4664 2.893 0.00385 **
## Ospedaleosp2 -14.8410 14.9831 -0.991 0.32202
## Ospedaleosp3 28.8595 15.0424 1.919 0.05516 .
## SessoM 87.2050 12.4551 7.002 3.24e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 305.5 on 2492 degrees of freedom
## Multiple R-squared: 0.6625, Adjusted R-squared: 0.6615
## F-statistic: 698.7 on 7 and 2492 DF, p-value: < 2.2e-16In questo caso, l’\(R^2\) aggiustato si è notevolmente abbassato, con una varianza spiegata del 66.1%. Pertanto, è opportuno mantenere la variabile Cranio nel modello.
Proviamo, invece, a rimuovere la variabile Lunghezza.
mod4 <- update(mod2, ~ . - Lunghezza)
summary(mod4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Cranio + Tipo.parto +
## Ospedale + Sesso, data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1340.2 -218.6 -11.9 211.6 1552.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.237e+03 1.641e+02 -38.006 <2e-16 ***
## N.gravidanze 4.903e+00 5.250e+00 0.934 0.350
## Gestazione 9.307e+01 4.059e+00 22.926 <2e-16 ***
## Cranio 1.710e+01 4.603e-01 37.151 <2e-16 ***
## Tipo.partoNat 8.442e+00 1.464e+01 0.577 0.564
## Ospedaleosp2 -8.718e-02 1.630e+01 -0.005 0.996
## Ospedaleosp3 2.683e+01 1.637e+01 1.639 0.101
## SessoM 1.179e+02 1.348e+01 8.745 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 332.3 on 2492 degrees of freedom
## Multiple R-squared: 0.6004, Adjusted R-squared: 0.5993
## F-statistic: 535 on 7 and 2492 DF, p-value: < 2.2e-16Si osserva immediatamente come la variabile Lunghezza sia essenziale nel modello. Infatti, rimuovendola, le variabili legate al numero di gravidanze e al tipo di parto perdono completamente la loro significatività, e l’\(R^2\) aggiustato scende addirittura a 0.60.
Proviamo, ora, ad approfondire il leggero effetto non lineare che risulta esserci per la variabile Gestazione osservando il diagramma di dispersione. Aggiungiamo, quindi, una nuova variabile al modello, ovvero l’effetto quadratico delle settimane di gestazione.
mod5 <- update(mod2, ~ . + I(Gestazione^2))
summary(mod5)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + I(Gestazione^2), data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1144.82 -182.58 -13.38 163.47 2641.94
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4804.2152 896.9686 -5.356 9.28e-08 ***
## N.gravidanze 12.4546 4.3296 2.877 0.00405 **
## Gestazione -74.2441 49.6348 -1.496 0.13483
## Lunghezza 10.4056 0.3036 34.278 < 2e-16 ***
## Cranio 10.5832 0.4273 24.770 < 2e-16 ***
## Tipo.partoNat 28.8495 12.0747 2.389 0.01696 *
## Ospedaleosp2 -10.2764 13.4311 -0.765 0.44428
## Ospedaleosp3 28.3093 13.4797 2.100 0.03582 *
## SessoM 75.3932 11.2082 6.727 2.15e-11 ***
## I(Gestazione^2) 1.4183 0.6607 2.147 0.03192 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.7 on 2490 degrees of freedom
## Multiple R-squared: 0.7292, Adjusted R-squared: 0.7282
## F-statistic: 745 on 9 and 2490 DF, p-value: < 2.2e-16Questo modello presenta un \(R^2\) aggiustato di 0.7282, quindi impercettibilmente maggiore del modello 2. Tuttavia, ora l’effetto principale delle settimane di gestazione non sono più significative, mentre è significativo l’effetto quadratico. Per appurare se il leggerissimo aumento di variabilità spiegata giustifica l’aggiunta di una variabile al modello, eseguiamo un test ANOVA, un test statistico che rapporta le varianze spiegate dai due modelli e informa se ci sia stato un aumento o diminuzione significativo.
anova(mod5, mod2)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso + I(Gestazione^2)
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2490 186554769
## 2 2491 186899996 -1 -345227 4.6078 0.03192 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Con un p-value di 0.03, quindi minore di 0.05, si rileva che l’aumento di varianza spiegata a seguito dell’aggiunta dell’effetto quadratico delle settimane di gestazione è significativo. Pertanto, anche il modello 5 è un buon candidato per la selezione del modello ottimale. Successivamente utilizzeremo anche altri criteri di valutazione, ovvero il Criterio di Informazione di Akaike (AIC) e il Criterio di informazione Bayesiano (BIC).
Ora, ripartendo dal modello 2, proviamo ad aggiungere l’effetto quadratico della variabile Lunghezza.
mod6 <- update(mod2, ~ . + I(Lunghezza^2))
summary(mod6)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + I(Lunghezza^2), data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1135.60 -177.99 -9.98 160.49 1771.74
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 150.487495 722.264437 0.208 0.83497
## N.gravidanze 13.962670 4.257593 3.279 0.00105 **
## Gestazione 42.112604 3.865949 10.893 < 2e-16 ***
## Lunghezza -20.040737 3.155098 -6.352 2.52e-10 ***
## Cranio 10.605029 0.417940 25.375 < 2e-16 ***
## Tipo.partoNat 26.917987 11.866355 2.268 0.02339 *
## Ospedaleosp2 -9.622633 13.194789 -0.729 0.46590
## Ospedaleosp3 29.357873 13.245527 2.216 0.02675 *
## SessoM 69.520169 11.004595 6.317 3.14e-10 ***
## I(Lunghezza^2) 0.031482 0.003259 9.661 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 269 on 2490 degrees of freedom
## Multiple R-squared: 0.7385, Adjusted R-squared: 0.7376
## F-statistic: 781.3 on 9 and 2490 DF, p-value: < 2.2e-16Questo risultato è ancora più interessante, poiché l’effetto quadratico della lunghezza dei neonati è estremamente significativo, e l’effetto principale rimane significativo. La significatività delle altre variabili, invece, rimane la stessa rispetto al modello di partenza, mentre l’\(R^2\) aggiustato è leggermente maggiore, con una variabilità spiegata del 73.6%.
Eseguiamo un test ANOVA sui modelli 2 e 6.
anova(mod6, mod2)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso + I(Lunghezza^2)
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2490 180146731
## 2 2491 186899996 -1 -6753265 93.344 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Con un p-value molto vicino allo zero, si può confermare che l’aumento della variabilità spiegata nel modello 6 rispetto al modello 2 è estremamente significativa. Pertanto, ora è diventato il modello 6 quello da preferire.
Completiamo l’esplorazione degli effetti non lineari provando ad aggiungere l’effetto quadratico della variabile Cranio al modello 2.
mod7 <- update(mod2, ~ . + I(Cranio^2))
summary(mod7)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + I(Cranio^2), data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1104.10 -180.25 -12.96 160.32 2603.25
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.22382 1149.32164 0.018 0.98596
## N.gravidanze 12.64963 4.30373 2.939 0.00332 **
## Gestazione 38.47448 3.92165 9.811 < 2e-16 ***
## Lunghezza 10.54676 0.30111 35.026 < 2e-16 ***
## Cranio -31.59620 7.15289 -4.417 1.04e-05 ***
## Tipo.partoNat 27.50533 12.00451 2.291 0.02203 *
## Ospedaleosp2 -9.55407 13.34856 -0.716 0.47422
## Ospedaleosp3 30.20426 13.40078 2.254 0.02429 *
## SessoM 72.62010 11.13073 6.524 8.25e-11 ***
## I(Cranio^2) 0.06226 0.01056 5.894 4.27e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 272.1 on 2490 degrees of freedom
## Multiple R-squared: 0.7324, Adjusted R-squared: 0.7315
## F-statistic: 757.3 on 9 and 2490 DF, p-value: < 2.2e-16In questo caso, l’\(R^2\) aggiustato è di 0.7315, quindi si rileva anche qui un aumento della variabilità spiegata rispetto al modello di partenza. Tuttavia, dato che il modello 6 ha un \(R^2\) aggiustato di 0.7376, a parità di numero di parametri, rimane da preferire.
Proviamo ora, ad aggiungere l’effetto quadratico delle settimane di gestazione al modello 6.
mod8 <- update(mod6, ~ . + I(Gestazione^2))
summary(mod8)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + I(Lunghezza^2) + I(Gestazione^2),
## data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1156.03 -178.88 -8.63 155.46 1375.80
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.508e+03 9.037e+02 -2.775 0.005560 **
## N.gravidanze 1.433e+01 4.239e+00 3.381 0.000733 ***
## Gestazione 3.455e+02 6.259e+01 5.520 3.74e-08 ***
## Lunghezza -3.228e+01 4.028e+00 -8.016 1.67e-15 ***
## Cranio 1.039e+01 4.183e-01 24.844 < 2e-16 ***
## Tipo.partoNat 2.720e+01 1.181e+01 2.303 0.021386 *
## Ospedaleosp2 -1.117e+01 1.314e+01 -0.850 0.395188
## Ospedaleosp3 3.057e+01 1.319e+01 2.318 0.020511 *
## SessoM 7.217e+01 1.097e+01 6.579 5.74e-11 ***
## I(Lunghezza^2) 4.390e-02 4.130e-03 10.628 < 2e-16 ***
## I(Gestazione^2) -3.997e+00 8.230e-01 -4.856 1.27e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.8 on 2489 degrees of freedom
## Multiple R-squared: 0.741, Adjusted R-squared: 0.7399
## F-statistic: 711.9 on 10 and 2489 DF, p-value: < 2.2e-16Notiamo che, a differenza del modello 5, dove l’effetto principale di Gestazione non era significativo, in questo modello esso è estremamente significativo, così come l’effetto quadratico. Inoltre, l’\(R^2\) aggiustato ora è di 0.7399, quindi la variabilità spiegata è aumentata.
Eseguiamo il test ANOVA per il confronto tra le varianze.
anova(mod8, mod6)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso + I(Lunghezza^2) + I(Gestazione^2)
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Ospedale + Sesso + I(Lunghezza^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2489 178455910
## 2 2490 180146731 -1 -1690821 23.583 1.271e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Il p-value è molto vicino allo zero, quindi l’aumento di variabilità spiegata è significativo. Sembrerebbe, quindi, che il modello 8 sia da preferire al modello 6.
Dal momento che, precedentemente, abbiamo rilevato che le misure antropometriche sono significativamente diverse tra i due sessi, proviamo ad aggiungere al modello 8 delle interazioni tra queste variabili. In particolare, poiché la significatività era risultata maggiore nel caso della variabile Lunghezza, aggiungiamo al modello l’interazione tra questa variabile e la variabile Sesso.
mod9 <- update(mod8, ~ . + Lunghezza * Sesso)
summary(mod9)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + I(Lunghezza^2) + I(Gestazione^2) +
## Lunghezza:Sesso, data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1154.73 -178.67 -9.19 155.99 1377.87
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.493e+03 9.060e+02 -2.751 0.00598 **
## N.gravidanze 1.431e+01 4.241e+00 3.375 0.00075 ***
## Gestazione 3.460e+02 6.263e+01 5.524 3.66e-08 ***
## Lunghezza -3.244e+01 4.079e+00 -7.953 2.73e-15 ***
## Cranio 1.040e+01 4.187e-01 24.830 < 2e-16 ***
## Tipo.partoNat 2.714e+01 1.182e+01 2.297 0.02171 *
## Ospedaleosp2 -1.122e+01 1.314e+01 -0.854 0.39328
## Ospedaleosp3 3.053e+01 1.319e+01 2.314 0.02074 *
## SessoM 1.236e+02 2.138e+02 0.578 0.56310
## I(Lunghezza^2) 4.410e-02 4.218e-03 10.456 < 2e-16 ***
## I(Gestazione^2) -4.003e+00 8.236e-01 -4.860 1.24e-06 ***
## Lunghezza:SessoM -1.040e-01 4.313e-01 -0.241 0.80954
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.8 on 2488 degrees of freedom
## Multiple R-squared: 0.741, Adjusted R-squared: 0.7398
## F-statistic: 647 on 11 and 2488 DF, p-value: < 2.2e-16Poiché l’effetto di interazione non è significativo e, per giunta, la variabile Sesso non è più significativa, si può concludere che questa variabile non vada inserita nel modello.
Dal momento che uno dei focus di questo studio è l’impatto del fumo materno sul peso del neonato, proviamo ad inserire nuovamente la relativa variabile nel modello, per verificare se in questo nuovo modello è significativa.
mod10 <- update(mod8, ~ . + Fumatrici)
summary(mod10)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + I(Lunghezza^2) + I(Gestazione^2) +
## Fumatrici, data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1156.59 -178.76 -8.74 155.74 1374.88
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.527e+03 9.040e+02 -2.795 0.005224 **
## N.gravidanze 1.453e+01 4.245e+00 3.422 0.000631 ***
## Gestazione 3.459e+02 6.259e+01 5.527 3.60e-08 ***
## Lunghezza -3.223e+01 4.028e+00 -8.002 1.85e-15 ***
## Cranio 1.039e+01 4.184e-01 24.831 < 2e-16 ***
## Tipo.partoNat 2.741e+01 1.182e+01 2.320 0.020433 *
## Ospedaleosp2 -1.123e+01 1.314e+01 -0.854 0.393019
## Ospedaleosp3 3.035e+01 1.319e+01 2.301 0.021467 *
## SessoM 7.233e+01 1.097e+01 6.593 5.23e-11 ***
## I(Lunghezza^2) 4.384e-02 4.131e-03 10.611 < 2e-16 ***
## I(Gestazione^2) -4.000e+00 8.231e-01 -4.860 1.25e-06 ***
## Fumatrici1 -2.428e+01 2.693e+01 -0.902 0.367275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.8 on 2488 degrees of freedom
## Multiple R-squared: 0.741, Adjusted R-squared: 0.7399
## F-statistic: 647.2 on 11 and 2488 DF, p-value: < 2.2e-16La variabile relativa al fumo materno rimane non significativa. Tuttavia, si potrebbe ipotizzare un’interazione tra il fumo materno e le settimane di gestazione, dovuta al fatto che il fumo materno potrebbe portare a nascite premature.
mod11 <- update(mod8, ~ . + Gestazione * Fumatrici)
summary(mod11)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso + I(Lunghezza^2) + I(Gestazione^2) +
## Fumatrici + Gestazione:Fumatrici, data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1155.84 -179.19 -9.77 157.51 1380.25
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.465e+03 9.061e+02 -2.720 0.00658 **
## N.gravidanze 1.457e+01 4.245e+00 3.432 0.00061 ***
## Gestazione 3.401e+02 6.286e+01 5.410 6.90e-08 ***
## Lunghezza -3.208e+01 4.031e+00 -7.959 2.61e-15 ***
## Cranio 1.039e+01 4.184e-01 24.831 < 2e-16 ***
## Tipo.partoNat 2.757e+01 1.182e+01 2.333 0.01974 *
## Ospedaleosp2 -1.091e+01 1.314e+01 -0.830 0.40675
## Ospedaleosp3 3.062e+01 1.319e+01 2.321 0.02038 *
## SessoM 7.277e+01 1.098e+01 6.628 4.15e-11 ***
## I(Lunghezza^2) 4.368e-02 4.134e-03 10.565 < 2e-16 ***
## I(Gestazione^2) -3.917e+00 8.273e-01 -4.734 2.32e-06 ***
## Fumatrici1 7.167e+02 7.430e+02 0.965 0.33484
## Gestazione:Fumatrici1 -1.887e+01 1.891e+01 -0.998 0.31841
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.8 on 2487 degrees of freedom
## Multiple R-squared: 0.7411, Adjusted R-squared: 0.7399
## F-statistic: 593.4 on 12 and 2487 DF, p-value: < 2.2e-16Da quest’ultimo modello si può constatare che né l’effetto principale del fumo materno, né il suo effetto di interazione con le settimane di gestazione è significativo. Pertanto, il fumo materno non risulta avere un impatto sul peso dei neonati.
Finora, come criterio di valutazione dei modelli, abbiamo utilizzato solamente il test ANOVA. Ora adotteremo anche il Criterio di Informazione di Akaike (AIC) e il Criterio di Informazione Bayesiano (BIC). Secondo entrambi i criteri, il modello migliore è quello che riporta il valore per il criterio di informazione più basso; tuttavia, l’AIC tende a preferire i modelli con più parametri, mentre il BIC li penalizza, preferendo modelli più semplici.
model_list <- list(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8, mod9, mod10, mod11)
model_names <- paste0("mod", 1:11)
model_comparison <- data.frame(Modello = model_names,
AIC = sapply(model_list, AIC),
BIC = sapply(model_list, BIC))
model_comparison %>%
mutate(AIC = cell_spec(round(AIC, 2),
bold = ifelse(AIC == min(AIC), TRUE, FALSE),
background = ifelse(AIC == min(AIC), "yellow", "transparent")),
BIC = cell_spec(round(BIC, 2),
bold = ifelse(BIC == min(BIC), TRUE, FALSE),
background = ifelse(BIC == min(BIC), "yellow", "transparent"))) %>%
kable(digits = 2,
caption = "Confronto AIC e BIC dei modelli",
col.names = c("Modello", "AIC", "BIC"),
align = "c",
escape = FALSE,
format = "html") %>%
kable_styling(full_width = FALSE) %>%
column_spec(1:3, extra_css = "padding-left: 30px; padding-right: 30px;")| Modello | AIC | BIC |
|---|---|---|
| mod1 | 35171.95 | 35241.84 |
| mod2 | 35169.79 | 35228.03 |
| mod3 | 35713.93 | 35766.34 |
| mod4 | 36135.6 | 36188.02 |
| mod5 | 35167.17 | 35231.23 |
| mod6 | 35079.78 | 35143.85 |
| mod7 | 35137.15 | 35201.21 |
| mod8 | 35058.21 | 35128.1 |
| mod9 | 35060.15 | 35135.86 |
| mod10 | 35059.39 | 35135.1 |
| mod11 | 35060.39 | 35141.93 |
Entrambi i criteri confermano che il modello migliore è il modello 8.
Ora che abbiamo scelto il modello, verifichiamo che non ci siano problemi di multicollinearità calcolando i Variance Inflation Factor (VIF).
vif_8_output <- vif(mod8)
vif_8_data <- as.data.frame(vif_8_output) %>%
tibble::rownames_to_column(var = "Variabile")
vif_8_data %>%
kable(caption = "Analisi della multicollinearità (VIF)",
digits = 3,
col.names = c("Variabile", colnames(vif_8_data)[-1])) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = F,
position = "center")| Variabile | GVIF | Df | GVIF^(1/(2*Df)) |
|---|---|---|---|
| N.gravidanze | 1.027 | 1 | 1.013 |
| Gestazione | 476.738 | 1 | 21.834 |
| Lunghezza | 391.643 | 1 | 19.790 |
| Cranio | 1.646 | 1 | 1.283 |
| Tipo.parto | 1.004 | 1 | 1.002 |
| Ospedale | 1.005 | 2 | 1.001 |
| Sesso | 1.049 | 1 | 1.024 |
| I(Lunghezza^2) | 373.131 | 1 | 19.317 |
| I(Gestazione^2) | 454.420 | 1 | 21.317 |
Questi valori mostrano un evidente problema di multicollinearità nelle variabili Lunghezza e Gestazione e nei loro effetti quadratici, poiché i VIF non dovrebbero superare un valore di 5. Questo risultato è comprensibile, poiché l’effetto quadratico di una variabile è naturalmente molto correlato al suo effetto principale.
Una possibile soluzione in questi casi è di rimuovere una delle variabili in questione, ma avevamo già rilevato in precedenza che la qualità del modello ne risentirebbe. Pertanto, procederemo con il metodo della “centratura” delle variabili, che prevede che alle variabili venga sottratta la loro media, prima di elevarle al quadrato. Definiamo, quindi, due nuove variabili per le settimane di gestazione e la lunghezza del neonato, alle quali verrà sottratta la rispettiva media.
Gestazione_c <- as.numeric(scale(Gestazione, center = TRUE, scale = FALSE))
Lunghezza_c <- as.numeric(scale(Lunghezza, center = TRUE, scale = FALSE))A questo punto, creiamo un nuovo modello di regressione utilizzando queste due nuove variabili.
mod13 <- lm(Peso ~ N.gravidanze + Gestazione_c + Lunghezza_c + Cranio +
Tipo.parto + Ospedale + Sesso + I(Lunghezza_c^2) + I(Gestazione_c^2),
data = newborn_data)
summary(mod13)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione_c + Lunghezza_c +
## Cranio + Tipo.parto + Ospedale + Sesso + I(Lunghezza_c^2) +
## I(Gestazione_c^2), data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1156.03 -178.88 -8.63 155.46 1375.80
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -341.98222 142.25804 -2.404 0.016291 *
## N.gravidanze 14.33278 4.23910 3.381 0.000733 ***
## Gestazione_c 33.87857 4.20550 8.056 1.21e-15 ***
## Lunghezza_c 11.14883 0.30508 36.543 < 2e-16 ***
## Cranio 10.39307 0.41834 24.844 < 2e-16 ***
## Tipo.partoNat 27.20023 11.81305 2.303 0.021386 *
## Ospedaleosp2 -11.17354 13.13924 -0.850 0.395188
## Ospedaleosp3 30.57491 13.18825 2.318 0.020511 *
## SessoM 72.16715 10.96858 6.579 5.74e-11 ***
## I(Lunghezza_c^2) 0.04390 0.00413 10.628 < 2e-16 ***
## I(Gestazione_c^2) -3.99680 0.82303 -4.856 1.27e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.8 on 2489 degrees of freedom
## Multiple R-squared: 0.741, Adjusted R-squared: 0.7399
## F-statistic: 711.9 on 10 and 2489 DF, p-value: < 2.2e-16Calcoliamo i VIF su questo nuovo modello per verificare la multicollinearità.
vif_13_output <- vif(mod13)
vif_13_data <- as.data.frame(vif_13_output) %>%
tibble::rownames_to_column(var = "Variabile")
vif_13_data %>%
kable(caption = "Analisi della multicollinearità (VIF)",
digits = 3,
col.names = c("Variabile", colnames(vif_13_data)[-1])) %>%
kable_styling(bootstrap_options = c("striped", "condensed"),
full_width = F,
position = "center")| Variabile | GVIF | Df | GVIF^(1/(2*Df)) |
|---|---|---|---|
| N.gravidanze | 1.027 | 1 | 1.013 |
| Gestazione_c | 2.153 | 1 | 1.467 |
| Lunghezza_c | 2.247 | 1 | 1.499 |
| Cranio | 1.646 | 1 | 1.283 |
| Tipo.parto | 1.004 | 1 | 1.002 |
| Ospedale | 1.005 | 2 | 1.001 |
| Sesso | 1.049 | 1 | 1.024 |
| I(Lunghezza_c^2) | 2.420 | 1 | 1.556 |
| I(Gestazione_c^2) | 2.952 | 1 | 1.718 |
Ora tutti i valori sono abbondantemente inferiori alla soglia di 5, quindi il problema della multicollinearità è stato risolto. Il modello 13, quindi, risulta essere il modello migliore.
mod_final <- mod13
summary(mod_final)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione_c + Lunghezza_c +
## Cranio + Tipo.parto + Ospedale + Sesso + I(Lunghezza_c^2) +
## I(Gestazione_c^2), data = newborn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1156.03 -178.88 -8.63 155.46 1375.80
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -341.98222 142.25804 -2.404 0.016291 *
## N.gravidanze 14.33278 4.23910 3.381 0.000733 ***
## Gestazione_c 33.87857 4.20550 8.056 1.21e-15 ***
## Lunghezza_c 11.14883 0.30508 36.543 < 2e-16 ***
## Cranio 10.39307 0.41834 24.844 < 2e-16 ***
## Tipo.partoNat 27.20023 11.81305 2.303 0.021386 *
## Ospedaleosp2 -11.17354 13.13924 -0.850 0.395188
## Ospedaleosp3 30.57491 13.18825 2.318 0.020511 *
## SessoM 72.16715 10.96858 6.579 5.74e-11 ***
## I(Lunghezza_c^2) 0.04390 0.00413 10.628 < 2e-16 ***
## I(Gestazione_c^2) -3.99680 0.82303 -4.856 1.27e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 267.8 on 2489 degrees of freedom
## Multiple R-squared: 0.741, Adjusted R-squared: 0.7399
## F-statistic: 711.9 on 10 and 2489 DF, p-value: < 2.2e-16Analizziamo i coefficienti del modello finale e i relativi significati:
- Per ogni gravidanza in più della madre, il neonato pesa mediamente 14.3 grammi in più, mantenendo costanti le altre variabili;
- Per ogni millimetro in più nel diametro del cranio, il neonato pesa circa 10.4 grammi in più;
- All’aumentare delle settimane di gestazione, aumenta anche il peso del neonato, ma l’incremento non è costante. Infatti, esso tende a rallentare verso la fine della gestazione;
- All’aumentare della lunghezza del neonato, non solo il peso aumenta, ma tende ad aumentare più rapidamente;
- I bambini nati con parto naturale pesano mediamente 27.2 grammi in più rispetto ai bambini nati con parto cesareo, sempre mantenendo costanti le altre variabili;
- I neonati tendono a pesare circa 72.2 grammi in più rispetto alle neonate;
- I bambini nati nell’ospedale 3 pesano mediamente 30.6 grammi in più rispetto ai bambini nati nell’ospedale 1. Non ci sono differenze significative, invece, tra i bambini nati nell’ospedale 2 e quelli nati nell’ospedale 1.
Analisi della Qualità del Modello
Ora che abbiamo scelto il modello di regressione finale, ne valutiamo la capacità predittiva. Come abbiamo già visto, il valore per l’\(R^2\) aggiustato è di 0.7399, indicando quindi che il modello spiega il 73.99% della variabilità. Un’altra importante metrica, tuttavia, è il Root Mean Squared Error (RMSE), che misura la differenza media tra i valori stimati di un modello e i valori reali. Calcoliamo quindi questa metrica estraendo i residui del modello, elevandoli al quadrato, trovandone la media e calcolando la radice quadrata di quest’ultima.
res <- residuals(mod_final)
rmse <- sqrt(mean(res^2))
rmse## [1] 267.1748Con un valore del RMSE di circa 267.2, possiamo affermare che il modello, in media, sbaglia le stime del peso dei neonati di circa 267 grammi. Considerato che, secondo la letteratura medica, il peso di un neonato può variare tra i 2500 e i 4500 grammi, il modello può essere definito sufficientemente accurato.
Un’altra importante fase nell’analisi della qualità del modello è la diagnostica sui residui. Un modello affidabile, infatti, deve anche rispettare tutte le assunzioni sui residui:
- Linearità;
- Normalità della distribuzione;
- Omoschedasticità (varianza costante);
- Indipendenza tra loro e dal regressore.
par(mfrow = c(2, 2))
plot(mod_final)
Analizziamo ciascun quadrante, tenendo però a mente che ogni assunzione dovrà essere poi verificata tramite i relativi test statistici:
- Nel grafico in alto a sinistra, i residui risultano sparsi casualmente attorno alla media di zero, senza quindi seguire un pattern. Ciò significherebbe che l’informazione viene filtrata bene dai regressori e non si riversa, quindi, sui residui. Sembrerebbe che si possa confermare, quindi, l’assunzione di linearità.
- Nel grafico in alto a destra vengono messi in relazione i residui con i quantili di una distribuzione normale. Se i punti giacciono sulla bisettrice del grafico, significa che i residui seguono una distribuzione normale. In questo caso, i residui risultano correttamente disposti sopra la retta, ad eccezione di alcuni punti alle due estremità del grafico. Pertanto, l’assunzione di normalità della distribuzione dei residui sembrerebbe grossomodo rispettata.
- Anche nel grafico in basso a sinistra non dovremmo visualizzare alcun pattern, ed effettivamente vediamo che i residui sono sparsi in una nuvola di punti più o meno orizzontale attorno a un valore di y. Ciò indicherebbe una varianza costante, quindi anche l’assunzione di omoschedasticità risulterebbe rispettata.
- Nel grafico in basso a destra si possono visualizzare i potenziali valori influenti, ovvero i valori di leva e gli outliers. Le linee grigie tratteggiate rappresentano la soglia di avvertimento (0.5) o di allarme (1) relative alla distanza di Cook; quando un residuo supera queste soglie, significa che la relativa osservazione è potenzialmente influente sulle stime di regressione. A prima vista, risultano esserci un paio di residui nella parte inferiore del grafico e uno nella parte superiore che possono rivelarsi valori influenti.
Rappresentiamo graficamente la funzione di densità dei residui del modello.
res_mod <- data.frame(index = 1:length(residuals(mod_final)),
res = residuals(mod_final))
ggplot(res_mod, aes(x = res))+
geom_density(fill = "steelblue",
alpha = 0.5)+
geom_vline(xintercept = mean(res_mod$res),
color = "red",
linetype = "dashed")+
labs(title = "Funzione di densità dei residui del modello",
x = "Residui",
y = "Densità")+
theme_minimal()
La distribuzione dei residui assomiglia abbastanza a una normale, anche se risulta leggermente leptocurtica. Per poterlo affermare, però, dobbiamo condurre un test statistico.
Eseguiamo, quindi, la diagnostica sui residui mediante i relativi test statistici. In particolare, utilizzeremo il test di Shapiro-Wilk per verificare l’ipotesi di normalità della distribuzione dei residui, il test Breusch-Pagan per verificare l’ipotesi di omoschedasticità, e il test Durbin-Watson per verificare l’ipotesi di incorrelazione.
shapiro.test(residuals(mod_final))##
## Shapiro-Wilk normality test
##
## data: residuals(mod_final)
## W = 0.98897, p-value = 5.697e-13Con un p-value molto vicino allo zero, si rifiuta l’ipotesi nulla, quindi la distribuzione dei residui non è completamente normale.
bptest(mod_final)##
## studentized Breusch-Pagan test
##
## data: mod_final
## BP = 92.92, df = 10, p-value = 1.409e-15Con un p-value molto vicino allo zero, si rifiuta l’ipotesi nulla di omoschedasticità. Pertanto, i residui non hanno varianza costante.
dwtest(mod_final)##
## Durbin-Watson test
##
## data: mod_final
## DW = 1.9475, p-value = 0.09469
## alternative hypothesis: true autocorrelation is greater than 0Con un p-value di 0.095, non si rifiuta l’ipotesi nulla di incorrelazione tra i residui.
Ora approfondiamo i valori di leva e gli outliers con i metodi numerici.
Cominciando dai valori di leva, identifichiamo il valore soglia e verifichiamo se esistono valori oltre quella soglia.
lev <- hatvalues(mod_final)
p <- sum(lev)
lev_threshold <- 2 * p/N
df_lev <- data.frame(index = 1:N, leverage = lev)
ggplot(df_lev, aes(x = index,
y = leverage))+
geom_point(shape = 21,
color = "black",
fill = "steelblue",
size = 3,
alpha = 0.6)+
geom_hline(yintercept = lev_threshold,
color = "red",
linetype = "solid",
size = 1)+
labs(title = "Valori di leva",
x = "Indice dell'osservazione",
y = "Leverage") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
lev[lev > lev_threshold]## 15 67 89 101 106 131
## 0.010330890 0.009606717 0.014458901 0.018023391 0.025711706 0.008811916
## 151 155 161 204 206 304
## 0.012757398 0.013568105 0.024556221 0.015664240 0.014478432 0.009616180
## 305 310 312 378 442 445
## 0.009033682 0.093121884 0.039393547 0.059424850 0.009270375 0.009773936
## 471 492 516 582 587 592
## 0.010949196 0.021370846 0.014272708 0.013189868 0.021756826 0.013761086
## 666 684 702 729 748 750
## 0.009695816 0.010907700 0.008926462 0.010664364 0.015688230 0.013767392
## 757 765 805 928 947 956
## 0.009220407 0.008892924 0.035944420 0.115980534 0.013682184 0.009955285
## 985 1014 1067 1091 1130 1188
## 0.010050958 0.028679206 0.016714419 0.010616375 0.037400760 0.011190064
## 1219 1248 1273 1283 1311 1321
## 0.032193759 0.026871368 0.013727271 0.009297824 0.010625464 0.010416090
## 1323 1357 1385 1400 1411 1428
## 0.009587102 0.013988124 0.041221519 0.009193286 0.010045943 0.030267964
## 1429 1450 1505 1513 1551 1553
## 0.045359247 0.016132525 0.015106865 0.010311426 0.248708686 0.009802599
## 1610 1619 1628 1686 1701 1718
## 0.014402825 0.066122040 0.009730338 0.020799481 0.019515686 0.009755128
## 1727 1780 1781 1809 1920 1977
## 0.014653543 0.130521698 0.018368233 0.012300044 0.009041713 0.011061747
## 2040 2086 2087 2089 2114 2115
## 0.020991564 0.014491131 0.009347549 0.010901084 0.039938716 0.023367207
## 2120 2140 2148 2149 2175 2200
## 0.070296033 0.008948819 0.009981219 0.023184310 0.063298407 0.025701688
## 2216 2221 2257 2307 2317 2359
## 0.019559178 0.022777095 0.011792278 0.026168151 0.008994372 0.011573014
## 2391 2408 2422 2437 2452 2458
## 0.012955038 0.014855006 0.022339588 0.093905591 0.095867442 0.019044760
## 2471 2478
## 0.023079602 0.012558689table(lev > lev_threshold)##
## FALSE TRUE
## 2402 98I valori di leva risultano essere 98 su 2500. In particolare, l’osservazione n. 1551, evidente anche nel grafico, presenta un valore leva di circa 0.25.
Ora indaghiamo i valori outliers, impostando le soglie convenzionali di -2 e 2.
res_student <- rstudent(mod_final)
df_res <- data.frame(index = 1:length(res_student),
residuals = res_student)
ggplot(df_res, aes(x = index,
y = residuals))+
geom_point(shape = 21,
color = "black",
fill = "steelblue",
size = 3,
alpha = 0.6)+
geom_hline(yintercept = c(-2, 2),
color = "red",
linetype = "solid",
size = 1)+
geom_hline(yintercept = 0,
color = "black",
size = 0.5) +
labs(title = "Valori outlier",
x = "Indice dell'osservazione",
y = "Residui studentizzati")+
theme_minimal()
Dal grafico si evince la presenza di molti outliers. Vediamo quanti sono significativi.
outlier_results <- outlierTest(mod_final)
df_outliers <- data.frame("Osservazione" = as.numeric(names(outlier_results$rstudent)),
"Residuo studentizzato" = outlier_results$rstudent,
"p-value non aggiornato" = outlier_results$p,
"p-value di Bonferroni" = outlier_results$bonf) %>%
arrange(Osservazione)
df_outliers %>%
kable(caption = "Outliers significativi",
digits = 10,
align = "c",
col.names = c("N° osservazione", "Residuo studentizzato", "p-value", "p-value Bonferroni")) %>%
kable_styling(bootstrap_options = c("striped", "condensed"),
full_width = F,
position = "center") %>%
column_spec(1, bold = TRUE, border_right = TRUE)| N° osservazione | Residuo studentizzato | p-value | p-value Bonferroni | |
|---|---|---|---|---|
| 155 | 155 | 4.486976 | 7.55260e-06 | 0.0188815113 |
| 1306 | 1306 | 5.015951 | 5.64900e-07 | 0.0014123217 |
| 1399 | 1399 | -4.342232 | 1.46688e-05 | 0.0366720707 |
| 1551 | 1551 | 5.968964 | 2.70000e-09 | 0.0000068208 |
| 1694 | 1694 | 4.304530 | 1.73816e-05 | 0.0434540736 |
Da questo test si evince che le osservazioni effettivamente confermate come outliers sono 5, poiché hanno un p-value di Bonferroni molto vicino allo zero.
Creiamo ora un grafico che consideri sia i valori di leva che gli outliers, utilizzando la distanza di Cook.
cook_values <- cooks.distance(mod_final)
df_cook <- data.frame(obs = 1:length(cook_values),
cook_distance = cook_values)
cook_threshold <- 4/nrow(model.frame(mod_final))
ggplot(df_cook, aes(x = obs,
y = cook_distance))+
geom_segment(aes(xend = obs, yend = 0),
color = "black",
size = 0.5)+
geom_point(aes(color = cook_distance > cook_threshold), size = 2)+
geom_hline(yintercept = cook_threshold,
linetype = "solid",
color = "red")+
scale_color_manual(values = c("black", "red"))+
labs(title = "Distanza di Cook",
x = "Indice dell'osservazione",
y = "Distanza di Cook")+
theme_minimal()+
theme(legend.position = "none")
Notiamo come l’osservazione n. 1551, ovvero quella con il p-value Bonferroni più vicino allo zero e quindi significativamente outlier, arrivi ben oltre la soglia di 4/N = 0.0016. Nel caso di questa osservazione, infatti, la distanza di Cook arriva oltre 1. Dal grafico si rilevano altri valori oltre la soglia, ma in maniera nettamente più contenuta.
Vediamo i valori di questa osservazione.
newborn_data[1551, ] %>%
t() %>%
as.data.frame() %>%
kable(col.names = "Valore", align = "c") %>%
kable_styling(
full_width = FALSE,
bootstrap_options = "condensed",
position = "center") %>%
column_spec(1, bold = TRUE, background = "#f8f9fa")| Valore | |
|---|---|
| Anni.madre | 35 |
| N.gravidanze | 1 |
| Fumatrici | 0 |
| Gestazione | 38 |
| Peso | 4370 |
| Lunghezza | 315 |
| Cranio | 374 |
| Tipo.parto | Nat |
| Ospedale | osp3 |
| Sesso | F |
Questi valori confermano la rarità dell’osservazione, in quanto si tratta di una bambina nata a 38 settimane di gestazione con una lunghezza di 31.5 centimetri, quindi ben al di sotto della media, con un peso di 4.37 kg.
3. Previsioni e Risultati
Per mezzo del modello di regressione scelto, proviamo ad eseguire alcune previsioni.
Per esempio, stimiamo il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana. Poiché dovremo assegnare un valore a ciascuna variabile, alle variabili che vanno mantenute costanti assegneremo il valore medio. Pertanto, come prima cosa dobbiamo calcolare le medie delle variabili. Per le variabili dummy, invece, useremo nel modello la modalità baseline.
mean_grav <- round(mean(N.gravidanze), digits = 0)
mean_gest <- mean(Gestazione)
mean_lung <- mean(Lunghezza)
mean_cran <- mean(Cranio)Ora impostiamo la previsione, ricordando però che, nel caso delle variabili Gestazione e Lunghezza, stiamo usando la loro versione “centrata”, e che quindi vanno considerate come distanze dalla media.
baby_1 <- data.frame(Sesso = "F",
N.gravidanze = 2,
Gestazione_c = 39 - mean_gest,
Lunghezza_c = 0,
Cranio = mean_cran,
Tipo.parto = "Nat",
Ospedale = "osp1")
predict(mod_final, newdata = baby_1)## 1
## 3248.493Basandoci su questo modello, il peso di una neonata considerando una madre alla terza gravidanza (che quindi ha già avuto due gravidanze) che partorirà alla 39esima settimana è di circa 3248 grammi.
Ora proviamo a stimare il peso di un neonato lungo 52 centimetri (520 millimetri) e nato con parto cesareo alla 42esima settimana.
baby_2 <- data.frame(Sesso = "M",
N.gravidanze = mean_grav,
Gestazione_c = 42 - mean_gest,
Lunghezza_c = 520 - mean_lung,
Cranio = mean_cran,
Tipo.parto = "Ces",
Ospedale = "osp1")
predict(mod_final, newdata = baby_2)## 1
## 3654.594In base alla previsione di questo modello, un neonato lungo 52 centimetri (520 millimetri) e nato con parto cesareo alla 42esima settimana peserà circa 3655 grammi.
Infine, stimiamo il peso di una bambina nata con parto naturale alla 34esima settimana di gestazione nell’ospedale 3, e con un diametro craniale di 320 millimetri.
baby_3 <- data.frame(Sesso = "F",
N.gravidanze = mean_grav,
Gestazione_c = 34 - mean_gest,
Lunghezza_c = 0,
Cranio = 320,
Tipo.parto = "Nat",
Ospedale = "osp3")
predict(mod_final, newdata = baby_3)## 1
## 2788.041In base alla previsione del modello, la neonata con le caratteristiche sopra riportate peserà circa 2788 grammi.
4. Visualizzazioni
Ora, per mezzo di rappresentazioni grafiche, visualizzeremo le relazioni principali tra le variabili maggiormente significative del modello, in particolare la durata della gestazione, le misure antropometriche e il sesso del neonato, e il numero di gravidanze della madre.
Innanzitutto, rappresentiamo l’impatto delle settimane di gestazione sul peso del neonato, distinguendo anche per il sesso del neonato. Dato che nel modello è stata inclusa anche la variabile \(Gestazione^2\), il grafico mostrerà l’effetto quadratico della variabile esplicativa sulla variabile risposta.
ggplot(data = newborn_data)+
geom_point(aes(x = Gestazione_c,
y = Peso,
col = Sesso),
position = "jitter",
alpha = 0.5)+
geom_smooth(aes(x = Gestazione_c,
y = Peso,
col = Sesso),
formula = y ~ x + I(x^2),
method = "lm",
se = F)+
labs(title = "Relazione tra Gestazione e Peso per Sesso",
x = "Settimane di gestazione (variabile centrata)",
y = "Peso (g)",
color = "Sesso")+
theme_minimal()
Dato che la variabile Gestazione è stata centrata, notiamo che la nuvola di punti si concentra attorno al valore 0 dell’asse x, ovvero attorno al valore medio.
In generale, comunque, il grafico mostra chiaramente la relazione non lineare tra Gestazione e Peso. In particolare, nonostante il peso del neonato aumenti all’avanzare delle settimane di gestazione, questo incremento tende a rallentare verso il termine della gravidanza. Si può notare, inoltre, che la curva condizionata per sesso maschile è più inclinata rispetto alla curva condizionata per sesso femminile, significando quindi che la variazione di peso di un neonato nel corso delle settimane di gestazione è più marcata rispetto a quella di una neonata.
Rappresentiamo ora la relazione tra la lunghezza e il peso del neonato. Anche in questo caso, utilizzeremo la variabile di Lunghezza centrata e includeremo l’effetto quadratico presente nel modello.
ggplot(data = newborn_data)+
geom_point(aes(x = Lunghezza_c,
y = Peso,
col = Sesso),
position = "jitter",
alpha = 0.5)+
geom_smooth(aes(x = Lunghezza_c,
y = Peso,
col = Sesso),
formula = y ~ x + I(x^2),
method = "lm",
se = F)+
labs(title = "Relazione tra Lunghezza e Peso per Sesso",
x = "Lunghezza in mm (variabile centrata)",
y = "Peso (g)",
color = "Sesso")+
theme_minimal()
Da questo grafico risulta ancora più evidente il significato del coefficiente positivo per l’effetto quadratico della lunghezza sul peso. Infatti, all’aumentare della lunghezza del neonato, non solo il peso aumenta, ma tende ad aumentare più rapidamente. Questo effetto è maggiormente accentuato nelle neonate, mentre per i neonati la linea risulta grossomodo una retta.
Ora creiamo un grafico che mostri la relazione tra il diametro craniale e il peso del neonato, distinguendo anche in questo caso per sesso.
ggplot(data = newborn_data)+
geom_point(aes(x = Cranio,
y = Peso,
col = Sesso),
position = "jitter",
alpha = 0.5)+
geom_smooth(aes(x = Cranio,
y = Peso,
col = Sesso),
formula = y ~ x,
method = "lm",
se = F)+
labs(title = "Relazione tra diametro craniale e peso per sesso",
x = "Diametro craniale (mm)",
y = "Peso (g)",
color = "Sesso")+
theme_minimal()
Da questo grafico emerge la relazione lineare positiva tra il diametro craniale e il peso del neonato, con la retta condizionata per sesso maschile che rimane sempre al di sopra della retta condizionata per sesso femminile. Un aspetto interessante da evidenziare è che il punto di inizio della linea per i neonati è più a destra rispetto all’inizio del grafico. Il motivo è che le osservazioni più basse, in questo caso, sono tutte su neonate.
Riproponiamo ora le stesse relazioni tra la variabile risposta e le variabili esplicative appena viste, ma condizionate per il tipo di parto.
ggplot(data = newborn_data)+
geom_point(aes(x = Gestazione_c,
y = Peso,
col = Tipo.parto),
position = "jitter",
alpha = 0.5)+
geom_smooth(aes(x = Gestazione_c,
y = Peso,
col = Tipo.parto),
formula = y ~ x + I(x^2),
method = "lm",
se = F)+
labs(title = "Relazione tra Gestazione e Peso per Tipo di parto",
x = "Settimane di gestazione (variabile centrata)",
y = "Peso (g)",
color = "Tipo di parto")+
scale_color_brewer(palette = "Dark2")+
theme_minimal()
Dal grafico non risultano differenze evidenti tra la curva condizionata per parto naturale e la curva condizionata per parto cesareo, indicando quindi un’influenza marginale del tipo di parto sul peso del neonato. Verso gli estremi del grafico, sia inferiore che superiore, si nota una separazione leggermente maggiore tra le linee, con la curva condizionata per parto naturale che si posiziona più in alto rispetto alla curva condizionata per parto cesareo. Ciò significa, quindi, che nei casi di parto molto prematuro oppure oltre il termine della gestazione, i bambini nati con parto naturale tendono a pesare lievemente di più rispetto a quelli nati con parto cesareo.
ggplot(data = newborn_data)+
geom_point(aes(x = Lunghezza_c,
y = Peso,
col = Tipo.parto),
position = "jitter",
alpha = 0.5)+
geom_smooth(aes(x = Lunghezza_c,
y = Peso,
col = Tipo.parto),
formula = y ~ x + I(x^2),
method = "lm",
se = F)+
labs(title = "Relazione tra Lunghezza e Peso per Tipo di parto",
x = "Lunghezza in mm (variabile centrata)",
y = "Peso (g)",
color = "Tipo di parto")+
scale_color_brewer(palette = "Dark2")+
theme_minimal()
Anche considerando l’impatto della lunghezza sul peso del neonato, non risultano differenze sostanziali tra parto naturale e parto cesareo. Si può notare, tuttavia, una separazione maggiore tra le due curve nella parte sinistra del grafico, indicando una maggiore influenza del tipo di parto sul peso del neonato se la lunghezza è minore (solitamente, ciò si verifica in caso di parto prematuro).
ggplot(data = newborn_data)+
geom_point(aes(x = Cranio,
y = Peso,
col = Tipo.parto),
position = "jitter",
alpha = 0.5)+
geom_smooth(aes(x = Cranio,
y = Peso,
col = Tipo.parto),
formula = y ~ x,
method = "lm",
se = F)+
labs(title = "Relazione tra diametro craniale e peso per tipo di parto",
x = "Diametro craniale (mm)",
y = "Peso (g)",
color = "Tipo di parto")+
scale_color_brewer(palette = "Dark2")+
theme_minimal()
La relazione che emerge da questo grafico è molto interessante. Infatti, sebbene non risultino differenze marcate tra i due tipi di parto, in presenza di un diametro craniale minore, un bambino nato con parto cesareo risulta pesare mediamente di più, mentre vale l’opposto in presenza di un diametro craniale maggiore. In corrispondenza della nuvola di punti rappresentante la maggior parte delle osservazioni, invece, le due linee si incontrano, indicando quindi l’assenza di un impatto rilevante.
ggplot(newborn_data, aes(x = factor(N.gravidanze),
y = Peso))+
geom_boxplot(outlier.shape = NA,
alpha = 0.7,
fill = "steelblue")+
geom_point(position = "jitter",
alpha = 0.1) +
labs(title = "Relazione tra peso del neonato e numero di gravidanze della madre",
x = "Numero di gravidanze",
y = "Peso (g)")+
theme_minimal()
Come si evince dal grafico, la quasi totalità delle osservazioni si concentra nella parte sinistra della distribuzione. La maggior parte delle madri osservate in questo dataset è alla prima gravidanza, e il numero di osservazioni cala drasticamente dalla quarta gravidanza in avanti. Soffermandosi sulla parte più densa del grafico, si può notare un leggero incremento del peso del neonato all’aumentare del numero di gravidanze della madre. Questa tendenza è visibile tra i valori 0 e 4 dell’asse x.
5. Conclusioni
L’analisi condotta ha portato a dei risvolti molto interessanti. Innanzitutto, l’età della madre e il fumo materno non sono risultati impattanti sul peso del neonato, contrariamente alle aspettative.
Le misure antropometriche, invece, si sono rivelate dei parametri chiave per la stima. In combinazione con il sesso del neonato, infatti, la lunghezza e il diametro del cranio, i quali possono essere misurati con le ecografie, possono fornire una previsione sufficientemente accurata di quanto pesa il neonato, intercettando quindi per tempo eventuali situazioni rischiose per la salute e necessità di cure o attenzioni specifiche dopo la nascita.
La durata della gestazione si conferma una variabile di forte impatto sul peso dei neonati. L’aspetto più interessante e utile, in particolare, è il suo effetto quadratico, ovvero il rallentamento della crescita durante le ultime settimane di gravidanza. Questo per esempio può significare che, nel caso di gravidanze di durata oltre il termine, attendere che la madre entri in travaglio senza induzione apporta un contributo solamente marginale nell’ulteriore aumento di peso del feto.
Un’altra variabile risultata significativa è il numero di gravidanze precedenti della madre. Poiché la relazione con il peso del neonato è positiva, il rischio di un neonato con un peso troppo basso ricade maggiormente su donne alla prima gravidanza. Ciò può essere un ulteriore aspetto da considerare nella valutazione di eventuali cure pediatriche neonatali necessarie.
Infine, il tipo di parto è un ulteriore parametro che incide sul peso del neonato; in particolare, i bambini nati con parto naturale pesano mediamente di più rispetto ai bambini nati con parto cesareo. Qualora la salute della madre non sia a rischio, quindi, si potrebbe valutare di optare per il parto naturale. Tuttavia, non sono state rilevate differenze nell’incidenza di parti cesarei tra i tre ospedali coinvolti, quindi, non risultano essere applicate pratiche e protocolli diversi.