#install.packages("tidyr")
library(kableExtra)
library(ggplot2)
library(psych)
## 
## Caricamento pacchetto: 'psych'
## I seguenti oggetti sono mascherati da 'package:ggplot2':
## 
##     %+%, alpha
library(DescTools)
## 
## Caricamento pacchetto: 'DescTools'
## I seguenti oggetti sono mascherati da 'package:psych':
## 
##     AUC, ICC, SD
library(dplyr)
## 
## Caricamento pacchetto: 'dplyr'
## Il seguente oggetto è mascherato da 'package:kableExtra':
## 
##     group_rows
## I seguenti oggetti sono mascherati da 'package:stats':
## 
##     filter, lag
## I seguenti oggetti sono mascherati da 'package:base':
## 
##     intersect, setdiff, setequal, union
library(paletteer)
library(tidyr)
#install.packages("gridExtra")
library(gridExtra)
## 
## Caricamento pacchetto: 'gridExtra'
## Il seguente oggetto è mascherato da 'package:dplyr':
## 
##     combine
library(grid)

Introduzione

L’obiettivo del progetto è quello di creare un modello statistico in grado di prevedere con precisione il peso dei neonati alla nascita, basandosi su variabili cliniche raccolte da tre ospedali. Si mira a migliorare la gestione delle gravidanze ad alto rischio, ottimizzare le risorse ospedaliere e garantire migliori risultati per la salute neonatale.

In questo contesto, la possibilità di prevedere il peso alla nascita dei neonati rappresenta un’opportunità fondamentale per migliorare la pianificazione clinica e ridurre i rischi associati a nascite problematiche, come parti prematuri o neonati con basso peso. Di seguito, i principali benefici che questo progetto porterà all’azienda e al settore sanitario:

  1. Miglioramento delle previsioni cliniche: il peso del neonato è un indicatore chiave della sua salute. Avere un modello predittivo accurato consente al personale medico di intervenire tempestivamente in caso di anomalie, riducendo le complicazioni perinatali come le difficoltà respiratorie o l’ipoglicemia.

  2. Ottimizzazione delle risorse ospedaliere: sapere in anticipo quali neonati potrebbero avere bisogno di cure intensive aiuta a organizzare le risorse umane e tecnologiche degli ospedali in modo efficiente. Questo si traduce in una riduzione dei costi operativi e una migliore pianificazione dell’utilizzo delle unità di terapia intensiva neonatale (TIN).

  3. Prevenzione e identificazione dei fattori di rischio: il modello potrà evidenziare i fattori che maggiormente influenzano negativamente il peso del neonato (come il fumo materno, gravidanze multiple o età avanzata della madre). Queste informazioni sono preziose per la prevenzione e la gestione personalizzata delle gravidanze, permettendo interventi proattivi in caso di rischio elevato.

  4. Valutazione delle pratiche ospedaliere: attraverso un’analisi comparativa tra i tre ospedali coinvolti, l’azienda potrà identificare eventuali differenze nei risultati clinici, come una maggiore incidenza di parti cesarei in una determinata struttura. Ciò consente di monitorare la qualità delle pratiche e armonizzare i protocolli tra i diversi centri ospedalieri, migliorando la coerenza delle cure.

  5. Supporto alla pianificazione strategica: l’analisi dei dati e le previsioni possono essere utilizzate per prendere decisioni informate non solo a livello clinico ma anche strategico. L’azienda potrà sfruttare queste informazioni per implementare nuove politiche di salute pubblica, garantendo un impatto positivo sui tassi di mortalità e morbilità neonatale.

1. Raccolta dei dati e struttura del dataset

Per costruire il modello predittivo, sono stati raccolti dati su 2500 neonati provenienti da tre ospedali, disponibili nel file neonati.csv.

dati<-read.csv("C:/Users/lara.nichetti/Desktop/Lara/neonati.csv")

Le variabili disponibili sono presentate nella seguente tabella:

variable_names<-c("Età della madre",
                  "Numero di gravidanze",
                  "Fumo materno", 
                  "Durata della gravidanza",
                  "Peso del neonato", 
                  "Lunghezza del cranio",
                  "Diametro del cranio", 
                  "Tipo di parto", 
                  "Ospedale di nascita",
                  "Sesso del neonato")

variables_descr<-c("Misura dell'età in anni", 
                   "Quante gravidanze ha avuto la madre", 
                   "Un indicatore binario (0=non fumatrice, 1=fumatrice)", 
                   "Numero di settimane di gestazione", 
                   "Peso alla nascita in grammi", 
                   "Lunghezza del neonato, misurabile anche durante la gravidanza tramite ecografie",
  
                   "Diametro craniale, misurabile anche durante la gravidanza tramite ecografie", 
                   "Naturale o cesareo", 
                   "Ospedale 1, 2 o 3", 
                   "Maschio (M) o femmina (F)")
variables<-data.frame("Variabili" = variable_names,
                         "Descrizione" = variables_descr)
variables %>%
  kbl() %>%
  kable_styling(bootstrap_options = c("striped", "bordered"),
                position = "center",
                full_width = FALSE)
Variabili Descrizione
Età della madre Misura dell’età in anni
Numero di gravidanze Quante gravidanze ha avuto la madre
Fumo materno Un indicatore binario (0=non fumatrice, 1=fumatrice)
Durata della gravidanza Numero di settimane di gestazione
Peso del neonato Peso alla nascita in grammi
Lunghezza del cranio Lunghezza del neonato, misurabile anche durante la gravidanza tramite ecografie
Diametro del cranio Diametro craniale, misurabile anche durante la gravidanza tramite ecografie
Tipo di parto Naturale o cesareo
Ospedale di nascita Ospedale 1, 2 o 3
Sesso del neonato Maschio (M) o femmina (F)

L’obiettivo principale è identificare quali di queste variabili sono più predittive del peso alla nascita, con un focus particolare sull’impatto del fumo materno e delle settimane di gestazione, che potrebbero indicare nascite premature.

2. Analisi e modellazione

2.1 Analisi Preliminare:

Durante la fase preliminare è importante esplorare le variabili tramite le statistiche descrittive in modo da comprenderne la distribuzione e identificare eventuali outlier o anomalie:

stat_descr<-describe(dati[,-c(3,8:10)])[,-c(1,6,7,13)]
stat_descr %>%
  kbl(caption="<span style='color:black'>Statistiche descrittive per variabili quantitative</span>") %>%
  kable_styling(bootstrap_options=c("striped", "bordered"),
                position="center")
