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:

  1. 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.

  2. 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).

  3. 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.

  4. 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.

  5. 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.

1. Raccolta dei Dati e Struttura del Dataset

Per costruire il modello predittivo, abbiamo raccolto dati su 2500 neonati provenienti da tre ospedali.

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.

dati = read.csv("neonati.csv", stringsAsFactors = T)
n = nrow(dati)
quantitative_columns = c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")
qualitative_columns = c("Fumatrici", "Tipo.parto", "Ospedale", "Sesso")

2. Analisi e Modellizzazione

2.1 Analisi Preliminare

Le variabili raccolte nel dataset includono:

  • Anni.madre: quantitativa discreta, misura dell’età della madre in anni.
  • N.gravidanze: quantitativa discreta, quante gravidanze ha avuto la madre.
  • Fumatrici: qualitativa binaria (0=non fumatrice, 1=fumatrice).
  • Gestazione: quantitativa discreta, numero di settimane di gestazione.
  • Peso: quantitativa discreta, peso alla nascita in grammi (Target).
  • Lunghezza: quatitativa discreta, lunghezza del neonato in millimetri.
  • Cranio: quantitativa discreta, diametro craniale in millimetri.
  • Tipo.parto: qualitativa binaria (“Ces”=cesareo, “Nat”=naturale), tipologia di parto.
  • Ospedale: qualitativa binaria (“osp1”, “osp2” o “osp3”), ospedale di riferimento.
  • Sesso: qualitativa binaria (“M”=maschio, “F”=femmina), sesso del neonato.

Abbiamo quindi, 2500 osservazioni di 9 variabili esplicative e una variabile risposta (o target), ovvero il Peso.

Calcoliamo gli indici di posizione, forma e variabilità per le variabili quantitative: Anni.madre, N.gravidanze, Gestazione, Peso, Lunghezza, Cranio.

# Definiamo una funzione per calcolare tutti gli indici
compute_indexes = function(x) {
  min_x     = min(x, na.rm=TRUE)
  q1        = quantile(x, 0.25, na.rm=TRUE)
  median_x  = median(x, na.rm=TRUE)
  mean_x    = mean(x, na.rm=TRUE)
  q3        = quantile(x, 0.75, na.rm=TRUE)
  max_x     = max(x, na.rm=TRUE)
  iqr_x     = IQR(x, na.rm=TRUE)
  var_x     = var(x, na.rm=TRUE)
  sd_x      = sd(x, na.rm=TRUE)
  cv_x      = sd_x / mean_x
  skew_x    = skewness(x)
  kurt_x    = kurtosis(x)
  c(
    Min       = min_x,
    '1st Qu.' = q1,
    Median    = median_x,
    Mean      = mean_x,
    '3rd Qu.' = q3,
    Max       = max_x,
    IQR       = iqr_x,
    Variance  = var_x,
    Std       = sd_x,
    CV        = cv_x,
    Skewness  = skew_x,
    Kurtosis  = kurt_x
  )
}

data_summary = as.data.frame(
  t(round(sapply(dati[, quantitative_columns], compute_indexes), 2))
)
kable(rownames_to_column(data_summary, "Variable"),
      format = 'html', digits = 2) %>%
  kable_styling(full_width = T)
Variable Min 1st Qu..25% Median Mean 3rd Qu..75% Max IQR Variance Std CV Skewness Kurtosis
Anni.madre 0 25 28 28.16 32 46 7 27.81 5.27 0.19 0.04 3.38
N.gravidanze 0 0 1 0.98 1 12 1 1.64 1.28 1.31 2.51 13.99
Gestazione 25 38 39 38.98 40 43 2 3.49 1.87 0.05 -2.07 11.26
Peso 830 2990 3300 3284.08 3620 4930 630 275665.68 525.04 0.16 -0.65 5.03
Lunghezza 310 480 500 494.69 510 565 30 692.67 26.32 0.05 -1.51 9.49
Cranio 235 330 340 340.03 350 390 20 269.79 16.43 0.05 -0.79 5.95

Dalla tabella riassuntiva possiamo fare le seguenti considerazioni preliminari:

  • La variabile Anni.madre presenta una lievissima asimmetria positiva leptocurtica e un indice di variabilità del 19%, con un minimo decisamente anomalo di zero anni. Tale valore può essere interpretato come un errore di inserimento dell’informazione, quindi andrebbe rimosso dal dataset per non influenzare il modello.
  • La variabile N.gravidanze presenta la variabilità decisamente più alta, una distribuzione asimmetrica positiva leptocurtica.
  • La variabile Gestazione risulta avere una distribuzione asimmetrica negativa e leptocurtica con variabilità del 5% e un minimo di 25 settimane. Tale valore potrebbe indicare la presenza di parti prematuri nel dataset. I valori di Media e Mediana sono molto simili, quindi l’asimmetria può essere determinata da alcuni valori leverage di gestazioni particolarmente lunghe.
  • Il Peso risulta avere una lieve asimmetria negativa leptocurtica con variabilità del 16%. Il valore minimo di 830 grammi conferma la presenza di parti prematuri.
  • La Lunghezza e il diametro del Cranio sembrano seguire caratteristiche correlate, ovvero variabilità del 5% con una asimmetrica negativa leptocurtica.

