Modello Statistico per la Previsione del Peso Neonatale

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.

1. Caricamento delle librerie necessarie

library(dplyr)
library(ggplot2)
library(tidyr)
library(car)
library(MASS)
library(htmltools)
library(moments)

2. Caricamento e ispezione dei dati

data <- read.csv("neonati.csv")
str(data)
## 'data.frame':    2500 obs. of  10 variables:
##  $ Anni.madre  : int  26 21 34 28 20 32 26 25 22 23 ...
##  $ N.gravidanze: int  0 2 3 1 0 0 1 0 1 0 ...
##  $ Fumatrici   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Gestazione  : int  42 39 38 41 38 40 39 40 40 41 ...
##  $ Peso        : int  3380 3150 3640 3690 3700 3200 3100 3580 3670 3700 ...
##  $ Lunghezza   : int  490 490 500 515 480 495 480 510 500 510 ...
##  $ Cranio      : int  325 345 375 365 335 340 345 349 335 362 ...
##  $ Tipo.parto  : chr  "Nat" "Nat" "Nat" "Nat" ...
##  $ Ospedale    : chr  "osp3" "osp1" "osp2" "osp2" ...
##  $ Sesso       : chr  "M" "F" "M" "M" ...
head(data)

Il dataset comprende dati su 2500 neonati provenienti da tre ospedali. Le variabili raccolte includono:

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.

3. Statistiche descrittive e visualizzazione dei dati

3.1 Variabiliti quantitative

data_quant <- dplyr::select(data, -Sesso, -Ospedale, -`Tipo.parto`, -Fumatrici)

calculate_statistics <- function(x) {
  c(
    coeff_var = sd(x) / mean(x),
    skewness  = skewness(x),
    kurtosis  = kurtosis(x) - 3  # Normalizzazione indice di curtosi
  )
}

statistics_table <- 
  data_quant%>%
  summarise(across(everything(), calculate_statistics)) %>%
  t() %>%
  as.data.frame()

#rinominare le colonne della tabella
colnames(statistics_table) <- c("coeff_var", "skewness", "kurtosis")
statistics_table <- round(statistics_table, 2)
statistics_table <- t(statistics_table) #inverto righe e colonne

#formattare i valori per rendere più leggibili senza alterare il formato originario
knitr::kable(
  summary(data_quant),
  caption = "Statistiche descrittive",
)
Statistiche descrittive
Anni.madre N.gravidanze Gestazione Peso Lunghezza Cranio
Min. : 0.00 Min. : 0.0000 Min. :25.00 Min. : 830 Min. :310.0 Min. :235
1st Qu.:25.00 1st Qu.: 0.0000 1st Qu.:38.00 1st Qu.:2990 1st Qu.:480.0 1st Qu.:330
Median :28.00 Median : 1.0000 Median :39.00 Median :3300 Median :500.0 Median :340
Mean :28.16 Mean : 0.9812 Mean :38.98 Mean :3284 Mean :494.7 Mean :340
3rd Qu.:32.00 3rd Qu.: 1.0000 3rd Qu.:40.00 3rd Qu.:3620 3rd Qu.:510.0 3rd Qu.:350
Max. :46.00 Max. :12.0000 Max. :43.00 Max. :4930 Max. :565.0 Max. :390 
knitr::kable(
  statistics_table,
  caption = "Statistiche descrittive",
  digits = 2,
  format.args = list(big.mark = ".", decimal.mark = ",")
)
Statistiche descrittive
Anni.madre N.gravidanze Gestazione Peso Lunghezza Cranio
coeff_var 0,19 1,31 0,05 0,16 0,05 0,05
skewness 0,04 2,51 -2,07 -0,65 -1,51 -0,79
kurtosis 0,38 10,99 8,26 2,03 6,49 2,95 
par(mfrow = c(2, 3))
invisible(lapply(names(data_quant), function(var) {
  x <- data_quant[[var]]
  dens <- density(x, na.rm = TRUE)  
  
  h <- hist(x, plot = FALSE)
  
  y_max <- max(c(h$density, dens$y)) * 1.05
  
  hist(x,
       main = var, xlab = var,
       col = "lightblue", border = "white",
       prob = TRUE,
       xlim = range(dens$x),
       ylim = c(0, y_max))
  lines(dens, col = "darkgreen", lwd = 2)
}))

Tutte le variabili hanno una distribuzione leptocurtica, con possibili valori di outlier molto distanti dalla media. Questo è immediatamente visibile dalla forma delle distribuzioni di alcune variaibili.

Anni.madre: La distribuzione è quasi simmetrica con età media di 28 anni circa. Il range va da 0 a 46 anni, con la maggior parte delle madri concentrate tra 25-32 anni di età. Un età pari a zero non è compatibile con lo status di madre, pertanto si andrà ad apportare un aggiustamento al dataset. La variabilità rispetto alla media è bassa.

N.gravidanze: ha variabilità molto alta, forte asimmetria positiva e un indice di curtosi molto alto indice di possibili outliers molto distanti dalla media, infatti il numero massimo di gravidanze è pari a 12, molto distante dalla mediana pari a 1. La variabilità rispetto alla media è piuttosto elevata.

Gestazione: asimmetria negativa e curtosi elevate. La media e la mediana sono a pari a 39 settimane e il minimo è di 25 settimane con un massimo di 43. I valori sembrano avere senso rispetto alla fisiologica durata della gravidanza (di nove mesi, 40 settimane) e indicano la presenza di nascite premature (intorno ai 6 mesi) come evidente dall’osservazione grafica. La variabilità rispetto alla media è bassa.

Peso: distribuzione quasi simmetrica, non si rilevano valori anomali dalla tabella delle statistiche. Il peso minimo sembra essere coerente rispetto ai parti con settimane di gestazione inferiori, infatti si nota una coda lunga a sinistra verso valori più bassi. La variabilità rispetto alla media è bassa.

Lunghezza: asimmetria negativa e curtosi abbastanza elevata. La variabilità rispetto alla media è bassa.

Cranio: lieve asimmetria negativa, la variabilità rispetto alla media è bassa.

Aggiustamento dataset:

controlliamo i valori anomali dell’età della madre e verifichiamo anche se vi sono colonne con valori mancanti

data[data$Anni.madre < 13,]
colSums(is.na(data))
##   Anni.madre N.gravidanze    Fumatrici   Gestazione         Peso    Lunghezza 
##            0            0            0            0            0            0 
##       Cranio   Tipo.parto     Ospedale        Sesso 
##            0            0            0            0

Eliminiamo dal dataset i valori anomali di Anni.Madre.

Dalla verifica non sono emersi invce valori mancanti nel dataset.

data_orig <- data  # backup
data <- data[data$Anni.madre >= 13, ]

3.2 Variabili qualitative

data_qual <- dplyr::select(data, Sesso, Ospedale, `Tipo.parto`, Fumatrici)

# Funzione per indice di Gini
gini_index <- function(x){
  freq_abs = table(x)
  freq_rel = freq_abs / length(x)
  freq_rel_squared = freq_rel^2
  J = length(freq_abs)
  gini = 1 - sum(freq_rel_squared)
  gini_norm = gini / ((J - 1) / J)
  return(round(gini_norm, 3))
}

# Funzione per tabella frequenze (senza Gini/Moda)
frequency <- function(x) {
  freq_abs <- table(x)
  freq_rel <- prop.table(freq_abs)
  
  df <- data.frame(
    Categoria = names(freq_abs),
    `Frequenza assoluta` = as.integer(freq_abs),
    `Frequenza relativa` = round(as.numeric(freq_rel), 2),
    check.names = FALSE
  )
  return(df)
}

# Lista variabili qualitative
variabili_qualitative <- names(data_qual)

# Ciclo per stampa tabella + commento con Gini e Moda
for (variabile in variabili_qualitative) {
  x <- data_qual[[variabile]]
  tabella <- frequency(x)
  gini <- gini_index(x)
  moda <- names(table(x))[which.max(table(x))]

  print(knitr::kable(
    tabella,
    align = "r",
    digits = 2,
    format.args = list(big.mark = ".", decimal.mark = ","),
    caption = paste("Distribuzione di", variabile)
  ))

    cat(
    sprintf('<p style="font-style: italic; color: #444; margin-top: -10px;">
              L’indice di Gini per <strong>%s</strong> è pari a <strong>%.3f</strong> 
              e la moda è <strong>%s</strong>.
            </p>', variabile, gini, moda)
  )
}
Distribuzione di Sesso
Categoria Frequenza assoluta Frequenza relativa
F 1.255 0,5
M 1.243 0,5

L’indice di Gini per Sesso è pari a 1.000 e la moda è F.

Distribuzione di Ospedale
Categoria Frequenza assoluta Frequenza relativa
osp1 816 0,33
osp2 848 0,34
osp3 834 0,33

L’indice di Gini per Ospedale è pari a 1.000 e la moda è osp2.

Distribuzione di Tipo.parto
Categoria Frequenza assoluta Frequenza relativa
Ces 728 0,29
Nat 1.770 0,71

L’indice di Gini per Tipo.parto è pari a 0.826 e la moda è Nat.

Distribuzione di Fumatrici
Categoria Frequenza assoluta Frequenza relativa
0 2.394 0,96
1 104 0,04

L’indice di Gini per Fumatrici è pari a 0.160 e la moda è 0. 

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))  