Statistiche descrittive per variabili quantitative
n mean sd median min max range skew kurtosis
Anni.madre 2500 28.1640 5.273578 28 0 46 46 0.0427858 0.3777127
N.gravidanze 2500 0.9812 1.280587 1 0 12 12 2.5127457 10.9782171
Gestazione 2500 38.9804 1.868639 39 25 43 18 -2.0640742 8.2491448
Peso 2500 3284.0808 525.038744 3300 830 4930 4100 -0.6466427 2.0275074
Lunghezza 2500 494.6920 26.318644 500 310 565 255 -1.5137904 6.4795857
Cranio 2500 340.0292 16.425330 340 235 390 155 -0.7845817 2.9414502

Considerando queste statistiche descrittive possono essere fatte le seguenti osservazioni:

  • Anni.madre: l’età media del campione considerato è circa 28 anni, con deviazione standard di 5,3 anni; tuttavia, il valore minimo del range uguale a 0, quindi si tratta di un errore. La distribuzione è abbastanza simmetrica (asimmetria ≈ 0,04).
which(dati$Anni.madre == 0)
## [1] 1380
dati[which(dati$Anni.madre == 0),]
##      Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio
## 1380          0            0         0         39 3060       490    330
##      Tipo.parto Ospedale Sesso
## 1380        Nat     osp3     M

Cosiderando meglio il record errato, solo la variabile Anni Madre sembra essere errata, le altre sono realstiche, quindi non è il caso di eliminarla dal dataset.

  • N.gravidanze: in media le donne intervistate nel campione hanno avuto una gravidanza, in questo caso la varianza è contenuta (deviazione standard ≈ 1,3). Il range va da un minimo di 0 ad un massimo di 12 gravidanze, con asimmetria positiva (asimmetria ≈ 2,5), si può concludere che la maggior parte ha poche gravidanze, ma c’è una coda lunga di donne con molte. La curtosi di questa variabile è maggiore di tutte le altre variabili considerate, infatti è molto alta ( curtosi ≈ 11), quindi nella distibuzioni ci sono tante osservazioni con valori estremi.

  • Gestazione: la media della durata della gestioni è di 39 settimane, con una deviazione standard di circa 1,9 settimana, quinidi i valori sono molto concentrati. Il range è tra 25 e 43 settimane, questo indica la presenza di alcune gestazioni molto premature. L’asimmetria, in questo caso, è negativa (asimmetria ≈ -2,06); il risultato indica che vi sono più valori bassi rispetto alla media, cioè più prematuri rispetto a molto post-termine. La curtosi è elevata (curtosi ≈ 8,2), molti casi sono vicini alla media, ma sono presenti anche outlier importanti.

  • Peso: il peso medio alla nascita è di circa 3284 grammi, con una deviazione standard pari a 525 grammi. Il range si estende da 830 a 4930 grammi, includendo quindi sia casi di peso molto basso sia valori eccezionalmente alti. La distribuzione mostra una lieve asimmetria negativa (asimmetria ≈ -0,65), segnalando una maggiore concentrazione di casi verso i pesi più bassi.

  • Lunghezza del neonato: la lunghezza media è pari a circa 495 millimetri, con una deviazione standard di 26. Il range varia da 310 a 565 millimetri, includendo alcuni valori particolarmente bassi probabilmente riferibili a nascite premature. La distribuzione è fortemente asimmetrica a sinistra (asimmetria ≈ -1,5), con diversi casi al di sotto della media. La curtosi (≈ 6,5) indica la presenza di outlier e una distribuzione con code piuttosto pesanti.

  • Diametro del cranio: la circonferenza cranica presenta un valore medio di circa 340 millimetri, con una deviazione standard di 16. I valori variano da 235 a 390 millimetri. La distribuzione tende leggermente verso valori inferiori alla media (asimmetria ≈ -0,78), evidenziando la presenza di alcuni casi con cranio particolarmente piccolo. La curtosi è circa 2,9, molto vicina a quella di una distribuzione normale (3), indicando che la variabile ha una forma abbastanza regolare, senza evidenza di code particolarmente pronunciate.

Per valutare graficamente le variabili quantitative è possibile utilizzare i boxplots:

dati_long <- dati %>% 
  pivot_longer(
    cols = c(Anni.madre, N.gravidanze, Gestazione, Peso, Lunghezza, Cranio),
    names_to = "Variabile",
    values_to = "Valore"
  )

ggplot(dati_long, aes(x = 1, y = Valore, fill = Variabile)) +
  geom_boxplot(width = 0.5) +
  scale_fill_paletteer_d("nationalparkcolors::Arches") +
  labs(
    x = "",
    y = "Valori",
    title = "Distribuzione delle variabili quantitative"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    axis.text.y = element_text(size = 8),
    axis.title.y = element_text(size = 10),
    legend.text = element_text(size = 8),
    legend.title = element_text(size = 9),
    plot.title = element_text(hjust = 0.5) 
  ) +
  facet_wrap(~Variabile, scales = "free_y") +
  scale_x_continuous(expand = expansion(mult = c(0.5, 0.5)))

I boxplot confermano quanto emerso dalle statistiche descrittive. In particolare, è evidente il gran numero di outlier presenti in molte delle variabili. Inoltre, per le variabili Gestazione e N.gravidanze, i boxplot mettono in luce una forte asimmetria: la prima mostra una concentrazione di casi verso un numero basso di settimane di gestazione, mentre la seconda evidenzia alcune madri con un numero di gravidanze eccezionalmente elevato rispetto alla maggioranza del campione.

Al fine di analizzare la distribuzione delle variabili quantitative possiamo utilizzare delle tabelle di frequenza e dei barplots:

make_table_grob<-function(var, varname, labels=NULL, scale=0.6) {
  df<-data.frame(
    "Categoria"=names(table(var)),
    "Frequenza Assoluta"=as.vector(table(var)),
    "Frequenza Relativa"=round(as.vector(table(var)/length(var)), 3)
  )
  if (!is.null(labels)) {
    df$Categoria<-labels[as.character(df$Categoria)]
  }
    theme<-ttheme_default(
    core=list(fg_params=list(cex=scale)),
    colhead=list(fg_params=list(cex=scale))
  )
  tableGrob(df, rows=NULL, theme=theme)
}

make_barplot<-function(var, varname, colors=NULL, labels=NULL) {
  df<-data.frame(var=factor(var))
  plot<-ggplot(df, aes(x=var, fill=var)) +
    geom_bar() +
    theme_minimal() +
    labs(title=paste("Barplot:", varname), x=NULL, y="Conteggio") +
    theme(legend.title=element_blank())
  if (!is.null(colors)) {
    if (!is.null(labels)) {
      plot<-plot + scale_fill_manual(values=colors, labels=labels)
    } else {
      plot<-plot + scale_fill_manual(values=colors)
    }
  }
  plot
}