Graficando la variabile Anni.madre confermiamo la distribuzione normale, nonostante la presenza di due osservazioni vicine allo zero (una è pari a 0 anni), mentre le restanti superano abbondantemente i 10 anni.

par(mfrow=c(1,2))
boxplot(dati$Anni.madre)
hist(dati$Anni.madre, col = 'orange', main = "", xlab = 'Peso', breaks = 18)

Anddiamo, quindi, a rimuovere dal dataset le osservazioni per le quali Anni.madre è inferiore a 10 anni.

dati = dati %>%
  filter(dati$Anni.madre >= 10)
n = nrow(dati)

Per quanto riguarda la variabile Gestazione, seguendo le indicazioni dell’Organizzazione Mondiale della Sanità (clicca qui per consultare i dati), si definisce parto pretermine quando questo avviene prima della 37° settimana. In particolare, moderatamente pretermine tra la 32° e 37°, molto pretermine tra la 28° e 32°, ed estremamente pretermine prima della 28° settimana. Visualizziamo allora la variabile tramite boxplot e scatterplot.

par(mfrow=c(1,2))
boxplot(dati$Gestazione)
hist(dati$Gestazione, col = 'orange', main = "", xlab = 'Gestazione', breaks = 18)

Osservando i grafici, possiamo confermare la distribuzione asimmetrica negativa per via di outlier che suggeriscono la presenza nel campione di parti pretermine.

Possiamo definire una nuova variabile Tipo.gestazione, con valore “regolare” se la gestazione ha raggiunto la 37° settimana, “molt. pretermine” se ha superato la 28° settimane e “estr. pretermine” negli altri casi, ed osservare come si distribuisce.

dati$Tipo.gestazione = factor(ifelse(dati$Gestazione >= 37,
                                     "regolare",
                                     ifelse(dati$Gestazione >= 28, 
                                            "molt. pretermine",
                                            "estr. pretermine")), 
                              levels = c("regolare","molt. pretermine","estr. pretermine"))
qualitative_columns = append(qualitative_columns, "Tipo.gestazione")

ni = table(dati$Tipo.gestazione)
fi = table(dati$Tipo.gestazione)/n
Ni = cumsum(ni)
Fi = cumsum(fi)

freq_distr = as.data.frame(cbind(ni, fi, Ni, Fi))
kable(freq_distr, format = "html", digits = 3) %>%
  kable_styling(full_width = T)
ni fi Ni Fi
regolare 2336 0.935 2336 0.935
molt. pretermine 158 0.063 2494 0.998
estr. pretermine 4 0.002 2498 1.000

Le gestazioni pretermine costituiscono il 6% del campione in esame. Vediamo, quindi, come la variabile Peso è influenzata dal tipo di gestazione. Mostriamo, quindi, il boxplot condizionato per ottenere le distribuzioni del Peso per ogni modalità di Tipo.gestazione.

par(mfrow=c(1,2))
boxplot(dati$Peso ~ dati$Tipo.gestazione, xlab = 'Tipo Gestazione', ylab = 'Peso')
plot(dati$Gestazione, dati$Peso, xlab = 'Gestazione', ylab = 'Peso')
abline(v=28, col=2)
text(27, 4900, "28°", col=2)
abline(v=32, col=7)
text(31, 4900, "32°", col=7)
abline(v=37, col=4)
text(36, 4900, "37°", col=4)

Il peso dei neonati tramite gestazioni regolari sembra essere mediamente più alto rispetto a quello pretermine (come è lecito pensare) e presenta diversi outlier. Inoltre, la variabile Peso sembra avere un andamento crescente non lineare con le settimane di gestazione.

Grafichiamo adesso le variabili Lunghezza e diametro del Cranio.

par(mfrow=c(1,2))
boxplot(dati$Lunghezza)
hist(dati$Lunghezza, col = 'orange', main = "", xlab = 'Lunghezza', breaks = 18)

par(mfrow=c(1,2))
boxplot(dati$Cranio)
hist(dati$Cranio, col = 'orange', main = "", xlab = 'Diam. Cranio', breaks = 18)

Dai grafici possiamo osservare diversi valori anomali che confermano l’asimmetria negativa della distribuzione. Se analizziamo gli scatterplot includendo anche l’informazione sul tipo di gestazione, possiamo notare un pattern chiaro dell’andamento del peso rispetto alla lunghezza e al diametro del cranio.

par(mfrow=c(1,2))
plot(dati$Lunghezza, dati$Peso, col=dati$Tipo.gestazione,
     xlab = "Lunghezza", ylab = "Peso")
plot(dati$Cranio, dati$Peso, col=dati$Tipo.gestazione,
     xlab = "Diam. Cranio", ylab = "Peso")

Si può notare una regione di addensamento delle osservazioni relativa alle gestazioni regolari (in nero), ed una coda di valori più bassi relativi alle gestazioni pretermine (in rosso e verde). L’andamento del peso è decisamente crescente all’aumentare della lunghezza e del diametro del cranio, quindi potremmo dire che le nascite pretermine mediamente determinano misure antropometriche inferiori rispetto a nascite regolari.

