Modello Statistico per la Previsione del Peso Neonatale

Contesto Aziendale

Azienda: Neonatal Health Solutions Obiettivo: Creare un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali. Il progetto mira a migliorare la gestione delle gravidanze ad alto rischio, ottimizzare le risorse ospedaliere e garantire migliori risultati per la salute neonatale.

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.

Dettagli del Progetto

1. Raccolta dei Dati e Struttura del Dataset

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.

2. Analisi e Modellizzazione

Analisi Preliminare

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:

    1. in alcuni ospedali si fanno più parti cesarei
    1. La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione
    1. Le misure antropometriche sono significativamente diverse tra i due sessi

Creazione del Modello di Regressione

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.

Selezione del Modello Ottimale

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.

Analisi della Qualità del Modello

Una volta ottenuto il modello finale, valuteremo la sua capacità predittiva utilizzando metriche come R2 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.

3. Previsioni e Risultati

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.

4. Visualizzazioni

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.

Conclusioni

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.

REGRESSIONE LINEARE MULTIPLA

import librerie

suppressMessages(library(dplyr))
suppressMessages(library(ggplot2))
suppressMessages(library(moments))
suppressMessages(library(patchwork))
suppressMessages(library(MASS))
suppressMessages(library(knitr))

Load data

setwd("//wsl.localhost/Ubuntu/home/fchiricotto/ProfessionAI/AI-courses/03-statistica-inferenziale/datasets/")
dati <- read.csv("./neonati.csv")
head(dati)
##   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
## 6         32            0         0         40 3200       495    340        Nat
##   Ospedale Sesso
## 1     osp3     M
## 2     osp1     F
## 3     osp2     M
## 4     osp2     M
## 5     osp3     F
## 6     osp2     F

1. Analisi Preliminare

A. Analisi descrittiva

Nella prima fase, esploreremo le variabili attraverso un’analisi descrittiva per comprenderne la distribuzione e identificare eventuali outlier o anomalie.

Variabili:

  • Anni.madre: variabile quantitativa discreta
  • N.gravidanze: variabile quantitativa discreta
  • Fumatrici: variabile qualitativa nominale dicotomica
  • Gestazione: variabile quantitativa discreta
  • Peso: variabile quantitativa continua
  • Lunghezza: variabile quantitativa continua
  • Cranio: variabile quantitativa continua
  • Tipo.parto: variabile qualitativa nominale
  • Ospedale: variabile qualitativa nominale
  • Sesso: variabile qualitativa nominale

Numero di osservazioni:

n <- nrow(dati)
n
## [1] 2500

Distribuzione di frequenze

Fumatrici
N = length(dati$Fumatrici)
# assolute
ni <- table(dati$Fumatrici)
# relative
fi <- ni / N
# trasformo la tabella in un dataframe
# trasformo la tabella in un dataframe
distr_freq_fumatrici <- as.data.frame(cbind(ni, fi))
distr_freq_fumatrici
##     ni     fi
## 0 2396 0.9584
## 1  104 0.0416

Nel campione di dati sono presenti per il 95% donne non fumatrici.

Tipo di parto
N = length(dati$Tipo.parto)
# assolute
ni <- table(dati$Tipo.parto)
# relative
fi <- ni / N
# trasformo la tabella in un dataframe
# trasformo la tabella in un dataframe
distr_freq_tipo_parto <- as.data.frame(cbind(ni, fi))
distr_freq_tipo_parto
##       ni     fi
## Ces  728 0.2912
## Nat 1772 0.7088

Il 70% delle nascite del campione è avvenuto tramite parto naturale mentre il restante 30% con parto cesareo.

Ospedale
N = length(dati$Ospedale)
# assolute
ni <- table(dati$Ospedale)
# relative
fi <- ni / N
# trasformo la tabella in un dataframe
# trasformo la tabella in un dataframe
distr_freq_ospedale <- as.data.frame(cbind(ni, fi))
distr_freq_ospedale
##       ni     fi
## osp1 816 0.3264
## osp2 849 0.3396
## osp3 835 0.3340

Le nascite registrate all’interno del campione provengono da 3 ospedali. Le 2500 nascite del campione sono distribuite tra i 3 ospedali come riportato in tabella precedente; l’ospedale 2 è quello che ha registrato il maggior numero di nascite.

Sesso
N = length(dati$Sesso)
# assolute
ni <- table(dati$Sesso)
# relative
fi <- ni / N
# trasformo la tabella in un dataframe
# trasformo la tabella in un dataframe
distr_freq_tipo_sesso <- as.data.frame(cbind(ni, fi))
distr_freq_tipo_sesso
##     ni     fi
## F 1256 0.5024
## M 1244 0.4976

I nuovi nascituri si dividono all’interno del campione quasi equamente fra maschi e femmine con una lieve maggioranza di femmine.

Distribuzioni

Peso
plot(density(dati$Peso),         
    xlab = "Peso")

skewness(dati$Peso)
## [1] -0.6470308
kurtosis(dati$Peso) - 3
## [1] 2.031532

La variabile Peso oggetto di studio presenta una distribuzione con asimmetria negativa (valori di skewness negativi) e leptocurtica (valori di curtosi positivi) ovvero più allungata rispetto alla distribuzione normale.

Vediamo anche come sono distribuite le altre variabili quantitative continue:

par(mfrow=c(1,2))
plot(density(dati$Lunghezza),         
    xlab = "Lunghezza")
plot(density(dati$Cranio),         
    xlab = "Diametro craniale")

skewness(dati$Lunghezza)
## [1] -1.514699
kurtosis(dati$Lunghezza) - 3
## [1] 6.487174

La variabile Lunghezza presenta una distribuzione con asimmetria negativa (valori di skewness negativi) e leptocurtica (valori di curtosi positivi) ovvero più allungata rispetto alla distribuzione normale.

summary(dati$Lunghezza)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   310.0   480.0   500.0   494.7   510.0   565.0
skewness(dati$Cranio)
## [1] -0.7850527
kurtosis(dati$Cranio) - 3
## [1] 2.946206

La variabile Cranio presenta una distribuzione con asimmetria negativa (valori di skewness negativi) e leptocurtica (valori di curtosi positivi) ovvero più allungata rispetto alla distribuzione normale.

summary(dati$Cranio)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     235     330     340     340     350     390

Guardiamo ora alle distribuzione delle variabili quantitative discrete:

ggplot(data = dati, aes(x = Gestazione)) +
  geom_histogram(binwidth = 1, 
                 col = "black", 
                 fill = "#69b3a2", 
                 alpha = 0.7) + 
  geom_text(stat = "count", aes(label = after_stat(count)), 
            vjust = -0.5, size = 3) + 
  labs(title = "Distribuzione del Numero di settimane di Gestazione",
       x = "Numero di settimane",
       y = "Frequenza Assoluta") +
  theme_minimal() + 
  scale_x_continuous(breaks = seq(min(dati$Gestazione), max(dati$Gestazione), by = 1)) + 
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) 

summary(dati$Gestazione)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   25.00   38.00   39.00   38.98   40.00   43.00

Il numero di settimane di gestazione media è di 40 settimane con un minimo di 25 ed un massimo di 43.

ggplot(data = dati, aes(x = N.gravidanze)) +
  geom_histogram(binwidth = 1, 
                 col = "black", 
                 fill = "#69b3a2", 
                 alpha = 0.7) + 
  geom_text(stat = "count", aes(label = after_stat(count)), 
            vjust = -0.5, size = 3) + 
  labs(title = "Distribuzione del Numero di Gravidanze Passate",
       x = "Numero di Gravidanze",
       y = "Frequenza Assoluta") +
  theme_minimal() + 
  scale_x_continuous(breaks = seq(min(dati$N.gravidanze), max(dati$N.gravidanze), by = 1)) + 
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) 

summary(dati$N.gravidanze)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  1.0000  0.9812  1.0000 12.0000

Su un totale di 2500 parti la maggior parte, il 43%, è stato per mettere alla luce il primo figlio. Il numero massimo di gravidanze passate è 12.

ggplot(data = dati, aes(x = Anni.madre)) +
  geom_histogram(binwidth = 2, 
                 col = "black", 
                 fill = "#69b3a2", 
                 alpha = 0.7) + 
  labs(title = "Distribuzione degli anni della madre",
       x = "Eta della madre",
       y = "Frequenza Assoluta") +
  theme_minimal() + 
  scale_x_continuous(breaks = seq(min(dati$Anni.madre), max(dati$Anni.madre), by = 2))

summary(dati$Anni.madre)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   25.00   28.00   28.16   32.00   46.00

Oserviamo che l’età media delle madri è di 28 anni. Si riscontra un minimo in 0 cosa che non è possibile quindi ci deve essere qualche anomalia nei dati. L’eta massima registrata è di 46 anni.

Esaminiamo in piu attenzione i dati anomali; andiamo a vedere quante osservazioni presentano Anni.madre minore di un valore di 10 anni.

table(dati$Anni.madre[dati$Anni.madre<10])
## 
## 0 1 
## 1 1

Abbiamo 2 osservazioni in cui l’eta della madre risulta essere addirittura 0 o 1. Essendo due soli dati decido di rimuoverli dal dataset.

dati <- dati[dati$Anni.madre > 10, ]
summary(dati$Anni.madre)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   13.00   25.00   28.00   28.19   32.00   46.00

Analisi della variabile risposta Peso

summary(dati$Peso)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     830    2990    3300    3284    3620    4930

Il peso medio del campione risulta essere di 3.284 Kg.

boxplot(dati$Peso ~ dati$Sesso, 
        xlab = "Genere", 
        ylab = "Peso (g)", 
        main = "Distribuzione del Peso per Sesso")

Graficamente osserviamo che il peso medio del neonato risulta essere più alto nei maschi che nelle femmine. Quantifichiamo questa differenza:

summary(dati$Peso[dati$Sesso == "F"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     830    2900    3160    3161    3470    4930
summary(dati$Peso[dati$Sesso == "M"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     980    3150    3430    3408    3720    4810

Le femmine sembrano pesare 250 gr in meno dei maschi.

Vediamo come si relazione invece il Peso alla variabile fumatrici:

#boxplot
boxplot(dati$Peso~dati$Fumatrici, 
        xlab = "Fumatrici", 
        ylab = "Peso (g)", 
        main = "Distribuzione del Peso per Fumatrici")

Nel caso di donne fumatrici il peso del neonato sembra essere minore.

Vediamo le statistiche di riepilogo per i 3 ospedali:

summary(dati$Peso[dati$Ospedale == "osp1"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     830    2980    3295    3270    3602    4600
summary(dati$Peso[dati$Ospedale == "osp2"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     990    2970    3280    3271    3600    4930
summary(dati$Peso[dati$Ospedale == "osp3"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     900    3030    3300    3312    3640    4810

Nell’ospedale 3 sembrano esserci nascite di bambini con in media un peso maggiore.

Vediamolo anche graficamente esaminando la distribuzione di Peso dei neonati per ospedale:

boxplot(Peso~Ospedale, data=dati, 
        xlab = "Ospedale", 
        ylab = "Peso (g)", 
        main = "Distribuzione del Peso per Ospedale")

Graficamente si conferma che le distribuzioni nei 3 ospedali sono molto simili.

Consideriamo ora il numero di parti cesarei e naturali ci sono in ciascun ospedale:

dati_raggruppati <- dati %>%
  group_by(Ospedale, Tipo.parto) %>%
  summarise(Totale = n())
## `summarise()` has grouped output by 'Ospedale'. You can override using the
## `.groups` argument.
ggplot(dati_raggruppati, aes(x = Ospedale, y = Totale, fill = Tipo.parto)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Confronto Parti Cesarei e Non Cesarei per Ospedale",
       x = "Ospedale",
       y = "Numero di Parti") +
  geom_text(aes(label = Totale), position = position_dodge(0.9), vjust = -0.5, size = 3.5)+
  theme_minimal()

L’ospedale 3 è quello con il maggior numero di parti naturali ed il minor numero di parti naturali.

Confrontiamo le distribuzioni di peso per capire se si vedono evidenti differenze tra parti cesarei e naturali:

#boxplot
boxplot(Peso~Tipo.parto, data = dati, 
        xlab = "Tipo.parto", 
        ylab = "Peso (g)", 
        main = "Distribuzione del Peso per Tipo.parto")

Da una analisi grafica non sembra esservi differenza evidente di peso medio dei neonati nati con parto naturale o cesario.

Quantifichiamo questa differenza in media:

summary(dati$Peso[dati$Tipo.parto == "Nat"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     830    3000    3300    3285    3620    4900
summary(dati$Peso[dati$Tipo.parto == "Ces"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     930    2980    3290    3282    3610    4930

In media non si osserva una differenza di peso per tipologie di parto differenti.

TEST IPOTESI 1: in alcuni ospedali si fanno più parti cesarei

Possiamo utilizzare un test Chi-Quadro per confrontare la distribuzione delle modalita di una variabile categorica (Tipo.parto) tra i livelli di una seconda variabile categorica (Ospedale).

Costruisco innanzitutto la tabella di contingenza:

# VECCHIO CODICE ERRATO:
#tabella_contingenza <- table(dati$Ospedale[dati$Tipo.parto == "Ces"], dati$Tipo.parto[dati$Tipo.parto == "Ces"])
#tabella_contingenza

# Creazione della tabella di contingenza completa (Ospedale vs Tipo.parto)
tabella_contingenza <- table(dati$Ospedale, dati$Tipo.parto)
tabella_contingenza
##       
##        Ces Nat
##   osp1 242 574
##   osp2 254 594
##   osp3 232 602

Eseguo il test del Chi-Quadro per verificare se le frequenze dei parti cesarei e naturali differiscono significativamente tra gli ospedali:

test_chi <- chisq.test(tabella_contingenza)
test_chi
## 
##  Pearson's Chi-squared test
## 
## data:  tabella_contingenza
## X-squared = 1.083, df = 2, p-value = 0.5819

Con un p-value di 0,6, quindi molto superiore alla soglia di significatività dello 0.05, possiamo accettare l’ipotesi nulla e quindi concludere che non si osservano differenze statisticamente significative nella distribuzione del tipo di parto tra ospedali.

TEST IPOTESI 2: La media del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione

Calcolo il valore medio di Peso e Lunghezza:

mu_peso_pop <- 3300
mu_lunghezza_pop <- 50

saggiamo l’ipotesi tramite un test-t:

t.test(dati$Peso, mu = mu_peso_pop)
## 
##  One Sample t-test
## 
## data:  dati$Peso
## t = -1.505, df = 2497, p-value = 0.1324
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3263.577 3304.791
## sample estimates:
## mean of x 
##  3284.184

il p-value è maggiore della soglia di significatività del 5% pertanto non si rifiuta l’ipotesi nulla, ovvero la media del Peso per il campione è significativamente uguale a quella della popolazione.

TEST IPOTESI 3: Le misure antropometriche sono significativamente diverse tra i due sessi
PESO
#boxplot
par(mfrow=c(1,2))
boxplot(dati$Peso, ylab = "Peso (g)")
boxplot(Peso~Sesso, data = dati, ylab = "Peso (g)")

Graficamente si osserva una differenza fra le medie di Peso tra i due sessi.

summary(dati$Peso[dati$Sesso == "M"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     980    3150    3430    3408    3720    4810
summary(dati$Peso[dati$Sesso == "F"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     830    2900    3160    3161    3470    4930

Se guardiamo al calcolo della medie osserviamo, in accordo con quanto osservato graficamente, come il peso medio dei maschi sia di 250 gr.

Tramite l’utilizzo di un t-test saggiamo quindi se la differenza del peso medio dei neonati di sesso maschile e femminile è uguale a 0.

Sistema di Ipotesi:

  • H0: diff=0
  • H1: diff!=0
t.test(
  dati$Peso[dati$Sesso == "M"], # maschi
  dati$Peso[dati$Sesso == "F"], # femmine
)
## 
##  Welch Two Sample t-test
## 
## data:  dati$Peso[dati$Sesso == "M"] and dati$Peso[dati$Sesso == "F"]
## t = 12.115, df = 2488.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  207.3844 287.4841
## sample estimates:
## mean of x mean of y 
##  3408.496  3161.061

Il p-value è praticamente 0 e quindi inferiore al valore di soglia del 5% e pertanto l’ipotesi nulla è rifiutata. La differenza tra il peso medio delle femmine ed il peso medio dei neonati maschi è da considerarsi significativamente differente da 0.

LUNGHEZZA
#boxplot
par(mfrow=c(1,2))
boxplot(dati$Lunghezza, ylab = "Lunghezza (cm)")
boxplot(Lunghezza~Sesso, data = dati, ylab = "Lunghezza (cm)")

Graficamente si osserva una differenza fra le medie di Lunghezza tra i due sessi.

summary(dati$Lunghezza[dati$Sesso == "M"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   320.0   490.0   500.0   499.7   515.0   560.0
summary(dati$Lunghezza[dati$Sesso == "F"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   310.0   480.0   490.0   489.8   505.0   565.0

In accordo con quanto visto graficamente, otteniamo che in media i maschi sono 10 cm piu lunghi delle femmine.

Tramite l’utilizzo di un t-test saggiamo quindi se la differenza della lunghezza media dei neonati di sesso maschile e femminile è uguale a 0.

Sistema di Ipotesi:

  • H0: diff=0
  • H1: diff!=0
t.test(
  dati$Lunghezza[dati$Sesso == "M"], # maschi
  dati$Lunghezza[dati$Sesso == "F"], # femmine
)
## 
##  Welch Two Sample t-test
## 
## data:  dati$Lunghezza[dati$Sesso == "M"] and dati$Lunghezza[dati$Sesso == "F"]
## t = 9.5823, df = 2457.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   7.882672 11.939001
## sample estimates:
## mean of x mean of y 
##  499.6750  489.7641

Il p-value è praticamente 0 e quindi inferiore al valore di soglia del 5% e pertanto l’ipotesi nulla è rifiutata e quindi la differenza tra la lunghezza media dei neonati di sesso maschile e femminile è da considerarsi significativamente differente da 0.

DIAMETRO CRANIALE
#boxplot
par(mfrow=c(1,2))
boxplot(dati$Cranio, ylab = "Cranio (cm)")
boxplot(Cranio~Sesso, data=dati, ylab = "Cranio (cm)")

Graficamente si osserva una lieve differenza fra le medie del diametro craniale tra i due sessi. Quantifichiamo questa differenza calcolando le medie per i due sessi:

summary(dati$Cranio[dati$Sesso == "M"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   265.0   334.0   343.0   342.5   352.0   390.0
summary(dati$Cranio[dati$Sesso == "F"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   235.0   330.0   340.0   337.6   348.0   390.0

In media osserviamo che i maschi hanno un diametri craniale piu grande di 5 cm.

Tramite l’utilizzo di un t-test saggiamo se la differenza del diametro craniale medio tra maschi e femmine dei neonati è uguale a 0.

Sistema di Ipotesi:

  • H0: diff=0
  • H1: diff!=0
t.test(
  dati$Cranio[dati$Sesso == "M"], # maschi
  dati$Cranio[dati$Sesso == "F"], # femmine
)
## 
##  Welch Two Sample t-test
## 
## data:  dati$Cranio[dati$Sesso == "M"] and dati$Cranio[dati$Sesso == "F"]
## t = 7.4366, df = 2489.4, p-value = 1.414e-13
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  3.560417 6.110504
## sample estimates:
## mean of x mean of y 
##  342.4586  337.6231

Il p-value è praticamente 0 e quindi inferiore al valore di soglia del 5% e pertanto l’ipotesi nulla è rifiutata e quindi le differenze del diametro craniale medio tra neonati di sesso maschile e femminile sono da considerarsi significativamente differenti da 0.

2. Creazione del Modello di Regressione

Iniziamo sviluppando un modello di regressione lineare multipla che includa tutte le variabili rilevanti. Possiamo poi 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.

A. NORMALITA DELLA DISTRIBUZIONE DI Y

Verifico innanzitutto che la variabile risposta sia approssimativamente distribuita NORMALMENTE con un test di shapiro e gli indici di forma.

Utilizziamo la libreria moments per stimare gli indici di forma:

moments::skewness(dati$Peso) 
## [1] -0.6474036
moments::kurtosis(dati$Peso)-3 
## [1] 2.028753

La Skewness negativa indica che la distribuzione è leggermente asimmetrica verso sinistra o negativa. La Kurtosis di 2 suggerisce che la distribuzione ha code più pesanti rispetto alla normale, con probabilità aumentata di valori estremi.

Utilizziamo il test di shapiro per saggiare la normalità della distribuzione dei pesi:

shapiro.test(dati$Peso)
## 
##  Shapiro-Wilk normality test
## 
## data:  dati$Peso
## W = 0.97068, p-value < 2.2e-16

Con un p-value molto piccolo (minore della soglia 0.05) si rifiuta l’ipotesi di normalita (per un livello di significatività alfa del 5%).

Analisi grafica

par(mfrow=c(1,2))
plot(density(dati$Peso), xlab = "Peso (g)")
boxplot(dati$Peso, main = "Boxplot", col = "lightblue", ylab = "Peso (g)")

La distribuzione è più appuntita e con code più lunghe rispetto alla normale. Come mostrato anche dal boxplot potrebbero esserci valori estremi (outlier), soprattutto nella coda sinistra (considerando la skewness negativa).

B. analisi delle correlazioni

Sfruttiamo il grafico pairs per il calcolo delle correlazioni:

# Selezioniamo solo le colonne numeriche
dati_num <- dati[, sapply(dati, is.numeric)]

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- (cor(x, y))
  txt <- format(c(r, 1), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = 1.5)
}
#correlazioni

suppressWarnings(pairs(dati_num,lower.panel=panel.cor, upper.panel=panel.smooth))

Sfruttando il grafico guardiamo alle correlazioni tra le variabili quantitative:

Correlazioni con la variabile risposta Peso:

  • si osserva una alta correlazione positiva con le variabili Lunghezza e Cranio
  • si ha una moderata correlazione positiva con la variabile Gestazione
  • Le altre variabili hanno una correlazione con il Peso prossima allo 0

Vediamo inoltre che per Gestazione sembra esserci una piccola non linearità dei dati.

Correlazioni tra regressori:

  • lieve correlazione positiva tra la Anni.madre e N.gravidanze;
  • correlazione media positiva tra Gestazione e le variabili Lunghezza e Cranio; un periodo piu lungo di gestazione fa pensare ad una conseguente crescita maggiore del neonato e quindi ci si attendeva una correlazione positiva del periodo di gestazione con le misure antropometriche.
  • lieve correlazione negativa (-0.1 ) tra la variabile Gestazione e le variaibli Anni.madre e N.gravidanze; le donne con maggiore età o maggiore numero di gravidanze hanno un periodo di gestazione decrescente.

Mentre le variabili esplicative con una alta correlazione con la variabile risposta sono quelle che portano piu informazione sulla risposta y riuscendo a spiegare gran parte della sua variabilita, dobbiamo porre attenzione ai regressori molto correlati fra loro in quanto possono portare a problemi di multicollinearita.

Approfondiremo questo ultimo aspetto con altre metriche.

Per quanto riguarda le variabili qualitative (Fumatrici, Sesso, Tipo.parto, Ospedale), per queste non sono utili gli scatterplot e la correlazione ma si usano altri strumenti statistici quali i BOXPLOT.

FUMATRICI

#boxplot
par(mfrow=c(1,2))
boxplot(dati$Peso, ylab = "Peso (g)")
boxplot(Peso~Fumatrici, data=dati, ylab = "Peso (g)")

Guardando alla distribuzione condizionata si osserva come i valori di peso per donne fumatrici e non fumatrici siano diversi; il peso dei neonati figli di donne fumatrici hanno un peso inferiore.

Calcolo il valore medio del peso nei due casi:

mean(dati$Peso[dati$Fumatrici==0]) 
## [1] 3286.262
mean(dati$Peso[dati$Fumatrici==1]) 
## [1] 3236.346

Il Peso medio nel campione dei neonati nel caso di madri fumatrici è minore (50 grammi in meno).

Eseguo un t-test relativamente alla variabile Fumatrici per saggiare l’ipotesi nulla per cui la variabile peso medio sia la stessa sia per Fumatrici che non.

t.test(Peso~Fumatrici, data = dati)
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Fumatrici
## t = 1.0362, df = 114.12, p-value = 0.3023
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -45.5076 145.3399
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.262        3236.346

Il T-test ci dice che con un p-value maggiore della soglia di significativita del 5% non si rifiuta l’ipotesi nulla Il peso medio dei neonati per Fumatrici e non fumatrici non è significativamente diverso.

Ci aspettiamo quindi un coefficiente beta di regressione per questa variabile Fumatrici non significativo; potrebbe non avere quindi senso inserire questa variabile nel modello.

SESSO

#boxplot
par(mfrow=c(1,2))
boxplot(dati$Peso, ylab = "Peso (g)")
boxplot(Peso~Sesso, data = dati, ylab = "Peso (g)")

Guardando alla distribuzione condizionata si osserva come i valori di peso per M e F siano diversi (M pesano piu delle F).

# PESO MEDIO PER M E F
mean(dati$Peso[dati$Sesso=="M"]) 
## [1] 3408.496
mean(dati$Peso[dati$Sesso=="F"]) 
## [1] 3161.061

Il peso medio per i neonati maschi è superiore (250 grammi) al peso medio delle femmine.

Eseguo un t-test relativamente alla variabile sesso: saggiamo l’ipotesi nulla per cui la variabile peso medio sia la stessa per entrambe i sessi.

t.test(dati$Peso~dati$Sesso)
## 
##  Welch Two Sample t-test
## 
## data:  dati$Peso by dati$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

Il T-test ci dice che: - con un p-value molto minore della soglia di significativita del 5% si rifiuta l’ipotesi nulla. Il peso medio dei due sessi è significativamente diverso.

Ci aspettiamo un coefficiente beta di regressione per questa variabile sesso significativo; per questo potrebbe avere senso inserire questa variabile nel modello.

OSPEDALE

#boxplot
par(mfrow=c(1,2))
boxplot(dati$Peso, ylab = "Peso (g)")
boxplot(Peso~Ospedale, data=dati, ylab = "Peso (g)")

Guardando alla distribuzione condizionata si osserva come i valori di peso si distribuiscono per i 3 ospedali; non sembra esserci una differenza importante tra i 3 ospedali in termine di peso dei neonati.

mean(dati$Peso[dati$Ospedale=="osp1"]) 
## [1] 3270.266
mean(dati$Peso[dati$Ospedale=="osp2"]) 
## [1] 3270.507
mean(dati$Peso[dati$Ospedale=="osp3"]) 
## [1] 3311.709

Il peso medio dei neonati per ospedale sembra essere simile. Saggiamo l’ipotesi nulla per cui la variabile peso medio sia la stessa per tutti gli ospedali. Per esguire il t-test contemporaneamente per i 3 ospedali sfruttiamo la funzione pairs:

pairwise.t.test(dati$Peso, 
                dati$Ospedale, 
                paired = F, # campioni indipendenti
                pool.sd = T, 
                p.adjust.method = "bonferroni") 
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  dati$Peso and dati$Ospedale 
## 
##      osp1 osp2
## osp2 1.00 -   
## osp3 0.33 0.32
## 
## P value adjustment method: bonferroni

Tutti i p-value sono elevati e sopra la soglia del 5%. Tutte le differenze in media del peso dei neonati possono essere considerate significativamente nulle.

Ci aspettiamo un coefficiente beta di regressione per questa variabile Ospedale non significativo; per questo potrebbe avere senso non inserire questa variabile nel modello.

D. Creazione del Modello di Regressione Multipla

Costruisco un primo modello includendo tutte le variabili:

mod<-lm(Peso~. , data = dati)
summary(mod)
## 
## Call:
## lm(formula = Peso ~ ., 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
  • Gestazione : coeff positivo e molto significativo (ogni settimana di gestazione in piu fa si che il peso del neonato aumenti in media di 30 grammi)
  • Lunghezza: coeff positivo e molto significativo (ogni centimetro di lunghezza in piu fa si che il peso del neonato aumenti in media di 10 grammi)
  • Cranio: coeff positivo e molto significativo (ogni centimetro di diametro craniale in piu fa si che il peso del neonato aumenti in media di 10 grammi)
  • N.gravidanze: coeff positivo e significativito (ad ogni gravidanza in piu si osserva un peso del neonato superiore in media di 11 grammi)

Per quanto riguarda le variabili qualitative R le trasforma in dummy e stima il parametro di una delle due modalita usando l’altra come baseline:

  • SessoM: tenendo fissate le altre variabili nei maschi si rileva significativamente un peso medio di 77.5 grammi in piu che nelle femmine.
  • Tipo.partoNat : tenendo fissate le altre variabili nei parti naturali si rileva significativamente un peso medio in media di 30 grammi in piu che nei parti cesarei.

La variabile Ospedale nelle sue 3 modalita mostra che: - Ospedaleosp3: tenendo fissate le altre variabili nell’ospedale 3 si rileva un peso medio di 28 grammi superiore agli altri ospedali. - Ospedaleosp2 non risulta essere significativo.

Le variabili Fumatrici e Anni.madre sembrano non essere significative.

Il coefficiente di Determinazione R2 ha un valore di 0.73. Nel caso di regressione lineare multipla è sempre meglio guardare al coefficiente di Determinazione aggiustato R2-ADJ che è comunque di 0.73: indica un buon adattamento, ma c’è ancora il 27% della varianza nel peso non spiegata dal modello.

Selezione dei regressori

rimozione Anni.madre

Per iniziare proviamo a rimuovere le variabili non significative quali Ospedale, Anni.madre,

mod2<-update(mod, ~.-Anni.madre)
summary(mod2)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1113.82  -180.30   -16.22   160.66  2616.32 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6708.6189   136.0211 -49.320  < 2e-16 ***
## N.gravidanze     12.5833     4.3400   2.899  0.00377 ** 
## Fumatrici       -30.4268    27.5455  -1.105  0.26944    
## Gestazione       32.2996     3.7997   8.501  < 2e-16 ***
## Lunghezza        10.2916     0.3008  34.209  < 2e-16 ***
## Cranio           10.4874     0.4257  24.638  < 2e-16 ***
## Tipo.partoNat    29.6654    12.0892   2.454  0.01420 *  
## Ospedaleosp2    -10.9509    13.4442  -0.815  0.41541    
## Ospedaleosp3     28.5171    13.4986   2.113  0.03474 *  
## SessoM           77.6452    11.1849   6.942 4.91e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274 on 2488 degrees of freedom
## Multiple R-squared:  0.7288, Adjusted R-squared:  0.7279 
## F-statistic: 743.1 on 9 and 2488 DF,  p-value: < 2.2e-16
  • Non si osserva un grande cambiamento
  • R2 resta lo stesso
  • N-gravidanze è piu significativo di prima

rimozione Ospedale

Proviamo a rimuovere anche la variabile Ospedale:

mod3<-update(mod2, ~.-Ospedale)
summary(mod3)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1129.88  -181.36   -16.24   160.63  2634.62 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6708.8092   136.0640 -49.306  < 2e-16 ***
## N.gravidanze     12.9927     4.3439   2.991  0.00281 ** 
## Fumatrici       -31.8823    27.5803  -1.156  0.24780    
## Gestazione       32.5970     3.8039   8.569  < 2e-16 ***
## Lunghezza        10.2684     0.3011  34.098  < 2e-16 ***
## Cranio           10.5015     0.4262  24.637  < 2e-16 ***
## Tipo.partoNat    30.4244    12.1041   2.514  0.01201 *  
## SessoM           78.1031    11.1998   6.974 3.94e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.4 on 2490 degrees of freedom
## Multiple R-squared:  0.7278, Adjusted R-squared:  0.7271 
## F-statistic: 951.3 on 7 and 2490 DF,  p-value: < 2.2e-16
  • Non si osserva un grande cambiamento
  • R2 resta lo stesso

rimozione Fumatrici

Proviamo a rimuovere anche la variabile Fumatrici

mod4<-update(mod3, ~.-Fumatrici)
summary(mod4)
## 
## 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
  • Non si osserva un grande cambiamento
  • R2 resta lo stesso

rimozione Tipo.parto

mod5<-update(mod4, ~.-Tipo.parto)
summary(mod5)
## 
## 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
  • Non si osserva un grande cambiamento rispetto al modello precedente.
  • R2 resta lo stesso 0.73.

Le variabili rimanenti sembrano essere tutti significative quindi mi fermerei qui nel togliere le variabili dal modello.

Aggiungo termine non lineare

Guardiamo il grafico di dispersione per capire meglio la relazione tra Peso e Gestazione:

plot(dati$Peso, dati$Gestazione, pch=20)

Il grafico potrebbe suggerire una lieve dipendenza non lineare (quadratica tra le due variabili). Proviamo ad inserire un termine quadratico.

mod6<-update(mod5, ~.+I(Gestazione^2))
summary(mod6)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2), 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    
## 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 ***
## I(Gestazione^2)     1.5168     0.6621   2.291  0.02206 *  
## ---
## 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

L’aggiunta di un termine quadratico nella Gestazione non sembra portare benefici e complicherebbe il modello ulteriormente, pertanto tenderei a scartare questo modello.

INTERAZIONE GESTAZIONE * TIPO PARTO

#CORREZIONE:
# mod_fumo_parto <- lm(Peso~.+Gestazione*Tipo.parto,data=dati)

mod_fumo_parto <- update(mod5, ~.+Gestazione*Tipo.parto,data=dati)
summary(mod_fumo_parto)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Tipo.parto + Gestazione:Tipo.parto, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1135.50  -181.98   -15.56   161.44  2639.12 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              -6451.5734   242.5207 -26.602  < 2e-16 ***
## N.gravidanze                12.6428     4.3386   2.914   0.0036 ** 
## Gestazione                  25.9500     6.2709   4.138 3.62e-05 ***
## Lunghezza                   10.2731     0.3010  34.135  < 2e-16 ***
## Cranio                      10.4981     0.4262  24.629  < 2e-16 ***
## SessoM                      78.5007    11.2073   7.004 3.18e-12 ***
## Tipo.partoNat             -307.4415   264.5754  -1.162   0.2453    
## Gestazione:Tipo.partoNat     8.6542     6.7751   1.277   0.2016    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.4 on 2490 degrees of freedom
## Multiple R-squared:  0.7279, Adjusted R-squared:  0.7271 
## F-statistic: 951.4 on 7 and 2490 DF,  p-value: < 2.2e-16

FUMO * GESTAZIONE

#CORREZIONE
#mod_fumo_gestazione <- lm(Peso~.+Fumatrici*Gestazione,data=dati)

mod_fumo_gestazione <- update(mod5, ~.+Fumatrici*Gestazione,data=dati)
summary(mod_fumo_gestazione)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Fumatrici + Gestazione:Fumatrici, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.83  -181.58   -16.95   163.76  2635.68 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -6699.6947   136.7290 -49.000  < 2e-16 ***
## N.gravidanze            12.7452     4.3470   2.932   0.0034 ** 
## Gestazione              33.2007     3.8418   8.642  < 2e-16 ***
## Lunghezza               10.2252     0.3011  33.957  < 2e-16 ***
## Cranio                  10.5313     0.4265  24.693  < 2e-16 ***
## SessoM                  78.7436    11.2241   7.016 2.94e-12 ***
## Fumatrici              795.7016   757.5315   1.050   0.2936    
## Gestazione:Fumatrici   -21.0469    19.2830  -1.091   0.2752    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2490 degrees of freedom
## Multiple R-squared:  0.7273, Adjusted R-squared:  0.7265 
## F-statistic: 948.6 on 7 and 2490 DF,  p-value: < 2.2e-16

ETA MADRE * NUMERO GRAVIDANZE

#CORREZIONE
#mod_eta_gravidanze <- lm(Peso~.+Anni.madre*N.gravidanze,data=dati)

mod_eta_gravidanze <- update(mod5, ~.+Anni.madre*N.gravidanze, data=dati)
summary(mod_eta_gravidanze)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Anni.madre + N.gravidanze:Anni.madre, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1159.62  -181.89   -15.04   164.14  2634.02 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             -6.713e+03  1.428e+02 -47.012  < 2e-16 ***
## N.gravidanze             1.185e+01  2.694e+01   0.440    0.660    
## Gestazione               3.269e+01  3.823e+00   8.552  < 2e-16 ***
## Lunghezza                1.025e+01  3.010e-01  34.038  < 2e-16 ***
## Cranio                   1.052e+01  4.274e-01  24.624  < 2e-16 ***
## SessoM                   7.790e+01  1.122e+01   6.946 4.77e-12 ***
## Anni.madre               9.108e-01  1.357e+00   0.671    0.502    
## N.gravidanze:Anni.madre -2.307e-02  8.382e-01  -0.028    0.978    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.8 on 2490 degrees of freedom
## Multiple R-squared:  0.7271, Adjusted R-squared:  0.7263 
## F-statistic: 947.7 on 7 and 2490 DF,  p-value: < 2.2e-16

L’aggiunta dei termini di interazione non ha portato alcun beneficio ed ha elevato la complessita del modello, pertanto tenderei ad escludere questi modelli.

3. Selezione del modello migliore

Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC), selezioniamo il modello più parsimonioso, eliminando le variabili non significative. Verranno considerati anche modelli con interazioni tra le variabili e possibili effetti non lineari.

Criterio AIC

AIC(mod, mod2, mod3, mod4, mod5, mod6, mod_eta_gravidanze, mod_fumo_gestazione, mod_fumo_parto)
##                     df      AIC
## mod                 12 35145.57
## mod2                11 35144.06
## mod3                 9 35149.33
## mod4                 8 35148.67
## mod5                 7 35152.89
## mod6                 8 35149.63
## mod_eta_gravidanze   9 35156.28
## mod_fumo_gestazione  9 35154.46
## mod_fumo_parto       9 35149.03

Secondo il criterio AIC il modello migliore è il 2 ovvero il modello che utlizza come regressori le varaibili:

  • N.gravidanze
  • Gestazione
  • Lunghezza
  • Cranio
  • Tipo.parto
  • Sesso
  • Fumatrici
  • Ospedale

Criterio BIC

BIC(mod, mod2, mod3, mod4, mod5, mod6, mod_eta_gravidanze, mod_fumo_gestazione, mod_fumo_parto)
##                     df      BIC
## mod                 12 35215.45
## mod2                11 35208.12
## mod3                 9 35201.73
## mod4                 8 35195.25
## mod5                 7 35193.65
## mod6                 8 35196.21
## mod_eta_gravidanze   9 35208.69
## mod_fumo_gestazione  9 35206.87
## mod_fumo_parto       9 35201.44

Secondo il criterio BIC, che tende a penalizzare i modelli over-parametrizzati a favore dei modelli piu semplici, il modello migliore è il 5 ovvero il modello che utlizza come regressori le varaibili:

  • N.gravidanze
  • Gestazione
  • Lunghezza
  • Cranio
  • Sesso

Confrontiamo la scelta con la procedura stepwise automatizzata

stepwise.mod <- MASS::stepAIC(mod, # modello con tutte le variabili (FULL MODEL)
              direction = "both", # direzione della procedura (forward, backwars, both)
              k=log(n)) # criterio da usare (AIC:2, BIC:log(n), n osservazioni)
## Start:  AIC=28118.62
## Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Anni.madre    1     36710 186779904 28111
## - Fumatrici     1     90677 186833870 28112
## - Ospedale      2    687555 187430749 28112
## - N.gravidanze  1    446244 187189438 28117
## - Tipo.parto    1    451073 187194266 28117
## <none>                      186743194 28119
## - Sesso         1   3610705 190353899 28159
## - Gestazione    1   5458852 192202046 28183
## - Cranio        1  45318506 232061700 28654
## - Lunghezza     1  87861708 274604902 29074
## 
## Step:  AIC=28111.29
## Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Fumatrici     1     91599 186871503 28105
## - Ospedale      2    693914 187473818 28105
## - Tipo.parto    1    452049 187231953 28110
## <none>                      186779904 28111
## - N.gravidanze  1    631082 187410986 28112
## + Anni.madre    1     36710 186743194 28119
## - Sesso         1   3617809 190397713 28151
## - Gestazione    1   5424800 192204704 28175
## - Cranio        1  45569477 232349381 28649
## - Lunghezza     1  87852027 274631931 29066
## 
## Step:  AIC=28104.69
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Ospedale + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Ospedale      2    702925 187574428 28098
## - Tipo.parto    1    444404 187315907 28103
## <none>                      186871503 28105
## - N.gravidanze  1    608136 187479640 28105
## + Fumatrici     1     91599 186779904 28111
## + Anni.madre    1     37633 186833870 28112
## - Sesso         1   3601860 190473363 28145
## - Gestazione    1   5358199 192229702 28168
## - Cranio        1  45613331 232484834 28642
## - Lunghezza     1  88259386 275130889 29063
## 
## Step:  AIC=28098.42
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Sesso
## 
##                Df Sum of Sq       RSS   AIC
## - Tipo.parto    1    467626 188042054 28097
## <none>                      187574428 28098
## - N.gravidanze  1    648873 188223301 28099
## + Ospedale      2    702925 186871503 28105
## + Fumatrici     1    100610 187473818 28105
## + Anni.madre    1     44184 187530244 28106
## - Sesso         1   3644818 191219246 28139
## - Gestazione    1   5457887 193032315 28162
## - Cranio        1  45747094 233321522 28636
## - Lunghezza     1  87955701 275530129 29051
## 
## Step:  AIC=28096.81
## Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## 
##                Df Sum of Sq       RSS   AIC
## <none>                      188042054 28097
## - N.gravidanze  1    621053 188663107 28097
## + Tipo.parto    1    467626 187574428 28098
## + Ospedale      2    726146 187315907 28103
## + Fumatrici     1     92548 187949505 28103
## + Anni.madre    1     45366 187996688 28104
## - Sesso         1   3650790 191692844 28137
## - Gestazione    1   5477493 193519547 28161
## - Cranio        1  46098547 234140601 28637
## - Lunghezza     1  87532691 275574744 29044
summary(stepwise.mod)
## 
## 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

Secondo la procedura automatizzata stepwise il modello migliore è il modello con regressori:

  • N.gravidanze
  • Gestazione
  • Lunghezza
  • Cranio
  • Sesso

che equivale al modello 5.

Controlliamo multicollinearita:

car::vif(mod5)
## N.gravidanze   Gestazione    Lunghezza       Cranio        Sesso 
##     1.023462     1.669779     2.075747     1.624568     1.040184

Tutti i VIF sono sotto la soglia indicativa di 5 quindi non abbiamo problemi di multicollinearità.

4. Analisi della Qualità del Modello

Una volta ottenuto il modello finale, valutiamo la sua capacità predittiva utilizzando metriche come R2 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 trovati i regressori che rappresentano il miglior compromesso tra sintesi e perdita di informazione resta l’analisi dei residui (ovvero della parte erratica).

A. ANALISI GRAFICA

Dividiamo la finestra in 4 parti e osserviamo i grafici dei residui

par(mfrow=c(2,2))
plot(mod5)

Osservando i 4 grafici notiamo che:

    1. in alto a sinistra: i punti non sembrano essere sparsi casualmente intorno alla media 0
    1. in alto a destra: rappresenta i residui in relazione con i quantili di una distro normale; i punti si trovano prevalentemente sulla bisettrice del grafico quindi hanno una distribuzione normale
    1. in basso a sinistra: non sembrano avere una distribuzione casuale intorno ad un valore di Y.
    1. in basso a destra: si valutano i valori anomali (leverage o outliers); nel grafico sono riportate tratteggiate le soglie da non oltrepassare (0.5 di avvertimento, 1 di allarme). Nessun punto supera le soglie. Non sembrano esserci problemi rilevanti di valori influenti.

B. METODI NUMERICI

Valori di leva

#leverage
lev<-hatvalues(mod5)

Grafico con le soglie dei valori di leva:

plot(lev)
p<-sum(lev) # somma dei leverage
n<-length(lev) # numero osservazioni
soglia=2*p/n
abline(h=soglia,col=2) # soglia dei leverage

kable(sum(lev>soglia), col.names = c("Numero Valori di Leva"), format = "html")
Numero Valori di Leva
152

Come vediamo ci sono molti punti oltre la soglia. Sono osservazioni che si trovano lontano rispetto al resto delle osservazioni nello spazio dei regressori.

Outliers

plot(rstudent(mod5))
abline(h=c(-2,2)) 

car::outlierTest(mod5) 
##       rstudent unadjusted p-value Bonferroni p
## 1551 10.046230         2.6345e-23   6.5810e-20
## 155   5.025345         5.3818e-07   1.3444e-03
## 1306  4.824963         1.4848e-06   3.7092e-03

Sono stati trovati 3 valori outlier.

dati[c(1551, 155, 1306), ]
##      Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1553         30            4         0         35 4520       520    360
## 155          30            0         0         36 3610       410    330
## 1307         29            0         0         42 3560       510    355
##      Tipo.parto Ospedale Sesso
## 1553        Nat     osp2     F
## 155         Nat     osp1     M
## 1307        Nat     osp1     F

Il primo valore risulta essere molto anomalo; analizzando nel dettaglio il dato osservo che si tratta di una nenonata alla 35 settimana di gestazione con un peso di 4.5 Kg ed una lunghezza oltre i 520 cm che sono valori alti rispetto alla media. Decido comunque di tenere questo dato nel dataset.

Distanza di Cook

Possiamo utilizzare la distanza di Cook per valutarli entrambe

cook<-cooks.distance(mod5)
plot(cook,ylim = c(0,1)) 

# LA DISTANZA MASSIMA DI COOK
max(cook)
## [1] 0.8297645

Una delle osservazioni supera il valore della soglia di avvertimento dello 0.5.

Omoschedasticità

lmtest::bptest(mod5)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod5
## BP = 90.297, df = 5, p-value < 2.2e-16

Si rifiuta l’ipotesi nulla; la varianza dei residui non è costante! Si ha eteroschedasticità significativa.

Autocorrelazione dei residui

lmtest::dwtest(mod5) 
## 
##  Durbin-Watson test
## 
## data:  mod5
## DW = 1.9532, p-value = 0.1209
## alternative hypothesis: true autocorrelation is greater than 0

Non si rifiuta l’ipotesi nulla: i residui non sono correlati!

Normalità

shapiro.test(mod5$residuals) 
## 
##  Shapiro-Wilk normality test
## 
## data:  mod5$residuals
## W = 0.97414, p-value < 2.2e-16

Si rifiuta l’ipotesi nulla di Normalità: i residui non sono normalmente distribuiti.

Vediamo la densita

plot(density(residuals(mod5)))

# si hanno delle code:

La distribuzione a campana sembra avere una deformazione sulla punta, riusltando cosi più appuntita e la coda destra è più accentuata della sinistra.

Il modello di previsione che abbiamo stimato non supera alcune delle assunzioni sui residui, ma nonostante questo proviamo ad utilizzarlo per fare delle previsioni e valutiamone i risultati.

4. Previsioni e Risultati

Una volta validato il modello, lo useremo per fare previsioni pratiche. Ad esempio, potremmo stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.

Per eseguire la predizione di esempio costruisco un modello di regressione multipla con i soli regressori Sesso, N.gravidanze e Gestazione:

mod_test<-lm(Peso~N.gravidanze+Gestazione+Sesso , data = dati)
summary(mod_test)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1494.51  -272.83   -14.09   267.03  1892.58 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3138.961    175.495 -17.886  < 2e-16 ***
## N.gravidanze    23.337      6.511   3.584 0.000344 ***
## Gestazione     162.087      4.501  36.009  < 2e-16 ***
## SessoM         165.125     16.741   9.864  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 414.3 on 2494 degrees of freedom
## Multiple R-squared:  0.3786, Adjusted R-squared:  0.3778 
## F-statistic: 506.5 on 3 and 2494 DF,  p-value: < 2.2e-16
esempio_test <- data.frame(
  Gestazione = c(39),
  Sesso = c("F"),
  N.gravidanze =c(3)
)

previsioni <- predict(mod_test, newdata = esempio_test)
output_df <- data.frame(Previsione = previsioni)
kable(output_df, col.names = c("Valore Predetto"), format = "html")
Valore Predetto
3252.431

Il modello ha stimato un peso medio della bambina di 3.252 Kg. Sembra essere un valore ragionevole considerando che il peso medio del campione era stato stimato essere di 3.161 KG.

Possiamo provare a fare la stessa stima con il modello 5, da me scelto, prendendo ad esempio i valori di Lunghezza e Cranio medi di neonate femmine del campione:

  • Lunghezza 490 cm
  • Cranio 338 cm
esempio_test_2 <- data.frame(
  Gestazione = c(39),
  Sesso = c("F"),
  N.gravidanze =c(3),
  Lunghezza = c(490),
  Cranio = c(338) 
)

previsioni <- predict(mod5, newdata = esempio_test_2)
output_df <- data.frame(Previsione = previsioni)
kable(output_df, col.names = c("Valore Predetto"), format = "html")
Valore Predetto
3201.682

Otteniamo con il modello 5 una stima di 3.201 Kg.

5. Visualizzazioni

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.

Vediamo con il seguente grafico come impatta il sesso sul modello di predizione

ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
               y=Peso,
               col=Sesso),
             position = "jitter")+
  geom_smooth(aes(x=Gestazione, 
                  y=Peso,
                  col=Sesso),se=F,method = "lm")
## `geom_smooth()` using formula = 'y ~ x'

Vediamo nel grafico come è molto influente il numero di settimane di gestazione nel determinare il peso del neonato. Inoltre si osserva una differenza tra bambini maschi e femmine con i maschi che mostrano nel corso delle settimane avere un peso superiore alle femmine.

dati$FumatriciSINO <- factor(dati$Fumatrici, levels = c(0, 1), labels = c("NO", "SI"))

ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
               y=Peso,
               col=FumatriciSINO))+
  geom_smooth(aes(x=Gestazione, 
                  y=Peso,
                  col=FumatriciSINO),se=F,method = "lm")
## `geom_smooth()` using formula = 'y ~ x'

Osservando la stima dell’andamento del peso rispetto alle settimane di gestazione si osserva come il peso stimato per i figli di madri fumatrici sia inferiore a quello di figli di madri non fumatrici.

CORREZIONI

analisi

  • anni.madre e settimane di gestazione sono continue (TEMPO)
  • è corretto da fare il chisq test, ma non capisco perchè lo hai fatto nella tabella solo con le frequene del cesareo. FAcendo così il test è distorto perhcè assumi che in tutti gli ospedali ci siano lo stesso identico numero di nascite, ma sebbene sia simile, non è così.

regressione:

  • non è detto che le anomalie ricadano dei residui
  • tipo di parto e ospedale potevano essere escluse a priori “logicamente” dal modello in quanto non hanno molto sneso in ottica previsionale
  • i commenti delle stime vanno fatti IN MEDIA non “DI CIRCA”
  • perchè quando hai inserito le interazioni hai reinserito tutte le variabili scartate prima?
  • per la previsione non si deve fare un modello ridotto, se no si perde tutta la potenza predittiva, devi usare la media per le altre variabili

Qualità dell’analisi

A parte le sbavature segnalate l’analisi è tutto sommato fatta bene, segue tutti i punti e la teoria. Però ti consiglio vivamente di sistemare le cose segnalate ;)

Qualità del documento

  • Ci sono dei print di console grezzi che non stanno benissimo… sempre meglio inserire i dati almeno in una tabella semplice per sistemare l’output, tipo kable()
  • i messaggi e i warnings della console, come per l’importazione di librerie e altro, si possono anche omettere per fare tutto più pulito
  • non stampare trafile di dati (tipo gli outliers)