1. Struttura del dataset

library(moments)
library(knitr)
library(kableExtra)
library(patchwork)
library(ggplot2)
library(dplyr)
library(glue)
library(tidyr)
library(ggpubr)
library(car)
library(Metrics)
library(broom)
library(lmtest)

dati_neonati <- read.csv("neonati.csv", stringsAsFactors = TRUE)
kable(head(dati_neonati,10))
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
26 0 0 42 3380 490 325 Nat osp3 M
21 2 0 39 3150 490 345 Nat osp1 F
34 3 0 38 3640 500 375 Nat osp2 M
28 1 0 41 3690 515 365 Nat osp2 M
20 0 0 38 3700 480 335 Nat osp3 F
32 0 0 40 3200 495 340 Nat osp2 F
26 1 0 39 3100 480 345 Nat osp3 F
25 0 0 40 3580 510 349 Nat osp1 M
22 1 0 40 3670 500 335 Ces osp2 F
23 0 0 41 3700 510 362 Ces osp2 F 
dimensioni <- data.frame(
  num_entries = nrow(dati_neonati), 
  num_variabili = ncol(dati_neonati)
)

kable(dimensioni, col.names = c("num_entries", "num_variabili"))%>%
  kable_styling(full_width = FALSE) %>% 
  column_spec(1, width = "5cm") %>%      
  column_spec(2, width = "5cm")
num_entries num_variabili
2500 10

Il dataset contiene 10 variabili:

  1. Anni.madre in anni: variabile quantitativa discreta.
  2. N.gravidanze: variabile quantitativa discreta.
  3. Fumatrici: variabile qualitativa (Non Fum = 0, Fum = 1).
  4. Gestazione: durata in settimane della gestazione, variabile quantitativa discreta.
  5. Peso: variabile quantitativa discreta, sarà la nostra variabile target.
  6. Lunghezza: variabile quantitativa discreta.
  7. Cranio: diametro del cranio, variabile quantitativa discreta.
  8. Tipo.parto: variabile qualitativa (Naturale = Nat, Cesareo = Ces).
  9. Ospedale: variabile qualitatva (osp1, osp2, osp3)
  10. Sesso: variabile qualitativa (maschio = M, femmina = F)

L’obiettivo di questa analisi è individuare quali di queste variabili sono più predittive del peso della nascita, con particolare attenzione all’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature.

Per coerenza con le altre variabili qualitative trasformiamo anche Fumatrici in un Factor:

dati_neonati$Fumatrici <- factor(dati_neonati$Fumatrici, 
                                 levels = c(0, 1), 
                                 labels = c("NF", "F"))

2. Analisi e Modellizzazione

Cerchiamo eventuali dati anomali (errori d’inserimento) nell’età della madre, facciamo un test vedendo se ci sono età della madre al di sotto dei 10 anni e nel caso eliminiamo quella serie di dati:

kable(dati_neonati[dati_neonati$Anni.madre<10, ])
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
1152 1 1 NF 41 3250 490 350 Nat osp2 F
1380 0 0 NF 39 3060 490 330 Nat osp3 M 
dati_neonati <- subset(dati_neonati, Anni.madre > 9)

Abbiamo eliminato dal dataset le due serie di dati anomali

2.1 Analisi Prelimininare

Procediamo con un’analisi descrittiva delle nostre variabile per comprenderne le loro distribuzioni. Procediamo determinando gli indici di posizione, di variabilità e di forma per le nostre variabile: Anni.madre, N.gravidanze, Gestazione, Peso, Lunghezza, Cranio

# Definizione funzione per il calcolo del coefficiente di variazione
CV <- function(x){
  return (sd(x)/mean(x))*100
}

# Lista delle colonne da analizzare
columns <- c("Anni.madre", "N.gravidanze", "Gestazione", "Peso", "Lunghezza", "Cranio")

# Creazione di una lista per salvare i risultati
results <- list()

# Ciclo sulle colonne
for (col in columns) {
  
  data <- dati_neonati[[col]] 
    
  # Calcola le statistiche base
  stats <- round(summary(data),2)  
    
  # Calcola ulteriori statistiche
  iqr <- round(IQR(data),2)  
  varianza <- round(var(data),2) 
  sd <- round(sd(data),2)  
  cv <- round(CV(data),2)
  skew <- round(skewness(data),2)  
  kurt <- round(kurtosis(data),2)-3  
    
    # Combina tutte le statistiche in una riga
  stats_row <- c(stats, IQR = iqr, Varianza = varianza, SD = sd, CV = cv, Skewness = skew, Kurtosis = kurt)
    
   
  results[[col]] <- as.data.frame(t(stats_row))

}

# Combina i risultati in un unico data.frame
results_table <- do.call(rbind, results)
results_table <- cbind(Variable = rownames(results_table), results_table)  
rownames(results_table) <- NULL  

# Aggiorna i nomi delle colonne per includere le nuove statistiche
colnames(results_table) <- c("Variable", "Min", "1st Qu.", "Median", "Mean", "3rd Qu.", "Max", 
                             "IQR", "Variance", "SD", "CV", "Skewness", "Kurtosis")

kable(results_table)
Variable Min 1st Qu. Median Mean 3rd Qu. Max IQR Variance SD CV Skewness Kurtosis
Anni.madre 13 25 28 28.19 32 46 7 27.22 5.22 0.19 0.15 -0.11
N.gravidanze 0 0 1 0.98 1 12 1 1.64 1.28 1.30 2.51 10.98
Gestazione 25 38 39 38.98 40 43 2 3.49 1.87 0.05 -2.07 8.26
Peso 830 2990 3300 3284.18 3620 4930 630 275865.90 525.23 0.16 -0.65 2.03
Lunghezza 310 480 500 494.70 510 565 30 693.21 26.33 0.05 -1.51 6.48
Cranio 235 330 340 340.03 350 390 20 269.93 16.43 0.05 -0.79 2.94

Guardando questa tabella possiamo fare alcune osservazioni iniziali:

  1. Anni.madre: presenta una leggera asimmetria positiva e leggermente platicurtica.
  2. N.gravidanze: come ci si poteva aspettare presenta un’assimmetria positiva ed è fortemente leptocurtica.
  3. Gestazione: tra i dati sono presenti alcuni parti prematuri, inoltre la distribuzione presenta una asimmetria negativa, seppur la media è molto simile alla mediana. Questo indica che è più probabile un parto prematuro di un parto extra lungo. La distribuzione è leptocurtica.
  4. Peso, Lunghezza, Cranio: presentano caratteristiche simili: asimmetria negativa e distribuzione leptocurtica.

Osserviamo queste distribuzioni anche attraverso istogrammi e boxplot:

# Definizione per il calcolo del numero di bin secondo il criterio di Rice
rice_bins <- function(data) {
  n <- length(data) # Numero di osservazioni
  return(2 * n^(1/3))
}

#Usiamo attach per accedere più facilmente alla variabili, dovendo costruire molteplici grafici

attach(dati_neonati)

# Lista per salvare i grafici
plots <- list()