Per quanto riguarda le variabili qualitative Fumatrici, Tipo.parto, Tipo.gestazione, Ospedale e Sesso, calcoliamo innanzitutto le frequenze percentuali, per capire se sono bilanciate rispetto alle modalità assunte. Convertiamo, preliminarmente, la variabile Fumatrici in tipo factor.

dati$Fumatrici = as.factor(dati$Fumatrici)
table_freq = list()
for (var in qualitative_columns) {
  values = dati[[var]]
  ni = table(values)
  fi = round(ni / length(values), 2)
  
  df_tmp = data.frame(
    Variabile = var,
    Categoria = names(fi),
    Frequeza_assoluta = as.integer(ni),
    Frequenza_relativa = as.numeric(fi)
  )
  
  table_freq[[var]] = df_tmp
}


frequenze = do.call(rbind, table_freq)
kable(frequenze, format = "html", digits = 2) %>%
  kable_styling(full_width = T)
Variabile Categoria Frequeza_assoluta Frequenza_relativa
Fumatrici.1 Fumatrici 0 2394 0.96
Fumatrici.2 Fumatrici 1 104 0.04
Tipo.parto.1 Tipo.parto Ces 728 0.29
Tipo.parto.2 Tipo.parto Nat 1770 0.71
Ospedale.1 Ospedale osp1 816 0.33
Ospedale.2 Ospedale osp2 848 0.34
Ospedale.3 Ospedale osp3 834 0.33
Sesso.1 Sesso F 1255 0.50
Sesso.2 Sesso M 1243 0.50
Tipo.gestazione.1 Tipo.gestazione regolare 2336 0.94
Tipo.gestazione.2 Tipo.gestazione molt. pretermine 158 0.06
Tipo.gestazione.3 Tipo.gestazione estr. pretermine 4 0.00

La variabile Sesso risulta perfettamente bilanciata, come anche la variabile Ospedale. La varibile Fumatrici, invece, risulta fortemente sbilanciata, infatti solo il 4% delle madri del campione risulta fumatrice. Anche la variabile Tipo.parto risulta sbilanciata, con il 30% dei parti di tipo cesareo. Infine, dalla variabile Tipo.gestazione, notiamo che il 94% delle osservazioni si riferiscono a gestazioni regolari (dalla 37° settimana).

Essendo dei conteggi, possono essere combinate con la variabile target Peso, e osservare come questa si distribuisce rispetto alle modalità assunte dalle variabili suddette. Creiamo dei boxplot condizionati.

attach(dati)
par(mfrow=c(1,2))
boxplot(Peso~Fumatrici)
boxplot(Peso~Sesso)

par(mfrow=c(1,2))
boxplot(Peso~Ospedale)
boxplot(Peso~Tipo.parto)

Dai boxplot possiamo notare diversi outlier per tutte le modalità assunte dalle variabili qualitative. Inoltre, il peso dei neonati da madri non fumatrici sembra essere mediamente più alto rispetto a quello da madri fumatrici, come anche il peso dei neonati maschi. Infine, le variabili Ospedale e Tipo.parto non sembrano avere una modalità particolarmente influente rispetto alle altre, quindi possiamo affermare che la variabilità del Peso dei neonati non sembra essere spiegata né dall’ospedale di nascita né dalla tipologia di parto. Tuttavia, non abbiamo informazioni sufficienti per stabilire come sono distribuite le tipologie di parto tra gli ospedali, infatti potrebbe esserci un ospedale che esegue un numero maggiore di parti cesarei.

2.1.1 Test 1

Per quantificare una possibile dipendenza tra le variabili Ospedale e Tipo.parto, andiamo a saggiare l’ipotesi che in alcuni ospedali si fanno più parti cesarei. Per osservare graficamente questa eventuale associazioni utilizziamo il balloonplot.

tabella_test1 = table(Ospedale, Tipo.parto)
ggballoonplot(data=as.data.frame(tabella_test1), fill="blue")

Possiamo osservare che i parti naturali siano più frequenti rispetto a quelli cesarei, ma che non c’è una sostanziale differenze tra gli ospedali nel campione. Questo, quindi, suggerisce una indipendenza tra Ospedale e Tipo.parto, che può essere quantificata calcolando la statistica Chi-quadro di Pearson.

pvalue_digits = function(x){
  value = ifelse(x<0.0001,"0.0001",x)
  return(value)
}

chi_square = chisq.test(tabella_test1)

tabella_test1_viz = data.frame(
  Statistica = chi_square$statistic,
  DF = chi_square$parameter,
  'p-value' = pvalue_digits(chi_square$p.value)
)
kable(tabella_test1_viz, format = "html", digits = 2) %>%
  kable_styling(full_width = T)
Statistica DF p.value
X-squared 1.08 2 0.58

La statistica Chi-quadro ha valore di 1.08 con un numero di gradi di libertà uguale a 2 e con un p-value di 0.58, pertanto non ci sono sufficienti evidenze per rifiutiare l’ipotesi nulla di indipendenza tra gli ospedali e il numero di parti cesarei. Grafichiamo la statistica test.