tab1<-make_table_grob(dati$Fumatrici, "Fumatrici", labels=c("0"="No", "1"="Sì"))
plt1<-make_barplot(dati$Fumatrici, "Fumatrici", 
                     colors=c("0"="antiquewhite3", "1"="azure3"), 
                     labels=c("0"="No", "1"="Sì"))

tab2<-make_table_grob(dati$Ospedale, "Ospedale", 
                        labels=c("osp1"="Ospedale 1", "osp2"="Ospedale 2", "osp3"="Ospedale 3"))
plt2<-make_barplot(dati$Ospedale, "Ospedale",
                     colors=c("osp1"="darksalmon", "osp2"="darkseagreen2", "osp3"="darkslategray2"),
                     labels=c("osp1"="Ospedale 1", "osp2"="Ospedale 2", "osp3"="Ospedale 3"))

tab3<-make_table_grob(dati$Sesso, "Sesso", labels=c("F"="Donna", "M"="Uomo"))
plt3<-make_barplot(dati$Sesso, "Sesso",
                     colors=c("F"="pink", "M"="lightblue"),
                     labels=c("F"="Donna", "M"="Uomo"))

grid.arrange(tab1, plt1, ncol=2, widths=c(0.5, 0.3))

  • Fumatrici: La tabella di frequenza e il barplot della variabile mostrano una forte asimmetria tra le osservazioni. Infatti, su un campione di 2500 intervistati, solo 104 donne sono fumatrici, corrispondenti al 4,2% del campione.
grid.arrange(tab2, plt2, ncol=2, widths=c(0.5, 0.3))

  • Ospedale: La distribuzione delle osservazioni è quasi uniforme tra i tre ospedali. L’ospedale che registra il numero maggiore di parti è il numero 2, con il 34% delle osservazioni; segue l’ospedale 3 con il 33,4%, mentre l’ospedale 1 registra il 32,6% delle osservazioni.
grid.arrange(tab3, plt3, ncol=2, widths=c(0.5, 0.3))

  • Sesso: Anche in questo caso, la distribuzione delle osservazioni è quasi uniforme tra i due sessi. Il sesso femminile ha solo 12 osservazioni in più rispetto al sesso maschile.

Inoltre, possiamo affermare che:

  • utilizzando un test chi-quadro per saggiare l’ipotesi di differenza tra il tipo di parto (Naturale o Cesareo) nei tre diversi ospedali:

    • H0: La proporzione di parti naturali e cesarei è uguale nei tre ospedali.

    • H1: Almeno un ospedale ha una proporzione di tipo di parto diversa.

tab <- table(dati$Ospedale, dati$Tipo.parto)
chisq.test(tab)
## 
##  Pearson's Chi-squared test
## 
## data:  tab
## X-squared = 1.0972, df = 2, p-value = 0.5778

Il p-value ottenuto è di 0.57, poiché il p-value è maggiore del livello di significatività standard (𝛼= 0.05), non ci sono evidenze sufficienti per rifiutare l’ipotesi nulla. Concludiamo quindi che la distribuzione dei tipi di parto non differisce significativamente tra i tre ospedali.

  • Utilizzando un test non parametrico con approccio frequentistico per determinare si può determinare il range in cui ricade la media della popolazione delle variabili Peso e Lunghezza con un intrevallo di confidenza del 95%. Risulta preferibile utilizzare un test non parametrico in quanto le due variabili considerati non rispettano una distribuzione normale:
shapiro.test(dati$Peso)
## 
##  Shapiro-Wilk normality test
## 
## data:  dati$Peso
## W = 0.97066, p-value < 2.2e-16
shapiro.test(dati$Lunghezza)
## 
##  Shapiro-Wilk normality test
## 
## data:  dati$Lunghezza
## W = 0.90941, p-value < 2.2e-16

Possiamo concludere che:

  • la media del peso della popolazione di un neonato è compreso tra 3263.089 e 3304.502 grammi, con una confidenza del 95%.

  • la media della lunghezza della popolazione di un neonato è compreso tra 493.6659 e 495.7260 centimetri, con una confidenza del 95%.

set.seed(123) 
B<-10000
n<-length(dati$Peso)

boot_means_peso<-numeric(B)
for(i in 1:B){
  sample_indices<-sample(1:n, size = n, replace = TRUE)
  boot_means_peso[i]<-mean(dati$Peso[sample_indices])
}

media_boot_peso<-mean(boot_means_peso)
IC_peso<-quantile(boot_means_peso, c(0.025, 0.975))
cat("IC Peso 95%:", IC_peso[1], "-", IC_peso[2], "\n\n")
## IC Peso 95%: 3263.089 - 3304.502
boot_means_lung<-numeric(B)
for(i in 1:B){
  sample_indices<-sample(1:n, size = n, replace = TRUE)
  boot_means_lung[i]<-mean(dati$Lunghezza[sample_indices])
}

media_boot_lung<-mean(boot_means_lung)
IC_lung<-quantile(boot_means_lung, c(0.025, 0.975))
cat("IC Lunghezza 95%:", IC_lung[1], "-", IC_lung[2], "\n\n")
## IC Lunghezza 95%: 493.6659 - 495.726
  • Utilizzando nuovamente è un test non parametrico: Wilcox test, è possibile saggiare l’ipotesi che le misure antropometriche sono significativamente diverse tra i due sessi.
shapiro.test(dati$Peso[dati$Sesso == "F"])
## 
##  Shapiro-Wilk normality test
## 
## data:  dati$Peso[dati$Sesso == "F"]
## W = 0.96285, p-value < 2.2e-16
shapiro.test(dati$Peso[dati$Sesso == "M"])
## 
##  Shapiro-Wilk normality test
## 
## data:  dati$Peso[dati$Sesso == "M"]
## W = 0.96647, p-value = 2.321e-16

  • il p-value ottenuto dal test applicato sulla variabile Peso è minore di 0.05, quindi non accettiamo H0 e concludiamo che le distribuzioni differiscono significativamente.

  • il p-value ottenuto dal test applicato sulla variabile Lunghezza è minore di 0.05, quindi non accettiamo H0 e concludiamo che le distribuzioni differiscono significativamente.

wilcox.test(Peso ~ Sesso, data = dati)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Peso by Sesso
## W = 538641, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(Lunghezza ~ Sesso, data = dati)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Lunghezza by Sesso
## W = 594455, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

2.2 Creazione del modello di Regressione:

In questa sezione, viene sviluppato un modello di regressione lineare multipla che includa tutte le variabili rilevanti. In questo modo, potremo quantificare l’impatto di ciascuna variabile indipendente sul peso del neonato ed eventuali interazioni.

In primo luogo, viene visualizzata la correlazione e gli scatter plot delle variabili quantitative:

dati$Fumatrici<-as.factor(dati$Fumatrici)
dati$Tipo.parto<-as.factor(dati$Tipo.parto)
dati$Ospedale<-as.factor(dati$Ospedale)
dati$Sesso<-as.factor(dati$Sesso)
attach(dati)
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(dati[,c(1,2,4:7)], upper.panel = panel.smooth, lower.panel = panel.cor)

Il peso è frotemente correlato con le settimane di gestazione, con una relazione direttamente proporzionale, ma che sembra non essere lineare, poiché all’aumentare delle settimane di gestazione il peso aumenta meno che proporzionalmente: la curva verso la fine sembra appiattirsi.

Inoltre, il peso è correlato positivamente con la lunghezza del neonato e il diametro del cranio, in questo caso la relazione sembra abbastanza lineare.

Date le prime evidenze, viene creato un modello di regressione lineare multipla nel quale vengono include tutte le variabili e le interazioni possibili:

model_1<-lm(Peso ~ .^2, data = dati)
summary(model_1)
## 
## Call:
## lm(formula = Peso ~ .^2, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1051.96  -181.59    -9.19   160.73  2138.57 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 3.216e+03  1.482e+03   2.169 0.030146 *  
## Anni.madre                 -6.717e+01  2.703e+01  -2.485 0.013025 *  
## N.gravidanze               -1.159e+02  9.906e+01  -1.170 0.242311    
## Fumatrici1                 -4.059e+02  9.861e+02  -0.412 0.680674    
## Gestazione                 -2.673e+02  6.826e+01  -3.917 9.22e-05 ***
## Lunghezza                   1.847e+01  5.531e+00   3.339 0.000854 ***
## Cranio                     -1.961e+01  8.456e+00  -2.319 0.020469 *  
## Tipo.partoNat              -5.431e+02  3.381e+02  -1.606 0.108378    
## Ospedaleosp2               -2.652e+02  3.503e+02  -0.757 0.449141    
## Ospedaleosp3               -2.109e+02  3.428e+02  -0.615 0.538572    
## SessoM                      4.629e+01  2.932e+02   0.158 0.874557    
## Anni.madre:N.gravidanze    -1.780e-01  8.584e-01  -0.207 0.835790    
## Anni.madre:Fumatrici1       3.385e+00  5.840e+00   0.580 0.562219    
## Anni.madre:Gestazione       2.772e+00  7.510e-01   3.691 0.000228 ***
## Anni.madre:Lunghezza       -2.280e-01  6.083e-02  -3.748 0.000182 ***
## Anni.madre:Cranio           2.150e-01  8.573e-02   2.508 0.012223 *  
## Anni.madre:Tipo.partoNat    1.076e+00  2.520e+00   0.427 0.669505    
## Anni.madre:Ospedaleosp2    -1.053e+00  2.754e+00  -0.383 0.702076    
## Anni.madre:Ospedaleosp3     1.153e+00  2.828e+00   0.407 0.683689    
## Anni.madre:SessoM          -1.709e+00  2.299e+00  -0.743 0.457453    
## N.gravidanze:Fumatrici1    -1.084e+01  2.424e+01  -0.447 0.654718    
## N.gravidanze:Gestazione    -5.832e+00  2.989e+00  -1.951 0.051149 .  
## N.gravidanze:Lunghezza      4.614e-01  2.485e-01   1.857 0.063496 .  
## N.gravidanze:Cranio         3.571e-01  3.827e-01   0.933 0.350775    
## N.gravidanze:Tipo.partoNat  1.335e+01  1.057e+01   1.263 0.206817    
## N.gravidanze:Ospedaleosp2  -2.937e+00  1.138e+01  -0.258 0.796384    
## N.gravidanze:Ospedaleosp3   3.065e+00  1.175e+01   0.261 0.794192    
## N.gravidanze:SessoM         7.673e+00  9.875e+00   0.777 0.437234    
## Fumatrici1:Gestazione      -4.557e+01  2.294e+01  -1.987 0.047072 *  
## Fumatrici1:Lunghezza        5.977e+00  1.752e+00   3.412 0.000656 ***
## Fumatrici1:Cranio          -2.126e+00  2.413e+00  -0.881 0.378365    
## Fumatrici1:Tipo.partoNat   -7.760e+01  6.454e+01  -1.202 0.229337    
## Fumatrici1:Ospedaleosp2    -1.229e+01  6.657e+01  -0.185 0.853485    
## Fumatrici1:Ospedaleosp3    -1.342e+02  6.874e+01  -1.953 0.050985 .  
## Fumatrici1:SessoM          -7.085e+01  6.167e+01  -1.149 0.250719    
## Gestazione:Lunghezza       -9.161e-03  1.158e-01  -0.079 0.936940    
## Gestazione:Cranio           6.921e-01  2.185e-01   3.167 0.001557 ** 
## Gestazione:Tipo.partoNat    2.063e+01  8.662e+00   2.382 0.017311 *  
## Gestazione:Ospedaleosp2    -1.653e+01  9.646e+00  -1.714 0.086740 .  
## Gestazione:Ospedaleosp3    -6.587e-01  9.329e+00  -0.071 0.943711    
## Gestazione:SessoM          -3.818e+00  7.799e+00  -0.490 0.624515    
## Lunghezza:Cranio           -6.489e-03  1.445e-02  -0.449 0.653413    
## Lunghezza:Tipo.partoNat    -7.103e-01  6.808e-01  -1.043 0.296884    
## Lunghezza:Ospedaleosp2      1.388e+00  7.419e-01   1.870 0.061579 .  
## Lunghezza:Ospedaleosp3      2.574e-01  7.392e-01   0.348 0.727671    
## Lunghezza:SessoM            9.997e-01  6.126e-01   1.632 0.102853    
## Cranio:Tipo.partoNat        2.623e-01  9.470e-01   0.277 0.781804    
## Cranio:Ospedaleosp2         7.318e-01  1.044e+00   0.701 0.483399    
## Cranio:Ospedaleosp3         3.534e-01  1.046e+00   0.338 0.735508    
## Cranio:SessoM              -7.910e-01  8.782e-01  -0.901 0.367860    
## Tipo.partoNat:Ospedaleosp2  4.481e+00  2.929e+01   0.153 0.878439    
## Tipo.partoNat:Ospedaleosp3 -2.801e+01  2.976e+01  -0.941 0.346734    
## Tipo.partoNat:SessoM       -9.861e+00  2.464e+01  -0.400 0.689057    
## Ospedaleosp2:SessoM        -6.733e+00  2.738e+01  -0.246 0.805770    
## Ospedaleosp3:SessoM         8.361e+00  2.740e+01   0.305 0.760249    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 270.4 on 2445 degrees of freedom
## Multiple R-squared:  0.7405, Adjusted R-squared:  0.7348 
## F-statistic: 129.2 on 54 and 2445 DF,  p-value: < 2.2e-16