# Ciclo per generare gli istogrammi
for (col in columns) {
  
  data <- dati_neonati[[col]]
  
  # Calcola il numero di bin e il binwidth
  if(!(col %in% c("Anni.madre", "Gestazione","N.gravidanze"))){
    bins <- rice_bins(data)
    binwidth <- diff(range(data, na.rm = TRUE)) / bins
  
  # Crea l'istogramma
    plot <- ggplot(dati_neonati) +
    geom_histogram(aes_string(x = col),
                   binwidth = binwidth,
                   col = "black",
                   fill = "lightblue") +
    labs(title = paste("HS per", col),
         x = col,
         y = "Count") +
    theme_minimal()
  }
  
  else{
    plot <- ggplot(dati_neonati) +
    geom_bar(aes_string(x = col),
                   col = "black",
                   fill = "lightblue") +
    labs(title = paste("HS per", col),
         x = col,
         y = "Count") +
    theme_minimal()
  }
  
  # Aggiungi il grafico alla lista
  plots[[col]] <- plot
}

# Disposizione dei grafici in una griglia 2x3
final_plot <- (plots[[1]] | plots[[2]] | plots[[3]]) /
              (plots[[4]] | plots[[5]] | plots[[6]])

# Mostra i grafici
print(final_plot)

plots <- list()

# Ciclo per generare i boxplot
for (col in columns) {
  
  data <- dati_neonati[[col]]
  
  plot <- ggplot(dati_neonati) +
    geom_boxplot(aes_string(y = col),
                   col = "black",
                   fill = "lightblue") +
    labs(title = paste("BP per", col),
         y = col) +
    theme_minimal()+
    theme(axis.text.x = element_blank(),  # Rimuove i numeri dall'asse X
          axis.ticks.x = element_blank())
  
  # Aggiungi il grafico alla lista
  plots[[col]] <- plot
}

# Disposizione dei grafici in una griglia 2x3
final_plot <- (plots[[1]] | plots[[2]] | plots[[3]]) /
              (plots[[4]] | plots[[5]] | plots[[6]])

# Mostra i grafici
print(final_plot)

Dai nostri grafici osserviamo la presenza di molti outlier: nella maggior parte dei casi con valori più piccoli (generano distribuzioni con asimmetria negativa) probabilmente legati a parti prematuri o di neonati particolarmente piccoli. Andiamo adesso ad analizzare le variabili qualitative: Fumatrici, Tipo.parto, Ospedale, Sesso per vedere se il nostro dataset è bilanciato rispetto a tutte le modalità oppure no.

colums_qual <- c("Fumatrici", "Tipo.parto", "Ospedale", "Sesso")

for (col in colums_qual) {
  # Calcolo delle frequenze relative
  tabella <- as.data.frame(prop.table(table(dati_neonati[[col]])) * 100)
  
  # Rinominazione delle colonne
  colnames(tabella) <- c(col, "Freq (%)")
  
  # Generazione della tabella con kable
  cat(
    kable(tabella, format = "html", digits = 2) %>%  # `digits = 2` per mostrare solo due decimali
      kable_styling(full_width = FALSE) %>%
      column_spec(1, width = "5cm") %>%
      column_spec(2, width = "5cm"),
    sep = "\n"
  )
}
Fumatrici Freq (%)
NF 95.84
F 4.16
Tipo.parto Freq (%)
Ces 29.14
Nat 70.86
Ospedale Freq (%)
osp1 32.67
osp2 33.95
osp3 33.39
Sesso Freq (%)
F 50.24
M 49.76

Possiamo osservare che le modalità di Fumatrici e Tipo.parto non sono bilanciate: le fumatrici nel dataset rappresentano solo il 4% e i parti cesarei circa il 30%. Per quanto riguarda gli ospedali coinvolti e il sesso il dataset risulta bilanciato.

2.1.A Parti Cesarei

Vogliamo verificare se in alcuni ospedali si eseguono più parti cesarei:

tabella <- dati_neonati %>%
  group_by(Ospedale) %>%
  summarise(
    Num.parti_cesarei = sum(Tipo.parto == 'Ces'),
    Num.parti_naturali = sum(Tipo.parto == 'Nat'),
    Perc.parti_cesarei = 100*sum(Tipo.parto == 'Ces')/n()
            )

kable(tabella, format = "html", digits = 2) %>%
  kable_styling(full_width = FALSE)
Ospedale Num.parti_cesarei Num.parti_naturali Perc.parti_cesarei
osp1 242 574 29.66
osp2 254 594 29.95
osp3 232 602 27.82 
# Creare una tabella di contingenza
contingenza <- table(dati_neonati$Ospedale, dati_neonati$Tipo.parto)