plot(density(rchisq(1000000, chi_square$parameter)), main = '')
abline(v=qchisq(0.95, chi_square$parameter), col=2)
points(chi_square$statistic, 0, cex=3, col=4, pch=20)

La statistrica ricade nella regione di accettazione.

2.1.2 Test 2

Come secondo test viene saggiata l’ipotesi che le medie del peso e della lunghezza di questo campione di neonati sono significativamente uguali a quelle della popolazione.

Secondo l’Organizzazione Mondiale della Sanità, il peso medio di un neonato regolare (non pretermine) è di circa 3200 g, mentre la lunghezza media è di 500 mm. Eseguiamo, allora, il test statistico t di Student sulle osservazioni relative ai parti regolari, avendo osservato che il 94% delle osservazioni sono di questo tipo.

mu0 = 3200
t_test_peso = t.test(Peso[dati$Tipo.gestazione == "regolare"],
       mu = mu0,
       conf.level = 0.95,
       alternative = "two.sided")

tabella_test_peso = data.frame(
  Statistica = t_test_peso$statistic,
  DF = t_test_peso$parameter,
  'p-value' = pvalue_digits(t_test_peso$p.value),
  Mean_x = t_test_peso$estimate[1],
  CI_lower = t_test_peso$conf.int[1],
  CI_upper = t_test_peso$conf.int[2]
)
kable(tabella_test_peso, format = "html", digits = 2) %>%
  kable_styling(full_width = T)
Statistica DF p.value Mean_x CI_lower CI_upper
t 16.29 2335 0.0001 3349.56 3331.55 3367.56

Per la variabile Peso, il p-value risulta inferiore a 1%, con un livello di singnificatività del 95% e con un intervallo di confidenza compreso fra 3331g e 3367g che non contiene la media della popolazione.

mu0 = 500
t_test_lun = t.test(Lunghezza,
       mu = mu0,
       conf.level = 0.95,
       alternative = "two.sided")

tabella_test_lun = data.frame(
  Statistica = t_test_lun$statistic,
  DF = t_test_lun$parameter,
  'p-value' = pvalue_digits(t_test_lun$p.value<0.0001),
  Mean_x = t_test_lun$estimate[1],
  CI_lower = t_test_lun$conf.int[1],
  CI_upper = t_test_lun$conf.int[2]
)
kable(tabella_test_lun, format = "html", digits = 2) %>%
  kable_styling(full_width = T)
Statistica DF p.value Mean_x CI_lower CI_upper
t -10.07 2497 TRUE 494.7 493.66 495.73

Anche nel caso della varabile Lunghezza, il p-value è inferiore a 1% con un livello di singnificatività del 95% e con un intervallo di confidenza compreso fra 497mm e 499mm che non contiene il valore di riferimento.

Pertanto rifiutiamo l’ipotesi che le medie sono significativamente uguali a quelle della popolazione.

2.1.3 Test 3

Saggiamo adesso l’ipotesi che le misure antropometriche (Peso, Lunghezza, Cranio) sono significativamente diverse tra i due sessi. Visualizziamo graficamente la variabilità delle misure antropometriche rispetto alle modalità assunte da Sesso tramite il boxplot condizionato.

par(mfrow=c(1,3))
boxplot(Peso~Sesso)
boxplot(Lunghezza~Sesso)
boxplot(Cranio~Sesso)

Tutti i boxplot mostrano misure antropometriche mediamente più alte per il sesso maschile. Eseguiamo adesso il test statistico t di Student per verificare una eventuale dipendenda dal sesso.

t_peso_sesso = t.test(Peso ~ Sesso)
t_lun_sesso = t.test(Lunghezza ~ Sesso)
t_cranio_sesso = t.test(Cranio ~ Sesso)

tabella_tests = data.frame(
  Misura = c("Peso","Lunghezza","Diametro Cranio"),
  Statistica = c(t_peso_sesso$statistic, t_lun_sesso$statistic, t_cranio_sesso$statistic),
  p_value = c(pvalue_digits(t_peso_sesso$p.value), pvalue_digits(t_lun_sesso$p.value), pvalue_digits(t_cranio_sesso$p.value)),
  CI_lower = c(t_peso_sesso$conf.int[1], t_lun_sesso$conf.int[1], t_cranio_sesso$conf.int[1]),
  CI_upper = c(t_peso_sesso$conf.int[2], t_lun_sesso$conf.int[2], t_cranio_sesso$conf.int[2])
)
kable(tabella_tests, format = "html", digits = 2) %>%
  kable_styling(full_width = T)
Misura Statistica p_value CI_lower CI_upper
Peso -12.11 0.0001 -287.48 -207.38
Lunghezza -9.58 0.0001 -11.94 -7.88
Diametro Cranio -7.44 0.0001 -6.11 -3.56

I test condotti per tutte le misure antropometriche riportano un p-value molto inferiore al 1%, quindi rifiutiamo l’ipotesi nulla di indipendenza tra le variabili e il sesso dei neonati.

2.2 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.