invisible(lapply(names(data_qual), function(var) {
  barplot(table(data_qual[[var]]),
          main = var,
          col = "lightblue",
          border = "white",
          las = 2)  
}))

  • Il sesso dei bambini è equamente distribuito tra maschi e femmine, come confermato dall’indice di Gini;
  • I dati sono equamenente distribuiti anche rispetto all’ospedale di nascita, infatti l’indice di Gini suggerisce un’alta eterogeneità;
  • Il tipo di parto è prevalentemente naturale, ma la numerosità di parti naturali e cesarei è abbastanza bilanciata nel dataset (alto indice di Gini);
  • Vi sono quasi esclusivamente madri non fumatrici nel campione, infatti l’indice di Gini è basso, indice che il comportamento rispetto al fumo è omogeneo, infatti la maggior parte delle madri sono non fumatrici.

4. Test di ipotesi

Test 1: in alcuni ospedali si fanno più parti cesarei

Utilizziamo il test del chi-quadro su tabelle di contingenza per verificare se esiste una relazione significativa tra le variaibli qualitative Ospedale e Tipo.Parto

# Tabella di contingenza Ospedale e Tipo.Parto

test1 <- prop.table(table(data_qual$Ospedale, data_qual$Tipo.parto), margin=1)
options(digits = 2)
test1
##       
##         Ces  Nat
##   osp1 0.30 0.70
##   osp2 0.30 0.70
##   osp3 0.28 0.72
mosaicplot(
  table(data$Ospedale, data$Tipo.parto),
  main = "Distribuzione dei tipi di parto per ospedale",
  color = TRUE, 
  xlab = "Ospedale",
  ylab = "Tipo di parto"
)

La tabella di contingenza e l’evidenza grafica suggeriscono che non c’è una forte evidenza di differenze tra ospedale e tipo di parto.

chisq.test(test1)
## 
##  Pearson's Chi-squared test
## 
## data:  test1
## X-squared = 0.001, df = 2, p-value = 1

Sotto l’ipotesi nulla le due variabili sono indipendenti.

Il valore del test è vicino allo zero e indica che le frequenze osservate sono quasi identiche a quelle attese sotto l’ipotesi di indipendenza.

Il p-value ci conferma che che c’è un’evidenza statistica sufficiente per non rifiutare l’ipotesi nulla.

qchisq(0.95, df = 2)
## [1] 6

Il risultato della statistica è inferiore al valore critico per 2 gradi di libertà e un livello di significatività del 5%.

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

Utilizziamo il test t di student a due code dove l’ipotesi nulla H0 sostiene che la media del campione è uguale alla media della popolazione. Valuteremo se vi è una differenza significativa nel peso medio e nella lunghezza media dei neonati rispetto a quelli della popolazione con un livello di significatività del 5%.

Dalle statistiche descrittive abbiamo visto che il peso ha una media di 3284 grammi e la lunghezza una media di 49,47 cm.

La media teorica sotto H0 (derivata dalla letteratura) è di 3300 grammi di peso e 50 cm di lunghezza

