library(ggplot2)
library(dplyr)
##
## Caricamento pacchetto: 'dplyr'
## I seguenti oggetti sono mascherati da 'package:stats':
##
## filter, lag
## I seguenti oggetti sono mascherati da 'package:base':
##
## intersect, setdiff, setequal, union
library(car)
## Warning: il pacchetto 'car' è stato creato con R versione 4.4.3
## Caricamento del pacchetto richiesto: carData
## Warning: il pacchetto 'carData' è stato creato con R versione 4.4.3
##
## Caricamento pacchetto: 'car'
## Il seguente oggetto è mascherato da 'package:dplyr':
##
## recode
knitr::opts_chunk$set(echo = TRUE)
setwd("C:/Users/g.morgese/Desktop/DATA SCIENCE & AI MASTER/3. statistica inferenziale/PROGETTO")
getwd()
## [1] "C:/Users/g.morgese/Desktop/DATA SCIENCE & AI MASTER/3. statistica inferenziale/PROGETTO"
dir()
## [1] "MORGESE-GIOVANNI---STATISTICA-INFERENZIALE.html"
## [2] "MORGESE-GIOVANNI---STATISTICA-INFERENZIALE.Rmd"
## [3] "MORGESE GIOVANNI - STATISTICA INFERENZIALE.docx"
## [4] "MORGESE GIOVANNI - STATISTICA INFERENZIALE.Rmd"
## [5] "neonati.csv"
dati = read.csv("neonati.csv",
stringsAsFactors = TRUE)
attach(dati)
n = nrow(dati)
Neonatal Health Solutions
Creare un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali. Il progetto mira a migliorare la gestione delle gravidanze ad alto rischio, ottimizzare le risorse ospedaliere e garantire migliori risultati per la salute neonatale. Il progetto si inserisce all’interno di un contesto di crescente attenzione verso la prevenzione delle complicazioni neonatali. La possibilità di prevedere il peso alla nascita dei neonati rappresenta un’opportunità fondamentale per migliorare la pianificazione clinica e ridurre i rischi associati a nascite problematiche, come parti prematuri o neonati con basso peso.
Di seguito, i principali benefici che questo progetto porterà all’azienda e al settore sanitario:
Miglioramento delle previsioni cliniche: Il peso del neonato è un indicatore chiave della sua salute. Avere un modello predittivo accurato consente al personale medico di intervenire tempestivamente in caso di anomalie, riducendo le complicazioni perinatali come le difficoltà respiratorie o l’ipoglicemia.
Ottimizzazione delle risorse ospedaliere: Sapere in anticipo quali neonati potrebbero avere bisogno di cure intensive aiuta a organizzare le risorse umane e tecnologiche degli ospedali in modo efficiente. Questo si traduce in una riduzione dei costi operativi e una migliore pianificazione dell’utilizzo delle unità di terapia intensiva neonatale (TIN).
Prevenzione e identificazione dei fattori di rischio: Il modello potrà evidenziare i fattori che maggiormente influenzano negativamente il peso del neonato (come il fumo materno, gravidanze multiple o età avanzata della madre). Queste informazioni sono preziose per la prevenzione e la gestione personalizzata delle gravidanze, permettendo interventi proattivi in caso di rischio elevato.
Valutazione delle pratiche ospedaliere: Attraverso un’analisi comparativa tra i tre ospedali coinvolti, l’azienda potrà identificare eventuali differenze nei risultati clinici, come una maggiore incidenza di parti cesarei in una determinata struttura. Ciò consente di monitorare la qualità delle pratiche e armonizzare i protocolli tra i diversi centri ospedalieri, migliorando la coerenza delle cure.
Supporto alla pianificazione strategica: L’analisi dei dati e le previsioni possono essere utilizzate per prendere decisioni informate non solo a livello clinico ma anche strategico. L’azienda potrà sfruttare queste informazioni per implementare nuove politiche di salute pubblica, garantendo un impatto positivo sui tassi di mortalità e morbilità neonatale.
Per costruire il modello predittivo, abbiamo raccolto dati su 2500 neonati provenienti da tre ospedali. Le variabili raccolte includono:
• Età della madre Misura dell’età in anni.
• Numero di gravidanze Quante gravidanze ha avuto la madre.
• Fumo materno Un indicatore binario (0=non fumatrice, 1=fumatrice).
• Durata della gravidanza Numero di settimane di gestazione.
• Peso del neonato Peso alla nascita in grammi.
• Lunghezza e diametro del cranio Lunghezza del neonato e diametro craniale, misurabili anche durante la gravidanza tramite ecografie.
• Tipo di parto Naturale o cesareo.
• Ospedale di nascita Ospedale 1, 2 o 3.
• Sesso del neonato Maschio (M) o femmina (F).
L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature.
Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie. Inoltre si saggeranno le seguenti ipotesi con i test adatti:
Verrà sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo, potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato ed eventuali interazioni. Ad esempio, ci aspettiamo che una maggiore durata della gestazione aumenterebbe in media il peso del neonato.
Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC), selezioneremo il modello più parsimonioso, eliminando le variabili non significative. Verranno considerati anche modelli con interazioni tra le variabili e possibili effetti non lineari.
Una volta ottenuto il modello finale, valuteremo la sua capacità predittiva utilizzando metriche come R² e il Root Mean Squared Error (RMSE). Un’attenzione particolare sarà rivolta all’analisi dei residui e alla presenza di valori influenti, che potrebbero distorcere le previsioni, indagando su di essi.
Una volta validato il modello, lo useremo per fare previsioni pratiche. Ad esempio, potremo stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.
Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili. Ad esempio, potremmo visualizzare l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.
Il progetto di previsione del peso neonatale è un’iniziativa fondamentale per Neonatal Health Solutions. Attraverso l’uso di dati clinici dettagliati e strumenti di analisi statistica avanzati, possiamo contribuire a migliorare la qualità della cura prenatale, ridurre i rischi per i neonati e ottimizzare l’efficienza delle risorse ospedaliere. Questo progetto rappresenta un punto di svolta per l’azienda, consentendo non solo un miglioramento della pratica clinica ma anche l’implementazione di politiche sanitarie più informate e proattive.
Il dataset fornito offre un elenco di nove variabili esplicative e una di risposta, che rappresenta proprio il peso del neonato. Procediamo al caricamento trasformando le stringhe in fattori in modo da poterle trattare correttamente nell’analisi e nella definizione del modello. Considerando la finalità dello studio, il campione e la biologia che governa la crescita del neonato, è lecito aspettarsi distribuzioni leptocurtiche per variabili come Lunghezza e Cranio, la cui crescita è direttamente proporzionale all’età del feto. Verifichiamo la struttura dei dati ed esploriamoli nel dettaglio per ottenere una prima statistica riassuntiva.
head(dati, 5)
## Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto
## 1 26 0 0 42 3380 490 325 Nat
## 2 21 2 0 39 3150 490 345 Nat
## 3 34 3 0 38 3640 500 375 Nat
## 4 28 1 0 41 3690 515 365 Nat
## 5 20 0 0 38 3700 480 335 Nat
## Ospedale Sesso
## 1 osp3 M
## 2 osp1 F
## 3 osp2 M
## 4 osp2 M
## 5 osp3 F
Analizzando il dataset possiamo classificare le variabili esplicative in due categorie:
Variabili quantitative:
Variabili qualitative:
La variabile risposta, oggetto del test, è la variabile Peso dalla quale iniziamo la nostra analisi.
Esploriamo i dati con un riassunto delle informazioni e con un’analisi visiva della distribuzione dei dati.
summary(Peso)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 830 2990 3300 3284 3620 4930
ggplot(dati, aes(x = Peso)) +
geom_histogram(bins = 30, fill = 3, col = 1) +
geom_density(aes(y = after_stat(density)), col = 2, linewidth = 1) +
labs(
title = "Distribuzione del Peso alla Nascita",
x = "Peso (grammi)",
y = "Frequenza"
) +
theme_minimal()
I dati risultano distribuiti secondo una normale leggermente irregolare che presenta valori più probabili centrati sul pero di circa 3300 e con una coda allungata a sinistra. Quantifichiamo di quanto la distribuzione si allontana dalla perfetta regolarità della distribuzione normale grazie all’indici di forma di Asimmetria e al coefficiente di Curtosi.
#indice di asimmetria
moments::skewness(Peso)
## [1] -0.6470308
La distribuzione ha un indice di asimmetria negativo, ovvero una frequenza maggiore di valori alti.
#indice di asimmetria
moments::kurtosis(Peso)-3
## [1] 2.031532
Un coefficiente di Curtosi di due evidenzia una distribuzione Leptocurtica, più allungata rispetto alla distribuzione normale, a sottolineare la presenza di valori outlier che mettiamo in evidenza grazie ad un grafico boxplot della variabile risposta
#indice di asimmetria
boxplot(Peso)
Il grafico mostra una presenza numerosa di outliers che potrebbe influenzare il risultato del test di Shapiro Wilk sulla normalità di questa distribuzione. E’ lecito aspettarsi un rifiuto dell’ipotesi nulla
#indice di asimmetria
shapiro.test(Peso)
##
## Shapiro-Wilk normality test
##
## data: Peso
## W = 0.97066, p-value < 2.2e-16
Come atteso, il test di Shapiro Wilk conferma, con un p-value praticamente prossimo allo zero, che si rifiuta l’ipotesi di normalità della distribuzione della variabile Peso.
Esploriamo la colonna del dataset contenente gli anni della madre:
summary(Anni.madre)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 25.00 28.00 28.16 32.00 46.00
E’ subito possibile notare quello che è a tutti gli effetti un errore di inserimento, che vede assegnare il valore di 0 agli anni della madre, valore biologicamente impossibile. Per verificare la presenza di altri errori analizziamo il dataset ordinando in ordine crescente la colonna Anni.madre e verifichiamo la presenza di un altro valore anomalo, che assegna il valore 1 alla variabile.
Si aggiorna il dataset escludendo le righe che presentano i valori errati di 0 e 1 relativamente agli anni della madre. Si prosegue l’analisi con un totale di 2498 osservazioni
dati = dati %>% filter(Anni.madre > 1)
attach(dati)
## I seguenti oggetti sono mascherati da dati (pos = 3):
##
## Anni.madre, Cranio, Fumatrici, Gestazione, Lunghezza, N.gravidanze,
## Ospedale, Peso, Sesso, Tipo.parto
E’ ora possibile procedere al ricalcolo dei valori di sintesi della variabile
summary(Anni.madre)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13.00 25.00 28.00 28.19 32.00 46.00
Verifichiamo la distribuzione dei dati
ggplot(dati, aes(x = Anni.madre)) +
geom_histogram(bins = 30, fill = 3, col = 1) +
labs(
title = "Distribuzione dell'Età della Madre",
x = "Età (anni)",
y = "Frequenza"
) +
theme_minimal()
I dati risultano distribuiti secondo una normale con le code allungate, che presenta il numero maggiore di misurazioni relative all’età di 28 anni, ma con ampi scarti tra il valore massimo di 46 anni e il valore minimo di 13. Questa forbice fa presagire una possibilità che il test di Shapiro Wilk possa rifiutare l’ipotesi nulla di distribuzione normale proprio a causa della presenza di outliers marcati.
Calcoliamo, prima, gli indici di forma.
#indice di asimmetria
moments::skewness(Anni.madre)
## [1] 0.1510624
La distribuzione ha un indice di asimmetria che si discosta leggermente dallo zero verso destra, ad evidenziare una lieve asimmetria positiva ed una coda leggermente più ma con un effetto molto contenuto.
#indice di asimmetria
moments::kurtosis(Anni.madre)-3
## [1] -0.1056061
Anche il coefficiente di Curtosi è prossimo allo zero, ma negativo. Questo evidenzia una curva lievemente platicurtica, più piatta di una distribuzione normale e con valori meno concentrati attorno la media, ma con un effetto molto blando.
Saggiamo, utilizzando il test di Shapiro Wilk, l’ipotesi di normalità della curva.
#indice di asimmetria
shapiro.test(Anni.madre)
##
## Shapiro-Wilk normality test
##
## data: Anni.madre
## W = 0.99491, p-value = 1.477e-07
Il test restituisce un p-value molto basso che conferma il rifiuto dell’ipotesi di normalità. Considerando i valori prossimi allo zero sia dell’indice di asimmetria, sia del coefficiente di Curtosi, possiamo asserire che sono gli outliers a spiegare questa non normalità.
Verifichiamo anche visivamente, utilizzando il boxplot, la loro distribuzione
#indice di asimmetria
boxplot(Anni.madre)
Esploriamo i dati con un riassunto delle informazioni e con un’analisi visiva della loro distribuzione.
summary(Lunghezza)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 310.0 480.0 500.0 494.7 510.0 565.0
L’analisi riassuntiva mette in evidenza una maggiore distanza tra valore minimo e primo quartile rispetto a valore massimo e terzo quartile. Questo evidenzia come la distribuzione dei dati possa avere, verosimilmente, un’accentuata asimmetria.
ggplot(dati, aes(x = Lunghezza)) +
geom_histogram(bins = 30, fill = 3, col = 1) +
geom_density(aes(y = after_stat(density)), col = 2, linewidth = 1) +
labs(
title = "Distribuzione della Lunghezza alla Nascita",
x = "Lunghezza (mm)",
y = "Frequenza"
) +
theme_minimal()
Anche l’istogramma delle frequenze conferma l’ipotesi secondo cui la curva possa essere asimmetrica negativa, con un’elevata presenza di valori a sinistra della curva.
Procediamo a quantificare l’indici di forma di Asimmetria e il coefficiente di Curtosi.
#indice di asimmetria
moments::skewness(Lunghezza)
## [1] -1.514575
La distribuzione ha un indice di asimmetria negativo, a conferma di una frequenza maggiore di nella parte destra della distribuzione.
#indice di asimmetria
moments::kurtosis(Lunghezza)-3
## [1] 6.48093
Un coefficiente di Curtosi di due evidenzia una distribuzione Leptocurtica accentuata, più allungata rispetto alla distribuzione normale, a sottolineare la presenza di valori outlier che mettiamo in evidenza grazie ad un grafico boxplot della variabile risposta
#indice di asimmetria
boxplot(Lunghezza)
Il grafico mostra una presenza numerosa di outliers che potrebbe influenzare il risultato del test di Shapiro Wilk sulla normalità di questa distribuzione. Anche per questa variabile ci aspettiamo il rifiuto dell’ipotesi di normalità.
#indice di asimmetria
shapiro.test(Lunghezza)
##
## Shapiro-Wilk normality test
##
## data: Lunghezza
## W = 0.90944, p-value < 2.2e-16
Anche i numeri confermano l’analisi visiva, con un p-value praticamente prossimo allo zero.
Riassumiamo i dati relativi alla circonferenza del cranio, per una prima analisi.
summary(Cranio)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 235 330 340 340 350 390
Proprio come era lecito aspettarsi, si evidenzia una maggiore distanza tra minimo e primo quartile rispetto a massimo e terzo quantile. Tuttavia, la mediana uguale alla media potrebbe far pensare ad una distribuzione simmetrica.
ggplot(dati, aes(x = Cranio)) +
geom_histogram(bins = 30, fill = 3, col = 1) +
geom_density(aes(y = after_stat(density)), col = 2, linewidth = 1) +
labs(
title = "Distribuzione della circonferenza del Cranio alla Nascita",
x = "Cranio (mm)",
y = "Frequenza"
) +
theme_minimal()
L’istogramma mostra sicuramente una coda allungata a sinistra, ma non risolve il dubbio relativo alla possibile simmetria della distribuzione. E’, dunque, necessario procedere a quantificare l’indice di forma di Asimmetria e il coefficiente di Curtosi.
#indice di asimmetria
moments::skewness(Cranio)
## [1] -0.7850906
La distribuzione ha un indice di asimmetria negativo, a conferma di una frequenza maggiore di valori concentrati nella parte alta della distibuzione.
#indice di asimmetria
moments::kurtosis(Cranio)-3
## [1] 2.94487
Decisamente più marcato è il coefficiente di Curtosi, che evidenzia una distribuzione Leptocurtica accentuata con numerosi outliers.
#indice di asimmetria
boxplot(Cranio)
Gli outliers sono evidenti anche nel boxplot rappresentato, e sono proprio questi outliers a farci immaginare un rifiuto dell’ipotesi di normalità come risultato del test di Shapiro Wilk.
#indice di asimmetria
shapiro.test(Cranio)
##
## Shapiro-Wilk normality test
##
## data: Cranio
## W = 0.96358, p-value < 2.2e-16
Come ipotizzato, il p-value ha un valore piccolissimo che ci porta a rifiutare l’ipotesi di distribuzione normale dei dati.
Per quanto riguarda la variabile gestazione, è lecito attendersi una curva Leptocurtica con Asimmetria Negativa, in quanto, biologicamente, il periodo del parto è sempre fissato a 40 settimane, con la possibilità di ritardarlo non oltre la settimana 43. A questo si contrappongono i parti anticipati, per cause diverse, che fissano il limite inferiore parecchio più in basso della media e fanno presagire la presenza di diversi outliers. Verifichiamo, con un riassunto delle informazioni, che non vi siano valori anomali.
summary(Gestazione)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 25.00 38.00 39.00 38.98 40.00 43.00
I valori sono in linea con quelli attesi, procediamo all’analisi.
ggplot(dati, aes(x = Gestazione)) +
geom_histogram(bins = 30, fill = 3, col = 0) +
labs(
title = "Distribuzione delle Settimane di Gestazione",
x = "Settimane",
y = "Frequenza"
) +
theme_minimal()
#indice di asimmetria
moments::skewness(Gestazione)
## [1] -2.065131
#indice di asimmetria
moments::kurtosis(Gestazione)-3
## [1] 8.255516
Anche i risultati numerici e l’analisi visiva confermano le ipotesi.
Nel campione sono state selezionate 2394 madri Non fumatrici e 104 Fumatrici.
table(Fumatrici)
## Fumatrici
## 0 1
## 2394 104
Nel campione è stata fatta una distinzione tra parti naturali e parti cesarei, con un’incidenza, per questi ultimi, di 728 su 2498 osservazioni.
table(Tipo.parto)
## Tipo.parto
## Ces Nat
## 728 1770
I campioni sono stati prelevati da 3 ospedali diversi, con una frequenza equilibrata. Infatti l’analisi del dataset mostra 816 osservazioni presso l’ospedale 1, 848 presso l’ospedale 2 e 834 presso l’ospedale 3.
table(Ospedale)
## Ospedale
## osp1 osp2 osp3
## 816 848 834
Tra i neonati sono stati osservati 1243 Maschi e 1255 Femmine.
table(Sesso)
## Sesso
## F M
## 1255 1243
Per la grande maggioranza delle osservazioni è possibile notare come più di un terzo (ben 1095) erano alla prima gravidanza, con casi osservati di madri fino a 12 gravidanze.
table(N.gravidanze)
## N.gravidanze
## 0 1 2 3 4 5 6 7 8 9 10 11 12
## 1095 817 340 150 48 21 11 1 8 2 3 1 1
La tabella mostra evidentemente una distribuzione asimmetrica positiva della variabile, che vede la maggior parte dei valori concentrati nella parte sinistra, infatti più del 75% delle osservazioni non supera le due gravidanze.
In alcuni ospedali si fanno più cesarei
Questa ipotesi mette in relazione la frequenza di una variabile categorica con un’altra variabile categorica. Per saggiare questa ipotesi settiamo come ipotesi nulla una frequenza di parti cesarei uguale tra i tre ospedali, valutando come indipendenti tra loro le due variabili.
Per l’analisi dell’associazione tra due variabili categoriche si costruisce il test di indipendenza Chi-Quadrato.
Iniziamo con la definizione di una tabella dei parti, che metta in relazione ospedale e tipologia di parto.
Si inizia con la creazione di una tabella contenente le informazioni da saggiare.
info_parti <- table(Ospedale, Tipo.parto)
info_parti
## Tipo.parto
## Ospedale Ces Nat
## osp1 242 574
## osp2 254 594
## osp3 232 602
Una volta ottenuta la tabella, si applica la funzione chisq.test ai dati contenuti in essa per ottenere il risultato del test.
chi_square <- chisq.test(info_parti)
chi_square
##
## Pearson's Chi-squared test
##
## data: info_parti
## X-squared = 1.083, df = 2, p-value = 0.5819
Il test mostra un p-value di gran lunga superiore alla soglia del 5%, che permette di non rifiutare l’ipotesi nulla e di confermare l’indipendenza tra le due varibili nel campione analizzato.
Il grafico a barre stacked mostra, anche graficamente, un’omogeneità nelle proporzioni della tipologia di parto per ciascun ospedale.
ggplot(dati, aes(x = Ospedale, fill = Tipo.parto)) +
geom_bar(position = "fill") + # proporzioni
labs(title = "Proporzione dei tipi di parto per ospedale",
x = "Ospedale",
y = "Proporzione",
fill = "Tipo di parto") +
theme_minimal()
La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione
L’ipotesi proposta può essere saggiata in due modi:
La media dei dati della popolazione può essere estratta dalle informazioni condivise dall’OMS, che definisce, per i neonati, una media peso di 3250g e una media altezza di 50cm.
# Peso: ipotetica media di popolazione = 3250g
t.test(dati$Peso, mu = 3250)
##
## One Sample t-test
##
## data: dati$Peso
## t = 3.2529, df = 2497, p-value = 0.001158
## alternative hypothesis: true mean is not equal to 3250
## 95 percent confidence interval:
## 3263.577 3304.791
## sample estimates:
## mean of x
## 3284.184
# Lunghezza: ipotetica media di popolazione = 49,5cm
t.test(dati$Lunghezza, mu = 495)
##
## One Sample t-test
##
## data: dati$Lunghezza
## t = -0.57754, df = 2497, p-value = 0.5636
## alternative hypothesis: true mean is not equal to 495
## 95 percent confidence interval:
## 493.6628 495.7287
## sample estimates:
## mean of x
## 494.6958
L’analisi dei risultati ottenuti ci porta ad analizzare separatamente peso e lunghezza.
Per quanto riguarda il peso otteniamo un p-value di 0.0012 che è inferiore alla soglia di accettazione del 5%, inoltre il valore medio è al di fuori dell’intervallo di confidenza, per cui possiamo rifiutare l’ipotesi nulla e asserire che la media del peso è significativamente diversa da quella della popolazione, risultando più alto del valore di riferimento.
Discorso diverso, invece, per la lunghezza. I dati comunicano che il valore medio è all’interno dell’intevallo di confidenza e il p-value di 0.56, nettamente maggiore della soglia del 5% rafforza la tesi secondo cui la lunghezza del campione è uguale a quella della media della popolazione.
Le misure antropometriche sono significativamente diverse tra i due sessi
Per il test di questa ipotesi confronteremo, usando il t-test e lavorando su campioni indipendenti, le seguenti misure antropometriche: - lunghezza - peso - cranio
Saremo in grado, così, di verificare l’eventuale differenza tra le medie dei due gruppi.
Assumiamo, come ipotesi nulla, che le medie di ciascuna misura antropometrica siano uguali per maschi e femmine
# Peso
t.test(Peso ~ Sesso, data = dati)
##
## Welch Two Sample t-test
##
## data: Peso by Sesso
## t = -12.115, df = 2488.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -287.4841 -207.3844
## sample estimates:
## mean in group F mean in group M
## 3161.061 3408.496
# Lunghezza
t.test(Lunghezza ~ Sesso, data = dati)
##
## Welch Two Sample t-test
##
## data: Lunghezza by Sesso
## t = -9.5823, df = 2457.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -11.939001 -7.882672
## sample estimates:
## mean in group F mean in group M
## 489.7641 499.6750
# Cranio
t.test(Cranio ~ Sesso, data = dati)
##
## Welch Two Sample t-test
##
## data: Cranio by Sesso
## t = -7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -6.110504 -3.560417
## sample estimates:
## mean in group F mean in group M
## 337.6231 342.4586
PESO La media del peso nelle femmine è di 3161g contro i 3404g dei maschi. La differenza di 247g e il p-value estremamente basso ci dimostrano la negazione dell’ipotesi nulla. Possiamo asserire che i maschi pesano significativamente di più delle femmine.
LUNGHEZZA La differenza delle medie delle lunghezze è di un centimetro a favore dei maschi(499,7mm). Anche in questo caso il p-value basso conferma che è necessario rifiutare l’ipotesi nulla: i maschi sono significativamente più lunghi delle femmine.
CRANIO La media della circonferenza del cranio nelle femmine è di 337,6mm, 4,8mm in meno dei maschi. Il valore basso del p-value ci permette di asserire che i maschi hanno una circonferenza cranica maggiore delle femmine.
Per tutte e tre le variabili antropometriche considerate il sesso ha avuto un’influenza rilevante, comunque supportata dalla letteratura medica. Questo dettaglio ci sottolinea l’importanza di considerare il sesso nel modello di regressione.
Il modello di regressione lineare multipla che include tutte le variabili esplicative e permette di stimare l’effetto marginale di ciascuna variabile sul peso del neonato, mantenendo costante le altre variabili. Il modello ha questa struttura
Peso=β0+β1Anni.madre+β2N.gravidanze+β3Fumatrici+β4Gestazione+β5Lunghezza+β6Cranio+β7Tipo.parto+β8Ospedale+β9Sesso+ε
Ogni coefficiente βi rappresenta la variazione media del peso associata all’incremento unitario della variabile esplicativa corrispondente.
Creiamo un modello totale mod_tot e calcoliamo la funzione summary per visualizzare i parametri chiave di tutte le variabili.
mod_tot <- lm(Peso ~ Anni.madre + N.gravidanze + Fumatrici +
Gestazione + Lunghezza + Cranio +
Tipo.parto + Ospedale + Sesso,
data = dati)
summary(mod_tot)
##
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione +
## Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1123.26 -181.53 -14.45 161.05 2611.89
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6735.7960 141.4790 -47.610 < 2e-16 ***
## Anni.madre 0.8018 1.1467 0.699 0.4845
## N.gravidanze 11.3812 4.6686 2.438 0.0148 *
## Fumatrici -30.2741 27.5492 -1.099 0.2719
## Gestazione 32.5773 3.8208 8.526 < 2e-16 ***
## Lunghezza 10.2922 0.3009 34.207 < 2e-16 ***
## Cranio 10.4722 0.4263 24.567 < 2e-16 ***
## Tipo.partoNat 29.6335 12.0905 2.451 0.0143 *
## Ospedaleosp2 -11.0912 13.4471 -0.825 0.4096
## Ospedaleosp3 28.2495 13.5054 2.092 0.0366 *
## SessoM 77.5723 11.1865 6.934 5.18e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274 on 2487 degrees of freedom
## Multiple R-squared: 0.7289, Adjusted R-squared: 0.7278
## F-statistic: 668.7 on 10 and 2487 DF, p-value: < 2.2e-16
Per valutare la bontà del modello consideriamo il parametro R2 aggiustato. Con il valore di circa 0.73 ci dice che il modello è in grado di spiegare circa il 73% della variabilità del peso dei neonati, a dimostrazione che le variabili considerate sono in grado di spiegare buona parte delle dinamiche che regolano il peso del neonato.
L’errore standard residuo è di 273,9 grammi, valore accettabile che permette di promuoverne l’uso operativo, dato supportato anche dal basso valore del p-value.
Valutiamo singolarmente le variabili e, per ciascuna, verifichiamo segno e stima del coefficiente. Partiamo dalle variabili più significative, individuate dal valore del p-value basso. - Gestazione ha un coefficiente di +32,57. Questo significa che per ogni settimana di gestazione il neonato incrementa il proprio peso di 32,6 grammi;
Lunghezza ha un coefficiente di 10,29, per cui ogni millimetro di altezza incrementa il peso del neonato di 10,3g;
Cranio ha un coefficiente di 10,47, quindi ogni millimetro in più di circonferenza cranica incrementa il peso di circa 10,5g;
SessoM ha un coefficiente di 77,54, che indica che i maschi pesano 78g in più delle femmine.
Scendendo di significatività troviamo: - numero di gravidanze che incrementano il peso di 11g (effetto trascurabile)
parto naturale influenza il peso di 30g (probabilmente legato al fatto che i cesarei vengono spesso eseguiti prima del termine)
ospedale3 vede la nascita di neonati 28g più pesanti, ma potrebbe trattarsi di un effetto legato al campione
Le altre variabili (anni della madre, fumatrici, ospedale 2) non risultano significative.
Per verificare se esiste multicollinearità usiamo il Variance Inflation Factor che non evidenzia problemi se inferiore a 5.
## GVIF Df GVIF^(1/(2*Df))
## Anni.madre 1.190241 1 1.090982
## N.gravidanze 1.189278 1 1.090540
## Fumatrici 1.007426 1 1.003706
## Gestazione 1.695675 1 1.302181
## Lunghezza 2.086879 1 1.444604
## Cranio 1.631049 1 1.277125
## Tipo.parto 1.004227 1 1.002111
## Ospedale 1.004267 2 1.001065
## Sesso 1.040743 1 1.020168
Il modello non evidenzia problemi di multicollinearità rilevanti, dando valore alle variazioni influenzate dalle variabili significative che, pur essendo correlate biologicamente, hanno coefficienti stabili.
Per la scelta del modello ottimale seguiremo la procedura Stepwise finalizzata alla rimozione delle variabili non significative e al miglioramento del modello mediante trasformazioni. Partiamo dal modello completo mod_tot e iniziamo a rimuovere le seguenti variabili deboli: - Anni.madre (p = 0,43) - Fumatrici (p = 0,27) - Ospedaleosp2 (p = 0,4)
mod_1 <- lm(Peso ~ N.gravidanze +
Gestazione + Lunghezza + Cranio +
Tipo.parto + Sesso,
data = dati)
summary(mod_1)
##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Tipo.parto + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1129.14 -181.97 -16.26 160.95 2638.18
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6708.0171 136.0715 -49.298 < 2e-16 ***
## N.gravidanze 12.7356 4.3385 2.935 0.00336 **
## Gestazione 32.3253 3.7969 8.514 < 2e-16 ***
## Lunghezza 10.2833 0.3009 34.177 < 2e-16 ***
## Cranio 10.5063 0.4263 24.648 < 2e-16 ***
## Tipo.partoNat 30.1601 12.1027 2.492 0.01277 *
## SessoM 77.9171 11.1994 6.957 4.42e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.4 on 2491 degrees of freedom
## Multiple R-squared: 0.7277, Adjusted R-squared: 0.727
## F-statistic: 1109 on 6 and 2491 DF, p-value: < 2.2e-16
Il risultato ottenuto restituisce un R quadro di 0.727, rispetto allo 0.7278 ottenuto dal modello completo. Questa pulizia ci ha permesso di eliminare variabili inutili mantenendo pressoché inalterata la capacità predittiva.
Proviamo ad affinare ancora il modello provando a tagliare la variabile Tipo.parto
mod_2 <- lm(Peso ~ N.gravidanze +
Gestazione + Lunghezza + Cranio + Sesso,
data = dati)
summary(mod_2)
##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio +
## Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1149.37 -180.98 -15.57 163.69 2639.09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6681.7251 135.8036 -49.201 < 2e-16 ***
## N.gravidanze 12.4554 4.3416 2.869 0.00415 **
## Gestazione 32.3827 3.8008 8.520 < 2e-16 ***
## Lunghezza 10.2455 0.3008 34.059 < 2e-16 ***
## Cranio 10.5410 0.4265 24.717 < 2e-16 ***
## SessoM 77.9807 11.2111 6.956 4.47e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared: 0.727, Adjusted R-squared: 0.7265
## F-statistic: 1327 on 5 and 2492 DF, p-value: < 2.2e-16
Questo affinamento ulteriore ci permette di gestire un modello ancor più semplice a parità di prestazioni.
La variabile N.gravidanze ha un p-value inferiore alla soglia del 5%, quindi il criterio statistico ci suggerirebbe di mantenerla, ma conviene comunque valutarne il taglio comparando sia il criterio di valutazione AIC, sia il criterio BIC.
mod_3 <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso,
data = dati)
summary(mod_3)
##
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso,
## data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1138.1 -184.4 -17.4 163.3 2626.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6651.6732 135.5952 -49.055 < 2e-16 ***
## Gestazione 31.3262 3.7884 8.269 < 2e-16 ***
## Lunghezza 10.2024 0.3009 33.909 < 2e-16 ***
## Cranio 10.6706 0.4247 25.126 < 2e-16 ***
## SessoM 79.1027 11.2205 7.050 2.31e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 275.1 on 2493 degrees of freedom
## Multiple R-squared: 0.7261, Adjusted R-squared: 0.7257
## F-statistic: 1652 on 4 and 2493 DF, p-value: < 2.2e-16
AIC(mod_2, mod_3)
## df AIC
## mod_2 7 35152.89
## mod_3 6 35159.12
BIC(mod_2, mod_3)
## df BIC
## mod_2 7 35193.65
## mod_3 6 35194.06
Alla luce dei risultati ottenuti possiamo scegliere il modello 2 come migliore in quanto ha AI più basso di 6.3. Per quanto riguarda il BIC, possiamo considerare i due modelli come praticamente identici.
Per il modello BIC non c’è differenza tra tenere o rimuovere la variabile N.gravidanze, per AIC invece è meglio tenerla per non perdere informazioni e garantirsi un R quadro leggermente migliore. Dato che N.gravidanze è una variabile incrementale che migliora, seppur leggermente, il modello, decidiamo di tenerla.
Analizziamo, ora, graficamente la non linearità sulle variabili significative.
Partiamo dalla gestazione che, si suppone, non abbia un comportamento lineare in quanto è lecito attendersi un appiattimento della curva dopo la quarantesima settimana di gestazione.
Per verificare la distribuzione paragoniamo la relazione lineare (colorata in blu nel grafico seguente) con la forma reale dei dati (in rosso)
ggplot(dati, aes(x = Gestazione, y = Peso)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "loess", color = "red", linewidth = 1.2) +
geom_smooth(method = "lm", color = "blue", linetype = "dashed") +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Come era lecito attendersi, la curva cresce regolare fino alla settimana 40, oltre la quale tende ad appiattirsi.
Per misurare la sovrastima del modello nelle settimane finali di gravidanza, aggiungiamo l’effetto quadratico e misuriamone l’impatto
mod_gquad <- lm(Peso ~ N.gravidanze +
Gestazione + I(Gestazione^2) +
Lunghezza + Cranio + Sesso,
data = dati)
summary(mod_gquad)
##
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + I(Gestazione^2) +
## Lunghezza + Cranio + Sesso, data = dati)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1144.0 -181.5 -12.9 165.8 2661.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4646.7158 898.6322 -5.171 2.52e-07 ***
## N.gravidanze 12.5489 4.3381 2.893 0.00385 **
## Gestazione -81.2309 49.7402 -1.633 0.10257
## I(Gestazione^2) 1.5168 0.6621 2.291 0.02206 *
## Lunghezza 10.3502 0.3040 34.045 < 2e-16 ***
## Cranio 10.6376 0.4282 24.843 < 2e-16 ***
## SessoM 75.7563 11.2435 6.738 1.99e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 274.5 on 2491 degrees of freedom
## Multiple R-squared: 0.7276, Adjusted R-squared: 0.7269
## F-statistic: 1109 on 6 and 2491 DF, p-value: < 2.2e-16
AIC(mod_2, mod_gquad)
## df AIC
## mod_2 7 35152.89
## mod_gquad 8 35149.63
BIC(mod_2, mod_gquad)
## df BIC
## mod_2 7 35193.65
## mod_gquad 8 35196.21
Il coefficiente quadratico ha un p-value significativo e dal confronto dei modelli emerge che AIC premia il modello quadratico (preferisce modelli sovraparametrizzati) mentre BIC premia leggermente il modello 2. A fronte di una prestazione quasi identica, preferiamo la semplicità di interpretazione ddel modello 2 essendo la differenza AIC comunque piccola e l’influenza quadratica oltre le classiche 40 settimane di gestazione.
Le altre variabili significative, per natura biologica, dovrebbero avere una relazione lineare con il peso. A supporto di questa tesi analizziamo la distribuzione grafica della relazione peso-lunghezza.
ggplot(dati, aes(x = Lunghezza, y = Peso)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "loess", color = "red") +
geom_smooth(method = "lm", color = "blue", linetype = "dashed")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Le curve sono pressoché sovrapposte per lunghezza maggiore di 400mm. Nella parte bassa la LOESS diverge ma possiamo notare un numero decisamente inferiore di misurazioni e una maggiore dispersione.
Il modello che utilizzeremo, quindi, sarà il modello 2.
Il modello scelto è in grado di spiegare circa il 73% dela variabilità del peso neonatale (dato da R quadro aggiustato).
Calcoliamo il RMSE
rmse = sqrt(mean(resid(mod_2)^2))
rmse
## [1] 274.3666
In media il modello sbaglia la previsione del peso di circa 274g.
Analizziamo, ora, i residui
par(mfrow = c(2,2))
plot(mod_2)
Il grafico dei residui (Residuals vs Fitted) mostra punti sparsi attorno alla media di 0, senza pattern forti, quindi rispetta la linearità dei residui. E’ notabile una leggera eteroschedasticità accettabile per N~2500 campioni.
Q-Q Residuals mostra punti allineati al centro con deviazioni nelle code (già evidenziate dai grafici di linearità precedenti). Vi sono outliers che evidenziano la presenza, nei campioni, di pesi estremi.
Scale Location mostra una varianza crescente con il peso ma che non influenza la stabilità del modello ed una eteroschedasticità lieve data dall’aumento della variabilità dei residui per valori più elevati della variabile risposta.
Residuals Vs Leverage non mostra alcuna osservazione fortemente influente e afferma la stabilità del modello.
Applichiamo, ora, il modello ad una caso realistiche per misurarne la performance. Vogliamo stimare il peso di una neonata considerando queste variabili: - Numero gravidanze: 3 - Gestazione: 39 settimane - Lunghezza: 500mm - Cranio: 250mm - Sesso: F
osservazione <- data.frame(
N.gravidanze = 3,
Gestazione = 39,
Lunghezza = 500,
Cranio = 330,
Sesso = "F"
)
predict(mod_2, newdata = osservazione, interval = "prediction")
## fit lwr upr
## 1 3219.809 2680.517 3759.102
Le condizioni considerate restituiscono un peso medio di 3220g con un intervallo di previsione che va da 2681g a 3759g.
In questo intervallo possiamo confermare che il feto è in una situazione normale e non presenta alcun rischio evidente.
Il modello si è dimostrato utilizzabile, in grado di produrre risultati plausibili. Rendiamo visibili questi risultati mostrando il grafico delle previsioni del peso. La linea blu rapresenta la relazione stimata del peso rispetto alle settimane di gestazione. I punti neri sono i dati reali mentre in rosso evidenziamo la previsione del caso precedente.
ggplot(dati, aes(x = Gestazione, y = Peso)) +
geom_point(alpha = 0.2) +
geom_smooth(method = "lm", color = "blue") +
geom_point(data = osservazione, aes(x = Gestazione, y = 3219.779),
color = "red", size = 4) +
labs(
title = "Previsione del peso neonatale",
subtitle = "Il punto rosso rappresenta la previsione del modello",
x = "Gestazione",
y = "Peso (g)"
) +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
Il punto si colloca lungo la retta stimata, coerentemente con i dati
osservati, segno di una previsione realistica.
I dati raccolti mostrano una elevata dispersione verticale, data dalla variabilità naturale del peso neonatale, a partire dalla settimana 37, che significa che la previsione è in linea con quello che è il comportamento medio dei dati.
A supporto di quanto costruito, mostriamo le relazioni tra le variabili più significative
ggplot(dati, aes(x = Gestazione, y = Peso)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", color = "blue") +
labs(
title = "Relazione tra settimane di gestazione e peso neonatale",
x = "Settimane di gestazione",
y = "Peso (g)"
) +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
ggplot(dati, aes(x = Sesso, y = Peso, fill = Sesso)) +
geom_boxplot() +
labs(
title = "Distribuzione del peso per sesso",
x = "Sesso",
y = "Peso (g)"
) +
theme_minimal()
ggplot(dati, aes(x = Lunghezza, y = Peso)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", color = "blue") +
labs(
title = "Relazione tra lunghezza e peso neonatale",
x = "Lunghezza (mm)",
y = "Peso (g)"
) +
theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'
ggplot(dati, aes(x = factor(Fumatrici), y = Peso, fill = factor(Fumatrici))) +
geom_boxplot() +
labs(
title = "Peso neonatale per abitudine al fumo",
x = "Fumatrici (0 = No, 1 = Sì)",
y = "Peso (g)"
) +
theme_minimal()
Il progetto ha permesso di sviluppare un modello predittivo solido per la previsione del peso neonatale, utilizzando un numero considerevole di campioni estratti da tre strutture ospedaliere distinte.
L’analisi preliminare ci ha permesso uno studio approfondito delle caratteristiche del dataset, permettendoci di individuare ed escludere valori e misurazioni anomale e comprendere la natura delle variabili considerate.
Il risultato è stato un modello in grado di spiegare il 73% delle variabili che influenzano il peso neonatale, tanto parsimonioso quanto efficace, con un errore medio di previsione di circa 274g, basato sul contributo delle seguenti variabili: - settimane di gestazione, determinante principale per la crescita - misure antropometriche, fortemente correlate al peso - sesso, con maschi mediamente più pesanti - numero di gravidanze della madre, che ha dato un contributo positivo, seppur contenuto.
L’applicazione del modello a un caso reale ha quindi dimostrato la sua coerenza e la capacità di fornire una stima interpretabile, sia numericamente che graficamente, in linea con i risultati raccolti, centrando quello che era l’obiettivo iniziale di Neonatal Health Solutions.