Il modello appena creato contiene un numero di regressori decisamente alto ed è di difficile interpretazione. Per andare a ridurre la complessità del modello, vengono manetenute solo le interazioni con una significatività maggiore del 95%.

model_2<-lm(Peso ~ .+Anni.madre:Gestazione+Anni.madre:Lunghezza+Anni.madre:Cranio+
              Fumatrici:Gestazione+Fumatrici:Lunghezza+Gestazione:Tipo.parto+
              Gestazione:Cranio, data = dati)
summary(model_2)
## 
## Call:
## lm(formula = Peso ~ . + Anni.madre:Gestazione + Anni.madre:Lunghezza + 
##     Anni.madre:Cranio + Fumatrici:Gestazione + Fumatrici:Lunghezza + 
##     Gestazione:Tipo.parto + Gestazione:Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1030.11  -181.09   -10.04   163.94  2351.12 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              2758.79269 1416.45662   1.948 0.051567 .  
## Anni.madre                -71.76158   24.54260  -2.924 0.003487 ** 
## N.gravidanze               12.59382    4.60336   2.736 0.006267 ** 
## Fumatrici1                -46.46733  804.57412  -0.058 0.953949    
## Gestazione               -231.66444   35.72398  -6.485 1.07e-10 ***
## Lunghezza                  16.46893    1.61047  10.226  < 2e-16 ***
## Cranio                    -18.66654    4.44916  -4.196 2.82e-05 ***
## Tipo.partoNat            -480.19641  262.60570  -1.829 0.067582 .  
## Ospedaleosp2               -7.14898   13.26979  -0.539 0.590114    
## Ospedaleosp3               27.98510   13.32863   2.100 0.035863 *  
## SessoM                     72.42418   11.08998   6.531 7.92e-11 ***
## Anni.madre:Gestazione       2.44574    0.67934   3.600 0.000324 ***
## Anni.madre:Lunghezza       -0.21560    0.05534  -3.896 0.000101 ***
## Anni.madre:Cranio           0.24755    0.07651   3.235 0.001231 ** 
## Fumatrici1:Gestazione     -63.22470   21.86548  -2.892 0.003867 ** 
## Fumatrici1:Lunghezza        5.07210    1.45361   3.489 0.000493 ***
## Gestazione:Tipo.partoNat   13.04078    6.72458   1.939 0.052581 .  
## Gestazione:Cranio           0.57672    0.09136   6.312 3.24e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 270.2 on 2482 degrees of freedom
## Multiple R-squared:  0.737,  Adjusted R-squared:  0.7352 
## F-statistic: 409.1 on 17 and 2482 DF,  p-value: < 2.2e-16

Il Modello 2 contiene un numero inferiore di variabili e l’\(R^2_{adj}\) è maggiore rispetto a quello del modello precedente. Tuttavia, possiamo escludere ancora altri regressore.

model_3 <- lm(Peso ~ Anni.madre+N.gravidanze+Fumatrici+
                Gestazione+Lunghezza+Cranio+
                Anni.madre:Gestazione + Anni.madre:Lunghezza+
                Anni.madre:Cranio + Fumatrici:Gestazione + Fumatrici:Lunghezza+
                Gestazione:Cranio, data = dati)
summary(model_3)
## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Anni.madre:Gestazione + Anni.madre:Lunghezza + 
##     Anni.madre:Cranio + Fumatrici:Gestazione + Fumatrici:Lunghezza + 
##     Gestazione:Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1029.36  -176.46   -13.64   168.51  2372.33 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           2865.85018 1398.13640   2.050 0.040492 *  
## Anni.madre             -71.06550   24.76198  -2.870 0.004140 ** 
## N.gravidanze            13.40988    4.64864   2.885 0.003952 ** 
## Fumatrici1            -226.44268  810.46697  -0.279 0.779963    
## Gestazione            -235.06782   35.23841  -6.671 3.12e-11 ***
## Lunghezza               16.56293    1.62404  10.199  < 2e-16 ***
## Cranio                 -20.49084    4.47380  -4.580 4.87e-06 ***
## Anni.madre:Gestazione    2.36684    0.68586   3.451 0.000568 ***
## Anni.madre:Lunghezza    -0.21249    0.05580  -3.808 0.000143 ***
## Anni.madre:Cranio        0.25073    0.07718   3.249 0.001174 ** 
## Fumatrici1:Gestazione  -59.32732   22.08621  -2.686 0.007276 ** 
## Fumatrici1:Lunghezza     5.13407    1.46703   3.500 0.000474 ***
## Gestazione:Cranio        0.62604    0.09185   6.816 1.17e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.2 on 2487 degrees of freedom
## Multiple R-squared:  0.7306, Adjusted R-squared:  0.7293 
## F-statistic:   562 on 12 and 2487 DF,  p-value: < 2.2e-16

A fronte di una riduzione del numero di regressori, l’\(R^2_{adj}\) è peggiorato di \(0,0059\), un risultato assolutamente accettabile a fronte di una maggiore semplicità del modello.

Successivamente, si considera il modello senza interazioni e con tutti i regressori:

model_4 <- lm(Peso ~ . , data=dati)
summary(model_4)
## 
## Call:
## lm(formula = Peso ~ ., data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1124.40  -181.66   -14.42   160.91  2611.89 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6738.4762   141.3087 -47.686  < 2e-16 ***
## Anni.madre        0.8921     1.1323   0.788   0.4308    
## N.gravidanze     11.2665     4.6608   2.417   0.0157 *  
## Fumatrici1      -30.1631    27.5386  -1.095   0.2735    
## Gestazione       32.5696     3.8187   8.529  < 2e-16 ***
## Lunghezza        10.2945     0.3007  34.236  < 2e-16 ***
## Cranio           10.4707     0.4260  24.578  < 2e-16 ***
## Tipo.partoNat    29.5254    12.0844   2.443   0.0146 *  
## Ospedaleosp2    -11.2095    13.4379  -0.834   0.4043    
## Ospedaleosp3     28.0958    13.4957   2.082   0.0375 *  
## SessoM           77.5409    11.1776   6.937 5.08e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2489 degrees of freedom
## Multiple R-squared:  0.7289, Adjusted R-squared:  0.7278 
## F-statistic: 669.2 on 10 and 2489 DF,  p-value: < 2.2e-16

Riducendo nuovamente il numero delle variabili, l’\(R^2_{adj}\) peggira leggermente, ma in modo accettabile rispetto alla semplificazione del modello.

Il modello viene ulteriormente semplificato andando a rimuovere i regressori non significativi:

model_5<-lm(Peso~ N.gravidanze+Gestazione+Lunghezza+Cranio+Tipo.parto+Ospedale+Sesso)
summary(model_5)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Tipo.parto + Ospedale + Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1113.18  -181.16   -16.58   161.01  2620.19 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -6707.4293   135.9438 -49.340  < 2e-16 ***
## N.gravidanze     12.3619     4.3325   2.853  0.00436 ** 
## Gestazione       31.9909     3.7896   8.442  < 2e-16 ***
## Lunghezza        10.3086     0.3004  34.316  < 2e-16 ***
## Cranio           10.4922     0.4254  24.661  < 2e-16 ***
## Tipo.partoNat    29.2803    12.0817   2.424  0.01544 *  
## Ospedaleosp2    -11.0227    13.4363  -0.820  0.41209    
## Ospedaleosp3     28.6408    13.4886   2.123  0.03382 *  
## SessoM           77.4412    11.1756   6.930 5.36e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.9 on 2491 degrees of freedom
## Multiple R-squared:  0.7287, Adjusted R-squared:  0.7278 
## F-statistic: 836.3 on 8 and 2491 DF,  p-value: < 2.2e-16

Riducendo il numero di variabili, nel Modello 5 l’l’\(R^2_{adj}\) rimase invariato.

Per scrupolo, vengono tolte anche le variabili con un livello di significatività inferiore al 99%.

model_6<-lm(Peso~ N.gravidanze+Gestazione+Lunghezza+Cranio+Sesso)
summary(model_6)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.44  -180.81   -15.58   163.64  2639.72 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.1445   135.7229 -49.226  < 2e-16 ***
## N.gravidanze    12.4750     4.3396   2.875  0.00408 ** 
## Gestazione      32.3321     3.7980   8.513  < 2e-16 ***
## Lunghezza       10.2486     0.3006  34.090  < 2e-16 ***
## Cranio          10.5402     0.4262  24.728  < 2e-16 ***
## SessoM          77.9927    11.2021   6.962 4.26e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2494 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1328 on 5 and 2494 DF,  p-value: < 2.2e-16

Il Modello 6 sembra offrire il miglior equilibrio tra capacità esplicativa e semplicità.

2.3 Selezione del Modello Ottimale:

Attraverso tecniche di selezione del modello, come la minimizzazione del criterio di informazione di Akaike (AIC) o di Bayes (BIC), si vuole confermare o smentire il risulato ottenuto nella sezione precedente al fine di ottenere il modello più parsimonioso, ma che abbia una buona capacità eplicativa dei dati. Verranno considerati anche modelli con interazioni tra le variabili e possibili effetti non lineari.

anova(model_1,model_2,model_3,model_4,model_5,model_6)
## Analysis of Variance Table
## 
## Model 1: Peso ~ (Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso)^2
## Model 2: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso + Anni.madre:Gestazione + 
##     Anni.madre:Lunghezza + Anni.madre:Cranio + Fumatrici:Gestazione + 
##     Fumatrici:Lunghezza + Gestazione:Tipo.parto + Gestazione:Cranio
## Model 3: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Anni.madre:Gestazione + Anni.madre:Lunghezza + Anni.madre:Cranio + 
##     Fumatrici:Gestazione + Fumatrici:Lunghezza + Gestazione:Cranio
## Model 4: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
## Model 5: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
##     Ospedale + Sesso
## Model 6: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##   Res.Df       RSS  Df Sum of Sq       F    Pr(>F)    
## 1   2445 178778861                                    
## 2   2482 181188760 -37  -2409899  0.8908 0.6584111    
## 3   2487 185591325  -5  -4402565 12.0420 1.529e-11 ***
## 4   2489 186762521  -2  -1171196  8.0087 0.0003414 ***
## 5   2491 186899996  -2   -137475  0.9401 0.3907445    
## 6   2494 188065546  -3  -1165550  5.3134 0.0011935 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Considerando il test ANOVA, i modelli 3, 4 e 6 risultano significativi almeno al livello del 95%. Ciò indica che, rispetto al modello di riferimento (Modello 1), l’esclusione delle variabili eliminate non compromette la capacità esplicativa, che rimane statisticamente rilevante. Poiché si tratta di modelli più semplici, essi rappresentano un buon compromesso tra capacità esplicativa e parsimonia.

Inoltre, cosiderando la minimizzazione del criterio di informazione di Bayes (BIC), si conferma che il Modello 6 è il migliore tra quelli considerati.

Infine, viene utilizzata la procedura stepwise per indagare se vi sono modelli migliori che non sono stati considerati.

## 
## Caricamento pacchetto: 'MASS'
## Il seguente oggetto è mascherato da 'package:dplyr':
## 
##     select
## 
## Call:
## lm(formula = Peso ~ Anni.madre + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso + Anni.madre:Gestazione + Anni.madre:Lunghezza + 
##     Anni.madre:Cranio + Fumatrici:Gestazione + Fumatrici:Lunghezza + 
##     Gestazione:Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1063.01  -181.69   -11.02   167.09  2358.05 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           2329.37310 1390.72886   1.675 0.094074 .  
## Anni.madre             -70.83130   24.59144  -2.880 0.004007 ** 
## Fumatrici1               7.94736  805.75589   0.010 0.992131    
## Gestazione            -218.06096   35.06761  -6.218 5.88e-10 ***
## Lunghezza               16.21047    1.61215  10.055  < 2e-16 ***
## Cranio                 -18.38141    4.45292  -4.128 3.78e-05 ***
## SessoM                  72.94950   11.12286   6.559 6.59e-11 ***
## Anni.madre:Gestazione    2.35358    0.68091   3.457 0.000556 ***
## Anni.madre:Lunghezza    -0.20897    0.05539  -3.773 0.000165 ***
## Anni.madre:Cranio        0.24946    0.07664   3.255 0.001150 ** 
## Fumatrici1:Gestazione  -61.25661   21.93338  -2.793 0.005265 ** 
## Fumatrici1:Lunghezza     4.81361    1.45687   3.304 0.000966 ***
## Gestazione:Cranio        0.57167    0.09156   6.244 5.00e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 271.3 on 2487 degrees of freedom
## Multiple R-squared:  0.7343, Adjusted R-squared:  0.733 
## F-statistic: 572.7 on 12 and 2487 DF,  p-value: < 2.2e-16

Questa procedura indica un modello simile al Modello 3, l’unica differenza è che in quest’ultimo la variabile N.gravidanze rientrava trai regressori, mentre nel modello ottenuto tramite stepwise questa non è inclusa.

Utilizziamo nuovamente l’ANOVA e il BIC per paragonare questo nuovo modello e il Modello 6:

anova(model_6,stepwise.mod)
## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ Anni.madre + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Sesso + Anni.madre:Gestazione + Anni.madre:Lunghezza + Anni.madre:Cranio + 
##     Fumatrici:Gestazione + Fumatrici:Lunghezza + Gestazione:Cranio
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1   2494 188065546                                  
## 2   2487 183046410  7   5019136 9.7419 5.078e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BIC(model_6,stepwise.mod)
##              df      BIC
## model_6       7 35220.10
## stepwise.mod 14 35207.24

Entrambi gli indicatori mostrano che è preferibile optare per il modello ottenuto tramite la procedura di stepwise. Tuttavia, data la semplicità del Modello 6 e la sua capacità di spiegare la variabilità della variabile risposta, si conclude che questo è il modello finale:

\[ \begin{aligned} Peso = \ & -6681.1445 + 12.4750 \, \cdot \, N.gravidanze + 32.3321 \, \cdot \, Gestazione \\ &+ 10.2486 \, \cdot\, Lunghezza + 10.5402 \, \cdot\, Cranio + 77.9927 \, \cdot\, Sesso \end{aligned} \]

2.4 Analisi della Qualità del Modello

Il modello ottenuto viene valutato qualitativamente in questa sezione, vengono considerati l’\(R^2\) e il RMSE, oltre alle analisi sui residui del modello.

L’\(R^2\) corretto misura quanta parte della variabilità della variabile risposta (Peso) è spiegata dal modello, tenendo conto anche del numero di predittori.Nel modello considerato circa il 72–73% della variabilità del peso alla nascita viene spiegata dai regressori inclusi.

Il RMSE è molto vicino a 275 g, questo significa che, in media, le previsioni del modello sul peso dei neonati si discostano di circa 275 grammi dal valore osservato.

Considerando i residue utilizziamo i seguneti test:

library(car)
## Caricamento del pacchetto richiesto: carData
## 
## Caricamento pacchetto: 'car'
## Il seguente oggetto è mascherato da 'package:dplyr':
## 
##     recode
## Il seguente oggetto è mascherato da 'package:DescTools':
## 
##     Recode
## Il seguente oggetto è mascherato da 'package:psych':
## 
##     logit
par(mfrow=c(2,2))
plot(model_6)

Il primo grafico Residuals vs Fitted serve a controllare la linearità e l’omoschedasticità, i residui del modello sembrano essere distribuiti in modo abbastanza casuale intorno allo zero, ma si nota una maggiore dispersione per valori fitted più alti.

Il Q-Q plot dei residui serve a valutare la normalità dei residui, nel modello la maggior parte dei punti segue la linea di riferimento, ma gli estremi deviano dalla linea, indicando che ci sono outlier o residui pesanti alle code.

Il plot Scale-Location serve a verificare l’omoschedasticità in modo più chiaro: la dispersione dei residui sembra aumentare leggermente con i fitted values più grandi, confermando quanto visto nel primo grafico.

Infine, il plot Residuals vs Leverage serve a identificare osservazioni influenti, nel modello ci sono alcuni punti chiaramente influenti. La Distanza di Cook conferma questi punti come influenti, tuttavia, nessun punto sembra superare chiaramente la soglia di 1, ma alcuni si avvicinano. Questo segnala che potrebbero avere un impatto non trascurabile sul modello.

par(mfrow=c(2,2))
#leverage
lev<-hatvalues(model_6)
plot(lev)
p=sum(lev)
soglia=2*p/n
abline(h=soglia,col=2)

#outliers
plot(rstudent(model_6))
abline(h=c(-2,2),col=2) # come fare un test t
outlierTest(model_6) # viene applicata la correzione di Bonferroni
##       rstudent unadjusted p-value Bonferroni p
## 1551 10.051908         2.4906e-23   6.2265e-20
## 155   5.027798         5.3138e-07   1.3285e-03
## 1306  4.827238         1.4681e-06   3.6702e-03
# distanza di cook
cook<-cooks.distance(model_6)
plot(cook)
max(cook)
## [1] 0.8300965
library(lmtest)
## Caricamento del pacchetto richiesto: zoo
## 
## Caricamento pacchetto: 'zoo'
## I seguenti oggetti sono mascherati da 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(model_6)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_6
## BP = 90.253, df = 5, p-value < 2.2e-16
dwtest(model_6)
## 
##  Durbin-Watson test
## 
## data:  model_6
## DW = 1.9535, p-value = 0.1224
## alternative hypothesis: true autocorrelation is greater than 0
shapiro.test(residuals(model_6))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_6)
## W = 0.97408, p-value < 2.2e-16
plot(density(residuals(model_6)))

Il grafico in alto a sinitra mostra i leverage dei singoli punti, la maggior parte dei punti ha leverage molto basso, vicini a 0. Ci sono pochi punti sopra la soglia, che potrebbero avere un’influenza più alta sul modello. Confermando sostanzialmente quanto concluso attravesto i grafici precedenti.

Il plot in alto a destra mostra i residui corretti per varianza non costante, la linea rossa a ±2 indica una soglia per considerare un punto un outlier. La maggior parte dei residui è compresa tra -2 e 2. Con l’outlierTest(stepwise.mod) andiamo a identificare i punti che sono significativamente outliers, applicando la correzione di Bonferroni; questi sono i punti 1551, 155 e 1306.

Il plot in basso a sinistra mostra la distanza di Cook, che misura l’influenza complessiva di ciascun punto sul modello. La maggior parte dei punti ha valori di Cook vicini a 0, quindi questi hanno una bassa influenza. Solo un punto ha valore più elevato, ma non estremi.

Infine, la distribuzione dei residui è visibile nell’ultimo plot, in basso a destra. Questo mostra la densità dei residui. La distribuzione sembra leggermente asimmetrica, con un picco centrato vicino a 0. Il test di Shapiro-Wilk verifica la normalità dei residui, poiché p-value < 0.05 allora i residui non sono normalmente distribuiti.

Il test di Breusch-Pagan verifica se la varianza dei residui è costante, il p-value < 0.05, quindi è presente l’eteroschedasticità.

Il test di Durbin-Watson verifica l’autocorrelazione dei residui, il p-value > 0.05, quindi l’autocorrelazione è pari a zero.

Per concludere, anche se vi sono evidenze dell’etereschedasticità e dell’asimmetria dei residui del modello, queste non sono particolarmente elvate. Anche gli outliers e i punti di leverage sono presenti, ma nessun punto sembra superare chiaramente la soglia di Cook. Quindi il modello di regressione lineare multipla può essere utilizzato come punto di partenza per prevedere il peso dei neonati. Tuttavia, questo potrebbe essere migliorato attraverso un GLM.

3. Previsioni e Risultati

I coefficienti del modello finale possono essere interpretati come segue:

Considerando le conclusioni tratte fino, è possibile stimare il peso di una neonata considerando una madre alla terza gravidanza che partorirà alla 39esima settimana, dove la lunghezza e diametro del hanno le misure medie del campione. In questo caso il peso della neonata sarà di 3271.09 grammi.

newdata <- data.frame(
  N.gravidanze = 3,
  Gestazione = 39,
  Lunghezza = mean(dati$Lunghezza),
  Cranio = mean(dati$Cranio),
  Sesso = factor("F", levels = levels(dati$Sesso))
)
predict(model_6, newdata = newdata)
##       1 
## 3271.09

4. Visualizzazioni

Infine, vengono utilizzati grafici e rappresentazioni visive per comunicare i risultati del modello e mostrare le relazioni più significative tra le variabili. Ad esempio, potremmo visualizzare l’impatto del numero di settimane di gestazione, della lunghezza del neonato e diametro del cranio.

Innanzitutto visualizzaimo un grafico in cui abbiamo sull’asse delle ascisse il peso predetto dal modello, mentre sull’asse delle ordinate il peso reale registrato nel dataset. Inoltre viene mostrata sia l’intercetta sia la retta di regrassione del modello. Da notare che nel grafico viene identificata, ancora una volta, una leggera eteroschedasticità dei residue, poiché all’aumentare del peso aumenta anche la dispersione dei residui intorno all’intercetta.

peso_predict <- predict(model_6, dati)
grafico_df <- data.frame(Peso_reale = dati$Peso, Peso_predetto = peso_predict)
ggplot(grafico_df, aes(x = Peso_predetto, y = Peso_reale)) +
  geom_point(alpha = 0.6) +  
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed", linewidth=1.5) +
  geom_smooth(method = "lm", se = TRUE, color = "blue", fill = "lightblue", linewidth=0.8) +
  annotate("text", x = 100, 
           y = -120 + 1, 
           label = "Linea y = x", color = "red", hjust = 0) +
  annotate("text", x = 500, 
           y = 250, 
           label = "Retta regressione", color = "blue", hjust = 0) +
  expand_limits(x = grafico_df$Peso_reale, y = grafico_df$Peso_reale) +
  labs(title = "Confronto tra peso predetto e peso reale",
       x = "Peso predetto",
       y = "Peso reale") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Inoltre, nei sueguenti grafici sono rappresentate le relazioni tra la variabile dipendente e i regressori considerati nel modello:

ggplot(dati, aes(x = N.gravidanze, y = Peso)) +
  geom_point(alpha = 0.6) +  
  geom_smooth(method = "lm", se = TRUE, color = "blue", fill = "lightblue", linewidth = 0.8) +
  annotate("text", x = max(dati$N.gravidanze) * 0.7, 
           y = max(dati$Peso) * 0.9, 
           label = "Retta regressione", color = "blue", hjust = 0) +
  labs(title = "Peso in funzione del numero di gravidanze",
       x = "Numero di gravidanze",
       y = "Peso") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Il primo grafico rappresenta il peso del neonato rispetto al numero di gravidanze (variabile discreta). Da questo si evinca che c’è una lieve relazione lineare positiva tra le due variabili: al crescere del numero di gravidante, il peso del neonato aumenta leggermete. In ogni caso, dal grafico si intuisce anche che vi sia una maggioranza di intervistate con un numero di gravidanze nullo o basso, infatti, la maggior parte delle osservazioni è concentrato a sinistra.

ggplot(dati, aes(x = Gestazione, y = Peso)) +
  geom_point(alpha = 0.6) +  
  geom_smooth(method = "lm", se = TRUE, color = "blue", fill = "lightblue", linewidth = 0.8) +
  annotate("text", x = max(dati$N.gravidanze) * 0.7, 
           y = max(dati$Peso) * 0.9, 
           label = "Retta regressione", color = "blue", hjust = 0) +
  labs(title = "Peso in funzione delle settimane di gestazione",
       x = "Settimane di gestazione",
       y = "Peso") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

La relazione tra il peso del neonato e le settimane di gestazione (variabile discreta), rappresentata nel grafico, è posiriva e particolarmente evidente: all’aumentare delle settimane di gestazione aumenta il peso del neonato.

ggplot(dati, aes(x = Lunghezza, y = Peso)) +
  geom_point(alpha = 0.6) +  
  geom_smooth(method = "lm", se = TRUE, color = "blue", fill = "lightblue", linewidth = 0.8) +
  annotate("text", x = max(dati$N.gravidanze) * 0.7, 
           y = max(dati$Peso) * 0.9, 
           label = "Retta regressione", color = "blue", hjust = 0) +
  labs(title = "Peso in funzione della lunghezza del neonato",
       x = "Lunghezza (mm) del neonato",
       y = "Peso") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

ggplot(dati, aes(x = Cranio, y = Peso)) +
  geom_point(alpha = 0.6) +  
  geom_smooth(method = "lm", se = TRUE, color = "blue", fill = "lightblue", linewidth = 0.8) +
  annotate("text", x = max(dati$N.gravidanze) * 0.7, 
           y = max(dati$Peso) * 0.9, 
           label = "Retta regressione", color = "blue", hjust = 0) +
  labs(title = "Peso in funzione del diametro del cranio del neonato",
       x = "Diametro (mm) del cranio del neonato",
       y = "Peso") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Anche in questo caso, la relazioni rappresentate nei grafici sono positive ed elevate: all’aumentare della lunghezza del neonato o del suo diametro del cranio, il peso aumenta.

ggplot(dati, aes(x = Sesso, y = Peso, fill = Sesso)) +
  geom_violin(alpha = 0.5, trim = FALSE) +
  geom_boxplot(width = 0.1, outlier.shape = NA, alpha = 0.7) +
  labs(
    title = "Distribuzione del peso per sesso del neonato",
    x = "Sesso",
    y = "Peso (gr)"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

Infine, considerando l’ultima variabile del modello è posiibile confrontgare i pesi dei neonati maschi e femmine, dal grafico a violino emerge che entrambi i gruppi presentano una distribuzione molto simile. Tuttavia, i maschi mostrino in media valori leggermente più elevati, in linea con il modello di regressione multipla.

Conclusioni

Il progetto di previsione del peso neonatale è un’iniziativa fondamentale per Neonatal Health Solutions. Attraverso l’uso di dati clinici dettagliati e strumenti di analisi statistica, possiamo contribuire a migliorare la qualità della cura prenatale, ridurre i rischi per i neonati e ottimizzare l’efficienza delle risorse ospedaliere. Questo progetto rappresenta un punto di svolta per l’azienda, consentendo non solo un miglioramento della pratica clinica ma anche l’implementazione di politiche sanitarie più informate e proattive.