Modello statistico per la previsione del peso dei neonati
Su richiesta dell’azienda Neonatal Health Solutions è stato sviluppato un modello statistico in grado di prevedere il peso dei neonati. Le fasi che hanno portato alla realizzazione del modello sono illustrate in questo documento.
Dopo un’analisi descrittiva dei dati forniti, si effettueranno test statistici al fine di verificare alcune ipotesi richieste dall’azienda committente. Infine, dopo aver costruito il modello e dopo averne valutata la qualità, si svolgeranno predizioni su dati nuovi e si analizzeranno i risultati ottenuti.
Importazione ed esplorazione dei dati
Questa sezione è dedicata all’importazione e a una prima esplorazione dei dati.
Iniziamo importando i dati da un file CSV fornito dall’azienda in un DataFrame R chiamato df.
df <- read.csv("neonati.csv", sep = ",", stringsAsFactors = T)La cella che segue mostra le dimensioni iniziali del dataset.
dim(df)## [1] 2500 10Il dataset consiste in 2500 osservazioni, una per ogni neonato, di cui sono state misurate 10 variabili. La cella che segue mostra una panoramica delle variabili in esame.
summary(df)## Anni.madre N.gravidanze Fumatrici Gestazione
## Min. : 0.00 Min. : 0.0000 Min. :0.0000 Min. :25.00
## 1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:0.0000 1st Qu.:38.00
## Median :28.00 Median : 1.0000 Median :0.0000 Median :39.00
## Mean :28.16 Mean : 0.9812 Mean :0.0416 Mean :38.98
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:0.0000 3rd Qu.:40.00
## Max. :46.00 Max. :12.0000 Max. :1.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. :390Notiamo subito la presenza di anomalie nelle colonne Anni.madre e Fumatrici. Questo e il sospetto della presenza di record anomali ci ha spinto a effettuare delle operazioni di pulizia del dataset, che approfondiremo nella prossima sezione.
Pulizia del dataset
In questa sezione ci occuperemo delle operazioni di pulizia del dataset fornito.
Iniziamo dalla conversione della colonna Fumatrici in una variabile categorica. Essa, infatti, è una variabile di controllo binaria che indica se la madre fuma o meno.
df$Fumatrici <- factor(df$Fumatrici, levels = c(0, 1), labels = c("No", "Sì"))Per quanto riguarda Anni.madre, osserviamo che il valore minimo osservato è 0. Dal momento che è impossibile per una neonata partorire, l’ipotesi è che si tratti di un errore di battitura nell’inserimento dei dati. Escludiamone altri filtrando questi valori e sostituendoli con la mediana dei valori validi, assumendo che il minimo valore valido sia 12 anni.
mediana_anni <- median(df$Anni.madre[df$Anni.madre >= 12], na.rm = TRUE)
df$Anni.madre <- ifelse(df$Anni.madre<12, mediana_anni, df$Anni.madre)La cella che segue mostra un nuovo riepilogo dei dati validi.
summary(df)## Anni.madre N.gravidanze Fumatrici Gestazione Peso
## Min. :13.00 Min. : 0.0000 No:2396 Min. :25.00 Min. : 830
## 1st Qu.:25.00 1st Qu.: 0.0000 Sì: 104 1st Qu.:38.00 1st Qu.:2990
## Median :28.00 Median : 1.0000 Median :39.00 Median :3300
## Mean :28.19 Mean : 0.9812 Mean :38.98 Mean :3284
## 3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:40.00 3rd Qu.:3620
## Max. :46.00 Max. :12.0000 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 osp3:835
## Mean :494.7 Mean :340
## 3rd Qu.:510.0 3rd Qu.:350
## Max. :565.0 Max. :390Malgrado questi valori siano plausibili dal punto di vista logico, non ce la sentiamo di escludere la presenza di record che in seguito possono inficiare negativamente il modello. Ci riferiamo, in particolare, alle variabili Gestazione, Lunghezza e Cranio in relazione al Peso.
Le celle che seguono mostrano in un grafico la relazione tra queste tre variabili e il peso del neonato.
plot(df$Gestazione, df$Peso, pch=20)
plot(df$Lunghezza, df$Peso, pch=20)
plot(df$Cranio, df$Peso, pch=20)
Notiamo la presenza di punti isolati che si staccano dal corpo principale dei dati. In particolare, nel grafico Gestazione-Peso, evidenziamo un neonato alla 35° settimana con peso intorno ai 4500 grammi, così come un neonato dello stesso peso lungo poco più di 30 centimetri nel grafico Lunghezza-Peso o anche i due neonati sotto i 2 kg di peso con una circonferenza cranica superiore ai 35 cm.
Per capire come comportarci, dobbiamo isolare questi valori. Innanzitutto, aggiungiamo le colonne Lunghezza_CL e Cranio_CL, in modo da dividere in fasce la lunghezza e il cranio dei neonati.
min_limit <- round(min(df$Lunghezza), digits = -2)
max_limit <- round(max(df$Lunghezza), digits = -1)
df$Lunghezza_CL <- cut(df$Lunghezza, breaks = seq(min_limit, max_limit+20, 20))min_limit <- round(min(df$Cranio), digits = -1)
max_limit <- round(max(df$Cranio), digits = -1)
df$Cranio_CL <- cut(df$Cranio, breaks = seq(min_limit-10, max_limit, 10))La funzione is_outlier, definita nella cella che segue, prende in ingresso un dataframe ed etichetta come outlier i pesi che si trovano al di sotto del primo percentile e al di sopra del 99° percentile.
is_outlier <- function(dati) {
if(nrow(dati) < 5) {
dati$is_outlier <- FALSE
return(dati)
}
s_bassa <- quantile(dati$Peso, 0.01)
s_alta <- quantile(dati$Peso, 0.99)
dati$is_outlier <- dati$Peso < s_bassa | dati$Peso > s_alta
return(dati)
}Per le future fasi di pulizia dei dati, creiamo la colonna ID che identifica con un numero ciascun record.
df$ID <- 1:nrow(df)Dividiamo i dati in base ai valori di Gestazione, Lunghezza_CL e Cranio_CL e applichiamo la funzione is_outlier, per poi salvare i valori anomali nel DataFrame df_outlier.
Gestazione_list <- split(df, df$Gestazione)
lista_processata <- lapply(Gestazione_list, is_outlier)
df <- do.call(rbind, lista_processata)
df_outlier <- df[df$is_outlier==T, ]Lunghezza_list <- split(df, df$Lunghezza_CL)
lista_processata <- lapply(Lunghezza_list, is_outlier)
df <- do.call(rbind, lista_processata)
df_outlier <- rbind(df_outlier, df[df$is_outlier==T,])Cranio_list <- split(df, df$Cranio_CL)
lista_processata <- lapply(Cranio_list, is_outlier)
df <- do.call(rbind, lista_processata)
df_outlier <- rbind(df_outlier, df[df$is_outlier==T,])Rimuoviamo le colonne create finora tranne ID ed estraiamo i valori unici.
df_outlier$is_outlier <- NULL
df_outlier$Lunghezza_CL <- NULL
df_outlier$Cranio_CL <- NULL
df_outlier <- unique(df_outlier)I tre grafici che seguono mostrano gli outlier (in rosso) all’interno del dataset completo.
plot(df$Gestazione, df$Peso, pch=20, col=rgb(0,0,0,0.2))
points(df_outlier$Gestazione, df_outlier$Peso, col="red", pch=19, cex=1.2)
plot(df$Lunghezza, df$Peso, pch=20, col=rgb(0,0,0,0.2))
points(df_outlier$Lunghezza, df_outlier$Peso, col="red", pch=19, cex=1.2)
plot(df$Cranio, df$Peso, pch=20, col=rgb(0,0,0,0.2))
points(df_outlier$Cranio, df_outlier$Peso, col="red", pch=19, cex=1.2)
La pulizia dei dati non solo ha individuato i punti anomali già segnalati nelle prime visualizzazioni, ma ha riguardato tutti quei record che si trovano intorno alla massa principale di punti. Quantifichiamo la percentuale di questi record rispetto al dataset intero nella cella che segue.
nrow(df_outlier); (nrow(df_outlier)/nrow(df))*100## [1] 150## [1] 6Gli outlier così individuati sono 150, che costituiscono il 6% dell’intero dataset. Vediamo un riepilogo di questi dati.
summary(df_outlier)## Anni.madre N.gravidanze Fumatrici Gestazione Peso
## Min. :15.00 Min. :0.000 No:147 Min. :27.00 Min. : 990
## 1st Qu.:24.00 1st Qu.:0.000 Sì: 3 1st Qu.:37.00 1st Qu.:2265
## Median :29.00 Median :1.000 Median :39.00 Median :2775
## Mean :28.45 Mean :1.033 Mean :37.86 Mean :3085
## 3rd Qu.:32.00 3rd Qu.:1.000 3rd Qu.:40.00 3rd Qu.:4068
## Max. :42.00 Max. :7.000 Max. :42.00 Max. :4930
## Lunghezza Cranio Tipo.parto Ospedale Sesso ID
## Min. :315.0 Min. :267.0 Ces: 40 osp1:49 F:80 Min. : 29.0
## 1st Qu.:456.2 1st Qu.:315.0 Nat:110 osp2:59 M:70 1st Qu.: 727.8
## Median :480.0 Median :332.5 osp3:42 Median :1374.5
## Mean :480.2 Mean :332.6 Mean :1316.5
## 3rd Qu.:515.0 3rd Qu.:351.0 3rd Qu.:1868.8
## Max. :560.0 Max. :390.0 Max. :2478.0Questo sottoinsieme è bilanciato e sensato dal punto di vista delle variabili di controllo e, per il criterio con cui è stato estratto dal dataset, mostra una variabilità nelle colonne quantitative. La nostra soluzione a beneficio del futuro modello è di eliminare questi dati, dal momento che l’informazione necessaria è stata comunque mantenuta nel dataset principale.
df$is_outlier <- NULL
df$Cranio_CL <- NULL
df$Lunghezza_CL <- NULL
df <- df[!(df$ID %in% df_outlier$ID), ]
df$ID <- NULLOra che abbiamo un dataset quanto più pulito e informativo possibile, possiamo procedere con le prossime fasi.
Analisi esplorativa del dataset
Questa sezione, dedicata a un’esplorazione più approfondita del dataset, sarà divisa in tre sezioni.
Nella prima ci occuperemo della variabile target del modello, nella seconda esploreremo la relazione delle variabili quantitative con il target, mentre nella terza verificheremo la relazione tra le variabili qualitative e il target.
La variabile target
Iniziamo questa sezione dedicata alla colonna Peso calcolando la media e la deviazione standard del target.
mu <- mean(df$Peso, na.rm = TRUE)
sigma <- sd(df$Peso, na.rm = TRUE)La cella che segue confronta in un grafico la distribuzione del peso dei neonati con una distribuzione normale avente stessa media e stessa deviazione standard.
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.4.3ggplot(df, aes(x = Peso)) +
geom_density(fill = "grey", alpha = 0.3, color = "black", linewidth = 0.8) +
stat_function(fun = dnorm,
args = list(mean = mu, sd = sigma),
color = "red",
linewidth = 1) +
theme_minimal() +
labs(title = "Confronto tra Densità del Peso e Distribuzione Normale",
x = "Peso (grammi)",
y = "Densità") +
annotate("text", x = mu + 2*sigma, y = 0.0001, label = "Normale Teorica", color = "red", hjust = 0)
Osserviamo una lieve asimmetria verso i valori più alti di peso e una campana più alta rispetto alla distribuzione normale. Inoltre, la distribuzione del peso presenta una coda a sinistra del grafico rappresentata dai neonati più piccoli.
Verifichiamo l’asimmetria e la curtosi della distribuzione nella cella che segue.
library(moments)
skewness(df$Peso); kurtosis(df$Peso)-3## [1] -0.7183713## [1] 2.506875Come anticipato, la distribuzione è asimmetrica negativa e leptocurtica rispetto alla normale. Con questi valori, sospettiamo che la distribuzione della popolazione del peso non sia una normale come richiederebbero le assunzioni del modello che utilizzeremo.
Verifichiamolo con un test di Shapiro-Wilk, che ha come ipotesi nulla che la distribuzione del peso sia normale.
shapiro.test(df$Peso)##
## Shapiro-Wilk normality test
##
## data: df$Peso
## W = 0.96676, p-value < 2.2e-16Il p-value del test è pressoché nullo, il che ci porta a rifiutare l’ipotesi di normalità. Ciononostante, valuteremo comunque un modello statistico di regressione multipla prima di optare per un’eventuale trasformazione della colonna target.
Il target e le variabili quantitative
Questa sottosezione è dedicata alle variabili quantitative e alla loro relazione con il target.
Nella cella che segue, raccogliamo in un elenco le variabili di nostro interesse.
cols_num <- c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")Costruiamo la matrice di correlazione tra le variabili quantitative e mostriamola con una mappa termica.
library(reshape2)
corr_matrix <- cor(df[, cols_num], use = "complete.obs")
melted_corr <- melt(corr_matrix)
ggplot(data = melted_corr, aes(x = Var1, y = Var2, fill = value)) +
geom_tile(color = "white") +
scale_fill_gradient2(low = "red", high = "blue", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Correlazione") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
coord_fixed() +
geom_text(aes(label = round(value, 2)), color = "black", size = 3) +
labs(title = "Matrice di Correlazione Variabili Cliniche", x = "", y = "")
Concentrandoci su Peso, notiamo una scarsa correlazione con l’età della madre e le gravidanze pregresse, mentre una correlazione più forte le abbiamo con le misure antropometriche e la durata della gestazione. Inoltre, osserviamo una correlazione tra lunghezza e cranio che ci fa sospettare che le due colonne indichino in realtà la stessa caratteristica, ovvero la taglia del bambino. Valuteremo se e come eliminare una delle due colonne nella sezione dedicata alla costruzione del modello.
Il target e le variabili qualitative
Questa sottosezione è dedicata alle variabili di controllo e alla loro relazione con il target.
Iniziamo dal fumo materno, come anche richiesto dall’azienda. La cella che segue mostra un boxplot della variabile target in base ai valori della colonna Fumatrici.
boxplot(df$Peso~df$Fumatrici)
Al di là della presenza di outlier nel boxplot delle non fumatrici, osserviamo un range di valori più ristretto tra le fumatrici, oltre che un range interquartile più spostato verso il basso. Verifichiamo se il fumo materno incide sul peso tramite un T test, che assume come ipotesi nulla che le medie di peso dei due sottocampioni siano significativamente uguali.
t.test(df$Peso~df$Fumatrici)##
## Welch Two Sample t-test
##
## data: df$Peso by df$Fumatrici
## t = 1.7693, df = 110.67, p-value = 0.07959
## alternative hypothesis: true difference in means between group No and group Sì is not equal to 0
## 95 percent confidence interval:
## -9.466406 167.249641
## sample estimates:
## mean in group No mean in group Sì
## 3300.179 3221.287Il valore del p-value ci porta a rifiutare l’ipotesi nulla, ma non il senso comune. Il sospetto è che il fumo influenzi altre variabili quantitative che, a loro volta, influenzano il peso, senza tuttavia che queste due variabili siano correlate tra loro.
Per verificarlo, creiamo una tabella che raccoglie i p-value del T test eseguito sulle variabili quantitative influenzate dal fumo e ordiniamo questi valori in modo crescente.
p_values <- data.frame(cols = c(), p_vals = c())
for(i in 1:length(cols_num)){
test_num <- t.test(df[[cols_num[i]]]~df$Fumatrici)
p_values <- rbind(p_values, list(cols=cols_num[i], p_vals=test_num$p.value))
}
p_values[order(p_values$p_vals),]## cols p_vals
## 2 N.gravidanze 0.02239669
## 5 Lunghezza 0.03949287
## 4 Peso 0.07959402
## 3 Gestazione 0.08330045
## 6 Cranio 0.39581625
## 1 Anni.madre 0.82515572Se escludiamo le colonne con p-value superiore a 0.05 (e che di conseguenza non rifiutano l’ipotesi nulla), notiamo che il fumo materno potrebbe essere correlato a Lunghezza e a N.gravidanze. Sull’influenza del fumo nella sezione dedicata alla costruzione del modello.
Proseguiamo con le variabili qualitative tracciando i boxplot del peso in base al tipo di parto eseguito.
boxplot(df$Peso~df$Tipo.parto)
I due grafici sono pressoché identici, il che ci porta a concludere che non ci sia differenza significativa nel tipo di parto.
In modo analogo, tracciamo il boxplot del peso in base all’ospedale del circuito in cui è avvenuta la nascita.
boxplot(df$Peso~df$Ospedale)
Anche qui, non rileviamo alcuna differenza significativa tra i boxplot, il che ci porta a supporre che l’ospedale non abbia alcun ruolo nel peso del neonato.
Concludiamo l’analisi delle variabili qualitative con il sesso del bambino.
boxplot(df$Peso~df$Sesso)
Notiamo che il boxplot dei maschi è spostato verso i valori più elevati e che la sua mediana coincide con il terzo quartile del boxplot delle femmine. Questo suggerisce quindi una differenza significativa tra le medie dei pesi che possiamo verificare tramite il T test.
t.test(df$Peso~df$Sesso)##
## Welch Two Sample t-test
##
## data: df$Peso by df$Sesso
## t = -11.919, df = 2345.8, 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:
## -260.9943 -187.2470
## sample estimates:
## mean in group F mean in group M
## 3184.823 3408.944Il valore del p-value ci porta ad accettare come ipotesi alternativa una significativa differenza di peso in media tra maschi e femmine e, quindi, un’influenza del sesso del bambino sul peso.
Riassumendo, delle variabili di controllo esaminate è più probabile che manterremo Sesso, mentre di Fumatrici si potrebbe valutare un’interazione con altre variabili.
I test di ipotesi
Questa sezione è dedicata all’esecuzione di test statistici per la verifica delle seguenti ipotesi richieste dall’azienda committente:
- l’associazione tra ospedale e tipo di parto;
- la bontà della stima della media della popolazione di lunghezza e peso dei neonati da parte del campione;
- l’associazione tra sesso e misure di lunghezza, cranio e peso.
A ciascuna di queste ipotesi sarà dedicata una sottosezione di questo documento.
Associazione tra ospedale e tipo di parto
Questa sottosezione verifica l’ipotesi per cui in certi ospedali si preferisca praticare un determinato tipo di parto. Questa ipotesi equivale a quella per cui le variabili Ospedale e Tipo.parto siano dipendenti tra loro.
Il primo passo è stato ottenere una tabella di contingenza dalle colonne Ospedale e Tipo.parto. I risultati sono riepilogati nel balloon plot della cella seguente.
osservazioni <- table(df$Ospedale,df$Tipo.parto)
ggpubr::ggballoonplot(data = as.data.frame(osservazioni))
Osserviamo che la prevalenza del parto naturale rispetto al cesareo si mantiene per tutti e tre gli ospedali e che, anzi, nell’ospedale 3 i parti cesarei sono addirittura meno frequenti rispetto agli altri due ospedali. Graficamente, quindi, possiamo supporre che non esiste associazione tra le due variabili in esame.
Verifichiamo numericamente questa supposizione costruendo la tabella delle osservazioni attese.
attese <- osservazioni
for(i in 1:nrow(osservazioni)){
for(j in 1:ncol(osservazioni)){
attese[i,j] <- table(df$Ospedale)[i]*table(df$Tipo.parto)[j]/nrow(df)
}
}La cella che segue calcola il valore della statistica \(X^2\) di Pearson sulla tabella di contingenza e sulle frequenze attese.
X2 <- sum((osservazioni-attese)^2/attese)Il grafico seguente mostra dove si posiziona questo valore (in verde) nella distribuzione \(\chi^2\) con i gradi di libertà ottenuti dal numero di righe e di colonne della tabella di contingenza. La retta tratteggiata di rosso delimita la zona del grafico oltre la quale rifiutiamo l’ipotesi nulla di indipendenza.
plot(density(rchisq(nrow(df), df=(ncol(osservazioni)-1)*(nrow(osservazioni)-1))),main = "Relazione tra ospedale e tipo di parto")
abline(v=qchisq(0.95, df=(ncol(osservazioni)-1)*(nrow(osservazioni)-1)),col="red",lty=2)
abline(v=X2, col="green")
Osserviamo che il valore ottenuto ricade nella zona di accettazione, pertanto possiamo accettare l’ipotesi nulla per cui Ospedale e Tipo.parto sono due variabili indipendenti. I risultati ottenuti sono confermati dal test \(\chi^2\) nella cella seguente, dove il p-value supera la significatività al \(5\%\).
chisq.test(osservazioni)##
## Pearson's Chi-squared test
##
## data: osservazioni
## X-squared = 1.7949, df = 2, p-value = 0.4076Media di lunghezza e peso
In questa sottosezione verificheremo se la media di lunghezza e peso del campione siano una buona stima delle medie della popolazione. Queste medie, impostate nella cella che segue, sono state ottenute dal WHO Child Growth Standards (2006) e dal Intergrowth-21st Project (University of Oxford).
LUNGHEZZA.POP <- 495
PESO.POP <- 3300Dal momento che non conosciamo la varianza delle due variabili nella popolazione, opteremo per un T test a due code con una significatività del \(5\%\).
La cella seguente effettua il test su Lunghezza e ne memorizza i risultati.
t.lunghezza <- t.test(df$Lunghezza, mu = LUNGHEZZA.POP, conf.level = 0.95, alternative = "two.sided")Il grafico che segue mostra la posizione dei risultati ottenuti nel test (in verde) all’interno della distribuzione T di Student.
plot(density(rt(nrow(df), df = nrow(df)-1)), main="Differenza tra media della popolazione e media campionaria")
abline(v=c(qt(0.025,df=nrow(df)-1), qt(0.975,df=nrow(df)-1)),col="red",lty=2)
abline(v=qt(t.lunghezza$p.value, df=nrow(df)-1), col="green")
Il valore ottenuto dal test ricade nella zona di accettazione delimitata dalle due rette rosse, pertanto non rifiutiamo l’ipotesi nulla per cui la media del campione è una buona stima della media della popolazione.
In modo analogo, eseguiamo lo stesso tipo di test sulla colonna Peso.
t.peso <- t.test(df$Peso, mu = PESO.POP, conf.level = 0.95, alternative = "two.sided")Il grafico nella cella seguente mostra graficamente i risultati del test.
plot(density(rt(nrow(df), df = nrow(df)-1)), main = "Differenza tra media della popolazione e media campionaria")
abline(v=c(qt(0.025,df=nrow(df)-1), qt(0.975,df=nrow(df)-1)),col="red",lty=2)
abline(v=qt(t.peso$p.value, df=nrow(df)-1), col="green")
Anche qui, il valore ottenuto dal test ricade nella zona di accettazione, perciò non possiamo rifiutare l’ipotesi nulla per cui la media stimata dal campione è significativamente uguale a quella della popolazione.
Sesso e misure antropometriche
In questa sottosezione, testiamo l’associazione tra sesso e misure antropometriche dei neonati.
Riportiamo nella cella che segue le prime cinque righe delle colonne Sesso, Peso, Lunghezza e Cranio.
head(df[c("Sesso","Peso","Lunghezza","Cranio")],5)## Sesso Peso Lunghezza Cranio
## (230,240] F 930 355 235
## (240,250] F 930 345 245
## (250,260].(300,320].28.928 F 830 310 254
## (250,260].(320,340].25 F 900 325 253
## (260,270].(300,320].27.2437 M 980 320 265Ripetiamo il test T eseguito sulla variabile Peso nella sezione precedente. Ricordiamo che il p-value di questo test ci aveva spinto a rifiutare l’ipotesi nulla a favore dell’alternativa, per cui la media del peso era influenzata dal sesso.
t.test(df$Peso~df$Sesso)##
## Welch Two Sample t-test
##
## data: df$Peso by df$Sesso
## t = -11.919, df = 2345.8, 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:
## -260.9943 -187.2470
## sample estimates:
## mean in group F mean in group M
## 3184.823 3408.944In modo analogo, eseguiamo lo stesso tipo di test su Lunghezza e Cranio. I risultati sono riportati nelle due celle seguenti.
t.test(df$Lunghezza~df$Sesso)##
## Welch Two Sample t-test
##
## data: df$Lunghezza by df$Sesso
## t = -9.179, df = 2329.4, 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:
## -10.927112 -7.080099
## sample estimates:
## mean in group F mean in group M
## 491.1165 500.1201t.test(df$Cranio~df$Sesso)##
## Welch Two Sample t-test
##
## data: df$Cranio by df$Sesso
## t = -6.6581, df = 2344.7, p-value = 3.442e-11
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -5.474464 -2.983421
## sample estimates:
## mean in group F mean in group M
## 338.3895 342.6184I due test, pur essendo eseguiti su colonne diverse, ci portano a concludere che le misure antropometriche dei neonati sono condizionate dal loro sesso e pertanto terremo in considerazione l’interazione tra queste variabili nel modello.
Costruzione del modello
In questa sezione ci occuperemo della costruzione del modello seguendo una procedura di tipo stepwise backward.
Iniziamo con un modello baseline creato in modo automatico nella cella che segue e visualizziamone i risultati principali.
model_base <- MASS::stepAIC(lm(Peso~.,data=df,direction="both",k=2))## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
## extra arguments 'direction', 'k' will be disregarded## Start: AIC=25726.77
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## - Anni.madre 1 23932 132296117 25725
## <none> 132272185 25727
## - Fumatrici 1 127286 132399471 25727
## - Tipo.parto 1 244850 132517035 25729
## - Ospedale 2 419589 132691774 25730
## - N.gravidanze 1 436395 132708581 25732
## - Sesso 1 3305421 135577606 25783
## - Gestazione 1 3854058 136126243 25792
## - Cranio 1 31542477 163814662 26227
## - Lunghezza 1 70691331 202963516 26731## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
## extra arguments 'direction', 'k' will be disregarded##
## Step: AIC=25725.2
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Ospedale + Sesso
##
## Df Sum of Sq RSS AIC
## <none> 132296117 25725
## - Fumatrici 1 128284 132424401 25726
## - Tipo.parto 1 246422 132542539 25728
## - Ospedale 2 420749 132716866 25729
## - N.gravidanze 1 595674 132891791 25734
## - Sesso 1 3314745 135610863 25781
## - Gestazione 1 3832204 136128322 25790
## - Cranio 1 31721415 164017532 26228
## - Lunghezza 1 70672846 202968964 26729summary(model_base)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso, data = df, direction = "both",
## k = 2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -694.07 -168.15 -12.77 155.52 875.51
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6340.7298 129.7062 -48.885 < 2e-16 ***
## N.gravidanze 12.6507 3.8974 3.246 0.00119 **
## FumatriciSì -36.5991 24.2968 -1.506 0.13212
## Gestazione 29.6410 3.6003 8.233 3.00e-16 ***
## Lunghezza 10.3783 0.2935 35.356 < 2e-16 ***
## Cranio 9.6135 0.4059 23.687 < 2e-16 ***
## Tipo.partoNat 22.5496 10.8010 2.088 0.03693 *
## Ospedaleosp2 -19.3756 12.0646 -1.606 0.10841
## Ospedaleosp3 13.0990 12.0467 1.087 0.27699
## SessoM 76.5557 9.9981 7.657 2.76e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 237.8 on 2340 degrees of freedom
## Multiple R-squared: 0.7442, Adjusted R-squared: 0.7433
## F-statistic: 756.6 on 9 and 2340 DF, p-value: < 2.2e-16Questo modello non può essere preso in considerazione come candidato. L’intercetta è impossibile e le variabili generate da Ospedale hanno una significatività diversa l’una dall’altra, per non parlare delle ipotesi formulate in precedenza contraddette dai regressori individuati dal modello.
Procediamo quindi manualmente e iniziamo definendo la funzione collect_results, che dato un modello in ingresso raccoglie i dati di BIC, percentuale di varianza spiegata e radice dell’errore quadratico medio.
collect_results <- function(model_name,model){
return(list(
name=model_name,
bic = BIC(model),
r_adj = summary(model)$adj.r.squared,
rmse = sqrt(mean(model$residuals^2))
))
}Iniziamo dal primo modello, che contiene tutte le variabili. Il dataframe results salverà i risultati ottenuti dai modelli per le nostre considerazioni finali.
model_1 <- lm(Peso~., data=df)
results <- data.frame(collect_results("model_1",model_1))
summary(model_1)##
## Call:
## lm(formula = Peso ~ ., data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -695.52 -167.80 -12.15 156.15 873.62
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6363.9873 134.5585 -47.295 < 2e-16 ***
## Anni.madre 0.6705 1.0307 0.651 0.51541
## N.gravidanze 11.6464 4.1925 2.778 0.00551 **
## FumatriciSì -36.4578 24.3008 -1.500 0.13368
## Gestazione 29.8809 3.6195 8.255 2.50e-16 ***
## Lunghezza 10.3802 0.2936 35.356 < 2e-16 ***
## Cranio 9.5996 0.4065 23.617 < 2e-16 ***
## Tipo.partoNat 22.4787 10.8029 2.081 0.03756 *
## Ospedaleosp2 -19.4991 12.0676 -1.616 0.10627
## Ospedaleosp3 12.9107 12.0516 1.071 0.28415
## SessoM 76.4568 10.0005 7.645 3.02e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 237.8 on 2339 degrees of freedom
## Multiple R-squared: 0.7443, Adjusted R-squared: 0.7432
## F-statistic: 680.8 on 10 and 2339 DF, p-value: < 2.2e-16Come atteso anche dalle considerazioni iniziali, le variabili meno significative per questo modello sono l’età della madre, il fumo e l’ospedale dove è avvenuto il parto, mentre il numero di gravidanze pregresse e il tipo di parto hanno una significatività superiore a quella prevista.
Creiamo un secondo modello, ottenuto dal primo togliendo Anni.madre.
model_2 <- update(model_1, ~.-Anni.madre)
results <- rbind(results, collect_results("model_2",model_2))
summary(model_2)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Ospedale + Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -694.07 -168.15 -12.77 155.52 875.51
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6340.7298 129.7062 -48.885 < 2e-16 ***
## N.gravidanze 12.6507 3.8974 3.246 0.00119 **
## FumatriciSì -36.5991 24.2968 -1.506 0.13212
## Gestazione 29.6410 3.6003 8.233 3.00e-16 ***
## Lunghezza 10.3783 0.2935 35.356 < 2e-16 ***
## Cranio 9.6135 0.4059 23.687 < 2e-16 ***
## Tipo.partoNat 22.5496 10.8010 2.088 0.03693 *
## Ospedaleosp2 -19.3756 12.0646 -1.606 0.10841
## Ospedaleosp3 13.0990 12.0467 1.087 0.27699
## SessoM 76.5557 9.9981 7.657 2.76e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 237.8 on 2340 degrees of freedom
## Multiple R-squared: 0.7442, Adjusted R-squared: 0.7433
## F-statistic: 756.6 on 9 and 2340 DF, p-value: < 2.2e-16In questo secondo modello è lievemente aumentata la percentuale di varianza spiegata, ma mantiene la scarsa significatività di Fumatrici e di Ospedale. Quest’ultima variabile viene rimossa dal terzo modello.
model_3 <- update(model_2, ~.-Ospedale)
results <- rbind(results, collect_results("model_3",model_3))
summary(model_3)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza +
## Cranio + Tipo.parto + Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -691.8 -169.7 -11.4 156.8 876.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6342.3275 129.7323 -48.888 < 2e-16 ***
## N.gravidanze 13.0113 3.8996 3.337 0.000862 ***
## FumatriciSì -37.3597 24.3212 -1.536 0.124649
## Gestazione 29.8648 3.6035 8.288 < 2e-16 ***
## Lunghezza 10.3469 0.2936 35.237 < 2e-16 ***
## Cranio 9.6290 0.4063 23.700 < 2e-16 ***
## Tipo.partoNat 23.3233 10.8094 2.158 0.031054 *
## SessoM 76.9900 10.0084 7.693 2.11e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 238.1 on 2342 degrees of freedom
## Multiple R-squared: 0.7434, Adjusted R-squared: 0.7427
## F-statistic: 969.5 on 7 and 2342 DF, p-value: < 2.2e-16A fronte della diminuzione della percentuale di varianza spiegata, abbiamo ottenuto un aumento della significatività della variabile N.gravidanze.
Nel prossimo modello, rimuoviamo la variabile Fumatrici.
model_4 <- update(model_3, ~.-Fumatrici)
results <- rbind(results, collect_results("model_4",model_4))
summary(model_4)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -690.27 -169.93 -10.53 157.93 879.58
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6342.8834 129.7694 -48.878 < 2e-16 ***
## N.gravidanze 12.7180 3.8961 3.264 0.00111 **
## Gestazione 29.4824 3.5959 8.199 3.95e-16 ***
## Lunghezza 10.3733 0.2932 35.377 < 2e-16 ***
## Cranio 9.6333 0.4064 23.704 < 2e-16 ***
## Tipo.partoNat 22.9745 10.8102 2.125 0.03367 *
## SessoM 76.6880 10.0093 7.662 2.67e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 238.1 on 2343 degrees of freedom
## Multiple R-squared: 0.7432, Adjusted R-squared: 0.7425
## F-statistic: 1130 on 6 and 2343 DF, p-value: < 2.2e-16In questo modello abbiamo ottenuto almeno una stella in tutte le variabili rimaste, pur con un’ulteriore diminuzione della variabilità spiegata.
Aggiorniamo questo modello aggiungendo l’iterazione del fumo materno con la lunghezza del bambino e le gravidanze pregresse, ossia le variabili più correlate a Fumatrici nel dataset.
model_5 <- update(model_4, ~.+ Lunghezza*Fumatrici + N.gravidanze*Fumatrici)
results <- rbind(results, collect_results("model_5",model_5))
summary(model_5)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso + Fumatrici + Lunghezza:Fumatrici + N.gravidanze:Fumatrici,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -691.75 -170.24 -11.32 157.53 873.66
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6322.4557 130.7150 -48.368 < 2e-16 ***
## N.gravidanze 13.2794 3.9804 3.336 0.000862 ***
## Gestazione 29.9602 3.6049 8.311 < 2e-16 ***
## Lunghezza 10.2972 0.2960 34.785 < 2e-16 ***
## Cranio 9.6321 0.4063 23.706 < 2e-16 ***
## Tipo.partoNat 23.3844 10.8101 2.163 0.030627 *
## SessoM 76.3527 10.0231 7.618 3.72e-14 ***
## FumatriciSì -787.8833 589.8946 -1.336 0.181799
## Lunghezza:FumatriciSì 1.5400 1.1925 1.291 0.196687
## N.gravidanze:FumatriciSì -5.2490 19.7496 -0.266 0.790435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 238.1 on 2340 degrees of freedom
## Multiple R-squared: 0.7436, Adjusted R-squared: 0.7426
## F-statistic: 754.2 on 9 and 2340 DF, p-value: < 2.2e-16Pur essendo aumentata la variabilità spiegata e la significatività di N.gravidanze, notiamo che le nuove variabili inserite non sono per nulla significative per il modello.
Esaminiamo ulteriormente questo modello calcolando il VIF delle variabili coinvolte.
car::vif(model_5)## there are higher-order terms (interactions) in this model
## consider setting type = 'predictor'; see ?vif## N.gravidanze Gestazione Lunghezza
## 1.068715 1.650505 2.126198
## Cranio Tipo.parto Sesso
## 1.651854 1.003331 1.041461
## Fumatrici Lunghezza:Fumatrici N.gravidanze:Fumatrici
## 593.501629 586.827251 2.106994Le nuove variabili introdotte presentano un problema di estrema multicollinearità che ci portano a scartare sin da subito questo modello e a riprendere la costruzione dei nostri candidati dal modello 4.
Concludiamo derivando due ulteriori modelli dal modello 4:
il modello 6 aggiorna il modello rimuovendo la variabile Lunghezza;
il modello 7 aggiorna il modello rimuovendo la variabile Cranio.
Questo perché, come accennato nell’analisi esplorativa, le due variabili hanno una moderata correlazione positiva e indicano in modo analogo la taglia del bambino. Otteniamo quindi i modelli sopra descritti nelle celle che seguono.
model_6 <- update(model_4, ~.-Lunghezza)
results <- rbind(results, collect_results("model_6",model_6))
summary(model_6)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Cranio + Tipo.parto +
## Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -914.5 -208.5 -7.3 198.6 846.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5743.0161 159.3213 -36.047 <2e-16 ***
## N.gravidanze 3.4034 4.8137 0.707 0.480
## Gestazione 87.2764 3.9669 22.001 <2e-16 ***
## Cranio 16.3413 0.4451 36.712 <2e-16 ***
## Tipo.partoNat 7.3725 13.3756 0.551 0.582
## SessoM 117.5309 12.3123 9.546 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 294.9 on 2344 degrees of freedom
## Multiple R-squared: 0.606, Adjusted R-squared: 0.6052
## F-statistic: 721 on 5 and 2344 DF, p-value: < 2.2e-16model_7 <- update(model_4, ~.-Cranio)
results <- rbind(results, collect_results("model_7",model_7))
summary(model_7)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Tipo.parto +
## Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -804.56 -184.16 -16.77 180.99 1215.64
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5186.7651 133.8724 -38.744 < 2e-16 ***
## N.gravidanze 21.6354 4.3170 5.012 5.80e-07 ***
## Gestazione 42.2843 3.9577 10.684 < 2e-16 ***
## Lunghezza 13.6163 0.2887 47.161 < 2e-16 ***
## Tipo.partoNat 30.4926 12.0291 2.535 0.0113 *
## SessoM 82.2782 11.1396 7.386 2.09e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 265.1 on 2344 degrees of freedom
## Multiple R-squared: 0.6816, Adjusted R-squared: 0.6809
## F-statistic: 1003 on 5 and 2344 DF, p-value: < 2.2e-16Nel modello 6, la rimozione di Lunghezza ha portato a un drastico calo della variabilità spiegata e alla perdita di significatività delle gravidanze pregresse e del tipo di parto. Dall’altra parte, il modello 7 ha consolidato la significatività delle variabili rimaste al netto di una minore perdita di variabilità spiegata.
La selezione e l’analisi del modello
Questa sezione è dedicata alla selezione e all’analisi del modello migliore tra i candidati costruiti in precedenza.
Mostriamo allora la tabella results ottenuta nella sezione precedente.
results## name bic r_adj rmse
## 1 model_1 32466.93 0.7432005 237.2468
## 2 model_2 32459.59 0.7432638 237.2682
## 3 model_3 32451.53 0.7426672 237.6452
## 4 model_4 32446.14 0.7425179 237.7649
## 5 model_5 32465.21 0.7426491 237.5521
## 6 model_6 33444.12 0.6051527 294.4974
## 7 model_7 32943.54 0.6809049 264.7446Notiamo che dopo il quinto modello (quello dell’iterazione del fumo materno), l’errore medio dei candidati è aumentato insieme al BIC, mentre la variabilità spiegata è scesa drasticamente sotto il 70%.
Mostriamo questi valori sui grafici, contrassegnando i modelli in cui raggiungono l’ottimo (il minimo nel BIC e nell’errore medio, il massimo nell’R aggiustato).
indice_ottimo<- which.min(results$bic)
valore_ottimo <- min(results$bic)
plot(results$bic, type = "o", pch = 20, col = "darkgray", xlab = "Modelli", ylab = "BIC", main = "Selezione del Modello tramite BIC")
points(indice_ottimo, valore_ottimo, col = "green2", pch = 20, cex = 2)
indice_ottimo<- which.max(results$r_adj)
valore_ottimo <- max(results$r_adj)
plot(results$r_adj, type = "o", pch = 20, col = "darkgray", xlab = "Modelli", ylab = "Var. Spiegata", main = "Selezione del Modello tramite R Aggiustato")
points(indice_ottimo, valore_ottimo, col = "green2", pch = 20, cex = 2)
indice_ottimo<- which.min(results$rmse)
valore_ottimo <- min(results$rmse)
plot(results$rmse, type = "o", pch = 20, col = "darkgray", xlab = "Modelli", ylab = "RMSE", main = "Selezione del Modello tramite RMSE")
points(indice_ottimo, valore_ottimo, col = "green2", pch = 20, cex = 2)
Tra i tre modelli ottimali selezionati, scegliamo il modello 4 per il principio di parsimonia, dal momento che richiede meno variabili rispetto agli altri due ed è quindi meno complesso.
Impostiamo quindi il quarto modello come modello finale ed esaminiamolo.
final_model <- model_4
summary(final_model)##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -690.27 -169.93 -10.53 157.93 879.58
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6342.8834 129.7694 -48.878 < 2e-16 ***
## N.gravidanze 12.7180 3.8961 3.264 0.00111 **
## Gestazione 29.4824 3.5959 8.199 3.95e-16 ***
## Lunghezza 10.3733 0.2932 35.377 < 2e-16 ***
## Cranio 9.6333 0.4064 23.704 < 2e-16 ***
## Tipo.partoNat 22.9745 10.8102 2.125 0.03367 *
## SessoM 76.6880 10.0093 7.662 2.67e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 238.1 on 2343 degrees of freedom
## Multiple R-squared: 0.7432, Adjusted R-squared: 0.7425
## F-statistic: 1130 on 6 and 2343 DF, p-value: < 2.2e-16Il modello finale è riassumibile come segue:
L’intercetta ha il peso maggiore in termini di valore assoluto. Dal momento che è biologicamente impossibile per un bambino avere lunghezza nulla e peso negativo, ipotizziamo che lo scopo dell’intercetta sia quella di calibrare l’influenza degli altri regressori;
Il predittore continuo più potente è Lunghezza, che aggiunge 10.4 grammi di peso per ogni millimetro di lunghezza del bambino, mentre dal punto di vista discreto abbiamo Gestazione con circa 29.5 grammi in più per ogni settimana di gravidanza. Per quanto riguarda le variabili di controllo, come anche anticipato nelle sezioni precedenti, è il sesso del neonato a prevalere, dove a parità delle altre variabili un maschio pesa in media 76.7 grammi in più delle femmine;
Contrariamente a quanto ipotizzato in precedenza, N.gravidanze e Tipo.parto sono rientrate nel modello con una discreta significatività. Per ogni gravidanza pregressa, infatti, abbiamo 12.7 grammi in più di peso, mentre a parità di condizioni il parto naturale aumenta il peso fino a quasi 23 grammi.
Continuiamo analizzando i residui del modello. La cella che segue traccia la diagnostica dei residui in quattro grafici.
par(mfrow = c(2,2))
plot(final_model)
La diagnostica dei residui ci mette di fronte a un modello valido. Il grafico Q-Q, dove i punti sono allineati lungo la bisettrice, attesta la normalità degli errori, mentre il plot Scale-Location e il Residuals vs Fitted mostrano rispettivamente la minima eteroschedasticità e non-linearità. L’assenza di punti critici nel grafico del Leverage conferma che la strategia di rimozione degli outlier all’inizio del documento è stata efficace per ottenere un modello robusto e non influenzato da singole osservazioni anomale.
Verifichiamo ulteriormente questo sguardo di insieme tracciando la distribuzione dei residui nella cella che segue. La linea tratteggiata rossa mostra la distribuzione normale a parità di media e deviazione standard dei residui.
residui=final_model$residuals
plot(density(residui))
lines(density(rnorm(nrow(df), mean=mean(residui), sd=sd(residui))),col="red", lty=2)
Le due curve sembrano molto simili nella loro forma. Verifichiamo se i residui hanno distribuzione normale tramite il test di Shapiro-Wilk.
shapiro.test(residui)##
## Shapiro-Wilk normality test
##
## data: residui
## W = 0.9966, p-value = 4.07e-05Il test di Shapiro-Wilk ha restituito un valore statistico vicino a 1, indicando una quasi perfetta aderenza dei residui alla distribuzione normale. Il risultato ottenuto dal p-value è da attribuirsi all’elevata numerosità del campione, mentre già l’ispezione visiva del Q-Q Plot aveva confermato che i residui si distribuiscono lungo la diagonale, validando l’assunzione di normalità necessaria per l’inferenza statistica.
Eseguiamo ora il test di Breusch-Pagan per valutare l’omoschedasticità dei residui.
lmtest::bptest(final_model)##
## studentized Breusch-Pagan test
##
## data: final_model
## BP = 14.383, df = 6, p-value = 0.02564Il test di Breusch-Pagan ha rilevato una lieve eteroschedasticità. Tuttavia, data l’elevata numerosità campionaria e l’ispezione visiva del grafico Scale-Location, il modello può essere comunque ritenuto affidabile per l’inferenza e la stima dei coefficienti, dal momento che la stabilità della varianza è mantenuta per la quasi totalità del range osservato.
Valutiamo ora l’autocorrelazione con il test di Durbin-Watson.
lmtest::dwtest(final_model)##
## Durbin-Watson test
##
## data: final_model
## DW = 1.9916, p-value = 0.3985
## alternative hypothesis: true autocorrelation is greater than 0Il p-value ci permette di non rifiutare l’assenza di autocorrelazione tra le osservazioni e, di conseguenza, garantisce che le stime dei coefficienti siano corrette e non distorte da dipendenze interne al campione.
Concludiamo l’analisi del modello identificando i valori anomali. Procediamo con un’ispezione visiva nelle tre celle che seguono, che mostrano rispettivamente i punti di leva, gli outliers e le distanze di Cook.
plot(hatvalues(final_model))
abline(h=2*sum(hatvalues(final_model))/nrow(df), col="red")
plot(rstudent(final_model))
abline(h=c(-2,2), col="red")
plot(cooks.distance(final_model))
abline(h=0.5, col="red")
Tolte le distanze di Cook, dove nessun residuo supera la soglia critica di 0.5, osserviamo una forte quantità di valori anomali. Estraiamoli dal dataset nella cella che segue e calcoliamone la percentuale su tutto il dataset.
leverage_pts <- which(hatvalues(final_model) > 2*sum(hatvalues(final_model))/nrow(df))
outlier_pts <- which(rstudent(final_model)< -2 | rstudent(final_model)>2)
cook_pts <- which(cooks.distance(final_model)>0.5)
strange_vals <- unique(c(leverage_pts, outlier_pts, cook_pts))
df_strange <- df[strange_vals, c("Peso","N.gravidanze","Gestazione","Lunghezza","Cranio","Tipo.parto","Sesso")]
nrow(df_strange); (nrow(df_strange)/nrow(df))*100## [1] 203## [1] 8.638298La cella che segue riassume il comportamento delle variabili di questi record.
summary(df_strange)## Peso N.gravidanze Gestazione Lunghezza
## Min. : 830 Min. : 0.000 Min. :25.00 Min. :310.0
## 1st Qu.:2635 1st Qu.: 0.000 1st Qu.:36.00 1st Qu.:455.0
## Median :3100 Median : 1.000 Median :39.00 Median :490.0
## Mean :3037 Mean : 2.177 Mean :37.52 Mean :475.5
## 3rd Qu.:3750 3rd Qu.: 3.000 3rd Qu.:40.00 3rd Qu.:500.0
## Max. :4470 Max. :12.000 Max. :43.00 Max. :555.0
## Cranio Tipo.parto Sesso
## Min. :235.0 Ces: 64 F:112
## 1st Qu.:316.5 Nat:139 M: 91
## Median :337.0
## Mean :331.8
## 3rd Qu.:350.0
## Max. :390.0Le 203 osservazioni (8,6% rispetto al dataset completo) inserite in df_strange riguardano condizioni cliniche estreme, come nascite a 25 settimane o pesi inferiori al chilogrammo, che il modello fatica a individuare correttamente. Ciononostante, il modello proposto rimane comunque valido per via dei coefficienti altamente significativi e l’R-quadro elevato.
Previsioni su dati reali
L’ultima sezione di questo documento è dedicata a una previsione richiesta dall’azienda committente. È stato chiesto infatti di stimare il peso di una neonata alla trentanovesima settimana la cui madre ha già avuto due gravidanze precedenti.
Con queste specifiche, abbiamo estratto casi simili dal dataset e li abbiamo riassunti nella cella che segue.
df_to_predict <- df[df$Sesso=="F"&df$Gestazione==39&df$N.gravidanze==2,]
summary(df_to_predict)## Anni.madre N.gravidanze Fumatrici Gestazione Peso
## Min. :21.00 Min. :2 No:40 Min. :39 Min. :2630
## 1st Qu.:27.00 1st Qu.:2 Sì: 1 1st Qu.:39 1st Qu.:3090
## Median :32.00 Median :2 Median :39 Median :3220
## Mean :30.59 Mean :2 Mean :39 Mean :3259
## 3rd Qu.:34.00 3rd Qu.:2 3rd Qu.:39 3rd Qu.:3440
## Max. :42.00 Max. :2 Max. :39 Max. :4100
## Lunghezza Cranio Tipo.parto Ospedale Sesso
## Min. :445.0 Min. :305 Ces:14 osp1: 9 F:41
## 1st Qu.:480.0 1st Qu.:330 Nat:27 osp2:15 M: 0
## Median :495.0 Median :335 osp3:17
## Mean :492.8 Mean :336
## 3rd Qu.:500.0 3rd Qu.:345
## Max. :520.0 Max. :362Limitandoci ai regressori individuati nella costruzione del modello, osserviamo quanto segue:
La bambina avrà un peso stimato intorno ai 3260 grammi, che coincide con la media del peso delle bambine con le stesse caratteristiche;
Anche la lunghezza e la circonferenza cranica si aggirano sui 493 mm e i 336 mm rispettivamente;
C’è una prevalenza di parti naturali tra i precedenti, ma non possiamo escludere un parto cesareo.
Non avendo ulteriori informazioni oltre quelle fornite dall’azienda, abbiamo deciso di essere quanto più esaustivi possibile nelle previsioni. In particolare esamineremo i seguenti casi:
La bambina ha lunghezza e circonferenza cranica pari alle medie del campione e nasce con parto naturale;
La bambina nasce con parto cesareo e ha lunghezza e circonferenza cranica pari alle medie del sottocampione di bambine nate con cesareo;
La bambina nasce con parto naturale e ha lunghezza e circonferenza cranica pari alle medie del sottocampione di bambine nate con parto naturale;
Creiamo le nostre casistiche e mostriamole nella cella che segue.
df_to_evaluate = data.frame(
N.gravidanze = c(2, 2, 2),
Gestazione = c(39, 39, 39),
Tipo.parto = c("Nat", "Ces", "Nat"),
Sesso = c("F","F","F")
)
splitted_means <- split(df_to_predict, df_to_predict$Tipo.parto)
df_to_evaluate$Lunghezza <- c(mean(df_to_predict$Lunghezza),mean(splitted_means$Ces$Lunghezza),mean(splitted_means$Nat$Lunghezza))
df_to_evaluate$Cranio <- c(mean(df_to_predict$Cranio),mean(splitted_means$Ces$Cranio),mean(splitted_means$Nat$Cranio))
df_to_evaluate$Tipo.parto <- as.factor(df_to_evaluate$Tipo.parto)
df_to_evaluate$Sesso <- as.factor(df_to_evaluate$Sesso)
df_to_evaluate## N.gravidanze Gestazione Tipo.parto Sesso Lunghezza Cranio
## 1 2 39 Nat F 492.8049 336.0244
## 2 2 39 Ces F 493.5714 333.2143
## 3 2 39 Nat F 492.4074 337.4815Effettuiamo le predizioni con il modello che abbiamo costruito.
df_to_evaluate$Pred_Peso <- predict(final_model, df_to_evaluate)
df_to_evaluate## N.gravidanze Gestazione Tipo.parto Sesso Lunghezza Cranio Pred_Peso
## 1 2 39 Nat F 492.8049 336.0244 3204.391
## 2 2 39 Ces F 493.5714 333.2143 3162.298
## 3 2 39 Nat F 492.4074 337.4815 3214.305Il confronto tra il peso medio osservato e i valori predetti dal modello evidenzia un’elevata affidabilità della stima. In particolare, si osserva come il modello colga correttamente le variazioni dovute alla circonferenza cranica e alla tipologia di parto.
Validiamo ulteriormente i risultati ottenuti tramite il T test. Per farlo, abbiamo assunto che la predizione del peso sia equiparabile alla stima della media del peso della popolazione di bambine con le stesse caratteristiche e che il campione di riferimento equivalga a un campione estratto da questa popolazione.
Eseguiamo i T test per ognuno dei casi sopra esposti nelle celle che seguono.
t.test(df_to_predict$Peso, mu=df_to_evaluate$Pred_Peso[1], alternative = "two.sided")##
## One Sample t-test
##
## data: df_to_predict$Peso
## t = 1.1615, df = 40, p-value = 0.2523
## alternative hypothesis: true mean is not equal to 3204.391
## 95 percent confidence interval:
## 3163.782 3354.755
## sample estimates:
## mean of x
## 3259.268t.test(splitted_means$Ces$Peso, mu=df_to_evaluate$Pred_Peso[2], alternative = "two.sided")##
## One Sample t-test
##
## data: splitted_means$Ces$Peso
## t = 1.6247, df = 13, p-value = 0.1282
## alternative hypothesis: true mean is not equal to 3162.298
## 95 percent confidence interval:
## 3131.263 3381.594
## sample estimates:
## mean of x
## 3256.429t.test(splitted_means$Nat$Peso, mu=df_to_evaluate$Pred_Peso[3], alternative = "two.sided")##
## One Sample t-test
##
## data: splitted_means$Nat$Peso
## t = 0.70453, df = 26, p-value = 0.4874
## alternative hypothesis: true mean is not equal to 3214.305
## 95 percent confidence interval:
## 3125.261 3396.221
## sample estimates:
## mean of x
## 3260.741In tutti gli scenari analizzati, i p-value ci hanno permesso di non rifiutare l’ipotesi nulla e, di conseguenza, di validare la bontà della previsione del modello. Questi risultati ci portano a concludere che il modello possiede un’elevata accuratezza e capacità di generalizzazione sia per la specifica classe di neonati analizzata che per tutte le tipologie di record sottoposte in addestramento.
Conclusioni
In questo documento sono state mostrate le fasi di analisi e realizzazione di un modello statistico per la previsione del peso dei neonati, come richiesto da Neonatal Health Solutions.
Dopo aver importato i dati nella prima sezione, nella seconda abbiamo rimosso i valori anomali che avrebbero inficiato in modo negativo sulle prestazioni del modello. Il risultato è stato un dataset ridotto nelle dimensioni, ma con sufficiente informazione da ottenere un modello robusto.
Nella terza sezione abbiamo esplorato i dati con l’obiettivo di formulare le prime ipotesi sul futuro modello. Abbiamo così osservato che il target non rispettava l’assunzione di normalità e abbiamo individuato dei possibili regressori forti per il modello come il sesso e le misure antropometriche.
Nella quarta sezione abbiamo effettuato dei test statistici richiesti dall’azienda. Abbiamo in particolare smentito la dipendenza tra ospedale e tipo di parto, abbiamo confrontato le medie di lunghezza e peso con quelle stimate dal WHO e abbiamo visto una forte correlazione tra sesso e misure antropometriche.
Dopo aver costruito una serie di modelli nella quinta sezione, nella sesta sezione abbiamo scelto e analizzato il modello migliore. Malgrado i residui non rispettassero numericamente le assunzioni di normalità e di omoschedasticità, la diagnostica grafica unita a un’ispezione dei valori anomali ci ha permesso di stabilire che il modello è solido e che le mancate assunzioni sono perlopiù dovute alla numerosità del campione e a record estremi come nati prematuri o bambini molto pesanti.
Infine, abbiamo testato il modello effettuando una predizione richiesta dall’azienda. Esaminando tre possibili casistiche, abbiamo confermato la validità del modello tramite T test che equiparano la predizione a una stima della media del peso su una popolazione di bambine con caratteristiche simili.