Rimuoviamo la variabile Tipo.gestazione, che risulta ridondante per la costruzione del modello, e verifichiamo innanzitutto l’ipotesi di normalità della variabile target Peso tramite rappresentazione grafica e il test di Shapiro-Wilk.

dati = dati[, -11]
ggplot(dati, aes(x = Peso)) +
  geom_density(color = NaN, fill = "orange", alpha = 0.3) +
  labs(title = "",
       x = "Peso",
       y = "Densità") +
  theme_minimal()

shapiro_test_peso = shapiro.test(Peso)
tabella_shapiro_test = data.frame(
  Misura = "Peso",
  Statistica = shapiro_test_peso$statistic,
  p_value = pvalue_digits(shapiro_test_peso$p.value)
)
kable(tabella_shapiro_test, format = "html", digits = 2) %>%
  kable_styling(full_width = T)
Misura Statistica p_value
W Peso 0.97 0.0001

Purtroppo, dovremmo rifiutare l’ipotesi di normalità della distribuzione del peso, tuttavia proviamo comunque a considerare adeguato il modello di regressione lineare multipla per questo studio. Visualizziamo la matrice di covarianza tramite una particolare rappresentazione grafica delle misure esplicative quantitative.

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

pairs(dati[,quantitative_columns], upper.panel = panel.smooth, lower.panel = panel.cor)

Dalla riga corrispondente alla variabile target Peso si possono osservare le correlzioni con tutte le variabili esplicative, e possiamo vedere che le correlazioni più significative sono con le variabili Gestazione, Lunghezza e Cranio, le quali sembrano avere una forte correlazione positiva, come visto in precedenza. A seguire Anni.madre e N.gravidanze che non sembrano spiegarne pienamente la variabilità.

Le variabili qualitative possono essere esplorate tramite boxplot, come fatto nelle sezioni precedenti. Abbiamo già saggiato l’ipotesi di dipendenza della variabile Peso dal Sesso, quindi saggiamo le ipotesi di dipendenza con Fumatrici, Tipo.parto e Ospedale.

t_fumatrici = t.test(Peso ~ Fumatrici)
t_tipo_parto = t.test(Peso ~ Tipo.parto)

tabella_tests = data.frame(
  Misura = c("Peso-Fumatrici","Peso-Tipo.parto"),
  Statistica = c(t_fumatrici$statistic, t_tipo_parto$statistic),
  p_value = c(pvalue_digits(t_fumatrici$p.value), pvalue_digits(t_tipo_parto$p.value))
)
kable(tabella_tests, format = "html", digits = 2) %>%
  kable_styling(full_width = T)
Misura Statistica p_value
Peso-Fumatrici 1.04 0.30
Peso-Tipo.parto -0.14 0.89 
t_ospedali = pairwise.t.test(Peso, Ospedale, paired = F, pool.sd = T, p.adjust.method = 'bonferroni')
kable(as.data.frame(t_ospedali$p.value), format = "html", digits = 2) %>%
  kable_styling(full_width = T)
osp1 osp2
osp2 1.00 NA
osp3 0.33 0.32

Dai p-value ottenuti, possiamo affermare che non ci sono abbastanza evidenze per rifiutare l’ipotesi nulla di indipendenza del peso dalle madri fumatrici, dal tipo di parto e dal tipo di ospedale. Vedremo durante la costruzione del modello se questi regressori spiegano la variabilità.

Definiamo il primo modello di regressione lineare multipla, inserendo tutte le variabili esplicative in forma additiva.

mod1 = lm(Peso ~ . , data=dati)
summary(mod1)
## 
## 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 *  
## Fumatrici1      -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

Il modello mostra una variabilità spiegata di 0.73, che rappresenta una buona base di partenza. Osservando i p-value ottenti possiamo affermare che le variabili Gestazione, Lunghezza e Cranio mostrano stime positive molto significative, ovvero per ogni settimana in più di gestazione aumenta il peso di 32g circa, oppure per ogni mm in più di lunghezza il peso aumenta di 10g circa, come anche per il diamentro del Cranio. Meno significativa sembra essere la stima delle variabili N.gravidanze e Anni.madre.

Fissate tutte le variabili quantitative, possiamo inoltre osservare che il peso dei neonati da madri fumatrici è circa 30g in meno rispetto alle madri non fumatrici, come anche un parto naturale determina un incremento di 30g circa rispetto ad un parto cesareo. Ancora, il peso dei neonati maschi è mediamente 77g in più delle femmine, mentre i neonati presso l’Ospedale 2 hanno un peso mediamente di 11g in meno.

Per migliorare il modello cominciamo con il rimuovere Anni.madre, in quanto contribuisce meno a spiegare la variabilità.