# Test chi-quadrato
risultato<-chisq.test(contingenza)
cat(glue("Risultati del Test Chi-Quadrato: \n
- Statistiche del test: {round(risultato$statistic, 2)} \n
- Gradi di libertà: {risultato$parameter} \n
- P-value: {round(risultato$p.value, 2)}"))

Risultati del Test Chi-Quadrato:

  • Statistiche del test: 1.08

  • Gradi di libertà: 2

  • P-value: 0.58

Con un P-value di 0.58 non rifiutiamo l’ipotesi nulla e quindi non c’è evidenza statistica che in un ospedale si eseguano più parti cesarei che in un altro.

2.1.B Peso e lunghezza dei bambini del campione rispetto alla popolazione

Siamo interessati a verifica se il nostro campione è rappresentativo rispetto alla popolazione. Prendiamo i dati dal sito dell’University of Rochester Medical Center: Newborn Measurements

A baby’s birth weight is an important indicator of health. The average weight for full-term babies (born between 37 and 41 weeks gestation) is about 7 pounds (3.2 kg).

Per quanto riguarda il peso abbiamo un valore medio di 3200 grammi per bambini a termine, vediamo con un boxplot se osserviamo una differenza tra bambini nati a termine, prematuri e oltre il termine:

# Creiamo una nuova variabile Categoria_Gestazione con tre categorie
dati_neonati$Categoria_Gestazione <- ifelse(dati_neonati$Gestazione < 37, 
                                              "Prematuri (<37 settimane)", 
                                              ifelse(dati_neonati$Gestazione > 41, 
                                                     "Oltre il termine (>41 settimane)", 
                                                     "A termine (37-41 settimane)"))

# Creiamo il boxplot
ggplot(dati_neonati) +
  geom_boxplot(aes(y = Peso, color = Categoria_Gestazione), fill = "lightblue") +
  geom_hline(aes(yintercept = 3200, color = "Valore di riferimento (3200g)"), linetype = "dashed", size = 1) +
  labs(title = "Boxplot per Peso e Categoria di Gestazione",
       y = "Peso (g)",
       color = "Legenda") +
  scale_color_manual(values = c("Prematuri (<37 settimane)" = "#F8766D",
                                 "A termine (37-41 settimane)" = "#00BA38",
                                 "Oltre il termine (>41 settimane)" = "#619CFF",
                                 "Valore di riferimento (3200 g)" = "black")) +
  theme_minimal() +
  theme(legend.position = "right")+
  theme(axis.text.x = element_blank(),  # Rimuove i numeri dall'asse X
        axis.ticks.x = element_blank())

Anche dal nostro boxplot osserviamo una differenza tra le distribuzioni del peso nelle tre categorie, verifichiamo se il campione dei neonati a termine è rappresentativo della popolazione:

format_p_value <- function(p) {
  if (p < 0.00000001) {
    return("< 0.00000001")
  } else {
    return(format(p, digits = 8))
  }
}
mu0 = 3200
peso_neonati_37_41 <- Peso[dati_neonati$Categoria_Gestazione == "A termine (37-41 settimane)"]

result <- t.test(peso_neonati_37_41,
       mu = mu0,
       conf.level = 0.95,
       alternative = "two.sided")

cat(glue("
Risultati del Test t-Student:

- Media osservata: {round(result$estimate, 0)}
- Ipotesi nulla (media attesa): {round(mu0, 2)}
- Statistica t: {round(result$statistic, 2)}
- Gradi di libertà: {round(result$parameter, 2)}
- P-value: {format_p_value(result$p.value)}
- Intervallo di confidenza al 95%: ({round(result$conf.int[1], 0)}, {round(result$conf.int[2], 0)})
"))

Risultati del Test t-Student:

  • Media osservata: 3343
  • Ipotesi nulla (media attesa): 3200
  • Statistica t: 15.45
  • Gradi di libertà: 2277
  • P-value: < 0.00000001
  • Intervallo di confidenza al 95%: (3325, 3361)

La nostra distribuzione non è compatibile con il valore trovato in letteratura, infatti il valore di riferimento (3200 g) è al di fuori dell’intervallo di confidenza. Svolgiamo passaggi analoghi per la lunghezza. Dalla stessa fonte sappiamo che la lunghezza media di un neonato nato a termine è 50 cm.

ggplot(dati_neonati) +
  geom_boxplot(aes(y = Lunghezza,  color = Categoria_Gestazione), fill = "lightblue") +
  geom_hline(aes(yintercept = 500, color = "Valore di riferimento (3200g)"), linetype = "dashed", size = 1)+
  labs(title = "Boxplot per Lunghezza e Categoria di Gestazione",
       y = "Lunghezza (mm)",
       color = "Legenda") +
  scale_color_manual(values = c("Prematuri (<37 settimane)" = "#F8766D",
                                 "A termine (37-41 settimane)" = "#00BA38",
                                 "Oltre il termine (>41 settimane)" = "#619CFF",
                                 "Valore di riferimento (500 mm)" = "black")) +
  theme_minimal() +
  theme(legend.position = "right")+
  theme(axis.text.x = element_blank(),  # Rimuove i numeri dall'asse X
        axis.ticks.x = element_blank()) 

Lu0 = 500
Lunghezza_37_41 <- Lunghezza[dati_neonati$Categoria_Gestazione == "A termine (37-41 settimane)"]

result_L <- t.test(Lunghezza_37_41,
       mu = Lu0,
       conf.level = 0.95,
       alternative = "two.sided")

cat(glue("
Risultati del Test t-Student:

- Media osservata: {round(result_L$estimate, 0)}
- Ipotesi nulla (media attesa): {round(Lu0, 0)}
- Statistica t: {round(result_L$statistic, 2)}
- Gradi di libertà: {round(result_L$parameter, 2)}
- P-value: {format_p_value(result$p.value)}
- Intervallo di confidenza al 95%: ({round(result_L$conf.int[1], 0)}, {round(result_L$conf.int[2], 0)})
"))

Risultati del Test t-Student:

  • Media osservata: 498
  • Ipotesi nulla (media attesa): 500
  • Statistica t: -5.07
  • Gradi di libertà: 2277
  • P-value: < 0.00000001
  • Intervallo di confidenza al 95%: (497, 499)

Tutti e due i test danno esito negativo per quanto riguarda la rappresentanza del nostro campione rispetto alla popolazione.

2.1.C Misure antropometriche tra i due sessi

Vogliamo adesso valutare se le misure antropometriche (peso, lunghezza e diametro del cranio) mostrano delle differenze tra i due sessi, iniziamo a confrontare i loro boxplot per poi passare a svolgete il t test fra i due gruppi:

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

plots <- list()

# Ciclo per generare i boxplot
for (col in columns) {
  
  # Crea il boxplot
  plot <- ggplot(dati_neonati) +
    geom_boxplot(aes_string(x = "Sesso", y = col, color = "Sesso"),
                 fill = "lightblue") +
    labs(title = paste("Boxplot per", col),
         y = col,
         x = "Sesso",
         color = "Sesso") +
    theme_minimal() +
    theme(legend.position = "bottom",
          axis.text.x = element_text(angle = 45, hjust = 1))
  
  # Aggiungi il grafico alla lista
  plots[[col]] <- plot
}

final_plot <- (plots[[1]] | plots[[2]] | plots[[3]])

# Mostra i grafici
print(final_plot)

I boxplot mostrano delle differenze tra le due categorie ma procediamo adesso con il t.test per vedere se c’è una dipendenza dal sesso per quanto riguarda le grandezze antropometiche:

tabella_contingenza <- dati_neonati %>%
  group_by(Sesso) %>%
  summarise(Peso_medio = mean(Peso),
            Lun_media = mean(Lunghezza),
            Di_Cr_medio = mean(Cranio))

kable(tabella_contingenza, format = "html", digits = 2) %>%
  kable_styling(full_width = FALSE)
Sesso Peso_medio Lun_media Di_Cr_medio
F 3161.06 489.76 337.62
M 3408.50 499.67 342.46 
t_peso <- t.test(Peso ~ Sesso, data = dati_neonati)
t_lunghezza <- t.test(Lunghezza ~ Sesso, data = dati_neonati)
t_cranio <- t.test(Cranio ~ Sesso, data = dati_neonati)

# Estrai i risultati principali
risultati <- data.frame(
  Variabile = c("Peso", "Lunghezza", "Cranio"),
  `T-Statistic` = c(round(t_peso$statistic,2), round(t_lunghezza$statistic,2), round(t_cranio$statistic,2)),
  `Gradi di Libertà` = c(round(t_peso$parameter,0), round(t_lunghezza$parameter,0), round(t_cranio$parameter,0)),
  `P-Value` = sapply(c(t_peso$p.value, t_lunghezza$p.value, t_cranio$p.value), format_p_value)
)

# Mostra la tabella in un formato HTML con knitr::kable
kable(risultati, format = "html", digits = 4, align = "c") %>%
  kableExtra::kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))
Variabile T.Statistic Gradi.di.Libertà P.Value
Peso -12.11 2489 < 0.00000001
Lunghezza -9.58 2457 < 0.00000001
Cranio -7.44 2489 < 0.00000001

I risultati ottenuti per tutti i test mostrano che possiamo scartare in tutti i casi l’ipotesi nulla, quindi il peso, la lunghezza e il diametro del cranio mostrano una dipendenza dal sesso.

2.2 Creazione del Modello di Regressione

Completato lo studio delle diverse variabili, possiamo passare allo studio per la costruzione del modello che meglio possa descrivere i nostri dati. Abbiamo già visto che la nostra variabile target si presenta con una asimmetria negativa e con una tendenza leptocurtica, eseguiamo il test di Shapiro-Wilk per verificare l’ipotesi di normalità.

dati_neonati <- dati_neonati %>% select(-Categoria_Gestazione) #Eliminiamo la categoria ausiliaria che avevamo aggiunto precedentemente
n <- nrow(Peso)
ggplot(dati_neonati, aes(x = Peso)) +
  geom_density(fill = "lightblue", alpha = 0.6) +
  labs(title = "Grafico di Densità per Peso",
       x = "Peso (g)",
       y = "Densità") +
  theme_minimal()

# Esegui il test di Shapiro-Wilk
shapiro_result <- shapiro.test(Peso)

# Stampa dei risultati con format.pval
cat(glue::glue("Risultati del Test di Shapiro-Wilk: \n
- W: {round(shapiro_result$statistic, 3)} \n
- P-value: {format_p_value(shapiro_result$p.value)}"))

Risultati del Test di Shapiro-Wilk:

  • W: 0.971

  • P-value: < 0.00000001

Dato il risultato del test rifiutiamo l’ipotesi di distribuzione normale, vedremo se questo avrà un impatto sull’analisi dei nostri residui. Studiamo la correlazione tra i regressori tra di loro e in particolare con la variabile target.

dati_neonati_num <- dati_neonati %>% select_if(is.numeric)

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

Sia dai valori delle correlazioni che dagli scatterplot sembrerebbe che le correlazioni del Peso più importanti siano con Gestazione, Lunghezza e Cranio, si vedono anche dei possibili effetti non di primo grado.

Analizziamo attraverso boxplot la dipendenza del Peso dai regressori qualitativi Fumatrici, Tipo.parto e Ospedale, abbiamo già verificato come ci sia una dipendenza del peso dal sesso.

columns_qual <- c("Fumatrici", "Tipo.parto", "Ospedale")

plots<-list()

for (col in columns_qual){
  dati_neonati[[col]] <- as.factor(dati_neonati[[col]])
  plot <- ggplot(dati_neonati) +
    geom_boxplot(aes_string(x = col, y = Peso, color = col),
                 fill = "lightblue") +
    labs(title = paste("Boxplot Peso per", col),
         y = "Peso (g)",
         x = col,
         color = col) +
    theme_minimal() +
    theme(legend.position = "right")+
    theme(axis.text.x = element_blank(), 
        axis.ticks.x = element_blank()) 
  
  # Aggiungi il grafico alla lista
  plots[[col]] <- plot
}

final_plot <- (plots[[1]]) / 
              (plots[[2]]) /
              (plots[[3]])

# Mostra i grafici
print(final_plot)

I boxplot non mostrano particolari dipendenze ma svolgiamo anche i corrispondenti test statistici:

# Esegui i test
t_fumatrici <- t.test(Peso ~ Fumatrici)
t_tipoparto <- t.test(Peso ~ Tipo.parto)
pairwise_ospedale <- pairwise.t.test(Peso, Ospedale, paired = FALSE, 
                                     pool.sd = TRUE, 
                                     p.adjust.method = 'bonferroni')

# Estrapola i risultati principali
risultati_test <- data.frame(
  Test = c("Peso ~ Fumatrici", "Peso ~ Tipo.parto"),
  Statistica = c(
    round(t_fumatrici$statistic, 2),
    round(t_tipoparto$statistic, 2)
  ),
  P_Value = c(
    round(t_fumatrici$p.value,2),
    round(t_tipoparto$p.value,2)
  )
)

# Visualizza la tabella in RMarkdown
kable(risultati_test, format = "html", digits = 3) %>%
  kable_styling(full_width = FALSE, bootstrap_options = "striped")
Test Statistica P_Value
Peso ~ Fumatrici 1.04 0.30
Peso ~ Tipo.parto -0.14 0.89 
pairwise_results <- as.data.frame(pairwise_ospedale$p.value)
pairwise_results <- tibble::rownames_to_column(pairwise_results, "Ospedale")

kable(pairwise_results, format = "html", digits = 2) %>%
  kable_styling(full_width = FALSE, bootstrap_options = "striped")
Ospedale osp1 osp2
osp2 1.00 NA
osp3 0.33 0.32

I test statistici ci dicono che non è possibile scartare l’ipotesi nulla, non ci sono differenze significative nel peso medio tra fumatrici e non, tra i vari tipi di parto e tra i neonati dei tre ospedali.

Fatte queste analisi procediamo con la costruzione del modello di regressione, ricordiamo che la richiesta che ci è stata fatta è di porre un focus particolare sull’impatto del fumo materno (vedremo se il modello confermerà, come ci aspettiamo, l’indipendenza del peso dal fumo) e delle settimane di gestazione.

#Funzione per gestire la stampa formattata dei nostri modelli
print_summary <- function(mod) {
  # Estrazione dei coefficienti
  coefficients <- summary(mod)$coefficients
  
  # Separare i p-value e formattarli
  p_values <- coefficients[, 4]
  formatted_p_values <- sapply(p_values, format_p_value)
  
  # Arrotondare le altre quantità della tabella a 2 cifre decimali
  coeff_table <- data.frame(
    Estimate = round(coefficients[, 1], 2),
    Std.Error = round(coefficients[, 2], 2),
    t_value = round(coefficients[, 3], 2),
    `P-value` = formatted_p_values
  )
  
  # Stampare la tabella
  print(kable(coeff_table, format = "html") %>%
        kable_styling(full_width = FALSE, bootstrap_options = "striped"))
  
  # Calcolare l'Adjusted R-squared
  adj_r_squared <- summary(mod)$adj.r.squared
  
  # Stampare l'Adjusted R-squared con il messaggio
  cat("L'Adjusted R-squared per questo modello è:", format(adj_r_squared, digits = 2), "\n")
}
mod1 <- lm(Peso ~ . , data=dati_neonati)
# Calcolare l'Adjusted R-squared
print_summary(mod1)
Estimate Std.Error t_value P.value
(Intercept) -6735.80 141.48 -47.61 < 0.00000001
Anni.madre 0.80 1.15 0.70 0.48448612
N.gravidanze 11.38 4.67 2.44 0.01484575
FumatriciF -30.27 27.55 -1.10 0.27191274
Gestazione 32.58 3.82 8.53 < 0.00000001
Lunghezza 10.29 0.30 34.21 < 0.00000001
Cranio 10.47 0.43 24.57 < 0.00000001
Tipo.partoNat 29.63 12.09 2.45 0.014315418
Ospedaleosp2 -11.09 13.45 -0.82 0.40956172
Ospedaleosp3 28.25 13.51 2.09 0.03656524
SessoM 77.57 11.19 6.93 < 0.00000001

L’Adjusted R-squared per questo modello è: 0.73

In questo modello abbiamo considerato tutte le variabili del nostro dataset originario usando Peso come variabile target. I primi risultati ci dicono che Anni.madre è il regressore meno correlato con la variabile target, proviamo a toglierla tenendo come riferimento il risultato dell’Adjusted R-squared di 0.73.

mod2 <- update(mod1, ~.-Anni.madre)
print_summary(mod2)
Estimate Std.Error t_value P.value
(Intercept) -6708.62 136.02 -49.32 < 0.00000001
N.gravidanze 12.58 4.34 2.90 0.0037719398
FumatriciF -30.43 27.55 -1.10 0.2694382
Gestazione 32.30 3.80 8.50 < 0.00000001
Lunghezza 10.29 0.30 34.21 < 0.00000001
Cranio 10.49 0.43 24.64 < 0.00000001
Tipo.partoNat 29.67 12.09 2.45 0.014200587
Ospedaleosp2 -10.95 13.44 -0.81 0.41541031
Ospedaleosp3 28.52 13.50 2.11 0.034735333
SessoM 77.65 11.18 6.94 < 0.00000001

L’Adjusted R-squared per questo modello è: 0.73

L’Adjusted R-squared non è cambiato, quindi effettivamente la variabile Anni.madre non aggiungeva varianza spiegata al modello stesso, eliminiamo la variabile Ospedale che solo nel caso dell’ospedale 3 mostra un p-value inferiore al 5% ma ancora molto alto.

mod3 <- update(mod2, ~.-Ospedale)
print_summary(mod3)
Estimate Std.Error t_value P.value
(Intercept) -6708.81 136.06 -49.31 < 0.00000001
N.gravidanze 12.99 4.34 2.99 0.0028077677
FumatriciF -31.88 27.58 -1.16 0.24779972
Gestazione 32.60 3.80 8.57 < 0.00000001
Lunghezza 10.27 0.30 34.10 < 0.00000001
Cranio 10.50 0.43 24.64 < 0.00000001
Tipo.partoNat 30.42 12.10 2.51 0.012014373
SessoM 78.10 11.20 6.97 < 0.00000001

L’Adjusted R-squared per questo modello è: 0.73

L’Adjusted R-squared non è cambiato. Il risultato del t.test sul Peso e sulle mamme Fumatrici ci aveva mostrato che non era possibile rifiutare l’ipotesi nulla, cioè l’assenza di correlazione, proviamo, quindi, ad eliminare dal nostro modello questo regressore, anche se è una delle variabili oggetto di studio:

mod4 <- update(mod3, ~.-Fumatrici)
print_summary(mod4)
Estimate Std.Error t_value P.value
(Intercept) -6708.02 136.07 -49.30 < 0.00000001
N.gravidanze 12.74 4.34 2.94 0.0033607508
Gestazione 32.33 3.80 8.51 < 0.00000001
Lunghezza 10.28 0.30 34.18 < 0.00000001
Cranio 10.51 0.43 24.65 < 0.00000001
Tipo.partoNat 30.16 12.10 2.49 0.012766914
SessoM 77.92 11.20 6.96 < 0.00000001

L’Adjusted R-squared per questo modello è: 0.73

L’Adjusted R-squared è sempre 0.73, quindi concludiamo che dal nostro dataset l’essere fumatrice o meno da parte della madre non ha avuto influenze sul peso del bambino. Proviamo ad eliminare dal nostro modello il regressore Tipo.parto in quando è quello che presenta il p-value più alto, seppur minore del 5%. Riflettendo sul significato di questa variabile non sembra avere un impatto sulla formazione e lo sviluppo prenatale del neonato stesso.

mod5 <- update(mod4, ~.-Tipo.parto)
print_summary(mod5)
Estimate Std.Error t_value P.value
(Intercept) -6681.73 135.80 -49.20 < 0.00000001
N.gravidanze 12.46 4.34 2.87 0.0041540727
Gestazione 32.38 3.80 8.52 < 0.00000001
Lunghezza 10.25 0.30 34.06 < 0.00000001
Cranio 10.54 0.43 24.72 < 0.00000001
SessoM 77.98 11.21 6.96 < 0.00000001

L’Adjusted R-squared per questo modello è: 0.73

Anche togliendo questa variabile l’adjusted R-squared è rimasto pari a 0.73. Per motivi analoghi proviamo a eliminare anche il regressore N.gravidanze che mostra un p-value molto maggiore degli altri, applicando il principio per cui un modello più semplice permette una maggiore interpretabilità:

mod6 <- update(mod5, ~.-N.gravidanze)
print_summary(mod6)
Estimate Std.Error t_value P.value
(Intercept) -6651.67 135.60 -49.06 < 0.00000001
Gestazione 31.33 3.79 8.27 < 0.00000001
Lunghezza 10.20 0.30 33.91 < 0.00000001
Cranio 10.67 0.42 25.13 < 0.00000001
SessoM 79.10 11.22 7.05 < 0.00000001

L’Adjusted R-squared per questo modello è: 0.73

Notiamo come l’adjusted R-squared non ha subito modifiche.

Vale la pena testare se ci sono delle collinearità tra le variabile rimaste, calcolando il Variance Inflation Factor che ci darà un risultato minore di 5 in assenza di collinearità:

kable(vif(mod6), format = "html") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "condensed"))
x
Gestazione 1.654101
Lunghezza 2.070582
Cranio 1.606316
Sesso 1.038918

Il test mostra che tra i regressori selezionati non ci sono collinearità che possono influenzare in negativo il nostro modello. A partire da queste variabile verifichiamo se il nostro modello dev’essere arricchito con effetti non lineari e interazion tra le variabili.

2.3 Selezione del modello ottimale

Verifichiamo una prima ipotesi: dobbiamo introdurre effetti non lineari per i regressosi Gestazione, Lunghezza e Cranio? Gli scatterplot mostrano un possibile effetto quadratico

mod7L <- update(mod6, ~.+I(Lunghezza^2))
print_summary(mod7L)
Estimate Std.Error t_value P.value
(Intercept) 153.60 725.07 0.21 0.83225127
Gestazione 41.22 3.86 10.67 < 0.00000001
Lunghezza -19.91 3.17 -6.29 < 0.00000001
Cranio 10.80 0.42 25.87 < 0.00000001
SessoM 71.34 11.05 6.45 < 0.00000001
I(Lunghezza^2) 0.03 0.00 9.55 < 0.00000001

L’Adjusted R-squared per questo modello è: 0.74

mod7C <- update(mod6, ~.+I(Cranio^2))
print_summary(mod7C)
Estimate Std.Error t_value P.value
(Intercept) 74.00 1153.56 0.06 0.94885607
Gestazione 37.78 3.92 9.64 < 0.00000001
Lunghezza 10.44 0.30 34.62 < 0.00000001
Cranio -31.40 7.18 -4.37 1.269182e-05
SessoM 74.28 11.18 6.65 < 0.00000001
I(Cranio^2) 0.06 0.01 5.87 < 0.00000001

L’Adjusted R-squared per questo modello è: 0.73

mod7G <- update(mod6, ~.+I(Gestazione^2))
print_summary(mod7G)
Estimate Std.Error t_value P.value
(Intercept) -4640.60 899.96 -5.16 2.7142655e-07
Gestazione -80.95 49.81 -1.62 0.10429162
Lunghezza 10.31 0.30 33.89 < 0.00000001
Cranio 10.77 0.43 25.25 < 0.00000001
SessoM 76.91 11.25 6.83 < 0.00000001
I(Gestazione^2) 1.50 0.66 2.26 0.023882927

L’Adjusted R-squared per questo modello è: 0.73

Per decidere quale di questi tre modelli confrontare col modello 5, eseguiamo anche il test AIC e BIC per selezionare l’effetto più importante:

# Tabella AIC
kable(AIC(mod7C, mod7G, mod7L), format = "html") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "condensed"))
df AIC
mod7C 7 35126.81
mod7G 7 35156.01
mod7L 7 35071.37 
# Tabella BIC
kable(BIC(mod7C, mod7G, mod7L), format = "html") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "condensed"))
df BIC
mod7C 7 35167.58
mod7G 7 35196.77
mod7L 7 35112.13

