Modello statistico per la previsione del peso neonatale
1 - Raccolta dei dati e struttura del dataset
Il dataset è composto da 2500 osservazioni di 10 variabili. Osservando i dati mi accorgo che due osservazioni presentano per la variabile “Anni.madre” i valori 0 e 1. Le considero come errori nella raccolta dati e sostituisco i valori 0 e 1 con la media dell’età della madri (28 anni). Inoltre, siccome è codificata in numeri, converto la variabile “Fumatrici” in una variabile qualitativa. Facciamo poi una prima analisi sulle variabili del dataset.
dati <- read.csv("neonati.csv",
stringsAsFactors = T)
dati$Anni.madre <- ifelse(dati$Anni.madre %in% c(0, 1), 28, dati$Anni.madre)
dati$Fumatrici <- factor(dati$Fumatrici, levels = c(0, 1), labels = c("Non fumatrice", "Fumatrice"))
summary(dati)## Anni.madre N.gravidanze Fumatrici Gestazione
## Min. :13.00 Min. : 0.0000 Non fumatrice:2396 Min. :25.00
## 1st Qu.:25.00 1st Qu.: 0.0000 Fumatrice : 104 1st Qu.:38.00
## Median :28.00 Median : 1.0000 Median :39.00
## Mean :28.19 Mean : 0.9812 Mean :38.98
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:40.00
## Max. :46.00 Max. :12.0000 Max. :43.00
## Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
## Min. : 830 Min. :310.0 Min. :235 Ces: 728 osp1:816 F:1256
## 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330 Nat:1772 osp2:849 M:1244
## Median :3300 Median :500.0 Median :340 osp3:835
## Mean :3284 Mean :494.7 Mean :340
## 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
## Max. :4930 Max. :565.0 Max. :390str(dati)## 'data.frame': 2500 obs. of 10 variables:
## $ Anni.madre : num 26 21 34 28 20 32 26 25 22 23 ...
## $ N.gravidanze: int 0 2 3 1 0 0 1 0 1 0 ...
## $ Fumatrici : Factor w/ 2 levels "Non fumatrice",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Gestazione : int 42 39 38 41 38 40 39 40 40 41 ...
## $ Peso : int 3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
## $ Lunghezza : int 490 490 500 515 480 495 480 510 500 510 ...
## $ Cranio : int 325 345 375 365 335 340 345 349 335 362 ...
## $ Tipo.parto : Factor w/ 2 levels "Ces","Nat": 2 2 2 2 2 2 2 2 1 1 ...
## $ Ospedale : Factor w/ 3 levels "osp1","osp2",..: 3 1 2 2 3 2 3 1 2 2 ...
## $ Sesso : Factor w/ 2 levels "F","M": 2 1 2 2 1 1 1 2 1 1 ...Anni.madre: misura dell’età in anni. Quantitativa continua anche se i valori sono stati registrati in anni interi.
N.gravidanze: le gravidanze che ha avuto in passato la madre. Quantitativa discreta.
Fumatrici: indicatore binario. Qualitativa nominale.
Gestazione: Numero di settimane di gestazione. Quantitativa continua anche se i valori sono stati registrati in settimane intere.
Peso: peso alla nascita in grammi. Quantitativa continua.
Lunghezza: lunghezza del neonato alla nascita. Quantitativa continua.
Cranio: diametro del cranio del neonato. Quantitativa continua.
Tipo.parto: indicatore binario (naturale o cesareo). Qualitativa nominale.
Ospedale: ospedale 1,2 e 3. Qualitativa nominale.
Sesso: indicatore binario. Qualitativa nominale.
2 - Analisi e modellizzazione
Analisi preliminare
Per le variabili quantitative (Anni.madre, N.gravidanze, Gestazione, Peso, Lunghezza, Cranio) calcoliamo gli indici di posizione, variabilità e forma, e creiamo boxplot e grafici di densità. Per le variabili qualitative (Sesso, Ospedale; Fumatrici; Tipo.parto), invece, calcoliamo la distribuzione di frequenza.
library(moments)
library(ggplot2)
attach(dati)
funz_ind <- function(x) {
result <- list(
mean = mean(x),
min = min(x),
max = max(x),
median = median(x),
Q1 = quantile(x)[2],
Q3 = quantile(x)[4],
sd = sd(x),
cv = (sd(x)/mean(x))*100,
variance = var(x),
skewness = skewness(x),
kurtosi= kurtosis(x)-3
)
}
ind_anni_madre <- funz_ind(Anni.madre)
ind_gravidanze <- funz_ind(N.gravidanze)
ind_gestazione <- funz_ind(Gestazione)
ind_peso <- funz_ind(Peso)
ind_lunghezza <- funz_ind(Lunghezza)
ind_cranio <- funz_ind(Cranio)
ind_var_quant <- data.frame(
anni_madre = unlist(ind_anni_madre),
n_gravidanze = unlist(ind_gravidanze),
gestazione = unlist(ind_gestazione),
peso = unlist(ind_peso),
lunghezza = unlist(ind_lunghezza),
cranio = unlist(ind_cranio)
)
ind_var_quant## anni_madre n_gravidanze gestazione peso lunghezza
## mean 28.1860000 0.981200 38.980400 3.284081e+03 494.692000
## min 13.0000000 0.000000 25.000000 8.300000e+02 310.000000
## max 46.0000000 12.000000 43.000000 4.930000e+03 565.000000
## median 28.0000000 1.000000 39.000000 3.300000e+03 500.000000
## Q1.25% 25.0000000 0.000000 38.000000 2.990000e+03 480.000000
## Q3.75% 32.0000000 1.000000 40.000000 3.620000e+03 510.000000
## sd 5.2151206 1.280587 1.868639 5.250387e+02 26.318644
## cv 18.5025212 130.512310 4.793792 1.598739e+01 5.320208
## variance 27.1974830 1.639903 3.491813 2.756657e+05 692.671004
## skewness 0.1512083 2.514254 -2.065313 -6.470308e-01 -1.514699
## kurtosi -0.1032773 10.989406 8.258150 2.031532e+00 6.487174
## cranio
## mean 340.0292000
## min 235.0000000
## max 390.0000000
## median 340.0000000
## Q1.25% 330.0000000
## Q3.75% 350.0000000
## sd 16.4253299
## cv 4.8305645
## variance 269.7914639
## skewness -0.7850527
## kurtosi 2.9462063boxplot(Anni.madre,
main = "Boxplot dell'età della madre",
ylab = "Età (anni)",
col = "lightblue",
border = "darkblue",
outcol = "red",
outpch = 19)
plot(density(Anni.madre),
main = "Densità dell'età della madre",
xlab = "Anni madre",
ylab = "Densità",
col = "blue",
lwd = 2)
abline(v = mean(Anni.madre, na.rm = TRUE), col = "red", lty = 2)
abline(v = median(Anni.madre, na.rm = TRUE), col = "green", lty = 1)
boxplot(Peso,
main = "Boxplot del peso in grammi",
ylab = "Peso (grammi)",
col = "lightblue",
border = "darkblue",
outcol = "red",
outpch = 19)
plot(density(Peso),
main = "Densità del peso",
xlab = "Peso in grammi",
ylab = "Densità",
col = "blue",
lwd = 2)
abline(v = mean(Peso, na.rm = TRUE), col = "red", lty = 2)
abline(v = median(Peso, na.rm = TRUE), col = "green", lty = 1)
boxplot(N.gravidanze,
main = "Boxplot del numero di gravidanze",
ylab = "Numero di gravidanze",
col = "lightblue",
border = "darkblue",
outcol = "red",
outpch = 19)
plot(density(N.gravidanze),
main = "Densità del numero di gravidanze",
xlab = "Numero di gravidanze",
ylab = "Densità",
col = "blue",
lwd = 2)
abline(v = mean(N.gravidanze, na.rm = TRUE), col = "red", lty = 2)
abline(v = median(N.gravidanze, na.rm = TRUE), col = "green", lty = 1)
boxplot(Lunghezza,
main = "Boxplot della lunghezza",
ylab = "Lunghezza (mm)",
col = "lightblue",
border = "darkblue",
outcol = "red",
outpch = 19)
plot(density(Lunghezza),
main = "Densità della lunghezza",
xlab = "Lunghezza (mm)",
ylab = "Densità",
col = "blue",
lwd = 2)
abline(v = mean(Lunghezza, na.rm = TRUE), col = "red", lty = 2)
abline(v = median(Lunghezza, na.rm = TRUE), col = "green", lty = 1)
boxplot(dati$Cranio,
main = "Boxplot del diametro craniale",
ylab = "Diametro craniale (mm)",
col = "lightblue",
border = "darkblue",
outcol = "red",
outpch = 19)
plot(density(Cranio),
main = "Densità del diametro craniale",
xlab = "Diametro craniale",
ylab = "Densità",
col = "blue",
lwd = 2)
abline(v = mean(Cranio, na.rm = TRUE), col = "red", lty = 2)
abline(v = median(Cranio, na.rm = TRUE), col = "green", lty = 1)
boxplot(Gestazione,
main = "Boxplot delle settimane di gestazione",
ylab = "Numero di settimane di gestazione",
col = "lightblue",
border = "darkblue",
outcol = "red",
outpch = 19)
plot(density(Gestazione),
main = "Densità delle settimane di gestazione",
xlab = "Settimane di gestazione",
ylab = "Densità",
col = "blue",
lwd = 2)
abline(v = mean(Gestazione, na.rm = TRUE), col = "red", lty = 2)
abline(v = median(Gestazione, na.rm = TRUE), col = "green", lty = 1)
Per quanto riguarda l’età della madre, la distribuzione nel grafico di densità appare abbastanza normale, con la media e la mediana quasi sovrapposte, confermando l’assenza di particolari anomalie nella distribuzione. I valori di skewness calcolati confermano questa simmetria, con una distribuzione lievemente asimmetrica a destra ma non in modo significativo.
Nel caso del numero di gravidanze, la distribuzione è totalmente asimmetrica verso destra, come evidenziato sia dal boxplot che dal grafico di densità. La maggior parte delle donne ha poche gravidanze (0 o 1), ma ci sono alcune madri con numeri elevati di gravidanze che spostano la distribuzione.
Per quanto riguarda la gurata della gestazione i valori di skewness suggeriscono un’inclinazione verso sinistra, e ciò è evidente anche dal grafico di densità e dal boxplot. Il boxplot mostra diversi outliers verso il basso. Questa asimmetria verso sinistra è confermata dalla densità che mostra una coda prolungata verso i valori bassi, rappresentando gravidanze premature.
Il peso del neonato mostra, invece, una distribuzione praticamente normale. Tuttavia, il boxplot rivela alcuni outlier, soprattutto in basso, che rappresentano neonati con peso molto inferiore alla media. Anche il grafico di densità mostra una leggera coda verso sinistra, che corrisponde a questi outliers, coerentemente alla skewness negativa calcolata.
La lunghezza del neonato presenta anch’essa un leggero squilibrio, con la media e la mediana che sono abbastanza vicine, ma il boxplot rivela alcuni outliers verso il basso. Il grafico di densità mostra una coda simile verso sinistra, confermando la presenza di questi valori anomali.
Infine, anche per quanto riguarda il diametro craniale, il boxplot evidenzia alcuni outliers verso il basso, rappresentando neonati con cranio di dimensioni più piccole. Il grafico di densità mostra una coda verso sinistra, che corrisponde a questi valori anomali. Anche in questo caso, la distribuzione è leggermente inclinata a sinistra, ma la maggior parte dei dati si concentra vicino alla media e alla mediana.
freq_relative_parto <- prop.table(table(Tipo.parto)) * 100
freq_relative_sesso <- prop.table(table(Sesso)) * 100
freq_relative_fumatrici <- prop.table(table(Fumatrici)) * 100
freq_relative_ospedale <- prop.table(table(Ospedale)) * 100
freq_relative_parto## Tipo.parto
## Ces Nat
## 29.12 70.88freq_relative_sesso## Sesso
## F M
## 50.24 49.76freq_relative_fumatrici## Fumatrici
## Non fumatrice Fumatrice
## 95.84 4.16freq_relative_ospedale## Ospedale
## osp1 osp2 osp3
## 32.64 33.96 33.40Date le frequenze relative delle variabili qualitative, osserviamo che le proporzioni di madri non fumatrici (95.84%) e parti naturali (70.88%) indicano delle tendenze “tipiche” della popolazione. D’altra parte la distribuzione equilibrata del sesso e la distribuzione quasi equa delle nascite tra i tre ospedali suggeriscono che il campione è rappresentativo e bilanciato in termini di luogo di nascita e sesso dei neonati.
A questo punto andremo ad analizzare le relazioni che intercorrono tra le variabili rilevanti per la nostra ricerca. Ci focalizzeremo in particolare sulla relazione tra le variabili materne (età, fumo, gravidanze precedenti) con il peso del neonato e se o quanto queste influiscono sulla variabile. Inizieremo dalla variabile del fumo, mostrando un boxplot in cui si evidenziano visivamente le distribuzioni dei due campioni (madri fumatrici e madri non fumatrici). C’è da specificare che vi è una differenza sostanziale nella dimensione dei due gruppi: 104 madri fumatrici contro 2494 non fumatrici. Un campione più grande di madri fumatrici avrebbe sicuramente ridotto la varianza campionaria e aumentato la precisione delle stime, rendendo più facile identificare una differenza significativa se questa fosse esistita. Andremo comunque ad analizzare la relazione che intercorre tra i due gruppi.
boxplot(Peso ~ Fumatrici,
ylab = "Peso del neonato (grammi)",
col = c("lightblue4", "yellow"))
summary(Peso[Fumatrici == "Non fumatrice"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 3000 3300 3286 3620 4900summary(Peso[Fumatrici == "Fumatrice"])## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1780 2940 3175 3236 3440 4930Al netto delle premesse, i neonati delle madri fumatrici hanno una distribuzione del peso leggermente spostata verso valori più bassi rispetto ai neonati delle madri non fumatrici. Tale fatto è evidente dai valori della mediana e della media più bassi. Tutto questo lascerebbe intendere che il fumo potrebbe essere associato a una riduzione del peso medio dei neonati. Per confrontare le distribuzioni dei due gruppi utilizzeremo il test di Mann-Whitney in quanto diversamente dal t.test non richiede l’assunzione di normalità e, inoltre, visto che le distribuzioni che andremo ad analizzare presentano diversi outliers, un test non parametrico potrebbe presentarsi più robusto.
wilcox.test(Peso ~ Fumatrici)##
## Wilcoxon rank sum test with continuity correction
##
## data: Peso by Fumatrici
## W = 138162, p-value = 0.05971
## alternative hypothesis: true location shift is not equal to 0Come anticipato, quasi in controtendenza con i risultati visivi del boxplot, i risultati del test sono influenzati dalla disparità dei campioni. Con un p-value di 0.05971 (lievemente superiore a 0.05) non è possibile riufiutare l’ipotesi nulla, e non c’è quindi evidenza statistica della differenza di peso tra i neonati di madri fumatrici e quelli delle non fumatrici (con un livello di confidenza del 95%). Tuttavia, il p-value di 0.05971 riflette il fatto che, come anche accennato prima, a causa della piccola dimensione del campione delle madri fumatrici non c’è abbastanza potenza statistica per rilevare una differenza significativa. Se avessimo avuto un campione più ampio di madri fumatrici, il p-value sarebbe potuto risultare inferiore a 0.05, permettendoci di rifiutare l’ipotesi nulla e concludere che esiste una differenza significativa nel peso dei neonati tra i due gruppi.
Per ciò che concerne, invece, la relazione tra l’altra variabile relativa alla madre “Anni.madre” e la variabile “Peso”, osserviamo i relativi scatterplot e boxplot che permettono di vedere se sussiste una relazione lineare o non lineare tra le due variabili.
plot(jitter(Anni.madre), Peso,
xlab = "Anni della madre",
ylab = "Peso del neonato (grammi)",
pch = 20, col = "blue")
boxplot(Peso ~ Anni.madre,
xlab = "Anni della madre",
ylab = "Peso del neonato (grammi)",
col = "lightblue")
cor(Anni.madre, Peso)## [1] -0.02377346Già dallo scatterplot e dal boxplot si intuisce che tra le due variabili non vi è nessun tipo di relazione. A conferma di ciò il loro coefficiente di correlazione conferma che le due variabili sono praticamente del tutto indipendenti tra loro.
Infine, per analizzare la relazione tra l’ultima variabile relativa alla madre “N.gravidanze” e la variabile “Peso” utilizzeremo lo stesso procedimento di sopra.
boxplot(Peso ~ N.gravidanze,
xlab = "Numero di gravidanze",
ylab = "Peso del neonato (grammi)",
col = "lightblue")
plot(jitter(N.gravidanze), Peso,
xlab = "Numero di gravidanze",
ylab = "Peso del neonato (grammi)",
pch = 20, col = "blue")
cor(N.gravidanze, Peso)## [1] 0.0024073Dal boxplot osserviamo che mentre nei gruppi con gravidanze da 0 a 5, il peso medio del neonato è relativamente stabile con variazioni poco significative; per numeri di gravidanze superiori a 6, il peso mostra una maggiore variabilità, ma questo potrebbe essere dovuto al ridotto numero di casi piuttosto che a una relazione significativa (solo un osservazione per 7,11,12 gravidanze). Il coefficiente di correlazione conferma quanto osservato, ossia che la relazione tra le due varibili è pressocché inesistente.
Sulla base delle analisi condotte, possiamo affermare che tra il peso del neonato e l’età della madre non vi è correlazione. Stesso discorso per il numero di gravidanze precedenti. Per quanto riguarda il fumo materno, invece, il discorso è leggermente diverso. Come evidenziato prima, i nostri dati potrebbero non essere sufficienti per rilevare un effetto significativo tra i due gruppi.
Creazione del Modello di Regressione
Prima di sviluppare un modello di regressione lineare multipla, andremo ad analizzare le relazioni non osservate fino ad ora che interrocorrono tra la variabile risposta e le altre variabili indipendenti attraverso la visualizzazione di una matrice di correlazione e la sua trasposizione grafica. Per le variabili qualitative gli scatterplot non rendono bene l’idea delle relazioni e le stesse correlazioni perdono significato dal momento che una variabile può avere, nel caso nostro, solo due modalità mentre l’altra è una variabile continua. Per trarre evidenza della correlazione tra la variabile risposta e le variabili qualitative binarie “Sesso” e “Tipo.parto” utilizzeremo dei boxplot. La variabile qualitativa “Ospedale” non la terrò in considerazione in quanto è una variabile relativa al luogo del parto e logicamente non può avere un legame diretto o causale con il peso del neonato. L’unica considerazione che faremo in relazione alla variabile “Ospedale” sarà in riferimento una possibile incidenza di parti cesarei in una determinata struttura.
ggplot(dati, aes(x = Ospedale, fill = Tipo.parto)) +
geom_bar(position = "dodge") +
geom_text(stat = "count", aes(label = ..count..),
position = position_dodge(width = 0.9), vjust = -0.5) +
labs(x = "Ospedale",
y = "Numero di Neonati",
fill = "Tipo di Parto")## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
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(dati[, c("Peso", "Anni.madre", "Gestazione","Lunghezza","Cranio","N.gravidanze")],
lower.panel = panel.smooth, upper.panel = panel.cor,
gap=0)
boxplot(Peso~Sesso)
wilcox.test(Peso~Sesso)##
## Wilcoxon rank sum test with continuity correction
##
## data: Peso by Sesso
## W = 538641, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0boxplot(Peso~Tipo.parto)
wilcox.test(Peso~Tipo.parto)##
## Wilcoxon rank sum test with continuity correction
##
## data: Peso by Tipo.parto
## W = 634748, p-value = 0.5315
## alternative hypothesis: true location shift is not equal to 0Il grafico a barre sovrapposte che mostra il numero di neonati per ospedale distinguendo tra i due tipi di parto, ci evidenzia che la distribuzione dei parti naturali e cesarei è simile tra i tre ospedali, con delle differenze minime. Questo ci suggerisce che nessuno dei tre ospedali è specializzato nei parti cesarei. Come antropologicamente ci aspettavamo, il boxplot delle distribuzioni del peso dei neonati nei due sessi evidenzia come i neonati maschi pesano mediamente di più rispetto alle femmine. Il test di Mann-Whitney conferma quanto osservato. Il valore del p-value indica che la differenza tra le distribuzioni dei due gruppi è statisticamente significativa. D’altra parte, non ci sono prove sufficienti per affermare che i neonati nati con parto naturale e quelli nati con cesareo abbiano distribuzioni del peso significativamente differenti. Sulla variabile “Tipo.parto” faremo una considerazione più avanti.
Entriamo ora nel dettaglio della matrice di correlazione. In primis, anche ai fini del futuro modello che creeremo, le variabili esplicative con un’alta correlazione con la variabile risposta saranno quelle che spiegheranno perlomeno gran parte della variabilità del peso del neonato. Della correlazione tra la variabile “Anni.madre” e “Peso” ne abbiamo parlato quando abbiamo osservato l’influenza delle variabili materne sul peso del neonato, e abbiamo visto come l’età della madre presa singolarmente non ha un impatto significativo sul peso. Stesso discorso per quanto riguarda il numero di gravidanze passate. Diverso è la questione relativa alla variabile “Gestazione”. La durata della gravidanza, infatti, presenta una correlazione positiva e relativamente forte con il peso del neonato (0.59). Questo è un risultato significativo, poiché la durata della gestazione ha un impatto diretto sulla crescita fetale e, quindi, sul peso alla nascita. Questa variabile è probabilmente uno dei predittori più importanti per il modello che andremo a creare. Andando avanti con l’analisi, sia la lunghezza che il diametro craniale del neonato presentano una correlazione positiva e forte con il peso (rispettivamente 0.8 e 0.7), come logicamente ci aspettavamo. Come prima cosa sviluppo un modello di regressione lineare multipla che includa tutte le variabili indipendenti ad eccezione della variabile “Ospedale” in virtù di quanto detto prima.
mod1 <- lm(Peso ~. -Ospedale, data = dati)
summary(mod1)##
## Call:
## lm(formula = Peso ~ . - Ospedale, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1140.23 -181.78 -14.89 160.22 2630.40
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6737.4886 141.4866 -47.619 < 2e-16 ***
## Anni.madre 0.8647 1.1475 0.754 0.4512
## N.gravidanze 11.7148 4.6712 2.508 0.0122 *
## FumatriciFumatrice -31.5829 27.5739 -1.145 0.2522
## Gestazione 32.8391 3.8219 8.592 < 2e-16 ***
## Lunghezza 10.2723 0.3010 34.129 < 2e-16 ***
## Cranio 10.4844 0.4266 24.575 < 2e-16 ***
## Tipo.partoNat 30.2562 12.0994 2.501 0.0125 *
## SessoM 78.0338 11.1925 6.972 3.99e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2491 degrees of freedom
## Multiple R-squared: 0.7279, Adjusted R-squared: 0.727
## F-statistic: 832.9 on 8 and 2491 DF, p-value: < 2.2e-16Analizziamo i risultati del primo modello creato. Il valore di R quadro adattato ci dice che il modello spiega il 72.7% della variabilità del peso del neonato. Complessivamente un buon risultato. Anche l’errore standard residuo di 274,3 grammi è accettabile, indicando che in media i pesi osservati differiscono dai pesi previsti dal modello di quel valore. Vediamo ora le stime dei coefficienti, ossia gli effetti attesi di ciascun predittore sul peso del neonato mantenendo costanti tutte le altre variabili. I coefficienti ritenuti dal modello significativi sono:
N.gravidanze: Dai risultati emerge che ogni gravidanza precedente aumenta il peso di circa 11.7 g. Da specificare che precedentemente, quando abbiamo analizzato la relazione tra questa variabile e il peso del neonato abbiamo concluso che le due variabili fossero praticamente indipendenti tra loro. Nel modello, diversamente, la variabile “N.gravidanze” è emersa rilevante, seppur di poco. Questo perché il modello tiene conto anche delle altre variabili e del loro potenziale effetto indiretto. Se il numero di gravidanze è indirettamente legato al peso attraverso altre variabili, questo effetto viene ignorato nella correlazione.
Gestazione: Molto significativa. Ogni settimana di gestazione aumenta il peso di circa 32.83 g
Lunghezza e Cranio: Fortemente significative, con un aumento rispettivamente di 10.27 g e 10.48 g per ogni millimetro in più.
Tipo.parto: Il coefficiente ci dice che i neonati nati con parto naturale pesano circa 30.25 g in più rispetto a quelli nati con cesareo, tenendo costanti le altre variabili. Il p-value è leggermente inferiore al livello di signficatività.
Sesso: Neonati maschi pesano mediamente 78.03 g in più rispetto alle femmine. Variabile di controllo molto significativa fondamentale per il modello, soprattutto quando si tratta di studi medici.
I coefficienti non significativi:
Anni.madre: con un p-value maggiore di 0.05 l’età della madre non influisce sul peso del neonato.
Fumatrici: la variabile fumo non è emersa significativamente rilevante per il modello. Questo, come anticipato più volte, potrebbe anche essere dovuto al campione sbilanciato.
Adotteremo un metodo di selezione stepwise con un approccio però volto all’eliminazione iniziale (backward), per poi inserire eventuali interazioni e/o non linearità. Partiremo quindi da un modello completo di tutte le variabili indipendenti per eliminare quelle meno significative. Nel prossimo modello che svilupperemo andremo ad escludere la variabile “Anni.madre”. Per quanto riguarda la variabile “Fumatrice”, in virtù di quanto analizzato e in ottica di generalizzare il modello finale a dati non osservati, la includeremo nel modello anche se è risultata non significativa.
mod2 <- update(mod1, ~. -Anni.madre)
summary(mod2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1130.04 -181.29 -16.31 160.59 2635.32
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6708.074 135.984 -49.330 < 2e-16 ***
## N.gravidanze 13.012 4.342 2.997 0.00276 **
## FumatriciFumatrice -31.759 27.570 -1.152 0.24946
## Gestazione 32.541 3.801 8.561 < 2e-16 ***
## Lunghezza 10.272 0.301 34.129 < 2e-16 ***
## Cranio 10.501 0.426 24.648 < 2e-16 ***
## Tipo.partoNat 30.296 12.098 2.504 0.01234 *
## SessoM 78.114 11.191 6.980 3.77e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.3 on 2492 degrees of freedom
## Multiple R-squared: 0.7278, Adjusted R-squared: 0.7271
## F-statistic: 952 on 7 and 2492 DF, p-value: < 2.2e-16Come osserviamo, sia l’R quadro che l’errore standard residuo sono rimasti praticamente invariati rispetto al modello prima. Rispetto al modello precedente, inoltre, l’esclusione delle variabili non significative ha reso il predittore “N.gravidanze” più significativo (p-value = 0.0122 → 0.0027) e la stima del coefficiente è leggermente aumentata (da 11.71 a 13.01 g). Questo è dato dal fatto che la variabile “N.gravidanze” è correlata con un’altra variabile nel modello completo, ora esclusa, nel nostro caso “Anni.madre”. Di conseguenza, quando questa variabile è inclusa parte della variabilità spiegata da “N.gravidanze” viene attribuita ad “Anni.madre”, riducendo quindi la significatività di “N.gravidanze”.
Vorrei aprire una piccola parentesi per la variabile “Tipo.parto”. Dal punto di vista pratico e biologico il tipo di parto è più probabilmente una conseguenza del peso e non una causa. La variabile risulta statisticamente significativa, ma il suo contributo al modello potrebbe essere logicamente debole o poco significativo. Per tale motivo, ma anche perché la significatività risulta comunque leggermente inferiore al livello standard 0.05 e nell’ottica di un modello più parsimonioso, procederemo a sviluppare un nuovo modello senza la variabile in questione.
mod3 <- update(mod2, ~. -Tipo.parto)
summary(mod3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1150.3 -181.3 -15.7 163.0 2636.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.6714 135.7178 -49.232 < 2e-16 ***
## N.gravidanze 12.7185 4.3450 2.927 0.00345 **
## FumatriciFumatrice -30.4634 27.5948 -1.104 0.26972
## Gestazione 32.5914 3.8051 8.565 < 2e-16 ***
## Lunghezza 10.2341 0.3009 34.011 < 2e-16 ***
## Cranio 10.5359 0.4262 24.718 < 2e-16 ***
## SessoM 78.1713 11.2028 6.978 3.83e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2493 degrees of freedom
## Multiple R-squared: 0.7271, Adjusted R-squared: 0.7265
## F-statistic: 1107 on 6 and 2493 DF, p-value: < 2.2e-16Dai risultati emerge ancora una volta che l’ottimizzazione non ha penalizzato il modello. L’R quadro adattato ha subito una insignificante riduzione e l’errore residuo è rimasto invariato. Ad ora, possiamo dire che il modello 3 è il miglior modello, in quanto presenta le sole variabili indipendenti significative in aggiunta al predittore “Fumatrici”. Creeremo, tuttavia, ulteriori modelli provando a catturare anche eventuali interazioni o effetti non lineari.
Nella matrice di correlazione e dispersione precedentemente analizzata si intravedeva una leggera non linearità nello scatterplot del peso in relazione alle settimane di gestazione. Proviamo a catturare nel modello anche tali effetti non lineari, aggiungendo il quadrato della variabile “Gestazione”.
mod4 <- update(mod3, ~. + I(Gestazione^2))
summary(mod4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso + I(Gestazione^2), data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1145.00 -180.96 -13.12 165.00 2659.10
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4669.1027 898.3246 -5.198 2.18e-07 ***
## N.gravidanze 12.7990 4.3416 2.948 0.00323 **
## FumatriciFumatrice -29.2442 27.5772 -1.060 0.28904
## Gestazione -79.7741 49.7260 -1.604 0.10878
## Lunghezza 10.3386 0.3042 33.989 < 2e-16 ***
## Cranio 10.6312 0.4280 24.841 < 2e-16 ***
## SessoM 75.9814 11.2351 6.763 1.68e-11 ***
## I(Gestazione^2) 1.4999 0.6618 2.266 0.02352 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.4 on 2492 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.7269
## F-statistic: 951.4 on 7 and 2492 DF, p-value: < 2.2e-16Nel nuovo modello il termine “Gestazione^2” è significativo anche se di poco, suggerendo che la relazione tra gestazione e peso del neonato non è perfettamente lineare. D’altra parte, dopo l’aggiunta del nuovo termine l’effetto lineare di “Gestazione” non è più significativo. Inoltre, tale effetto è diventato negativo, in quanto si compensa con l’effetto del termine quadratico. Comunque, il modello non migliora significativamente in termini di R quadro ed errore standard residuo, facendoci intuire che l’aggiunta del termine quadratico ha un impatto minimo sulla validità del modello.
Un’altro modello che potremmo sviluppare, potrebbe includere le interazioni tra le variabili, in modo da modellare situazioni in cui l’effetto di una variabile dipende dal livello di un’altra variabile. Nel nostro caso, ad esempio, la relazione tra il peso e la lunghezza del neonato potrebbe variare a seconda della durata della gestazione. Aggiungeremo quindi al modello l’interazione tra la variabile “Gestazione” e “Lunghezza”.
mod5 <- update(mod3, ~. + Gestazione * Lunghezza)
summary(mod5)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso + Gestazione:Lunghezza, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1134.35 -179.93 -11.01 168.13 2650.24
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.011e+03 9.204e+02 -2.185 0.029015 *
## N.gravidanze 1.326e+01 4.324e+00 3.067 0.002189 **
## FumatriciFumatrice -2.742e+01 2.746e+01 -0.998 0.318155
## Gestazione -9.320e+01 2.481e+01 -3.757 0.000176 ***
## Lunghezza -5.090e-02 2.027e+00 -0.025 0.979969
## Cranio 1.075e+01 4.263e-01 25.232 < 2e-16 ***
## SessoM 7.246e+01 1.120e+01 6.469 1.18e-10 ***
## Gestazione:Lunghezza 2.717e-01 5.297e-02 5.130 3.12e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 273.2 on 2492 degrees of freedom
## Multiple R-squared: 0.73, Adjusted R-squared: 0.7292
## F-statistic: 962.5 on 7 and 2492 DF, p-value: < 2.2e-16Come osserviamo, il termine di interazione risulta altamente significativo e con un impatto importante sul peso del neonato. Il coefficiente positivo (0.26878) ci suggerisce un aumento del peso quando entrambe le variabili aumentano. Tuttavia, l’aggiunta dell’interazione ha fatto perdere di significatività la variabile Lunghezza. Infine, l’R quadro e l’errore standard residuo sono rimasti invariati. Possiamo concludere che l’interazione tra la variabile “Gestazione” e la variabile “Lunghezza” non va aggiunta al modello.
Facciamo lo stesso ragionamento per l’interazione tra la variabile “Gestazione” e “Cranio” per valutare se l’effetto della circonferenza cranica sul peso dipende dalla durata della gestazione.
mod6 <- update(mod3, ~. + Gestazione * Cranio)
summary(mod6)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Sesso + Gestazione:Cranio, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1138.0 -180.6 -11.7 167.1 2692.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -212.4233 1106.5877 -0.192 0.84779
## N.gravidanze 13.3615 4.3174 3.095 0.00199 **
## FumatriciFumatrice -27.7272 27.4141 -1.011 0.31191
## Gestazione -139.9398 29.5351 -4.738 2.28e-06 ***
## Lunghezza 10.4566 0.3013 34.708 < 2e-16 ***
## Cranio -9.7816 3.4754 -2.815 0.00492 **
## SessoM 72.2366 11.1734 6.465 1.21e-10 ***
## Gestazione:Cranio 0.5318 0.0903 5.890 4.38e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 272.8 on 2492 degrees of freedom
## Multiple R-squared: 0.7309, Adjusted R-squared: 0.7301
## F-statistic: 966.8 on 7 and 2492 DF, p-value: < 2.2e-16Il termine di interazione è altamente significativo, suggerendoci che esiste un effetto combinato reale tra queste due variabili sul peso del neonato. In particolare, anche qui il coefficiente positivo (0.5262) ci suggerisce un aumento del peso quando entrambe le variabili aumentano. Vediamo come gli effetti lineari negativi di “Gestazione” e “Cranio” si compensano con il termine di interazione. Diversamente da prima, l’aggiunta del termine di interazione non ha fatto perdere di significatività l’effetto delle variabili principali. Questo porta il modello 6 ad essere presumibilmente tra i modelli migliori secondo i criteri di valutazione che utilizzeremo. In generale, possiamo dire che la significatività di entrambi gli effetti combinati riflettono un fenomeno plausibile: gestazioni più lunghe permettono al cranio (e al corpo in generale) di crescere di più, portando quindi a un peso maggiore. In generale, Un metodo per capire se aggiungere o meno una variabile al modello è quello dell’ANOVA, ossia un test utile a confrontare due modelli, per valutare se il modello più complesso (con più termini) è significativamente migliore di uno più semplice.
anova(mod4, mod3)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Sesso + I(Gestazione^2)
## Model 2: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 187587020
## 2 2493 187973654 -1 -386634 5.1362 0.02352 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod5, mod3)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Sesso + Gestazione:Lunghezza
## Model 2: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 186009153
## 2 2493 187973654 -1 -1964501 26.319 3.117e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod6, mod3)## Analysis of Variance Table
##
## Model 1: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Sesso + Gestazione:Cranio
## Model 2: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Sesso
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 2492 185392735
## 2 2493 187973654 -1 -2580919 34.692 4.383e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Il confronto effettuato è con il modello 3, ossia il modello con le sole variabili principali significative. Dai risultati del test di ANOVA, osserviamo come nei modelli con l’interazioni Gestazione:Lunghezza, e Gestazione:Cranio, entrambi i termini di interazione sono significativi. Stesso discorso per il termine quadratico nel modello 4. Questo ci suggerisce che tutti e tre i termini portano un contributo significativo al modello. Tuttavia, mentre nel modello 4 l’effetto del termine è relativamente piccolo (l’RSS diminuisce di poco), nei modelli 5 e 6 gli effetti sono molto rilevanti (le differenze di RSS tra i modelli sono molto maggiori). Per concludere i test di ANOVA ci suggeriscono di lasciare i termini aggiunti, tuttavia, al fine di scegliere il modello migliore, andremo comunque ad utilizzare altri criteri di confronto, ossia l’AIC (Akaike Information Criterion) e il BIC (Bayesian Information Criterion).
** Selezione del Modello Ottimale
AIC(mod6, mod5, mod4, mod3, mod2, mod1)## df AIC
## mod6 9 35147.55
## mod5 9 35155.84
## mod4 9 35176.96
## mod3 8 35180.11
## mod2 9 35175.83
## mod1 10 35177.26BIC(mod6, mod5, mod4, mod3, mod2, mod1)## df BIC
## mod6 9 35199.96
## mod5 9 35208.26
## mod4 9 35229.38
## mod3 8 35226.70
## mod2 9 35228.24
## mod1 10 35235.50Come sappiamo, l’AIC valuta quanto bene un modello si adatta ai dati, penalizzando però la complessità del modello (ovvero il numero di parametri). BIC, d’altra parte, è simile all’AIC, ma introduce una penalità più forte che dipende dalla dimensione del campione. Il modello con AIC più basso risulta essere il modello 6, quindi secondo tale criterio il modello 6 è quello da preferire. Si tratta del modello che presenta il termine di interazione Gestazione:Cranio. Allo stesso modo, il modello migliore per il criterio BIC è sempre il modello 6. Anche il modello 3, quello con le sole variabili significative, è un buon modello nonostante i test ci dicono che ci sono modelli migliori. La variabilità spiegata della variabile risposta è molto simile agli altri modelli, ma questo a differenza degli altri è più parsimonioso presentando un regressore in meno. Di conseguenza, prenderemo come modelli migliori il modello 3 e 6. Per verificare che non ci siano problemi di multicollinearità utilizziamo la funzione VIF, che serve appunto per calcolare gli indicatori di multicollinearità, che per andare bene devono essere al di sotto di 5. Per il modello 6, che presenta un termine di interazione, utilizzo l’opzione “predictor” per calcolare il VIF in modo specifico per ciascun predittore individuale (senza interazioni), dal momento che inevitabilmnete i termini dell’interazione hanno forte correlazione con i loro termini principali. Per valutare la qualità dei due modelli utilizzeremo anche un’altra metrica, ossia la MSE, che rappresenta l’errore quadratico medio tra le osservazioni e le predizioni.
library(car)## Warning: il pacchetto 'car' è stato creato con R versione 4.4.2## Caricamento del pacchetto richiesto: carData## Warning: il pacchetto 'carData' è stato creato con R versione 4.4.2vif(mod3)## N.gravidanze Fumatrici Gestazione Lunghezza Cranio Sesso
## 1.026120 1.006607 1.675575 2.078644 1.624603 1.040271vif(mod6, type = "predictor")## GVIFs computed for predictors## GVIF Df GVIF^(1/(2*Df)) Interacts With
## N.gravidanze 1.026776 1 1.013300 --
## Fumatrici 1.006896 1 1.003442 --
## Gestazione 2.144164 3 1.135559 Cranio
## Lunghezza 2.111850 1 1.453221 --
## Cranio 2.144164 3 1.135559 Gestazione
## Sesso 1.048800 1 1.024109 --
## Other Predictors
## N.gravidanze Fumatrici, Gestazione, Lunghezza, Cranio, Sesso
## Fumatrici N.gravidanze, Gestazione, Lunghezza, Cranio, Sesso
## Gestazione N.gravidanze, Fumatrici, Lunghezza, Sesso
## Lunghezza N.gravidanze, Fumatrici, Gestazione, Cranio, Sesso
## Cranio N.gravidanze, Fumatrici, Lunghezza, Sesso
## Sesso N.gravidanze, Fumatrici, Gestazione, Lunghezza, Cranioresidui_mod3 <- residuals(mod3)
mse_mod3 <- mean(residui_mod3^2)
mse_mod3## [1] 75189.46residui_mod6 <- residuals(mod6)
mse_mod6 <- mean(residui_mod6^2)
mse_mod6## [1] 74157.09Dai risultati ottenuti emerge che entrambi i modelli sono bilanciati, dato che nel caso del modello 6 tutti i predittori hanno un GVIF^(1/(2*Df)) inferiore a 2, mentre nel modello 3 il VIF inferiore a 5. La multicollinearità nei modelli scelti, quindi, è bassa e non preoccupante. Per quanto rigurda l’errore quadratico medio, un valore più basso (nel nostro caso si presenta per il modello 6) indica che il modello ha una capacità predittiva migliore, con errori mediamente più piccoli. Tuttavia, la differenza tra i due valori è relativamente piccola, suggerendoci che il miglioramento apportato al modello 6 non è particolarmente significativo.
Utilizzeremo ora la funzione stepAIC, attraverso il quale R effettuerà una selezione stepwise delle variabili del modello completo, che faremo basare sul criterio di informazione di BIC.
library(MASS)
n <- nrow(dati)
stepwise_mod <- MASS::stepAIC(mod1,
direction = "both",
k=log(n))## Start: AIC=28132.98
## Peso ~ (Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso) - Ospedale
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 42734 187501837 28126
## - Fumatrici 1 98728 187557831 28127
## - Tipo.parto 1 470578 187929681 28131
## - N.gravidanze 1 473298 187932401 28132
## <none> 187459103 28133
## - Sesso 1 3658038 191117141 28174
## - Gestazione 1 5555847 193014950 28198
## - Cranio 1 45450123 232909226 28668
## - Lunghezza 1 87653102 275112206 29084
##
## Step: AIC=28125.73
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso
##
## Df Sum of Sq RSS AIC
## - Fumatrici 1 99840 187601677 28119
## - Tipo.parto 1 471817 187973654 28124
## <none> 187501837 28126
## - N.gravidanze 1 675718 188177555 28127
## + Anni.madre 1 42734 187459103 28133
## - Sesso 1 3665908 191167745 28166
## - Gestazione 1 5514506 193016343 28190
## - Cranio 1 45711714 233213551 28663
## - Lunghezza 1 87642114 275143951 29077
##
## Step: AIC=28119.23
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto +
## Sesso
##
## Df Sum of Sq RSS AIC
## - Tipo.parto 1 463870 188065546 28118
## <none> 187601677 28119
## - N.gravidanze 1 651066 188252743 28120
## + Fumatrici 1 99840 187501837 28126
## + Anni.madre 1 43845 187557831 28127
## - Sesso 1 3649259 191250936 28160
## - Gestazione 1 5444109 193045786 28183
## - Cranio 1 45758101 233359778 28657
## - Lunghezza 1 88054432 275656108 29074
##
## Step: AIC=28117.58
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 188065546 28118
## - N.gravidanze 1 623141 188688687 28118
## + Tipo.parto 1 463870 187601677 28119
## + Fumatrici 1 91892 187973654 28124
## + Anni.madre 1 45044 188020502 28125
## - Sesso 1 3655292 191720838 28158
## - Gestazione 1 5464853 193530399 28181
## - Cranio 1 46108583 234174130 28658
## - Lunghezza 1 87632762 275698308 29066summary(stepwise_mod)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.44 -180.81 -15.58 163.64 2639.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.1445 135.7229 -49.226 < 2e-16 ***
## N.gravidanze 12.4750 4.3396 2.875 0.00408 **
## Gestazione 32.3321 3.7980 8.513 < 2e-16 ***
## Lunghezza 10.2486 0.3006 34.090 < 2e-16 ***
## Cranio 10.5402 0.4262 24.728 < 2e-16 ***
## SessoM 77.9927 11.2021 6.962 4.26e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1328 on 5 and 2494 DF, p-value: < 2.2e-16Il modello scelto dalla funzione stepAIC coincide con il nostro modello 3, ad eccezione per il regressore “Fumatrici”. Questo Perché la funzione non considera automaticamente i termini di interazione o quadratici dal momento che non sono inclusi nel modello di partenza. Sulla base di quando detto, sceglieremo come migliore il modello 3, quello senza interazione, anche se i test fatti ci dicono diversamente, in quanto a mio parere è il modello che presenta il miglior compromesso tra sintesi e perdita di informazione. Inoltre, seguendo il principio del Rasoio di Occam, modelli più semplici sono preferiti a modelli complessi e non c’è necessità di usare parametri addizionali se essi non sono strettamente necessari.
Andremo ora ad analizzare i residui, e quindi la parte erratica, del modello 3.
** Analisi della Qualità del Modello
Come sappiamo, i residui del modello devono essere “puliti”, ossia devono rispettare le assunzioni di base: normalità, varianza costante (omoschedasticità), indipendenza e media di zero.
par(mfrow = c(2,2))
plot(mod3)
Il grafico in alto a sinistra, che mostra i residui standardizzati
rispetto i valori predefiniti, i residui sembrano distribuiti in modo
abbastanza casuale, ma ci sono alcuni outliers evidenti (i casi 1551,
155, 1306). Non sembra esserci un pattern sistematico (giusto una
leggerissima curvatura), quindi l’assunzione di linearità e
omoschedasticità è ragionevolmente soddisfatta, anche se gli outliers
visti potrebbero influire.
Nel Q-Q plot in alto a destra, invece, vengono messi in relazione i residui con i quantili di una distribuzione normale. Nel nostro caso sembra che i punti siano correttamente disposti sopra la retta, al netto di qualche discordanza sopratttutto nella coda superiore (ci sono le 3 osservazioni menzionate prima che risutano più lontane delle altre).
Nel Scale-Location plot in basso a sinistra, che mostra i residui standardizzati in funzione della radice quadrata dei valori predetti, si può notare una leggerissima curvatura e gli stessi tre casi estremi.
Nell’ultimo grafico si visualizzano i potenziali valori influenti, ovvero quei valori con alto leverage e/o molto outliers. Nel nostro caso i residui di osservazioni potenzialmente influenti sulle stime di regressione sono quelli delle osservazioni già menzionate. In particolare l’osservazione 1551 ha una distanza di Cook che supera la soglia di avvertimento (0.5) ma non la soglia di allarme (1).
Per andare ad indagare i valori di leva e gli outliers anche numericamente utilizzeremo la funzione hatvalues.
lev <- hatvalues(mod3)
plot(lev)
p <- sum(lev)
n <- length(lev)
soglia = 2*p/n
abline(h = soglia,col=2)
lev[lev > soglia]## 13 15 34 67 89 99
## 0.005701303 0.007063521 0.006754498 0.005892580 0.012910981 0.010439294
## 101 105 106 120 128 131
## 0.007547922 0.010619066 0.014502002 0.010038865 0.011383848 0.007234866
## 134 140 151 155 161 182
## 0.007594012 0.011383097 0.010954678 0.007236095 0.020448546 0.011314658
## 194 204 206 220 234 242
## 0.010838364 0.014576811 0.009483391 0.007440530 0.010852299 0.010235987
## 251 279 294 296 306 310
## 0.010893497 0.010509023 0.005974827 0.010171095 0.010856933 0.028812882
## 312 321 335 378 391 413
## 0.013203779 0.010712320 0.010921679 0.015934740 0.010945631 0.010550439
## 424 442 445 473 492 516
## 0.010778599 0.016100607 0.007513013 0.011311052 0.008306692 0.013188913
## 538 557 567 572 582 587
## 0.012114256 0.010675402 0.010349049 0.010612745 0.011720774 0.008412394
## 592 593 638 656 658 668
## 0.006384761 0.010423807 0.006695658 0.006018345 0.011305325 0.011515183
## 684 697 699 703 748 750
## 0.008836176 0.005872566 0.011104503 0.010765338 0.008580665 0.006965520
## 757 758 765 805 828 913
## 0.008193218 0.011588686 0.006075696 0.014369467 0.007252180 0.005643794
## 928 932 946 947 956 984
## 0.022756502 0.010471866 0.006929519 0.008436064 0.007853821 0.010404872
## 985 1014 1017 1026 1037 1051
## 0.007132127 0.008519880 0.011236061 0.011627171 0.010357179 0.010765729
## 1067 1091 1106 1110 1118 1130
## 0.008473411 0.008952316 0.006006209 0.010413324 0.010363532 0.032004825
## 1170 1175 1181 1188 1219 1227
## 0.010797117 0.010504719 0.005679157 0.006482056 0.030876924 0.011904008
## 1238 1248 1262 1271 1273 1282
## 0.005910463 0.014637569 0.012913716 0.010119368 0.007085833 0.010434740
## 1285 1291 1293 1311 1321 1326
## 0.012201635 0.006160032 0.006100177 0.009678501 0.009353884 0.011062394
## 1333 1357 1368 1379 1385 1397
## 0.011373022 0.006965052 0.011081209 0.010729161 0.012637438 0.011242355
## 1398 1400 1410 1411 1415 1425
## 0.010898342 0.005925138 0.012145967 0.008128147 0.010394855 0.010298475
## 1426 1428 1429 1443 1449 1450
## 0.013000422 0.008245853 0.021763751 0.011273205 0.011022883 0.015264031
## 1458 1473 1480 1505 1512 1525
## 0.010509023 0.010704585 0.011506315 0.013426687 0.011239772 0.010429984
## 1537 1551 1553 1556 1576 1583
## 0.011945191 0.048894280 0.008506884 0.005958271 0.010627832 0.012627876
## 1593 1610 1619 1626 1652 1660
## 0.005669713 0.008727229 0.015088232 0.011104701 0.011301902 0.011289554
## 1672 1686 1691 1701 1712 1718
## 0.010903644 0.009349482 0.010807413 0.010857075 0.006998227 0.007038575
## 1720 1727 1761 1763 1780 1781
## 0.010997269 0.013376581 0.011311923 0.010748635 0.025543096 0.016923009
## 1789 1809 1827 1902 1906 1920
## 0.010797999 0.008710061 0.006077091 0.010576575 0.010376577 0.014344029
## 1929 1933 1971 1977 2003 2016
## 0.012560840 0.010996875 0.012328192 0.006934103 0.011147851 0.013531116
## 2040 2046 2049 2086 2089 2101
## 0.011541669 0.014286949 0.010439294 0.013303759 0.015640622 0.011513947
## 2110 2114 2115 2120 2140 2145
## 0.010608705 0.013332484 0.011775621 0.018660016 0.006263186 0.010268351
## 2146 2148 2149 2157 2175 2200
## 0.005833949 0.007983119 0.013606612 0.005967649 0.032596366 0.011670531
## 2202 2216 2220 2221 2224 2237
## 0.010368796 0.008120263 0.013757017 0.021754315 0.005847734 0.010698549
## 2238 2244 2245 2256 2257 2270
## 0.010965046 0.006995300 0.013619106 0.010582603 0.006185756 0.011002949
## 2282 2285 2307 2317 2337 2353
## 0.010998766 0.010708229 0.013979507 0.007749993 0.014207152 0.012972794
## 2359 2361 2408 2412 2422 2437
## 0.010106830 0.010626804 0.009708738 0.010412911 0.021698506 0.023956769
## 2450 2452 2458 2459 2465 2471
## 0.010627186 0.023838748 0.008506091 0.010213071 0.011317924 0.021047568
## 2478
## 0.005775174sum(lev > soglia)## [1] 211plot(rstudent(mod3))
car::outlierTest(mod3)## rstudent unadjusted p-value Bonferroni p
## 1551 10.039719 2.8060e-23 7.0149e-20
## 155 5.022108 5.4723e-07 1.3681e-03
## 1306 4.823102 1.4986e-06 3.7465e-03Abbiamo trovato 211 osservazioni i cui residui presentano un alto leverage e 3 outliers (gli stessi visti prima). Per considerare entrambe le cose contemporaneamente utilizziamo la distanza di Cook.
cook <- cooks.distance(mod3)
plot(cook)
max(cook)## [1] 0.7117513Come notiamo la distanza di Cook massima registrata è di 0.71 che
come anche visto prima viene riportata dall’osservazione 1551. In ogni
caso, tale osservazione, anche se supera la soglia di avvertimento di
0.5, non supera la soglia di allarme di 1, quindi possiamo dire che
l’impatto non sia critico per il modello complessivo e non sia influente
sulle stime di regressione.
Andiamo ad osservare i valori di questa osservazione e proviamo a capire
perché ha una distanza di Cook così elevata. Il neonato in questione
presenta la seconda lunghezza più bassa di tutto il dataset, nonostante
il peso sia maggiore di circa 1 kg rispetto alla media.
osservazioni <- which(dati$Cranio > dati$Lunghezza)
dati[osservazioni, ]## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1551 35 1 Non fumatrice 38 4370 315 374
## Tipo.parto Ospedale Sesso
## 1551 Nat osp3 FTra l’altro, abbiamo appena visto che è l’unica osservazione del dataset in cui la lunghezza del neonato è minore del diametro craniale. Dal punto di vista biologico è difficile da pensare una situazione simile, tuttavia, non essendo un medico e in virtù di quanto detto prima, continuiamo con l’analisi del modello.
Osservati i grafici diagnostici, andremo ora ad effettuare sui residui i test di omoschedasticità di Breusch-Pagan, di autocorrelazione di Durbin-Watson e infine il test di Shapiro-Wilk.
library(lmtest)## Warning: il pacchetto 'lmtest' è stato creato con R versione 4.4.1## Caricamento del pacchetto richiesto: zoo##
## Caricamento pacchetto: 'zoo'## I seguenti oggetti sono mascherati da 'package:base':
##
## as.Date, as.Date.numericbptest(mod3)##
## studentized Breusch-Pagan test
##
## data: mod3
## BP = 89.798, df = 6, p-value < 2.2e-16dwtest(mod3)##
## Durbin-Watson test
##
## data: mod3
## DW = 1.9542, p-value = 0.126
## alternative hypothesis: true autocorrelation is greater than 0shapiro.test(residuals(mod3))##
## Shapiro-Wilk normality test
##
## data: residuals(mod3)
## W = 0.9741, p-value < 2.2e-16plot(density(residuals(mod3)))
Il test di Breusch-Pagan indica la presenza di eteroschedasticità.
Questo perché, come visto prima nel grafico Residuals vs Fitted ci sono
diversi outliers evidenti che potrebbero aver influito sui risultati del
test. Tuttavia, come anche affermato prima, graficamente, non sembra
esserci un pattern sistematico, e l’assunzione di omoschedasticità può
essere considerata ragionevolmente soddisfatta.
Dai risultati del test di Durbin-Watson, invece, risulta che i residui
sono indipendenti tra loro e non ci sono problemi di autocorrelazione.
Infine, dall’ultimo test di normalità effettuato, risulta che i residui
non seguono una distribuzione normale. Dal grafico di densità si evince
che il problema sta, come sappiamo, nelle code della distribuzione.
Questo ci interessa relativamente poco in quanto il resto della
distribuzione assomiglia proprio ad una distribuzione normale.
Nonostante i due test citati non sono andati a buon fine per le
problematiche viste, il modello 3 rimane il migliore modello per
descrivere i dati e fare previsioni.
3 Previsioni e Risultati
Ad esempio, come richiesto, possiamo stimare il peso di una neonata considerando una madre alla terza gravidanza, non fumatrice e che partorirà alla 39esima settimana. Per “Lunghezza” e “Cranio” useremo la media dei valori osservati, ossia rispettivamente 495 mm e 340 mm.
new_data <- data.frame(
N.gravidanze = 3,
Fumatrici = "Non fumatrice",
Gestazione = 39,
Lunghezza = 495,
Cranio = 340,
Sesso = "F"
)
predict(mod3, newdata = new_data)## 1
## 3275.609Il peso previsto nel caso citato risulta molto vicino alla media. Questo perché abbiamo assodato che il neonato sia alto e abbia il diametro craniale uguale alla media dei valori nel dataset.
4 Visualizzazioni
Genereremo una serie di grafici con diverse combinazioni di variabili, rappresentando le relazioni tridimensionali tra la variabile dipendente “Peso” e le variabili indipendenti “Gestazione”, “Lunghezza”, “Cranio” con raggruppamenti per “Fumatrici” e “Sesso”.
library(rgl)## Warning: il pacchetto 'rgl' è stato creato con R versione 4.4.2library(mgcv)## Caricamento del pacchetto richiesto: nlme## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.scatter3d(
x = dati$Gestazione,
y = dati$Cranio,
z = dati$Peso,
xlab = "Settimane di gestazione",
ylab = "Diametro craniale (mm)",
zlab = "Peso (grammi)")
rglwidget()
scatter3d(
x = dati$Gestazione,
y = dati$Lunghezza,
z = dati$Peso,
xlab = "Settimane di gestazione",
ylab = "Lunghezza (mm)",
zlab = "Peso (grammi)")
rglwidget()
scatter3d(
x = dati$Gestazione,
y = dati$Cranio,
z = dati$Peso,
groups = dati$Fumatrici,
xlab = "Settimane di gestazione",
ylab = "Diametro craniale (mm)",
zlab = "Peso (grammi)")
rglwidget()
scatter3d(
x = dati$Gestazione,
y = dati$Lunghezza,
z = dati$Peso,
groups = dati$Sesso,
xlab = "Settimane di gestazione",
ylab = "Lunghezza (mm)",
zlab = "Peso (grammi)")
rglwidget()
scatter3d(
x = dati$Gestazione,
y = dati$Cranio,
z = dati$Peso,
groups = dati$Sesso,
xlab = "Settimane di gestazione",
ylab = "Diametro craniale (mm)",
zlab = "Peso (grammi)")
rglwidget()
scatter3d(
x = dati$Gestazione,
y = dati$Lunghezza,
z = dati$Peso,
groups = dati$Fumatrici,
xlab = "Settimane di gestazione",
ylab = "Lunghezza (mm)",
zlab = "Peso (grammi)")
rglwidget()
I primi due grafici ci permettono di osservare come due tra le variabili indipendenti influenzano la variabile “Peso”, mentre negli altri quattro sono stati separati i due gruppi per le variabili binarie.