mod2 = update(mod1, ~. -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 ** 
## Fumatrici1      -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

Il modello mantiene lo stesso livello di variabilità spiegata di 0.73, e le variabili Gestazione, Lunghezza, Cranio si confermano stimatori significativi positivi. Tale risultato conferma che il primo regressore rimosso non contribuiva a spiegare la variabilità del peso dei neonati. Rimuoviamo la variabile Ospedale, avendo un p-value superiore ad 1%.

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 ** 
## Fumatrici1      -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

Il livello di variabilità spiegata resta invariato a 0.73, quindi anche il regressore Ospedale non contribuisce significativamente a spiegare la variabilità del peso, confermando il risultato del test sulla indipendenza. Procediamo con lo step successivo, rimuovendo il regressore Fumatrici, per il quale il test suggeriva, come per gli ospedali, l’assenza di differenze significative nel peso dei neonati da madri fumatrici e non 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

Il modello si mantiene stabile su un R-quadro adjusted di 0.73, quindi confermndo il risultato del test, nel campione preso in esame l’essere fumatrice non determina una variazione del peso significativo nei neonati. Rimuoviamo adesso il regressore Tipo.parto, poiché risulta il regressore con il p-value maggiore.

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

Il modello mantiene una variabilità spiegata di 0.73, quindi anche in questo caso confermiamo una indipendenza del peso dal tipo di parto. Rimuoviamo adesso il regressore N.gravideze.

mod6 = update(mod5, ~. -N.gravidanze)
summary(mod6)
## 
## 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

Anche il numero di gravidanze sembra non essere significativo per spiegare la variabilità del peso, infatti il modello mantiene un R-quadro di 0.73.

Dovendo costruire un modello applicabile per l’intera popolazione, ritengo che la variabile N.gravidanze debba essere inclusa nel modello, poiché, non peggiorando la variabilità spiegata, potrebbe risultare nuovamente significativa nel caso di un altro campione.

Verifichiamo la presenza di collinearità tra le variabili rimaste, tramite il calcolo del VIF.

kable(vif(mod5), format = 'html', digits = 2) %>%
  kable_styling(full_width = T)
x
N.gravidanze 1.02
Gestazione 1.67
Lunghezza 2.08
Cranio 1.62
Sesso 1.04

I valori del VIF sono tutti inferiori a 5, quindi possiamo escludere collinearità tra i regressori del modello. Infine, cerchiamo eventuali effetti di interazione o quadratici tra i regressori, come sembrano mostrare gli scatterplot di Lunghezza, diametro Cranio e Gestazione nella trasposizione grafica della matrice di covarianza.

mod6 = update(mod5,~.+I(Lunghezza^2))
summary(mod6)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Lunghezza^2), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1169.62  -181.77   -12.79   163.77  1786.03 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    212.288548 723.852095   0.293 0.769336    
## N.gravidanze    14.085464   4.266175   3.302 0.000975 ***
## Gestazione      42.551398   3.876629  10.976  < 2e-16 ***
## Lunghezza      -20.267001   3.162718  -6.408 1.76e-10 ***
## Cranio          10.651783   0.418894  25.428  < 2e-16 ***
## SessoM          69.968733  11.038797   6.338 2.75e-10 ***
## I(Lunghezza^2)   0.031655   0.003267   9.690  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.7 on 2491 degrees of freedom
## Multiple R-squared:  0.7369, Adjusted R-squared:  0.7363 
## F-statistic:  1163 on 6 and 2491 DF,  p-value: < 2.2e-16

L’effetto quadratico del regressore Lunghezza risulta singificativo, e la variabilità spiegata del modello è aumentata a 0.74.

Verifichiamo la presenza dell’effetto quadratico del regerssore Cranio.

mod7 = update(mod5,~.+I(Cranio^2))
summary(mod7)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Cranio^2), data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.6  -179.4   -14.8   163.4  2622.6 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    84.10118 1151.77280   0.073  0.94180    
## N.gravidanze   12.76356    4.31259   2.960  0.00311 ** 
## Gestazione     38.90540    3.93291   9.892  < 2e-16 ***
## Lunghezza      10.48745    0.30157  34.776  < 2e-16 ***
## Cranio        -31.79371    7.16973  -4.434 9.63e-06 ***
## SessoM         73.10236   11.16590   6.547 7.11e-11 ***
## I(Cranio^2)     0.06262    0.01059   5.915 3.77e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 272.8 on 2491 degrees of freedom
## Multiple R-squared:  0.7308, Adjusted R-squared:  0.7301 
## F-statistic:  1127 on 6 and 2491 DF,  p-value: < 2.2e-16

Si riduce il valore di Adjusted R-squared a 0.73, quindi scartiamo questo effetto e verifichiamo quello sulla Gestazione.

mod8 = update(mod5,~.+I(Gestazione^2))
summary(mod8)
## 
## 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

Anche in questo caso non registriamo un miglioramento un peggioramento della spiegabilità della varianza rispetto al modello 5 (senza effetti quadratici), quindi escludiamo anche effetti quadratici dovuti al regressore Gestazione. Verifichiamo, infine, effetti di interazione tra Lunghezza e Cranio.