Per entrambi i test il mod7L dove abbiamo introdotto la lunghezza al quadrato è quello col valore minore, quindi quello da preferire in ancora con i risulati dell’adjusted R-squared, facciamo anche il test ANOVA per la significatività della varianza spiegata rispetto al mod6.

# Esegui ANOVA e salva i risultati
anova_results <- anova(mod7L, mod6)

# Estrai il p-value del confronto tra i modelli
p_value <- anova_results$`Pr(>F)`[2] 
formatted_p_value <- format_p_value(p_value)

# Stampa il risultato con un messaggio
cat("Il p-value del confronto tra i modelli è:", formatted_p_value, "\n")

Il p-value del confronto tra i modelli è: < 0.00000001

Anche questo test conferma che il mod7L è in grado di spiegare maggiormente i dati, prima di decidere se sceglieremo o meno il mod7L, verifichiamo se è necessario introdurre l’interazione tra le variabile Lunghezza e Cranio nel nostro modello. Considerando le variabili rimaste queste sono le uniche per cui è ragionevole considerare un effetto combinato:

mod8A <-update(mod7L, ~.+Lunghezza*Cranio)
print_summary(mod8A)
Estimate Std.Error t_value P.value
(Intercept) -2299.20 1005.07 -2.29 0.022244732
Gestazione 39.14 3.90 10.04 < 0.00000001
Lunghezza -20.99 3.17 -6.61 < 0.00000001
Cranio 27.38 4.74 5.78 < 0.00000001
SessoM 73.20 11.04 6.63 < 0.00000001
I(Lunghezza^2) 0.04 0.00 8.98 < 0.00000001
Lunghezza:Cranio -0.03 0.01 -3.52 0.00044696936