(fonte: +https://www.ospedalebambinogesu.it/da-0-a-30-giorni-come-si-presenta-e-come-cresce-80012/#)

Test per la media del peso

Test Shapiro-Wilk: saggiare la normalità della variabile campionaria Peso

Il test t di Student che applicheremo per confrontare la media del nostro campione con quella della popolazione, richiede di verificare che i dati osservati siano distribuiti normalmente.

shapiro.test(data$Peso)
## 
##  Shapiro-Wilk normality test
## 
## data:  data$Peso
## W = 1, p-value <2e-16
qqnorm(data$Peso)
qqline(data$Peso, col = "steelblue", lwd = 2)

Il valore della statistica è prossimo a 1 è il valore del p-value è molto basso, inferiore al livello di significatività del 5% pertanto si rifiuta l’ipotesi nulla di distribuzione normale del campione di dati osservato. Il Q-Q Plot conferma una deviazione dalla normalità.

Nonostante il test di Saphiro-Wilk non ci permetta di accettare l’ipotesi nulla, graficamente possiamo evincere che i punti seguono approssimativamente la linea diagonale nella parte centrale della distribuzione (tra -2 e +2) deviando prevalentemente per valori bassi, inoltre considerata la numerosità campionaria (2500 osservazioni) possiamo considerare l’assunzione di normalità approssimativamente soddisfatta.

t.test(data$Peso, mu=3300)
## 
##  One Sample t-test
## 
## data:  data$Peso
## t = -2, df = 2497, p-value = 0.1
## alternative hypothesis: true mean is not equal to 3300
## 95 percent confidence interval:
##  3264 3305
## sample estimates:
## mean of x 
##      3284

Possiamo concludere che la media campionaria del peso dei neonati non è significativamente diversa da quella della popolazione. Il valore medio osservato cade all’interno dell’intervallo di confidenza.

Test per la media della lunghezza

Adoperiamo lo stesso processo di analisi: saggiamo la normalità della distribuzione campionaria della variabile lunghezza ed effettuiamo il test t di student per confrontare la media campionaria rispetto a quella della popolazione.

shapiro.test(data$Lunghezza)
## 
##  Shapiro-Wilk normality test
## 
## data:  data$Lunghezza
## W = 0.9, p-value <2e-16
qqnorm(data$Lunghezza)
qqline(data$Lunghezza, col = "steelblue", lwd = 2)

Il valore della statistica è prossimo a 1 è il valore del p-value è molto basso, inferiore al livello di significatività del 5% pertanto si rifiuta l’ipotesi nulla di distribuzione normale del campione di dati osservato.

Anche nel caso della lunghezza risulta evidente dal Q-Qplot che la distribuzione si discosta da una distribuzione normale sebbene i punti non aderiscano alla retta in particolar modo per valori bassi, pertanto, anche a fronte della numerosità campionaria, consideriamo l’assunzione di normalità valida e procediamo con il test t di Student.

t.test(data$Lunghezza, mu=500)
## 
##  One Sample t-test
## 
## data:  data$Lunghezza
## t = -10, df = 2497, p-value <2e-16
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  494 496
## sample estimates:
## mean of x 
##       495

Il risultato del test mostra un valore t molto elevato e un p-value molto basso pertanto la media campionaria osservata è significativamente diversa dal valore medio della popolazione.

La media della lunghezza dei neonati è significativamente diversa dalla media attesa. L’intervallo di confidenza per la lunghezza, infatti, non include il valore atteso di 50 cm.

Test 3: Le misure antropometriche sono significativamente diverse tra i due sessi

Utilizziamo il test t di student dove l’ipotesi nulla H0 sostiene che la media del peso e la media della lunghezza non sono significativamente differente tra maschi e femmine, ad un livello di significatività del 5%.

I requisiti che devono essere verificati per questo test sono:

- la distribuzione normale dei due campioni: abbiamo già indagato la forma della distribuzione per peso e lunghezza, concludendo che la distribuzione è approssimativamente normale posto che si nota una deviazione dalla normalità nelle code; elaboriamo lo stesso grafico anche per Cranio.

qqnorm(data$Cranio)
qqline(data$Cranio, col = "steelblue", lwd = 2)

Anche per Cranio deriviamo la stessa conclusione, inoltre, considerata la numerosità campionaria, possiamo nel complesso considerare valida l’ipotesi di normalità per le variabili antropometriche;

- l’indipendenza dei due campioni: ogni neonato è o maschio o femmina come naturalmente definito alla nascita.

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

variabili <- c("Peso", "Lunghezza", "Cranio")

for (var in variabili) {
  # Dati per femmine e maschi
  x_f <- data[[var]][data$Sesso == "F"]
  x_m <- data[[var]][data$Sesso == "M"]
  
  # Tabella riassuntiva
  summary_f <- summary(x_f)
  summary_m <- summary(x_m)
  
  tabella <- data.frame(
    Statistica = names(summary_f),
    Femmine = round(as.numeric(summary_f), 2),
    Maschi  = round(as.numeric(summary_m), 2)
  )
  
  # Stampa tabella
  print(
    knitr::kable(
      tabella,
      caption = paste("Statistiche descrittive di", var, "per sesso"),
      digits = 2,
      format.args = list(decimal.mark = ",", big.mark = ".")
    )
  )
  
  # Calcolo CV e IQR
  cv_f <- round(sd(x_f) / mean(x_f), 3)
  cv_m <- round(sd(x_m) / mean(x_m), 3)
  iqr_f <- round(IQR(x_f), 1)
  iqr_m <- round(IQR(x_m), 1)
  
  # Stampa testo
  cat(paste0(
    "\n",
    "CV femmine: ", cv_f, " | maschi: ", cv_m, "\n",
    "IQR femmine: ", iqr_f, " | maschi: ", iqr_m, "\n\n"
  ))
}
## 
## 
## Table: Statistiche descrittive di Peso per sesso
## 
## |Statistica | Femmine| Maschi|
## |:----------|-------:|------:|
## |Min.       |     830|    980|
## |1st Qu.    |   2.900|  3.150|
## |Median     |   3.160|  3.430|
## |Mean       |   3.161|  3.408|
## |3rd Qu.    |   3.470|  3.720|
## |Max.       |   4.930|  4.810|
## 
## CV femmine: 0.167 | maschi: 0.145
## IQR femmine: 570 | maschi: 570
## 
## 
## 
## Table: Statistiche descrittive di Lunghezza per sesso
## 
## |Statistica | Femmine| Maschi|
## |:----------|-------:|------:|
## |Min.       |     310|    320|
## |1st Qu.    |     480|    490|
## |Median     |     490|    500|
## |Mean       |     490|    500|
## |3rd Qu.    |     505|    515|
## |Max.       |     565|    560|
## 
## CV femmine: 0.056 | maschi: 0.048
## IQR femmine: 25 | maschi: 25
## 
## 
## 
## Table: Statistiche descrittive di Cranio per sesso
## 
## |Statistica | Femmine| Maschi|
## |:----------|-------:|------:|
## |Min.       |     235|    265|
## |1st Qu.    |     330|    334|
## |Median     |     340|    343|
## |Mean       |     338|    342|
## |3rd Qu.    |     348|    352|
## |Max.       |     390|    390|
## 
## CV femmine: 0.05 | maschi: 0.046
## IQR femmine: 18 | maschi: 18

Mediana e media di peso, lunghezza e diametro del cranio per i maschi sono più alte di quelle delle femmine.

La dispersione della parte centrale della distribuzione è la stessa.

C’è una variabilità relativa bassa per entrambi i sessi e leggermente più alta per le femmine, ma possiamo ritenerla una differenza trascurabile.

Sono presenti numerosi outliers sia tra i maschi che tra le femmine.

Procediamo con il test t:

Effettueremo il test t di Welch (più robusto) che non assume omogeneità della varianza tra gruppi.

Effettueremo anche il test non parametrico Wilcoxon che non risente dei valori anomali ed asimmetrie, infatti esso si basa sulla trasformazione dei valori osservati in ranghi e quindi può essere utilizzato qualora non sia consigliabile effettuare un confronto tra medie in assenza di normalità nelle distribuzioni.

t.test(data=data,
       Peso~Sesso,
       )
## 
##  Welch Two Sample t-test
## 
## data:  Peso by Sesso
## t = -12, df = 2489, p-value <2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287 -207
## sample estimates:
## mean in group F mean in group M 
##            3161            3408
t.test(data=data,
       Lunghezza~Sesso,
       )
## 
##  Welch Two Sample t-test
## 
## data:  Lunghezza by Sesso
## t = -10, df = 2457, p-value <2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -11.9  -7.9
## sample estimates:
## mean in group F mean in group M 
##             490             500
t.test(data=data,
       Cranio~Sesso,
       )
## 
##  Welch Two Sample t-test
## 
## data:  Cranio by Sesso
## t = -7, df = 2489, p-value = 1e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.1 -3.6
## sample estimates:
## mean in group F mean in group M 
##             338             342

il test t ci porta a rifiutare l’ipotesi nulla, pertanto la media del peso, della lunghezza e del diametro del cranio dei due gruppi sono statisticamente diversi.

Eseguiamo il test non-parametrico di Wilcoxon.

wilcox.test(Peso ~ Sesso, data = data)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Peso by Sesso
## W = 5e+05, p-value <2e-16
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(Lunghezza ~ Sesso, data = data)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Lunghezza by Sesso
## W = 6e+05, p-value <2e-16
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(Cranio~ Sesso, data = data)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Cranio by Sesso
## W = 6e+05, p-value = 7e-15
## alternative hypothesis: true location shift is not equal to 0

Si giunge alla stessa conclusione, ovvero si rifiuta l’ipotesi nulla pertanto le medie di peso, lunghezza e diametro del cranio tra maschi e femmine sono statisticamente differenti.

4.1 Conclusione sui test di ipotesi

  1. Abbiamo visto finora che non vi è dipendenza tra il tipo di parto e l’ospedale in cui nasce il neonato.

  2. Abbiamo appurato che la media del peso del campione non è statisticamente differente da quello della popolazione, mentre la media della lunghezza del campione è significativamente diversa da quella della popolazione.

  3. Si è verificato inoltre che la media delle misure antropometriche è significativamente diversa tra maschi e femmine, in particolare le misure medie risultano maggiori per i maschi.

Prossimi step

Lo scopo dell’analisi è quello di identificare quale variabile è più predittiva del peso alla nascita al fine di ridurre i rischi associati a nascite problematiche (es. neonati con basso peso).

I rischi sono identificabili come quei fattori che hanno un’influenza negativa sul peso alla nascita.

5. Analisi della correlazione

Verifichiamo se vi è una correlazione tra il peso e le seguenti variabili:

1) Peso vs madre fumatrice e Peso vs tipo parto: analisi grafica e t-test

#Relazione tra peso vs madre fumatrice/tipo parto
par(mfrow=c(1, 2)) 
boxplot(data$Peso~data$Fumatrici, data = data)
boxplot(data$Peso~data$Tipo.parto, data = data)

mean_fum <- t.test(data=data,
       Peso~Fumatrici,
       )
mean_parto <-t.test(data=data,
       Peso~Tipo.parto,
       )

knitr::kable(broom::tidy(mean_fum), caption = paste("T-Test Fumatrici"))
T-Test Fumatrici
estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high method alternative
50 3286 3236 1 0.3 114 -46 145 Welch Two Sample t-test two.sided 
knitr::kable(broom::tidy(mean_parto), caption = paste("T-Test Tipologia Parto"))
T-Test Tipologia Parto
estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high method alternative
-3 3282 3285 -0.14 0.89 1494 -46 40 Welch Two Sample t-test two.sided 
wilcox.test(data=data,
       Peso~Fumatrici,)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Peso by Fumatrici
## W = 1e+05, p-value = 0.06
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(data=data,
       Peso~Tipo.parto,)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Peso by Tipo.parto
## W = 6e+05, p-value = 0.5
## alternative hypothesis: true location shift is not equal to 0

Peso vs madri fumatrici: visivamente sembra esserci una differenza nel peso dei neonati di madri fumatrici (il valore medio è 50 grammi inferiore infatti nel caso di madre fumatrice), tuttavia non possiamo affermare che tale differenza sia significativa guardando al risultato della statistica test.

Nonostante anche il test Wilcoxon (meno influenzato dagli outliers) mostri un valore del p-value più vicino al livello di significatività del 5% non possiamo rifiutare l’ipotesi nulla.

Peso vs tipologia di parto: Non c’è differenza sul peso dei neonati nemmeno rispetto alla tipologia di parto, come evidente dai boxplot e dai risultati della statistica t (nonostante un p-value del test Wilcoxon pari alla soglia di significatività).

2) Correlazione tra il peso e le altre variabili quantitative: calcolo dell’indice di correlazione di Pearson (parametrico) e di Spearman (non-parametrico) e analisi grafica;

L’indice di Spearman non richiede ipotesi sulla distribuzione delle variabili ed è meno sensibile alla presenza di outliers.

L’indice di Pearson richiede invece sia l’ipotesi di normalità che di linearità tra le variabili.

#correlazione tra il peso e le altre variabili
dati_cor <- data[, c("Peso", "Anni.madre", "Gestazione", "N.gravidanze", "Lunghezza", "Cranio")]

# Calcolo della matrice di correlazione (Pearson)
correlazioni_P <- cor(dati_cor, method = "pearson", use = "complete.obs")
round(correlazioni_P, 2)
##               Peso Anni.madre Gestazione N.gravidanze Lunghezza Cranio
## Peso          1.00      -0.02       0.59         0.00      0.80   0.70
## Anni.madre   -0.02       1.00      -0.13         0.38     -0.06   0.02
## Gestazione    0.59      -0.13       1.00        -0.10      0.62   0.46
## N.gravidanze  0.00       0.38      -0.10         1.00     -0.06   0.04
## Lunghezza     0.80      -0.06       0.62        -0.06      1.00   0.60
## Cranio        0.70       0.02       0.46         0.04      0.60   1.00
# Calcolo della matrice di correlazione (Spearman)
correlazioni_S<- cor(dati_cor, method = "spearman", use = "complete.obs")
round(correlazioni_S, 2)
##               Peso Anni.madre Gestazione N.gravidanze Lunghezza Cranio
## Peso          1.00      -0.01       0.43         0.02      0.75   0.63
## Anni.madre   -0.01       1.00      -0.14         0.39     -0.05   0.03
## Gestazione    0.43      -0.14       1.00        -0.14      0.45   0.29
## N.gravidanze  0.02       0.39      -0.14         1.00     -0.07   0.06
## Lunghezza     0.75      -0.05       0.45        -0.07      1.00   0.52
## Cranio        0.63       0.03       0.29         0.06      0.52   1.00