mod9 = update(mod6,~.+Lunghezza*Cranio)
summary(mod9)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Lunghezza^2) + Lunghezza:Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1176.59  -179.44   -11.47   166.18  1306.63 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.263e+03  1.003e+03  -2.256 0.024128 *  
## N.gravidanze      1.424e+01  4.256e+00   3.345 0.000835 ***
## Gestazione        4.046e+01  3.912e+00  10.343  < 2e-16 ***
## Lunghezza        -2.136e+01  3.170e+00  -6.738 1.98e-11 ***
## Cranio            2.740e+01  4.727e+00   5.796 7.66e-09 ***
## SessoM            7.183e+01  1.103e+01   6.515 8.76e-11 ***
## I(Lunghezza^2)    4.474e-02  4.916e-03   9.102  < 2e-16 ***
## Lunghezza:Cranio -3.447e-02  9.692e-03  -3.556 0.000383 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.1 on 2490 degrees of freedom
## Multiple R-squared:  0.7383, Adjusted R-squared:  0.7375 
## F-statistic:  1003 on 7 and 2490 DF,  p-value: < 2.2e-16

L’effetto di interazione non migliora il coefficiente Adjusted R-squared, quindi per il principio di semplicità del modello, decidiamo di scartarlo.

2.3 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. Includeremo nella selezione anche i modelli contenenti l’effetto quadratico della Lunghezza.

kable(AIC(mod1,mod2,mod3,mod4,mod5,mod6),
      format = 'html') %>%
  kable_styling()
df AIC
mod1 12 35145.57
mod2 11 35144.06
mod3 9 35149.33
mod4 8 35148.67
mod5 7 35152.89
mod6 8 35062.46 
kable(BIC(mod1,mod2,mod3,mod4,mod5,mod6),
      format = 'html') %>%
  kable_styling()
df BIC
mod1 12 35215.45
mod2 11 35208.12
mod3 9 35201.73
mod4 8 35195.25
mod5 7 35193.65
mod6 8 35109.04

Sulla base dei valori ottenuti dai criteri AIC e BIC, scegliamo il modello con il valore BIC più basso, che favorisce i modelli più semplici e meno complessi, ovvero mod6:

summary(mod6)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Lunghezza^2), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1169.62  -181.77   -12.79   163.77  1786.03 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    212.288548 723.852095   0.293 0.769336    
## N.gravidanze    14.085464   4.266175   3.302 0.000975 ***
## Gestazione      42.551398   3.876629  10.976  < 2e-16 ***
## Lunghezza      -20.267001   3.162718  -6.408 1.76e-10 ***
## Cranio          10.651783   0.418894  25.428  < 2e-16 ***
## SessoM          69.968733  11.038797   6.338 2.75e-10 ***
## I(Lunghezza^2)   0.031655   0.003267   9.690  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269.7 on 2491 degrees of freedom
## Multiple R-squared:  0.7369, Adjusted R-squared:  0.7363 
## F-statistic:  1163 on 6 and 2491 DF,  p-value: < 2.2e-16

Possiamo trarre le seguenti considerazioni:

  • Il peso dei neonati aumenta di 43g circa per ogni settimana di gestazione in più.
  • Il peso dei neonati maschi è mediamente 70g in più rispetto alle femmine.
  • Il peso dei neonati cresce di 11g per ogni mm in più di diametro del cranio.
  • La lunghezza in valore assoluto ha un effetto maggiore sulla variabilità del peso. In particolare, all’aumentare delle lunghezza il peso inizialmente diminuisce fino ad un minimo, per poi iniziare ad aumentare. Il punto di minimo del peso è pari a 333g circa.
minimo = -20/(2*0.03)

2.4 Analisi della Qualità del Modello

Valutiamo, adesso, la capacità predittiva del modello utilizzando metriche come R² e il Root Mean Squared Error (RMSE). Per quanto riguarda R² abbiamo ottenuto dal modello finale un valore di 0.74. Procediamo con il calcolo di RMSE.

predetti = predict(mod6, dati)
osservati = dati
rmse = sqrt(mean((osservati$Peso - predetti)^2))
cat("RMSE:", round(rmse, 2))
## RMSE: 269.34

Le previsioni del nostro modello, quindi, si discostano di circa 269g dal valore reale del peso.

Effettuiamo adesso l’analisi dei residui e valutiamo la presenza di valori influenti, che potrebbero distorcere le previsioni, come osservazioni leverage o outlier.

Le assunzioni da verificare sono:

  1. Linearità
  2. Normalità
  3. Omoschedasticità
  4. Indipendenza
par(mfrow=c(2,2))
plot(mod6)

Graficando i residui del modello, osservare che:

  1. Dal grafico Residuals vs Fitted i punti si posizionano intorno alla media di zero, anche se non perfettamente schiacciati, soprattutto le osservazioni 155, 1549 e 1305.

  2. Dal grafico Q-Q Residuals, dove i residui vengono confronati con i quantili di una distribuzione normale, i punti giacciono sulla bisettricce, ad eccezione per le code ricurve dovute probabilmente a gestazioni pretermine o più lunghe (stesse osservazioni di prima). Quindi i residui seguono una distribuzione normale.

  3. Dal grafico Scale-Location non osserviamo pattern specifici, anche se la linea della varianza sembra avere una leggera pendenza.

  4. Dal grafico Residuals vs Leverage si visualizzano i potenziali valori influenti, ovvero outlier o leverage. Nessuno punto supera la soglia di avvertimo della linea di Cook a 0.5 e nessuna supera la soglia di allerme di 1, eccetto l’osservazione 1549 che risulta anomala.