L’Adjusted R-squared per questo modello è: 0.74

mod8B <- update(mod6, ~.+Lunghezza*Cranio)
print_summary(mod8B)
Estimate Std.Error t_value P.value
(Intercept) -1840.74 1019.67 -1.81 0.071160332
Gestazione 36.97 3.95 9.35 < 0.00000001
Lunghezza -0.20 2.21 -0.09 0.92716773
Cranio -4.40 3.20 -1.38 0.16820293
SessoM 74.47 11.21 6.64 < 0.00000001
Lunghezza:Cranio 0.03 0.01 4.76 2.0464178e-06

L’Adjusted R-squared per questo modello è: 0.73

Osserviamo che aver introdotto questa interazione non ha modificato il nostro adjusted R-squared (rimasto a 0.74) quindi decidiamo di scartarla. Selezioniamo il nostro modello tra quelli creati (da mod1 a mod6 più mod7L) attraverso il criterio di Informazione di Akaike (AIC) e quello di Bayes (BIC):

# Tabella AIC
kable(AIC(mod1,mod2,mod3,mod4,mod5,mod6,mod7L), format = "html") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "condensed"))
df AIC
mod1 12 35145.57
mod2 11 35144.06
mod3 9 35149.33
mod4 8 35148.67
mod5 7 35152.89
mod6 6 35159.12
mod7L 7 35071.37 
# Tabella BIC
kable(BIC(mod1,mod2,mod3,mod4,mod5,mod6,mod7L), format = "html") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "condensed"))
df BIC
mod1 12 35215.45
mod2 11 35208.12
mod3 9 35201.73
mod4 8 35195.25
mod5 7 35193.65
mod6 6 35194.06
mod7L 7 35112.13