I risultati dei due test (Person e Spearman) sono simili, pertanto rafforzano la conclusione che si può trarre su ogni variabile rispetto alla loro correlazione, come andiamo a documentare a seguire interprentanto la matrice di correlazione corredata di scatterplot.

Rappresentiamo un diagramma a matrice di dispersione e verifichiamo visivamente tutte le correlazioni tra variaibili quantitative:

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
data_corr<- data[, c("Anni.madre", "N.gravidanze", "Gestazione", "Lunghezza", "Cranio",  "Peso")]

pairs(data_corr,lower.panel=panel.cor, upper.panel=panel.smooth)

Gli scatterplot indicano che tra le variabili analizzate solo la lunghezza, la circonferenza del cranio e il numero di settimane di gestazione sembrano avere una correlazione positiva con il peso.

Gli indici di correlazione di Pearson e di Spearman indicano una relazione moderata-positiva (nel range 60%-80%), ma non perfettamente lineare in nessun caso.

Si può concludere che:

Sembra esserci inoltre una moderata correlazione tra Gestazione-Lunghezza e Lunghezza-Cranio; ad ogni modo, non è detto a priori che questo possa portare a problemi di multicollinearità nella nostra analisi di predizione tuttavia sarà importante verificarlo con gli opportuni test.

6. Modellazione: regressione lineare multipla

# Conversione variabili categoriche in fattori
data$Tipo.parto <- factor(data$Tipo.parto)
data$Ospedale <- factor(data$Ospedale)
data$Sesso <- factor(data$Sesso)

#modello di regressione lineare multipla:
mod1 = lm(Peso ~ Gestazione + Lunghezza + Ospedale + Cranio + Sesso + N.gravidanze + Tipo.parto + Fumatrici + Anni.madre, data=data)
summary(mod1)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Ospedale + Cranio + 
##     Sesso + N.gravidanze + Tipo.parto + Fumatrici + Anni.madre, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1123.3  -181.5   -14.5   161.0  2611.9 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6735.796    141.479  -47.61  < 2e-16 ***
## Gestazione       32.577      3.821    8.53  < 2e-16 ***
## Lunghezza        10.292      0.301   34.21  < 2e-16 ***
## Ospedaleosp2    -11.091     13.447   -0.82    0.410    
## Ospedaleosp3     28.250     13.505    2.09    0.037 *  
## Cranio           10.472      0.426   24.57  < 2e-16 ***
## SessoM           77.572     11.187    6.93  5.2e-12 ***
## N.gravidanze     11.381      4.669    2.44    0.015 *  
## Tipo.partoNat    29.634     12.091    2.45    0.014 *  
## Fumatrici       -30.274     27.549   -1.10    0.272    
## Anni.madre        0.802      1.147    0.70    0.484    
## ---
## 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.729,  Adjusted R-squared:  0.728 
## F-statistic:  669 on 10 and 2487 DF,  p-value: <2e-16

Osservando i p-value degli stimatori si deduce che le variabili altamente significative sono: Gestazione, Lunghezza, Cranio, Sesso. Il numero di gravidanze ha invece un significatività contenuta. Tipo di parto e Ospedale contribuiscono in parte alla spiegazione del peso alla nascita ma con un debole livello di significatività: il parto naturale si associa a un aumento medio del peso rispetto al parto cesareo, mentre Ospedale3 indica un possibile effetto positivo sul peso rispetto all’ospedale di riferimento.

Gestazione, Lunghezza, Cranio, Sesso, e Nr.gravidanze hanno un coefficiente positivo, all’aumentare di queste variabili aumenta il peso alla nascita.

Il modello spiega quasi il 73% della variabilità di peso di un neonato.

L’età della madre e il fatto che sia fumatrice non sembrano avere effetti significativi sul peso, proviamo a eliminare queste variabili dal modello.

Testiamo altri modelli:

mod2 <- update(mod1, ~ . - Fumatrici - Anni.madre)
summary(mod2)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Ospedale + Cranio + 
##     Sesso + N.gravidanze + Tipo.parto, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1113.1  -181.7   -16.7   161.1  2619.6 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6707.925    136.026  -49.31  < 2e-16 ***
## Gestazione       32.039      3.792    8.45  < 2e-16 ***
## Lunghezza        10.306      0.301   34.29  < 2e-16 ***
## Ospedaleosp2    -10.894     13.445   -0.81   0.4179    
## Ospedaleosp3     28.792     13.497    2.13   0.0330 *  
## Cranio           10.492      0.426   24.65  < 2e-16 ***
## SessoM           77.466     11.184    6.93  5.5e-12 ***
## N.gravidanze     12.336      4.334    2.85   0.0045 ** 
## Tipo.partoNat    29.408     12.087    2.43   0.0150 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274 on 2489 degrees of freedom
## Multiple R-squared:  0.729,  Adjusted R-squared:  0.728 
## F-statistic:  836 on 8 and 2489 DF,  p-value: <2e-16
anova(mod2, mod1)

L’analisi della varianza indica un p-value (0.43) non significativo, quindi l’aggiunta delle variabili non migliora il modello in modo statisticamente significativo.

L’R² ha lo stesso valore del modello completo (73%) pertanto l’esclusione dell’età della madre e dello status di fumatrici non ha alterato la significatività delle altre variaibli del modello sul peso; concludiamo quindi di poterle escludere dal modello a favore di un modello più semplice.

Proviamo a eliminare anche Tipo.Parto ed Ospedale:

mod3 <- update(mod2, ~ . - Tipo.parto - Ospedale)
summary(mod3)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     N.gravidanze, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1149.4  -181.0   -15.6   163.7  2639.1 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.725    135.804  -49.20  < 2e-16 ***
## Gestazione      32.383      3.801    8.52  < 2e-16 ***
## Lunghezza       10.245      0.301   34.06  < 2e-16 ***
## Cranio          10.541      0.426   24.72  < 2e-16 ***
## SessoM          77.981     11.211    6.96  4.5e-12 ***
## N.gravidanze    12.455      4.342    2.87   0.0042 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.726 
## F-statistic: 1.33e+03 on 5 and 2492 DF,  p-value: <2e-16
anova(mod3, mod2)

Il fatto che diminuendo il numero di variabili l’R² aggiustato rimanga stabile ci dice che quelle variabili non aggiungono vera capacità esplicativa.

Anche se l’anova tra il mod3 e il mod2 sembra suggerire che il tipo di ospedale e di parto aumentino la significatività del modello, notiamo che togliendo tali variaibli aumenta la significatività della variabile che esprime il numero di gravidanze, a fronte di un R² complessivo sostanzialmente invariato (72,7% vs 72,9%).

Pertanto scegliamo il modello più semplice con meno variaibili e tutte significative, ossia mod3.

Avendo osservato nei grafici di correlazione una relazione non lineare tra Gestazione e Peso, verifichiamo se l’aggiunta della relazione quadratica tra queste variabili ha effetto sul modello:

mod_quad <- update(mod3, ~ . + I(Gestazione^2))
summary(mod_quad)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     N.gravidanze + I(Gestazione^2), data = data)
## 
## 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.716    898.632   -5.17  2.5e-07 ***
## Gestazione        -81.231     49.740   -1.63   0.1026    
## Lunghezza          10.350      0.304   34.04  < 2e-16 ***
## Cranio             10.638      0.428   24.84  < 2e-16 ***
## SessoM             75.756     11.244    6.74  2.0e-11 ***
## N.gravidanze       12.549      4.338    2.89   0.0039 ** 
## I(Gestazione^2)     1.517      0.662    2.29   0.0221 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274 on 2491 degrees of freedom
## Multiple R-squared:  0.728,  Adjusted R-squared:  0.727 
## F-statistic: 1.11e+03 on 6 and 2491 DF,  p-value: <2e-16
anova(mod_quad, mod3)

Sebbene l’ANOVA restituisca un p-value che indica significatività (seppure debole), è altresì vero che la significatività del termine quadratico è debole e il termine lineare perde di significatività rispetto al modello 3.

Inoltre non c’è alcun miglioramento dell’R2.

Concluderei preferendo il modello 3 a favore della semplicità.

Metodo AIC/BIC

Confrontiamo i modelli analizzati finora:

AIC(mod1, mod2, mod3, mod_quad)
BIC(mod1, mod2, mod3, mod_quad)

Il BIC tende a preferire modelli più semplici, ed infatti suggerisce un valore più basso degli altri per il modello 3 (quindi preferibile). L’analisi AIC invece suggerirebbe il modello 2. Tra le due opzioni, scegliamo numavamente il modello con meno variabili.

Verifichiamo gli indicatori di multicollinearità del modello:

library(car)
vif(mod3)
##   Gestazione    Lunghezza       Cranio        Sesso N.gravidanze 
##          1.7          2.1          1.6          1.0          1.0

Nel modello i Vif sono minori di 5 quindi non è presente multicollinearità.

Modello definitivo:

scegliamo il modello 3 come modello definitivo.

mod_def = lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso + N.gravidanze, data=data)
summary(mod_def)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     N.gravidanze, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1149.4  -181.0   -15.6   163.7  2639.1 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.725    135.804  -49.20  < 2e-16 ***
## Gestazione      32.383      3.801    8.52  < 2e-16 ***
## Lunghezza       10.245      0.301   34.06  < 2e-16 ***
## Cranio          10.541      0.426   24.72  < 2e-16 ***
## SessoM          77.981     11.211    6.96  4.5e-12 ***
## N.gravidanze    12.455      4.342    2.87   0.0042 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.726 
## F-statistic: 1.33e+03 on 5 and 2492 DF,  p-value: <2e-16

7. Selezione del modello tramite procedura stepwise (con BIC)

Verifichiamo tramite la procedura automatica stepwise che il modello migliore sia il mod3.

mod_step <- stepAIC(mod1, direction = "both", k=log(nrow(data))) #log di n per BIC
## Start:  AIC=28119
## Peso ~ Gestazione + Lunghezza + Ospedale + Cranio + Sesso + N.gravidanze + 
##     Tipo.parto + Fumatrici + Anni.madre
## 
##                Df Sum of Sq      RSS   AIC
## - Anni.madre    1     36710 1.87e+08 28111
## - Fumatrici     1     90677 1.87e+08 28112
## - Ospedale      2    687555 1.87e+08 28112
## - N.gravidanze  1    446244 1.87e+08 28117
## - Tipo.parto    1    451073 1.87e+08 28117
## <none>                      1.87e+08 28119
## - Sesso         1   3610705 1.90e+08 28159
## - Gestazione    1   5458852 1.92e+08 28183
## - Cranio        1  45318506 2.32e+08 28654
## - Lunghezza     1  87861708 2.75e+08 29074
## 
## Step:  AIC=28111
## Peso ~ Gestazione + Lunghezza + Ospedale + Cranio + Sesso + N.gravidanze + 
##     Tipo.parto + Fumatrici
## 
##                Df Sum of Sq      RSS   AIC
## - Fumatrici     1     91599 1.87e+08 28105
## - Ospedale      2    693914 1.87e+08 28105
## - Tipo.parto    1    452049 1.87e+08 28109
## <none>                      1.87e+08 28111
## - N.gravidanze  1    631082 1.87e+08 28112
## + Anni.madre    1     36710 1.87e+08 28119
## - Sesso         1   3617809 1.90e+08 28151
## - Gestazione    1   5424800 1.92e+08 28175
## - Cranio        1  45569477 2.32e+08 28649
## - Lunghezza     1  87852027 2.75e+08 29066
## 
## Step:  AIC=28105
## Peso ~ Gestazione + Lunghezza + Ospedale + Cranio + Sesso + N.gravidanze + 
##     Tipo.parto
## 
##                Df Sum of Sq      RSS   AIC
## - Ospedale      2    702925 1.88e+08 28098
## - Tipo.parto    1    444404 1.87e+08 28103
## <none>                      1.87e+08 28105
## - N.gravidanze  1    608136 1.87e+08 28105
## + Fumatrici     1     91599 1.87e+08 28111
## + Anni.madre    1     37633 1.87e+08 28112
## - Sesso         1   3601860 1.90e+08 28145
## - Gestazione    1   5358199 1.92e+08 28167
## - Cranio        1  45613331 2.32e+08 28642
## - Lunghezza     1  88259386 2.75e+08 29063
## 
## Step:  AIC=28098
## Peso ~ Gestazione + Lunghezza + Cranio + Sesso + N.gravidanze + 
##     Tipo.parto
## 
##                Df Sum of Sq      RSS   AIC
## - Tipo.parto    1    467626 1.88e+08 28097
## <none>                      1.88e+08 28098
## - N.gravidanze  1    648873 1.88e+08 28099
## + Ospedale      2    702925 1.87e+08 28105
## + Fumatrici     1    100610 1.87e+08 28105
## + Anni.madre    1     44184 1.88e+08 28106
## - Sesso         1   3644818 1.91e+08 28139
## - Gestazione    1   5457887 1.93e+08 28162
## - Cranio        1  45747094 2.33e+08 28636
## - Lunghezza     1  87955701 2.76e+08 29051
## 
## Step:  AIC=28097
## Peso ~ Gestazione + Lunghezza + Cranio + Sesso + N.gravidanze
## 
##                Df Sum of Sq      RSS   AIC
## <none>                      1.88e+08 28097
## - N.gravidanze  1    621053 1.89e+08 28097
## + Tipo.parto    1    467626 1.88e+08 28098
## + Ospedale      2    726146 1.87e+08 28103
## + Fumatrici     1     92548 1.88e+08 28103
## + Anni.madre    1     45366 1.88e+08 28104
## - Sesso         1   3650790 1.92e+08 28137
## - Gestazione    1   5477493 1.94e+08 28161
## - Cranio        1  46098547 2.34e+08 28637
## - Lunghezza     1  87532691 2.76e+08 29044
summary(mod_step)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     N.gravidanze, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1149.4  -181.0   -15.6   163.7  2639.1 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.725    135.804  -49.20  < 2e-16 ***
## Gestazione      32.383      3.801    8.52  < 2e-16 ***
## Lunghezza       10.245      0.301   34.06  < 2e-16 ***
## Cranio          10.541      0.426   24.72  < 2e-16 ***
## SessoM          77.981     11.211    6.96  4.5e-12 ***
## N.gravidanze    12.455      4.342    2.87   0.0042 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.726 
## F-statistic: 1.33e+03 on 5 and 2492 DF,  p-value: <2e-16
BIC(mod_def, mod_step)

la procedura conferma tale scelta, il valore di R2 di 0.727 indica che le variabili sono in grado di spiegare quasi il 73% della variabilità del campione.

fitted <- fitted(mod_def)

plot(data$Peso, fitted, 
     main = "Osservati vs Predetti", 
     xlab = "Osservati", 
     ylab = "Predetti",
     pch = 20, col = "darkgreen")

abline(a = 0, b = 1, col = "red", lty = 4, lwd = 4) 

# Estrazione tabella coefficienti
coef_tab <- summary(mod_def)$coefficients


coef_df <- as.data.frame(coef_tab)
coef_df$Variabile <- rownames(coef_df)
rownames(coef_df) <- NULL


colnames(coef_df) <- c("Stima", "Errore Std", "t value", "p value", "Variabile")


coef_df <- coef_df[, c("Variabile", "Stima", "Errore Std", "t value", "p value")]


print(
    knitr::kable(
  coef_df,
  caption = "Coefficienti stimati del modello lineare",
  digits = 2,
  format.args = list(decimal.mark = ",", big.mark = ".")
)
)
## 
## 
## Table: Coefficienti stimati del modello lineare
## 
## |Variabile    |  Stima| Errore Std| t value| p value|
## |:------------|------:|----------:|-------:|-------:|
## |(Intercept)  | -6.682|     135,80|   -49,2|       0|
## |Gestazione   |     32|       3,80|     8,5|       0|
## |Lunghezza    |     10|       0,30|    34,1|       0|
## |Cranio       |     11|       0,43|    24,7|       0|
## |SessoM       |     78|      11,21|     7,0|       0|
## |N.gravidanze |     12|       4,34|     2,9|       0|

Il modello spiega quasi il 73% della variabilità nel peso dei neonati.

Tutti i predittori sono significativi (p < 0.05), con forte effetto di Lunghezza e Cranio stando al valore della stastica t (entrambe le variabili sono molto correlate a Peso); l’effetto è di +10 e +11 grammi in più per ogni ulteriore mm.

Gestazione ha un effetto positivo e significativo sul peso, ma più debole rispetto a Lunghezza e Cranio;

Sesso: i maschi pesano mediamente 78 grammi in più delle femmine;

Numero di gravidanze: effetto positivo moderato, 12 grammi in più sul peso alla nascita.

Graficamente notiamo che non tutti i punti del grafico si distribuiscono bene lungo la retta, indicando una certa non-linearità.