Quantifiamo quello che abbiamo visto dai grafici, tramite i test statistici di Shapiro Wilk (normalità dei residui), Breusch-Pagan (omoschedasticità) e Durbin-Watson (indipendenza)

shapirotest.res = shapiro.test(residuals(mod6))
bgtest.res = bgtest(mod6)
dwtest.res = dwtest(mod6)

tabella_tests = data.frame(
  Test = c("Normalità","Omoschedasticità","Indipendenza"),
  Statistica = c(shapirotest.res$statistic, bgtest.res$statistic, dwtest.res$statistic),
  p_value = c(pvalue_digits(shapirotest.res$p.value), pvalue_digits(bgtest.res$p.value), pvalue_digits(dwtest.res$p.value))
)
kable(tabella_tests, format = "html", digits = 2) %>%
  kable_styling(full_width = T)
Test Statistica p_value
W Normalità 0.99 0.0001
LM test Omoschedasticità 1.69 0.193858072465883
DW Indipendenza 1.95 0.0960113993694909

Dai p-value ottenuti, possiamo affermare che i residui non rispetttano l’ipotesi di normalità, rispettano l’ipotesi di varianza costante e non presentano autocorrelazione.

Analizziamo la presenza di eventuali valori influenti.

# leverage
lev = hatvalues(mod6)
plot(lev)
p = sum(lev)
soglia = 2 * p/n
abline(h=soglia, col=2)

leverage = length(lev[lev>soglia])

Abbiamo trovato 142 osservazioni leverage, ovvero osservazioni che si trovano lontano dalla media nello spazio dei regressori. Analizziamo adesso gli outlier:

# ouliers
plot(rstudent(mod6))
abline(h=c(-2,2), col=2)

outliers = outlierTest(mod6)

Otteniamo 5 osservazioni outlier stattisticamente significative (tre delle quali individuate dai quattro grafici del plot del modello), ovvero osservazioni lontane dalla media nello spazio della variabile target. Analizziamo la distanza di Cook

# distanza di cook
cook=cooks.distance(mod6)
plot(cook)

max_cook = c(which.max(cook),max(cook))

Abbiamo trovato l’osservazione 1549 anomale, quindi affetta probabilmente da errori di inserimento.

3. Previsioni e Risultati

Utilizziamo il modello per effettuare la seguente previsione:

Stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana.

Fissiamo i regressori impostando il sesso del neonato in ‘F’ , numero di gravidanze 2 (la terza sarebbe quella per la quale stiamo facendo la previsione), settimana di gestazione 39.

new_data = data.frame(
  Sesso = 'F',
  N.gravidanze = 2,
  Gestazione = 39,
  Lunghezza = mean(Lunghezza),
  Cranio = mean(Cranio)
)
previsione = predict(mod6, newdata = new_data)

La previsione è di 3243g.

4. Visualizzazioni

Abbiamo a visualizzare l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.

ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col=Fumatrici), position = 'jitter')+
  geom_smooth(aes(x=Gestazione,
                  y=Peso,
                  col=Fumatrici), se=F, method = 'lm')+
  geom_smooth(aes(x=Gestazione,
                  y=Peso), col='black', se=F, method = 'lm')+
  theme_minimal()

La relazione tra fumatrici e non fumatrici è molto simile e crescente con l’aumentare delle settimane di gestazione, anche se quella delle non fumatrici ha una pendenza maggiore, perché hanno neonati con peso maggiore. Dal grafico seguente, invece, confrontiamo i valori di peso osservati con quelli predetti, aspettandoci come situazione ideale una distribuzione di punti lungo la bisettrice.

dati$Peso.previsto = predict(mod6, newdata=dati)
ggplot(data=dati)+
  geom_point(aes(x=Peso,y=Peso.previsto), position = 'jitter')+
  geom_abline(intercept = 0, slope = 1, col = 'orange')+
  theme_minimal()

Il modello di previsione rispecchia bene i dati osservati. Nelle code abbiamo un leggero disallineamento, in particolare per valori bassi del peso il modello sembra sovrastimare la variabile, mentre per valori alti del peso il modello sembra sottostimarlo.

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.

Le evidenze trovate sono:

  1. Indipendenza tra tipologia di parto e ospedale, ovvero le strutture ospedalieri sembrano lavorare con lo stesso livello di offerta ed efficacia.
  2. Le medie campionarie di lunghezza e diametro del cranio sono diverse da quelle dichiarate dall’OMS, quindi potrebbe essere utile considerare valori di riferimento più localizzati rispetto alla regione di prelevamento del campione.
  3. Le misure antropometrice sono dipendenti dal sesso del neonato, e sembrano maggiori per i maschi.
  4. L’essere fumatrici sembra non influenzare il peso del neonato.
  5. Il modello trovato ha una variabilità spiegata dello 0.74, con un errore sulle previsioni del peso dell’8%.
  6. Il modello sembra non stimare bene i valori di peso nelle code, quindi in caso di parti pretermine o oltre le 40 settimane.