Per entrambi i criteri il mod7L è quello che descrive al meglio i nostri dati, ristampiamo il summary per poter commentare i coefficienti dei nostri regressori:

print_summary(mod7L)
Estimate Std.Error t_value P.value
(Intercept) 153.60 725.07 0.21 0.83225127
Gestazione 41.22 3.86 10.67 < 0.00000001
Lunghezza -19.91 3.17 -6.29 < 0.00000001
Cranio 10.80 0.42 25.87 < 0.00000001
SessoM 71.34 11.05 6.45 < 0.00000001
I(Lunghezza^2) 0.03 0.00 9.55 < 0.00000001

L’Adjusted R-squared per questo modello è: 0.74

  1. Il Peso aumenta di circa 41 grammi per ogni settimana di Gestazione.
  2. Per quanto riguarda la Lunghezza abbiamo un doppio effetto che va spiegato insieme: il primo regressore è -20 che indicherebbe una diminuizione di 20 grammi per ogni millimetro ma viene compensato dallo 0.03 per la lunghezza al quadrato che risulta dominante per valori oltre 333 mm. Quindi complessivamente il Peso cresce in modo quadratico rispetto alla Lunghezza.
  3. Per ogni millimetro in più nel diametro del Cranio il Peso cresce di circa 11 grammi.
  4. I neonati maschi avranno un Peso medio maggiore di circa 71 grammi rispetto alle femmine.

2.2.C Analisi della qualità del modello

Verifichiamo la presenza o meno di collinearità per il modello scelto:

kable(round(vif(mod7L),2), format = "html") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "condensed"))
x
Gestazione 1.78
Lunghezza 237.72
Cranio 1.61
Sesso 1.04
I(Lunghezza^2) 229.68

Il risultato non è preoccupante perché è normale che Lunghezza e Lunghezza al quadrato mostrino una collinearità, mentre le altre grandezze mostrano valori perfettamente nella norma.

Valutiamo la qualità del nostro modello attraverso alcune metrice:

  1. La prima che usiamo è quello dell’\[R^2_{adj}=0.74\].
  2. Valutiamo anche il Root Mean Squared Error (RMSE):
# Valori predetti
predicted <- predict(mod7L, dati_neonati)

# Valori osservati
observed <- dati_neonati$Peso

# RMSE
rmse_value <- rmse(observed, predicted)
cat("RMSE del modello mod7L:", round(rmse_value, 2))

RMSE del modello mod7L: 269.93

Il risultato del RMSE ci di che le nostre predizioni presenterano un’indeterminazione media di più o meno 270 g, circa l’8% rispetto al valore medio del peso del campione.

L’ultimo punto per valutare la qualità del nostro modello è l’analisi dei residui e la determinazione di valori influenti che potrebbero distorcere le nostre previsioni:

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

  1. Il grafico Residuals vs Fitted mostra che i dati sono sufficientemente equidistribuiti attorno allo zero con tre valori che si discostano più degli altri dalla media zero (155, 1306, 1551).
  2. Il grafico Normal Q-Q mostra delle discrepanze dalla nomale nelle code, in particolare per valori positivi (vengono segnalati gli stessi valori di prima).
  3. Nel grafico Scale-location la linea rossa che nel caso ideale dovrebbe essere orizzontale, mostra una leggera curvatura.
  4. Nel grafico Residual vs Leverage solo il punto 1551 mostra una distanza di Cook superiore alla distanza di allarme.

Dopo queste prime considerazioni vediamo più in dettaglio i nostri grafici con l’ausilio anche di test quantitativi:

Verifichiamo le caratteristiche della curva dei residui:

residui <- residuals(mod7L)

ggplot(data.frame(residui = residui), aes(x = residui)) +
  geom_density(fill = "blue", alpha = 0.3) +
  geom_vline(xintercept = 0, color = "red", linetype = "dashed") + 
  labs(title = "Distribuzione della densità dei residui",
       x = "Residui",
       y = "Densità") +
  theme_minimal()