Per alcuni valori bassi o alti di peso il modello potrebbe sotto- o sovrastimare il peso alla nascita.

Approfondiamo questo comportamento del modello con l’analisi dei residui.

8.Diagnostica sui residui

Analisi grafica

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

Verifichiamo le seguenti ipotesi:

  • I residui seguono una distribuzione normale;
  • hanno media pari a zero;
  • sono indipendenti e identicamente distribuiti (i.i.d.) con varianza costante (omoschedasticità) lungo il range dei valori predetti;
  • hanno una relazione lineare con i predittori.

Media costante zero (residual vs fitted): I residui tendono a distribuirsi più nella parte destra della distribuzione che nella parte sinistra, pertanto il modello potrebbe non adattarsi bene per neonati sottopeso; tuttavia per valori più alti il trend si stabilizza attorno allo zero, intorno al quale i valori sembrano sparsi in maniera casuale e abbastanza centrata.

Grafico Q-Q (Quantile-Quantile): I punti seguono approssimativamente la linea diagonale nella parte centrale della distribuzione deviando soprattutto per valori alti, sembrano pertanto distribuiti quasi normalmente.

Omogeneità della varianza dei residui (Scale-location): La linea rossa non è perfettamente costante (omoschedasticità) quindi potrebbe esserci una lieve eteroschedasticità.

Valori influent (residual vs leverage): l’oservazione 1551 è un valore influente (distanza di Cook oltre 0.5).

leva <-hatvalues(mod_def)
plot(leva)
p<-sum(leva)
n<-length(leva)
soglia=2*p/n
abline(h=soglia,col=2)

plot(rstudent(mod_def))
abline(h=c(-2,2))
abline(h=mean(residuals(mod_def)), col="red")

car::outlierTest(mod_def)
##      rstudent unadjusted p-value Bonferroni p
## 1551     10.0            2.6e-23      6.6e-20
## 155       5.0            5.4e-07      1.3e-03
## 1306      4.8            1.5e-06      3.7e-03
cook<-cooks.distance(mod_def)
plot(cook,ylim = c(0,1)) 

plot(mod_def, which = 4, id.n = 3)

  • La maggior parte delle osservazioni ha leverage basso, ma alcune superano la soglia e potrebbero essere influenti nella stima dei coefficienti.

  • I residui hanno media zero e sembrano sparsi casualmente attorno alla media.

  • L’outlier test ci mostra i tre residui estremi con p-value molto bassi, quindi sono outliers statisticamente significativi anche dopo la correzione di Bonferroni, ma l’osservazione 1551 è un outlier estremo.

  • Il cook.distance conferma che l’osservazione 1551 è molto influente (>0.5)

Test diagnostici sui residui

library(lmtest)
bptest(mod_def)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod_def
## BP = 90, df = 5, p-value <2e-16
dwtest(mod_def)
## 
##  Durbin-Watson test
## 
## data:  mod_def
## DW = 2, p-value = 0.1
## alternative hypothesis: true autocorrelation is greater than 0
shapiro.test(residuals(mod_def))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_def)
## W = 1, p-value <2e-16
plot(density(residuals(mod_def)))

Test Breusch-Pagan: si rifiuta l’ipotesi nulla, c’è eteroschedasticità (la varianza degli errori non è costante).

Test Durbin-Watson: non si rifiuta l’ipotesi nulla quindi i residui non sono autocorrelati (residui indipendenti)

Test Shapiro-Wilk: si rifiuta l’ipotesi nulla, quindi i residui non sono distribuiti normalmente.

La densità dei residui mostra ad ogni modo una distribuzione compatibile con una normale con una coda lunga a destra.

9. Gestione degli outliers

Abbiamo rilevato che l’osservazione 1551 è un outlier estremo. Procediamo a rimuovere l’outlier dal dataset e analizziamo nuovamente il nostro modello.

#Eliminiamo 1551 dal dataset e verifichiamo
rownames(data) <- NULL
data_wo_outlier <- data[-1551, ]
identical(data[1551, ], data_wo_outlier[1551, ])  # deve restituire FALSE
## [1] FALSE
data[1550:1552, ]
data_wo_outlier[1550:1552, ]

Eliminazione avvenuta con successo.

mod_no_outlier <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso + N.gravidanze, data=data_wo_outlier)
car::outlierTest(mod_no_outlier)
##      rstudent unadjusted p-value Bonferroni p
## 1549     10.0            2.6e-23      6.6e-20
## 155       5.0            5.3e-07      1.3e-03
## 1305      4.8            1.4e-06      3.6e-03
cook<-cooks.distance(mod_no_outlier)
plot(cook,ylim = c(0,1)) 

plot(mod_no_outlier, which = 4, id.n = 3)

Dopo la correzione del dataset risulta che vi siano ancora tre osservazioni potenzialmente influenti con p-value corretti molto bassi (1549, 155, 1305). Tuttavia la distanza di Cook per tali osservazioni è sotto la soglia, tranne per l’osservazione 1549 che potrebbe avere un’ influenza considerevole sul modello. Procediamo con la rimozione anche di tale valore.

#Eliminiamo 1549 dal dataset e verifichiamo
rownames(data) <- NULL
data_wo_outlier <- data_wo_outlier[-1549, ]
identical(data[1549, ], data_wo_outlier[1549, ])  # deve restituire FALSE
## [1] FALSE
data[1548:1552, ]
data_wo_outlier[1548:1552, ]

Eliminazione avvenuta con successo. Verifichiamo nuovamente eventuali valori influenti.

mod_no_outlier <- lm(Peso ~ Gestazione + Lunghezza + Cranio + Sesso + N.gravidanze, data=data_wo_outlier)
car::outlierTest(mod_no_outlier)
##      rstudent unadjusted p-value Bonferroni p
## 155       5.3            1.3e-07      0.00034
## 1305      4.9            8.0e-07      0.00199
## 1397     -4.3            1.4e-05      0.03594
cook<-cooks.distance(mod_no_outlier)
plot(cook,ylim = c(0,1)) 

plot(mod_no_outlier, which = 4, id.n = 3)

Il test sugli outliers presenta ancora delle osservazioni con p-value molto bassi tuttavia tali valori si possono non considerare influenti sul modello osservando la cook distance e pertanto non verranno rimossi.

Compariamo il modello senza outlier con quello precedente.

summary(mod_def)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     N.gravidanze, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1149.4  -181.0   -15.6   163.7  2639.1 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.725    135.804  -49.20  < 2e-16 ***
## Gestazione      32.383      3.801    8.52  < 2e-16 ***
## Lunghezza       10.245      0.301   34.06  < 2e-16 ***
## Cranio          10.541      0.426   24.72  < 2e-16 ***
## SessoM          77.981     11.211    6.96  4.5e-12 ***
## N.gravidanze    12.455      4.342    2.87   0.0042 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.726 
## F-statistic: 1.33e+03 on 5 and 2492 DF,  p-value: <2e-16
summary(mod_no_outlier)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     N.gravidanze, data = data_wo_outlier)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1162.8  -180.7   -12.7   162.9  1408.5 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6684.809    132.899  -50.30  < 2e-16 ***
## Gestazione      30.470      3.738    8.15  5.6e-16 ***
## Lunghezza       10.846      0.302   35.97  < 2e-16 ***
## Cranio           9.890      0.422   23.44  < 2e-16 ***
## SessoM          79.041     10.975    7.20  7.8e-13 ***
## N.gravidanze    12.566      4.253    2.95   0.0032 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269 on 2490 degrees of freedom
## Multiple R-squared:  0.738,  Adjusted R-squared:  0.737 
## F-statistic: 1.4e+03 on 5 and 2490 DF,  p-value: <2e-16

Dopo la rimozione degli outliers influenti (1551, 1549) il modello mostra un miglioramento generale.

L’R² sale da 72,7% a 73,8% e l’errore standard residuo scende da 275 a 269, indicando una maggiore precisione delle stime.

Tutti i predittori restano significativi, con i p-value di Sesso e Nr. Gravidanze migliorati.

Il coefficiente di Gestazione si riduce leggermente, suggerendo che gli outliers influivano probabilmente in particolare su questa variabile.

Il nuovo modello può quindi essere considerato più affidabile e rappresentativo.

BIC(mod_def, mod_no_outlier)
plot(mod_no_outlier)

library(lmtest)
bptest(mod_no_outlier)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod_no_outlier
## BP = 10, df = 5, p-value = 0.08
dwtest(mod_no_outlier)
## 
##  Durbin-Watson test
## 
## data:  mod_no_outlier
## DW = 2, p-value = 0.1
## alternative hypothesis: true autocorrelation is greater than 0
shapiro.test(residuals(mod_no_outlier))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_no_outlier)
## W = 1, p-value = 9e-13
plot(residuals(mod_no_outlier))
abline(h=mean(residuals(mod_no_outlier)),col="red")

plot(density(residuals(mod_no_outlier)))