shapiro <- shapiro.test(residuals(mod7L))
bp <- bptest(mod7L)
dw <- dwtest(mod7L)

# Creazione di un data frame per la tabella
test_results <- data.frame(
  Test = c("Shapiro-Wilk Test", "Breusch-Pagan Test", "Durbin-Watson Test"),
  Statistic = round(c(shapiro$statistic, bp$statistic, dw$statistic), 3),
  `P-Value` = sapply(c(shapiro$p.value, bp$p.value, dw$p.value), format_p_value)
)

# Stampa della tabella
kable(test_results, format = "html", col.names = c("Test", "Statistic", "P-Value")) %>%
  kable_styling(full_width = FALSE, bootstrap_options = "striped")
Test Statistic P-Value
W Shapiro-Wilk Test 0.986 < 0.00000001
BP Breusch-Pagan Test 126.080 < 0.00000001
DW Durbin-Watson Test 1.951 0.10831412

I nostri residui non soddisfano l’ipotesi di distribuzione normale (Shapito Test), mostrano eteroschedasticità, non hanno una varianza costante (Breusch-Pagan Test) ma non presentano autocorrelazione (Durbin-Watson test).

Cerchiamo i leverages e gli outliers:

lev <- hatvalues(mod7L)
p =sum(lev)
n=length(Peso)
soglia <- 2 * p / n
lev_data <- data.frame(Index = seq_along(lev), Leverage = lev)

ggplot(lev_data, aes(x = Index, y = Leverage)) +
  geom_point(color = "blue", size = 2) +  
  geom_hline(yintercept = soglia, color = "red", linetype = "dashed") +  
  labs(title = "Grafico dei valori di leverage", 
       x = "Indice delle osservazioni", 
       y = "Valori di leverage") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))  

length(lev[lev>soglia])

[1] 132

Tramite il grafico dei leverages e il filtro sulla soglia abbiamo individuato 132 valori leverages, cioè 132 dati che nello spazio dei regressori hanno valori molto distanti da quelli medi. Vediamo come si comportano gli outliers:

resid_student <- rstudent(mod7L)
resid_data <- data.frame(Index = seq_along(resid_student), Residuals = resid_student)

ggplot(resid_data, aes(x = Index, y = Residuals)) +
  geom_point(color = "blue", size = 2) +  # Punti per i residui studentizzati
  geom_hline(yintercept = c(-2, 2), color = "red", linetype = "dashed") +  # Linee soglia
  labs(title = "Grafico dei residui studentizzati", 
       x = "Indice delle osservazioni", 
       y = "Residui studentizzati") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))  # Centra il titolo

# Eseguire il test degli outlier
outliers <- outlierTest(mod7L)

# Controlla se ci sono risultati
if (!is.null(outliers)) {
  # Estrarre i dati rilevanti dall'oggetto outlierTest
  outliers_df <- data.frame(
    Row = names(outliers$rstudent), # Nomi delle righe (gli indici)
    `Bonferroni P` = sapply(outliers$bonf.p, format_p_value) # P-value di Bonferroni
  )
  
  # Stampare la tabella formattata senza duplicare i numeri
  kable(outliers_df, format = "html", digits = 3, row.names = FALSE) %>%
    kable_styling(full_width = FALSE, bootstrap_options = "striped")
} else {
  cat("Non sono stati rilevati outlier.")
}
Row Bonferroni.P
1551 < 0.00000001
1306 0.004660181
155 0.0082057635
1399 0.044339766
1694 0.045897464

Anche se il grafico dei residui studentizzati mostrerebbe più outliers, tramite l’outlierTest individuiamo quelli che sono statisticamente significativi e sono i seguenti: 1551, 1306, 155, 1399 e 1694 i primi 3 gli avevamo incontrati anche nei 4 scatterplot dei residui. Facciamo un ultimo test: quello della distanza di Cook per verifica se ci sono dati sperimentali che hanno un forte impatto sul nostro modello.

cook <- cooks.distance(mod7L)
plot(cook)

indice_max_cook <- which.max(cook)

indice_originale <- rownames(dati_neonati)[indice_max_cook]

# Estrai la riga corrispondente
riga <- dati_neonati[indice_max_cook, ]

# Crea un messaggio formattato con entrambi gli indici
messaggio <- glue(
  "Dati per la posizione corrente {indice_max_cook} (indice originale: {indice_originale}): \n
  Gestazione: {riga$Gestazione} \n
  Peso: {riga$Peso} \n
  Lunghezza: {riga$Lunghezza} \n
  Cranio: {riga$Cranio} \n
  Sesso: {riga$Sesso}"
)

# Stampa i dati
cat(messaggio)

Dati per la posizione corrente 1549 (indice originale: 1551):

Gestazione: 38

Peso: 4370

Lunghezza: 315

Cranio: 374

Sesso: F

Osserviamo che il dato che occupa la posizione 1549 è il dato 1551 (a causa dei due dati anomali cancellati all’inizio). Utilizzando la piattaforma presente nel sito Valutazione antropometrica neonatale che permette di calcolare i centili di appartenenza scopriamo che il neonato del dato 1551(numerazione del dataset originario) sarebbe al 100 centile per il Peso e allo 0 per la Lunghezza, addirittura non è previsto un Cranio con tale diametro, possiamo quindi concludere che questo dato è stato affetto da qualche errore durante la registrazione quindi può essere eliminato. Riaggiorniamo i parametri del modello una volta eliminata la misura, questo sarà il nostro modello definitivo:

dati_neonati <- dati_neonati[rownames(dati_neonati) != "1551", ]
mod7L <- update(mod7L, data = dati_neonati)
print_summary(mod7L)
Estimate Std.Error t_value P.value
(Intercept) -1666.39 760.26 -2.19 0.028480687
Gestazione 36.44 3.88 9.39 < 0.00000001
Lunghezza -11.37 3.35 -3.39 0.00069825441
Cranio 10.30 0.42 24.58 < 0.00000001
SessoM 73.58 10.94 6.72 < 0.00000001
I(Lunghezza^2) 0.02 0.00 6.66 < 0.00000001

L’Adjusted R-squared per questo modello è: 0.74

shapiro <- shapiro.test(residuals(mod7L))
bp <- bptest(mod7L)
dw <- dwtest(mod7L)

# Creazione di un data frame per la tabella
test_results <- data.frame(
  Test = c("Shapiro-Wilk Test", "Breusch-Pagan Test", "Durbin-Watson Test"),
  Statistic = round(c(shapiro$statistic, bp$statistic, dw$statistic), 3),
  `P-Value` = sapply(c(shapiro$p.value, bp$p.value, dw$p.value), format_p_value)
)

# Stampa della tabella
kable(test_results, format = "html", col.names = c("Test", "Statistic", "P-Value")) %>%
  kable_styling(full_width = FALSE, bootstrap_options = "striped")
Test Statistic P-Value
W Shapiro-Wilk Test 0.990 < 0.00000001
BP Breusch-Pagan Test 13.578 0.018527103
DW Durbin-Watson Test 1.952 0.11448895

Osserviamo che i parametri del nostro modello sono leggermente cambiati pur mantenendo lo stesso adjusted R-squared. L’altro cambiamento è legato al test di Breusch-Pagan che mostra un valore pari al 2% che seppur ancora ci porta ad un rifiuto dell’ipotesi di omoschedasticità, vuol dire che i nostri residui si discostano meno da questa ipotesi.

Visto l’impatto che ha avuto sui residui aver eliminato il datto 1551, per completezza, risvolgiamo una veloce rianalisi del modello, stavolta usando la funzione stepAIC(). Il modello risultante lo confronteremo con il mod7L e usaremo il test BIC per determinare se il nostro modello è ancora adeguato.

mod1B <- lm(Peso ~., data=dati_neonati)
suppress_output <- capture.output(
  stepwise.mod <- MASS::stepAIC(mod1B, direction = "both", k = log(n))
)  

print_summary(stepwise.mod)
Estimate Std.Error t_value P.value
(Intercept) -6683.83 133.16 -50.19 < 0.00000001
N.gravidanze 13.14 4.26 3.09 0.0020417163
Gestazione 29.63 3.74 7.93 < 0.00000001
Lunghezza 10.89 0.30 36.08 < 0.00000001
Cranio 9.92 0.42 23.46 < 0.00000001
SessoM 78.14 10.99 7.11 < 0.00000001

L’Adjusted R-squared per questo modello è: 0.74

kable(BIC(mod7L,stepwise.mod), format = "html") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "condensed"))
df BIC
mod7L 7 35046.88
stepwise.mod 7 35081.40

Il risultato del test BIC conferma che il modello mod7L è ancora il più adeguato.

3. Previsioni e Risultati

Usiamo il nostro modello per fare una previsione e vedere come si comporta:

nuovi_dati <- data.frame(
  Gestazione = 39,
  Sesso = 'F',  
  Lunghezza = mean(dati_neonati$Lunghezza),  
  Cranio = mean(dati_neonati$Cranio),       
  `I(Lunghezza^2)` = mean(dati_neonati$Lunghezza)^2 
)

previsione <- round(predict(mod7L, newdata = nuovi_dati),0)

cat("La previsione è:", previsione)

La previsione è: 3232

Con questo esempio abbiamo verificato che il nostro modello è in grado di fare una stima del peso dati i parametri richiesi.

4. Visualizzazioni

Mostreremo alcuni grafici che ci possono permettere di visualizzare la qualità del nostro grafico. Aggiungiamo al nostro dataset la colonna con le previsioni del nostro modello e costruiamo i grafici in funzione delle settimane di gestazione, della lunghezza e del diametro del cranio, oltre ad un grafico che fa il plot delle predizioni vs i dati dei pesi del dataset. Non mostriamo la dipendenza dal peso per non rendere i grafici poco chiari:

dati_neonati$predizione <- predict(mod7L, newdata = dati_neonati)

ggplot(dati_neonati, aes(x = Gestazione)) +
  geom_point(aes(y = Peso, color = "Dati reali"), alpha = 0.6, size = 3) +  # Dati reali
  geom_point(aes(y = predizione, color = "Predizioni"), alpha = 0.6, size = 3) +  # Predizioni
  labs(title = "Peso in funzione della Gestazione",
       x = "Gestazione (settimane)",
       y = "Peso (g)") +
  theme_minimal() +
  theme(legend.position = "right") +
  scale_color_manual(values = c("Dati reali" = "blue", "Predizioni" = "red"))  # Colori per dati e predizioni

ggplot(dati_neonati, aes(x = Lunghezza)) +
  geom_point(aes(y = Peso, color = "Dati reali"), alpha = 0.6, size = 3) +  # Dati reali
  geom_point(aes(y = predizione, color = "Predizioni"), alpha = 0.6, size = 3) +  # Predizioni
  labs(title = "Peso in funzione della Lunghezza",
       x = "Lunghezza (cm)",
       y = "Peso (g)") +
  theme_minimal() +
  theme(legend.position = "right") +
  scale_color_manual(values = c("Dati reali" = "blue", "Predizioni" = "red"))  # Colori per dati e predizioni

ggplot(dati_neonati, aes(x = Cranio)) +
  geom_point(aes(y = Peso, color = "Dati reali"), alpha = 0.6, size = 3) +  # Dati reali
  geom_point(aes(y = predizione, color = "Predizioni"), alpha = 0.6, size = 3) +  # Predizioni
  labs(title = "Peso in funzione del Cranio",
       x = "Cranio (cm)",
       y = "Peso (g)") +
  theme_minimal() +
  theme(legend.position = "right") +
  scale_color_manual(values = c("Dati reali" = "blue", "Predizioni" = "red"))  # Colori per dati e predizioni

ggplot(dati_neonati, aes(x = Peso, y = predizione)) +
  geom_point(color = "blue", alpha = 0.6, size = 3) +  # Punti dei dati
  geom_abline(intercept = 0, slope = 1, color = "red") +  # Linea ideale (peso osservato = peso predetto)
  labs(title = "Peso osservato vs Peso predetto",
       x = "Peso osservato (g)",
       y = "Peso predetto (g)") +
  theme_minimal()

Come possiamo vedere dai primi 3 grafici c’è un’ottima sovrapposizione tra i dati e le previsioni del nostro modello. Nell’ultimo grafico invece abbiamo creato uno scatterplot con le nostre previsioni in funzione dei valori della variabile Peso, in un modello ideale tutti i punti del grafico dovrebbero trovarsi sulla bisettrice del primo quadrante o comunque distribuirsi equamente attorno ad essa. Da questo grafico vediamo come il nostro modello nelle code tende a sottostimare i pesi dei neonati mentre nella zona centrale ha una buona risposta.

5. Conclusioni

Tramite l’analisi statistica del dataset sul peso dei neonati abbiamo ottenuto molteplci risultati:

  1. I parti cesarei sono equamente presenti nei tre ospedali, quindi probabilmente usano gli stessi protocolli per attivare questa procedura.
  2. Il nostro campione non è compatibile con i risultati trovati in letteratura per quanto riguarda il peso e la lunghezza di un neonato nato tra la 37-esima e la 41-esima settimana.
  3. Le misure antropometriche sono significativamente diverse tra i due sessi.
  4. La nostra analisi ci dice anche che in caso di madre tabagista questo non influenza il peso del neonato.

Considerazioni sul modello selezionato:

  1. Il modello presenta un adjusted R-squared pari a 0.74 quindi riesce a spiegare il 74% della varianza del dataset e un RMSE pari a circa l’8% del peso. Questi primi risultati sembrano soddisfacenti e sono corroborati anche dai grafici mostrati nell’ultima sezione in cui c’è una buona sovrapposizione tra i dati e le previsioni.
  2. Il modello mod7L non è perfetto infatti non soddisfa l’ipotesi di distribuzione normale e di omoschedasticità dei residui, con un p-value del 2%, (rispetta solo l’ipotesi d’indipendenza dei residui). Si presenta comunque come un buon compromesso tra un modello troppo semplice ed uno troppo complesso e quindi di difficile interpretazione.
  3. Il difetto di sottostimare il peso nelle code è un errore non grave per i bambini prematuri, in quanto porterebbe comunque ad attivare dei protocolli di cura per il neonato, è un problema più importante per le gravidanze extralunghe in cui le difficoltà possono essere legate alle eccessive dimensioni del nascituro.