L’indice BIC è più basso, quindi indica un modello migliore.

Il modello, come precedentemente osservato, sembra sottostimare i valori meno accuratamente per i neonati sottopeso.

Il test Breusch-Pagan indica che il problema di eteroschedasticità sembra risolto grazie all’eliminazione degli outliers.

Continua a non esserci autocorrelazione tra i residui.

I residui hanno media zero e sembrano sparsi casualmente intorno alla media. Il test di normalità indica che i residui non sono perfettamente normali, tuttarvia osservando graficamente potremmo concludere che i residui si distribuiscono quasi normalmente soprattutto nella massa centrale della distribuzione, con deviazioni sulle code.

In conclusione, il modello complessivamente sembra un buon strumento per condurre un’analisi predittiva sul peso dei neonati.

Come ultimo step, proviamo a considerare una trasformazione logarirmica dei per trattare la non-normalità.

10. Trasformazione logaritmica

Applichiamo ora una trasformazione logaritmica alla variabile target peso poi confrontiamo i modelli:

mod_log <- lm(log(Peso) ~ Gestazione + Lunghezza + Cranio + Sesso + N.gravidanze, data=data_wo_outlier)
summary(mod_no_outlier)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     N.gravidanze, data = data_wo_outlier)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1162.8  -180.7   -12.7   162.9  1408.5 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6684.809    132.899  -50.30  < 2e-16 ***
## Gestazione      30.470      3.738    8.15  5.6e-16 ***
## Lunghezza       10.846      0.302   35.97  < 2e-16 ***
## Cranio           9.890      0.422   23.44  < 2e-16 ***
## SessoM          79.041     10.975    7.20  7.8e-13 ***
## N.gravidanze    12.566      4.253    2.95   0.0032 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269 on 2490 degrees of freedom
## Multiple R-squared:  0.738,  Adjusted R-squared:  0.737 
## F-statistic: 1.4e+03 on 5 and 2490 DF,  p-value: <2e-16
summary(mod_log)
## 
## Call:
## lm(formula = log(Peso) ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     N.gravidanze, data = data_wo_outlier)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.4108 -0.0536  0.0006  0.0549  0.5069 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.40e+00   4.19e-02  105.02  < 2e-16 ***
## Gestazione   1.79e-02   1.18e-03   15.22  < 2e-16 ***
## Lunghezza    3.72e-03   9.51e-05   39.17  < 2e-16 ***
## Cranio       3.32e-03   1.33e-04   24.92  < 2e-16 ***
## SessoM       1.84e-02   3.46e-03    5.31  1.2e-07 ***
## N.gravidanze 4.07e-03   1.34e-03    3.04   0.0024 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.085 on 2490 degrees of freedom
## Multiple R-squared:  0.786,  Adjusted R-squared:  0.786 
## F-statistic: 1.83e+03 on 5 and 2490 DF,  p-value: <2e-16

R2 e adjusted R2 sono aumentati passando a 78,6% (entrambi) rispetto a 73,8% e 73,7%.

Tutti i predittori rimangono altamente significativi.

par(mfrow = c(2,2),
    mar = c(2, 4.1, 2, 2.1))

plot(mod_log)

BIC(mod_log, mod_no_outlier, mod_def)
bptest(mod_log)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod_log
## BP = 90, df = 5, p-value <2e-16
dwtest(mod_log)
## 
##  Durbin-Watson test
## 
## data:  mod_log
## DW = 2, p-value = 0.1
## alternative hypothesis: true autocorrelation is greater than 0
shapiro.test(mod_log$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod_log$residuals
## W = 1, p-value = 4e-12
fitted_log<- exp(fitted(mod_log))

# Grafico degli osservati vs predetti (sui dati originali)
plot(data_wo_outlier$Peso, fitted_log, 
     main = "Osservati vs Predetti", 
     xlab = "Osservati", 
     ylab = "Predetti",
     pch = 20, col = "blue")

# Aggiungi la linea ideale (linea di identità)
abline(a = 0, b = 1, col = "red", lty = 2, lwd = 2)  # Linea di identità

plot(residuals(mod_log))
abline(h=mean(residuals(mod_log)),col="red")

plot(density(residuals(mod_log)))

Nonostante Shapiro Test ci porti a rifiutare l’ipotesi nulla, il grafico di distribuzione dei residui sembra suggerire una distribuzione quasi normale.

Il grafico dei valori predetti ed osservati, l’R2 migliore e il BIC più basso ci indicano che il modello logaritmico sembrerebbe il migliore tra quelli testati.

I redisui hanno media zero e sembrano sparsi casualmente intorno alla media.

Tuttavia c’è eteroschedasticità, infatti rifiutiamo l’ipotesi nulla del test BP. Pertanto, rispetto al modello lineare semplice che soddisfava tutti i test e presentava anch’esso una distribuzione dei residui quasi normale, non sembra ragionevole passare ad una trasformazione logaritmica che peggiora la varianza dei residui.

Consideriamo pertanto di prediligere il modello lineare senza outliers che non presenta eteroschedasticità e nemmeno autocorrelazione tra i residui, sebbene rimanga violata l’ipotesi di normalità.

Calcoliamo i p-value robusti e verifichiamo se i coefficienti sono statisticamente significativi anche senza assumere normalità.

library(sandwich)
library(lmtest)


coeftest(mod_no_outlier, vcov. = vcovHC(mod_no_outlier, type = "HC1"))
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6684.809    167.917  -39.81  < 2e-16 ***
## Gestazione      30.470      4.395    6.93  5.2e-12 ***
## Lunghezza       10.846      0.372   29.19  < 2e-16 ***
## Cranio           9.890      0.511   19.35  < 2e-16 ***
## SessoM          79.041     11.039    7.16  1.1e-12 ***
## N.gravidanze    12.566      4.503    2.79   0.0053 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nessuna variabile perde significatività quindi il modello è solido.

11 Calcolo RMSE

Il Root Mean Squared Error (RMSE) è una metrica utilizzata per valutare la qualità di un modello di regressione.

residui <- data_wo_outlier$Peso - predict(mod_no_outlier, newdata = data_wo_outlier)
# Calcolo del Mean Squared Error (MSE)
mse <- mean(residui^2)
# Calcolo del Root Mean Squared Error (RMSE)
rmse <- sqrt(mse)
cat("Il Root Mean Squared Error (RMSE) del modello lineare semplice è:", rmse, "\n")
## Il Root Mean Squared Error (RMSE) del modello lineare semplice è: 268
residui_log <- data_wo_outlier$Peso - exp(predict(mod_log, newdata = data_wo_outlier))
# Calcolo del Mean Squared Error (MSE)
mse_log <- mean(residui_log^2)
# Calcolo del Root Mean Squared Error (RMSE)
rmse_log <- sqrt(mse_log)
cat("Il Root Mean Squared Error (RMSE) del modello con trasformazione logaritmica è:", rmse_log, "\n")
## Il Root Mean Squared Error (RMSE) del modello con trasformazione logaritmica è: 272

In media il modello lineare semplice sbaglia di circa 268g, quindi performa meglio (di pochi grammi) rispetto al modello con trasformazione logartmica.

cat("R²:", round(summary(mod_no_outlier)$r.squared, 2), "\n")
## R²: 0.74
cat("R²:", round(summary(mod_log)$r.squared, 2), "\n")
## R²: 0.79

Il modello lineare semplice spiega circa il 74% della variabilità osservata nel nel peso, mentre il modello con trasformazione logaritmica spiega il 79% della variabilità osservata.

12.Conclusioni - scelta del modello:

Useremo il modello lineare semplice “mod_wo_outlier” che si definisce come un una regressione lineare semplice che stima il peso neonatale in funzione delle seguenti variabili esplicative:

Il modello rappresenta un buon modello per fare predizioni del peso di un neonato alla nascita in quanto, tra i vari testati, rappresenta un buon compromesso tra efficienza e semplicità di interpretazione.

Tutti i predittori sono significativi, l’R² è inferiore rispetto al modello logaritmico ma non di molto, e si è dimostrato tra i più alti registrati tra i vari modelli testati;

l’R² ci suggerisce che il 74% della variabilità del peso è spiegata dal modello. Il BIC predilige questo modello ad altri e l’analisi dei residui non ha rilevato problemi di omoschedasticità e di autocorrelazione e la loro distribuzione è sufficientemente normale.

Il modello di regressione lineare “mod_wo_outlier” si considera pertanto un buon modello per predire il peso dei neonati.

summary(mod_no_outlier)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso + 
##     N.gravidanze, data = data_wo_outlier)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1162.8  -180.7   -12.7   162.9  1408.5 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6684.809    132.899  -50.30  < 2e-16 ***
## Gestazione      30.470      3.738    8.15  5.6e-16 ***
## Lunghezza       10.846      0.302   35.97  < 2e-16 ***
## Cranio           9.890      0.422   23.44  < 2e-16 ***
## SessoM          79.041     10.975    7.20  7.8e-13 ***
## N.gravidanze    12.566      4.253    2.95   0.0032 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 269 on 2490 degrees of freedom
## Multiple R-squared:  0.738,  Adjusted R-squared:  0.737 
## F-statistic: 1.4e+03 on 5 and 2490 DF,  p-value: <2e-16

Interpretazione dell’effetto marginale di ogni variabile:

13. Previsioni e Risultati

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

Creazione di un nuovo dataset:

Gestazione (settimane): 39;

Lunghezza (mm): si assume valore medio;

Cranio(mm): si assume valore medio;

Sesso: Femmina;

N.gravidanze: 3

lunghezza_f <- mean(data_wo_outlier$Lunghezza[data_wo_outlier$Sesso == "F"])
cranio_f <- mean(data_wo_outlier$Cranio[data_wo_outlier$Sesso == "F"])

data_test <- data.frame(
  Gestazione = 39,
  Lunghezza = lunghezza_f,
  Cranio = cranio_f,
  N.gravidanze = 3,
  Sesso = "F"
)

peso_test <- predict(mod_no_outlier, newdata = data_test, interval = "prediction")

cat(
  "Il peso stimato della neonata alla nascita è:",
  peso_test[1, "fit"], "gr.",
  "L'intervallo di previsione al 95% va da:",
  peso_test[1, "lwr"], "a:",
  peso_test[1, "upr"], "gr.\n"
)
## Il peso stimato della neonata alla nascita è: 3193 gr. L'intervallo di previsione al 95% va da: 2665 a: 3721 gr.

14.Visualizzazioni

Infine, utilizzeremo grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili.

Rappresentiamo la relazione tridimensionali tra la variabile risultato “Peso” e le variabili “Lunghezza” e “Cranio” in base alla variabile “Sesso”.

library(plotly)

plot_ly(
  data = data_wo_outlier,
  x = ~Lunghezza,
  y = ~Peso,
  z = ~Cranio,
  color = ~Sesso,
  colors = c("green", "blue"),
  size = ~Lunghezza,
  sizes = c(5, 30),
  type = "scatter3d",
  mode = "markers"
)

Questo grafico illustra la relazione positiva tra la lunghezza del neonato e il suo peso alla nascita, così come tra dimensione del cranio e peso alla nascita, differenziata per sesso.

Si osserva che a parità di lunghezza e a parità di dimensioni del cranio i neonati maschi tendono a pesare leggermente di più delle femmine e questa risultanza è in linea con quanto descritto dal nostro modello.


Visualizziamo questa relazione tra le variabili lunghezza e diametro del cranio sul peso anche attraverso un’altra tipologia di rappresentazione.

library(ggplot2)

# Variabili indipendenti da confrontare
variabili <- c("Cranio", "Lunghezza")

# Genera i grafici e stampali uno per uno
for (var in variabili) {
  p <- ggplot(data_wo_outlier, aes_string(x = var, y = "Peso", color = "Sesso")) +
    geom_point(position = position_jitter(width = 0.3, height = 0), alpha = 0.6) +
    geom_smooth(method = "lm", se = FALSE) +
    labs(
      title = paste("Relazione tra", var, "e Peso alla Nascita"),
      x = paste(var, "(mm)"),
      y = "Peso alla Nascita (grammi)",
      color = "Sesso"
    ) +
    theme_minimal()
  
  print(p)
}

Visualizziamo ora l’andamento del peso al variare delle settimane di gestazione per genere e a seguire l’andamento del peso al variare del numero di gravidanze della madre per genere.

ggplot(data = data_wo_outlier) +
  geom_point(aes(x = Gestazione,
                 y = Peso,
                 col = Sesso), position = "jitter") +
  geom_smooth(aes(x = Gestazione,
                  y = Peso,
                  col = Sesso), se = FALSE, method = "lm") +
  labs(title = "Relazione tra settimane di gestazione e peso alla nascita",
       x = "Gestazione (settimane)",
       y = "Peso alla Nascita (grammi)",
       color = "Sesso") +
  theme_minimal()

Le rette hanno la stessa pendenza tuttavia la retta dei maschi è più alta rispetto a quella delle femmine in quanto a parità di settimane di gestazione il peso dei maschi alla nascita è maggiore di quello delle femmine.

ggplot(data = data_wo_outlier) +
  geom_point(aes(x = N.gravidanze,
                 y = Peso,
                 col = Sesso), position = "jitter") +
  geom_smooth(aes(x = N.gravidanze,
                  y = Peso,
                  col = Sesso), se = FALSE, method = "lm") +
  labs(title = "Relazione tra numero di gravidanze e peso alla nascita",
       x = "Nr. Gravidanze",
       y = "Peso alla Nascita (grammi)",
       color = "Sesso") +
  theme_minimal()

Il peso dei maschi risulta essere maggiore di quelle delle femmine tuttavia le pendenza della retta suggerisce che il peso per entrambi i sessi rimane pressoché costante all’aumentare del numero delle gravidanze della madre, questo si può spiegare dal fatto che nella maggior parte dei casi il numero delle gravidanze per ogni madre è inferiore a 3.

Visualiziamo anche l’impatto del numero di settimane di gestazione e del fumo sul peso previsto.

ggplot(data = data_wo_outlier) +
  geom_point(aes(x = Gestazione, y = Peso, color = factor(Fumatrici)), alpha = 0.6) +
  geom_smooth(aes(x = Gestazione, y = Peso, color = factor(Fumatrici)), method = "lm", se = FALSE) +
  labs(
    title = "Relazione tra settimane di gestazione e peso alla nascita",
    x = "Settimane di Gestazione",
    y = "Peso alla Nascita (grammi)",
    color = "Fumatrici"
  ) +
  scale_color_manual(values = c("steelblue", "tomato")) +
  theme_minimal()

Ricordiamo che nel dataset le madri fumatrici sono poche, pertanto questa relazione risulta difficilmente confrontabile tra i due gruppi. Ad ogni modo, possiamo notare che verso la fine del periodo di gestazione la retta delle madri fumatrici si colloca al di sotto della retta delle madri non fumatrici, evidenziando una potenziale incidenza negativa dello status di madre fumatrice sul sesso alla nascita.

15.Conclusioni

L’analisi del dataset mirava ad indagare l’influenza di alcune variabili fisiologiche e di altri potenziali fattori di rischio sul peso di un neonato, al fine di consentire una più corretta gestione di tali fattori nella fase pre-natale.

Abbiamo osservato che la durata della gestazione è un importante predittore del peso alla nascita, infatti per ogni settimana in più il peso aumenta di ben 30gr. Dal punto di vista clinico, questa informazione suggerisce che la prevenzione dei parti prematuri può influire positivamente sulla salute dei nascituri.

Vi è una differenza significativa di peso tra i neonati maschi e femmine, a parità di altre condizioni, pertanto è importanto tenere presente il sesso nelle considerazioni cliniche di monitoraggio della gestazione.

Purtoppo il dataset presentava poche osservazioni relative a madri fumatrici. Sebbene i pochi dati disponibili sembrino indicare un’influenza negativa dell’abitudine a fumare delle madri sul peso alla nascita, ricordiamo che il t-test non ci permetteva di concludere che la differenza tra l’essere o meno fumatrice fosse significativa. Il p-value del regressore per le madri fumatrici, inoltre, non era significativo e la variabile non apportava alcun beneficio al modello di regressione.

Il fumo, in quanto fattore di rischio per la salute di qualsiasi persona, potrebbe essere un fattore di rischio da monitorare con attenzione qualora si potesse constatare l’effettiva influenza sul peso alla nascita. Si suggerisce pertanto di valutare tale potenziale effetto conducendo un’analisi ulteriore riperformando l’esperimento su un altro campione di dati.

Il modello di regressione lineare semplice si è dimostrato essere il modello migliore tra tutti i modelli testati, soprattutto dopo l’eliminazione di outliers problematici.

Nonostante graficamente il modello sembri non catturare a pieno la relazione non perfettamente lineare tra i valori osservati e predetti, l’aggiunta della relazione quadratica (che sembra esserci tra le settimane di gestazione e il peso) non ha portato benefici al modello.

L’analisi diagnostica dei residui non ha permesso di validare l’ipotesi di normalità, sebbene l’analisi grafica dimostri che la distribuzione dei residui si possa considerare quasi normale. Si è testato a tal fine un modello con trasformazione logaritmica che sebbene permettesse di aumentare l’R² del modello, ha fatto emergere la presenza di eteroschedasticità nei residui senza risolvere la non-normalità della loro distribuzione. Per questo motivo, si è nuovamente prediletto il modello lineare, più semplice e robusto, in grado di spiegare il 74% della varianza del peso alla nascita con un errore mediamente pari a ca. 268 grammi.