1. Introduzione

Questo progetto mira a sviluppare un modello statistico per prevedere il peso dei neonati alla nascita utilizzando variabili cliniche e demografiche raccolte in tre diversi ospedali. L’obiettivo è supportare decisioni cliniche e ottimizzare la gestione delle gravidanze a rischio.

Caricamento del dataset

dati <- read.csv("neonati.csv", sep=",", stringsAsFactors = T)
dati <- na.omit(dati)

dati <- dati %>%
  select(Peso, everything())

summary(dati)
##       Peso        Anni.madre     N.gravidanze       Fumatrici     
##  Min.   : 830   Min.   : 0.00   Min.   : 0.0000   Min.   :0.0000  
##  1st Qu.:2990   1st Qu.:25.00   1st Qu.: 0.0000   1st Qu.:0.0000  
##  Median :3300   Median :28.00   Median : 1.0000   Median :0.0000  
##  Mean   :3284   Mean   :28.16   Mean   : 0.9812   Mean   :0.0416  
##  3rd Qu.:3620   3rd Qu.:32.00   3rd Qu.: 1.0000   3rd Qu.:0.0000  
##  Max.   :4930   Max.   :46.00   Max.   :12.0000   Max.   :1.0000  
##    Gestazione      Lunghezza         Cranio    Tipo.parto Ospedale   Sesso   
##  Min.   :25.00   Min.   :310.0   Min.   :235   Ces: 728   osp1:816   F:1256  
##  1st Qu.:38.00   1st Qu.:480.0   1st Qu.:330   Nat:1772   osp2:849   M:1244  
##  Median :39.00   Median :500.0   Median :340              osp3:835           
##  Mean   :38.98   Mean   :494.7   Mean   :340                                 
##  3rd Qu.:40.00   3rd Qu.:510.0   3rd Qu.:350                                 
##  Max.   :43.00   Max.   :565.0   Max.   :390
head(dati, 5)
##   Peso Anni.madre N.gravidanze Fumatrici Gestazione Lunghezza Cranio Tipo.parto
## 1 3380         26            0         0         42       490    325        Nat
## 2 3150         21            2         0         39       490    345        Nat
## 3 3640         34            3         0         38       500    375        Nat
## 4 3690         28            1         0         41       515    365        Nat
## 5 3700         20            0         0         38       480    335        Nat
##   Ospedale Sesso
## 1     osp3     M
## 2     osp1     F
## 3     osp2     M
## 4     osp2     M
## 5     osp3     F
n <- nrow(dati)
cat(n)
## 2500

Il dataset contiene 2500 osservazioni totali e le seguenti 10 variabili:

  • Età della madre: misura dell’età in anni
  • Numero di gravidanze: numero di gravidanze avute dalla madre
  • Fumo materno: indicatore binario (0 = non fumatrice, 1 = fumatrice)
  • Durata della gravidanza: numero di settimane di gestazione
  • Peso del neonato: peso alla nascita in grammi
  • Lunghezza: lunghezza del neonato
  • Diametro del cranio: diametro craniale
  • Tipo di parto: naturale o cesareo
  • Ospedale di nascita: ospedale 1, 2 o 3
  • Sesso del neonato: maschio (M) o femmina (F)

2. Analisi descrittiva

In questa sezione si analizzerà ogni variabile in modo univariato, fornendo una descrizione generale di ciascuna.

# Funzione per analisi descrittiva di una variabile
descrizione_variabile <- function(x, nome) {
  stats <- data.frame(
    Statistica = c("Minimo", "1° Quartile", "Mediana", "3° Quartile", "Massimo",
                   "Media", "Deviazione standard", "Asimmetria", "Curtosi"),
    Valore = c(round(min(x),2),
               round(quantile(x, 0.25),2),
               round(median(x),2),
               round(quantile(x, 0.75),2),
               round(max(x),2),
               round(mean(x),2),
               round(sd(x),2),
               round(skewness(x),2),
               round(kurtosis(x)-3,2))
  )
  kable(stats, caption = paste("Statistiche descrittive per", nome), align = "l")
}

2.1 Analisi della variabile Peso

La variabile Peso è di tipo quantitativo continuo e rappresenta il peso dei neonati alla nascita (misurato in grammi).

descrizione_variabile(dati$Peso, "Peso")
Statistiche descrittive per Peso
Statistica Valore
Minimo 830.00
1° Quartile 2990.00
Mediana 3300.00
3° Quartile 3620.00
Massimo 4930.00
Media 3284.08
Deviazione standard 525.04
Asimmetria -0.65
Curtosi 2.03 
hist(dati$Peso, main="Distribuzione Peso Neonati", xlab="Peso (g)", col="slategray")

boxplot(dati$Peso, main="Boxplot Peso (g)", col="slategray")

La variabile Peso mostra una leggera asimmetria negativa ed è leptocurtica. L’indice di curtosi positivo indica che la distribuzione è più concentrata intorno alla media rispetto a una distribuzione normale. Dal boxplot si osservano diversi outlier, in particolare valori più bassi.

2.2 Analisi della variabile Anni.madre

La variabile Anni.madre è quantitativa continua e rappresenta l’età della madre al momento del parto.

descrizione_variabile(dati$Anni.madre, "Anni madre")
Statistiche descrittive per Anni madre
Statistica Valore
Minimo 0.00
1° Quartile 25.00
Mediana 28.00
3° Quartile 32.00
Massimo 46.00
Media 28.16
Deviazione standard 5.27
Asimmetria 0.04
Curtosi 0.38 
hist(dati$Anni.madre, main="Distribuzione Età Madre", xlab="Età", col="firebrick")

boxplot(dati$Anni.madre, main="Boxplot Età Madre", col="firebrick")

Questa variabile è sostanzialmente simmetrica e leggermente leptocurtica. Si osserva la presenza di pochi outliers. La media è di circa 28 anni, la deviazione standard di 5. L’età delle madri varia mediamente quindi tra i 23 e i 33 anni.

Dal boxplot e dai dati sommari si nota un valore minimo di 0 e controllando il dataset si osserva anche un altro valore anomalo (pari a 1). Questi valori sono sicuramente degli errori.

2.3 Analisi della variabile N.gravidanze

È una variabile quantitativa discreta che indica il numero di gravidanze già avute dalla madre.

descrizione_variabile(dati$N.gravidanze, "Numero gravidanze")
Statistiche descrittive per Numero gravidanze
Statistica Valore
Minimo 0.00
1° Quartile 0.00
Mediana 1.00
3° Quartile 1.00
Massimo 12.00
Media 0.98
Deviazione standard 1.28
Asimmetria 2.51
Curtosi 10.99 
hist(dati$N.gravidanze, main="Distribuzione Numero Gravidanze", xlab="Numero di gravidanze", col="goldenrod")

boxplot(dati$N.gravidanze, main="Boxplot Età Madre", col="goldenrod")

La mediana è 1 e la media si avvicina anch’essa a 1. Nonostante ciò, la distribuzione presenta una forte asimmetria positiva (coda lunga a destra) ed è nettamente leptocurtica. Questo indica che la maggior parte delle madri ha avuto poche gravidanze, mentre pochi casi presentano un numero elevato di gravidanze.

2.4 Analisi della variabile Fumatrici

È una variabile categorica dicotomica (dummy). Per questa variabile, di fatto qualitativa anche se rappresentata con numeri (0 e 1) ha senso calcolare la tabella delle frequenze.

kable(as.data.frame(prop.table(table(dati$Fumatrici))),
      col.names = c("Madre fumatrice", "Frequenza relativa"),
      digits = 2,
      align = "l",
      caption = "Frequenze relative delle madri fumatrici")
Frequenze relative delle madri fumatrici
Madre fumatrice Frequenza relativa
0 0.96
1 0.04 
barplot(table(dati$Fumatrici),
        names.arg = c("Non fumatrice","Fumatrice"),
        col = c("lightblue","salmon"),
        main = "Madri fumatrici")

La maggior parte delle madri del campione in oggetto non è fumatrice (quasi il 96%).

2.5 Analisi della variabile Gestazione

È una variabile quantitativa continua che indica il numero di settimane di gestazione.

descrizione_variabile(dati$Gestazione, "Gestazione")
Statistiche descrittive per Gestazione
Statistica Valore
Minimo 25.00
1° Quartile 38.00
Mediana 39.00
3° Quartile 40.00
Massimo 43.00
Media 38.98
Deviazione standard 1.87
Asimmetria -2.07
Curtosi 8.26 
hist(dati$Gestazione, main="Gestazione", xlab="Gestazione (settimane)", col="mediumseagreen")

boxplot(dati$Gestazione, main="Gestazione", col="mediumseagreen")

La gestazione media è di 39 settimane, si osserva però una netta asimmetria negativa e una distribuzione decisamente leptocurtica. L’asimmetria negativa indica la presenza di alcune gestazioni più brevi rispetto alla media, che generano una coda sinistra della distribuzione.

2.6 Analisi della variabile Lunghezza

È quantitativa continua e misura la lunghezza del neonato alla nascita (in millimetri).

descrizione_variabile(dati$Lunghezza, "Lunghezza")
Statistiche descrittive per Lunghezza
Statistica Valore
Minimo 310.00
1° Quartile 480.00
Mediana 500.00
3° Quartile 510.00
Massimo 565.00
Media 494.69
Deviazione standard 26.32
Asimmetria -1.51
Curtosi 6.49 
hist(dati$Lunghezza, main="Lunghezza", xlab="Lunghezza (mm)", col="steelblue")

boxplot(dati$Lunghezza, main="Lunghezza (mm)", col="steelblue")

La media della lunghezza dei neonati è di quasi 495 mm. La distribuzione è negativamente asimmetrica e l’indice di curtosi indica una distribuzione leptocurtica.

2.7 Analisi della variabile Cranio

È anche questa una variabile quantitiativa continua, misurata in millimetri, che indica la lunghezza del diametro del cranio del neonato.

descrizione_variabile(dati$Cranio, "Cranio")
Statistiche descrittive per Cranio
Statistica Valore
Minimo 235.00
1° Quartile 330.00
Mediana 340.00
3° Quartile 350.00
Massimo 390.00
Media 340.03
Deviazione standard 16.43
Asimmetria -0.79
Curtosi 2.95 
hist(dati$Cranio, main="Diametro del Cranio", xlab="Diametro (mm)", col="darkorange")

boxplot(dati$Cranio, main="Diametro del Cranio", col="darkorange")

Anche quesa distribuzione è asimmetrica (negativa) seppur più vicino alla simmetria rispetto alle altre distribuzioni. L’indice di curtosi è nuovamente positivo indicando una distribuzione leptocurtica.

2.8 Analisi della variabile Tipo.parto

Questa è una variabile qualitativa nominale (non ordinabile) e indica le due tipologie di parto: cesareo e naturale.

kable(as.data.frame(prop.table(table(dati$Tipo.parto))),
      col.names = c("Tipo di parto", "Frequenza relativa"),
      digits = 2,
      align = "l",
       caption = "Frequenze relative del tipo di parto")
Frequenze relative del tipo di parto
Tipo di parto Frequenza relativa
Ces 0.29
Nat 0.71 
barplot(table(dati$Tipo.parto),
        names.arg = c("Cesareo","Naturale"),
        col = c("lightblue","lightgreen"),
        main = "Tipologia di parto")

Dai dati si può osservare che più del 70% dei parti è di tipo naturale.

2.9 Analisi della variabile Ospedale

È una variabile qualitativa nominale che riassume i tre diversi ospedali in cui sono avvenuti i 2500 parti oggetto del dataset.

kable(as.data.frame(prop.table(table(dati$Ospedale))),
      col.names = c("Ospedale", "Frequenza relativa"),
      digits = 2,
      align = "l",
      caption = "Frequenze relative degli ospedali")
Frequenze relative degli ospedali
Ospedale Frequenza relativa
osp1 0.33
osp2 0.34
osp3 0.33 
barplot(table(dati$Ospedale),
        names.arg = c("Ospedale 1","Ospedale 2","Ospedale 3"),
        col = c("pink","pink1","pink2"),
        main = "Distribuzione per ospedale")

# Funzione per il calcolo dell'indice di eterogeneità di Gini
gini.index <- function(x){
  ni <- table(x)
  fi <- ni / length(x)
  fi2 <- fi^2
  J <- length(ni)
  
  gini <- 1 - sum(fi2)
  gini.normalizzato <- gini / ((J - 1) / J)
  
  return(gini.normalizzato)
}

# Indice normalizzato di Gini
gini_ospedale <- gini.index(dati$Ospedale)
cat(gini_ospedale)
## 0.9998683

L’indice di Gini normalizzato vicinissimo a 1 indica quasi la massima eterogeneità possibile,con le tre diverse modalità equamente rappresentate dal campione (circa il 33% ciascuna).

2.10 Analisi della variabile Sesso

È qualitativa nominale, indica il sesso dei neonati.

kable(as.data.frame(prop.table(table(dati$Sesso))),
      col.names = c("Sesso", "Frequenza relativa"),
      digits = 4,
      align = "l",
      caption = "Frequenze relative del genere")
Frequenze relative del genere
Sesso Frequenza relativa
F 0.5024
M 0.4976 
barplot(table(dati$Sesso),
        names.arg = c("Femmina", "Maschio"),
        col = c("pink", "lightblue"),
        main = "Distribuzione per genere")

Il dataset contiene un numero di neonati femmina appena superiore a quelli maschi (il 50,24% del campione è di sesso femminile).

3 Analisi bivariata

3.1 Relazione tra tipo di ospedale e tipo di parto

Per analizzare la relazione tra ospedale e tipo di parto, si crea prima la tabella delle frequenze per le due variabili categoriche:

tabella <- table(dati$Ospedale, dati$Tipo.parto)

kable(tabella, 
      digits = 2,
      align = "l",
      caption = "Frequenze assolute tra ospedale e tipo di parto")
Frequenze assolute tra ospedale e tipo di parto
Ces Nat
osp1 242 574
osp2 254 595
osp3 232 603 
tabella_rel <- prop.table(tabella)

kable(tabella_rel, 
      digits = 2,
      align = "l",
      caption = "Frequenze relative tra ospedale e tipo di parto")
Frequenze relative tra ospedale e tipo di parto
Ces Nat
osp1 0.10 0.23
osp2 0.10 0.24
osp3 0.09 0.24 
ggpubr::ggballoonplot(data=as.data.frame(tabella),
                      fill="blue")

Dato che entrambe le variabili sono categoriche, si effettua un test Chi-quadro per verificare l’indipendenza.

test.indipendenza <- chisq.test(tabella)
test.indipendenza
## 
##  Pearson's Chi-squared test
## 
## data:  tabella
## X-squared = 1.0972, df = 2, p-value = 0.5778

Il test restituisce un p-value di 0.58, quindi non si rifiuta l’ipotesi nulla di indipendenza tra ospedale e tipo di parto. Non emergono evidenze statistiche che la frequenza dei parti cesarei differisca tra i tre ospedali.

3.2 Confronto di peso e lunghezza con i valori noti della popolazione

Per confrontare il campione con valori di riferimento della popolazione, sono stati utilizzati gli standard di crescita dell’Organizzazione Mondiale della Sanità (OMS/WHO), che separa però i dati tra i due generi. I valori sono consultabili ai seguenti link:

Dato che il campione in analisi comprende sia maschi che femmine, è stata calcolata la media delle mediane maschile e femminile (per peso e lunghezza) dalle tabelle di riferimento dell’OMS allegate. Convertendo i valori nelle stesse unità del dataset, si è ottenuto:

  • Peso medio di riferimento: 3289 g
  • Lunghezza media di riferimento: 495 mm

Questi valori sono stati utilizzati come riferimento per i t-test sul peso e sulla lunghezza dei neonati.

# Calcolo della media dei pesi di riferimento
peso_medio_M <- 3.3464*1000 # g maschi
peso_medio_F <- 3.2322*1000 # g femmine
peso_medio <- round((peso_medio_M + peso_medio_F)/2)
cat(peso_medio)
## 3289
# t-test sul peso rispetto alla media di riferimento
t.test(dati$Peso,
       mu = peso_medio,
       conf.level = 0.95,
       alternative = "two.sided")
## 
##  One Sample t-test
## 
## data:  dati$Peso
## t = -0.46846, df = 2499, p-value = 0.6395
## alternative hypothesis: true mean is not equal to 3289
## 95 percent confidence interval:
##  3263.490 3304.672
## sample estimates:
## mean of x 
##  3284.081

Con un p-value di 0.64, non si rifiuta l’ipotesi nulla: il peso medio dei neonati del campione non è significativamente diverso dalla media della popolazione.

# Calcolo della media delle lunghezze di riferimento
lunghezza_media_M <- 49.8842*10 # mm maschi
lunghezza_media_F <- 49.1477*10 # mm femmine
lunghezza_media <- round((lunghezza_media_M + lunghezza_media_F)/2)
cat(lunghezza_media)
## 495
# t-test sulla lunghezza rispetto alla media di riferimento
t.test(dati$Lunghezza,
       mu = lunghezza_media,
       conf.level = 0.95,
       alternative = "two.sided")
## 
##  One Sample t-test
## 
## data:  dati$Lunghezza
## t = -0.58514, df = 2499, p-value = 0.5585
## alternative hypothesis: true mean is not equal to 495
## 95 percent confidence interval:
##  493.6598 495.7242
## sample estimates:
## mean of x 
##   494.692

Con un p-value di 0.56, non si rifiuta l’ipotesi nulla: il campione non si discosta significativamente dalla popolazione nemmeno per quanto riguarda la lunghezza media alla nascita.

3.3 Confronto tra le misure antropometriche condizionate al genere

Si effettua un t-test tra due campioni indipendenti all’interno del campione di riferimento, partendo dal Peso.

boxplot(dati$Peso~dati$Sesso, main = "Peso condizionato al genere", xlab = "Sesso", ylab = "Peso (g)")

t.test(dati$Peso~dati$Sesso, data=dati)
## 
##  Welch Two Sample t-test
## 
## data:  dati$Peso by dati$Sesso
## t = -12.106, df = 2490.7, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -287.1051 -207.0615
## sample estimates:
## mean in group F mean in group M 
##        3161.132        3408.215

Il t-test della variabile Peso condizionata al Sesso del neonato restituisce un p-value pari a zero e un intervallo di confidenza che non contiene 0, indicando di accettare l’ipotesi alternativa H1: la media del Peso tra i due generi è sensibilmente diversa.

Si procede analizzando la Lunghezza condizionata al genere.

boxplot(dati$Lunghezza~dati$Sesso, main = "Lunghezza condizionata al genere", xlab = "Sesso", ylab = "Lunghezza (mm)")

t.test(dati$Lunghezza~dati$Sesso, data=dati)
## 
##  Welch Two Sample t-test
## 
## data:  dati$Lunghezza by dati$Sesso
## t = -9.582, df = 2459.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -11.929470  -7.876273
## sample estimates:
## mean in group F mean in group M 
##        489.7643        499.6672

Anche in questo caso si arriva alla stessa conclusione: il p-value nullo impone di rifiutare l’ipotesi nulla di uguaglianza delle medie della Lunghezza tra i due generi.

Infine, si studia la variabile Cranio.

boxplot(dati$Cranio~dati$Sesso, main = "Diametro del cranio condizionato al genere", xlab = "Sesso", ylab = "Diametro (mm)")

t.test(dati$Cranio~dati$Sesso, data=dati)
## 
##  Welch Two Sample t-test
## 
## data:  dati$Cranio by dati$Sesso
## t = -7.4102, df = 2491.4, p-value = 1.718e-13
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -6.089912 -3.541270
## sample estimates:
## mean in group F mean in group M 
##        337.6330        342.4486

Il p-value molto basso indica che la differenza tra maschio e femmina è statisticamente significativa anche per quanto riguarda la misura del diametro del cranio con quasi 5 mm di differenza tra i due sessi.

4 Modello di regressione lineare

In questo modello lineare la variabile risposta è il peso. Si effettua preeliminarmente un test di normalità sulla variabile risposta.

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

Il test dà esito negativo, nel senso che la variabile Peso non si dispone normalmente. Come si è visto nella sezione sull’analisi descrittiva, questo probabilmente avviene a causa della lunga coda a sinistra mostrata dalla distribuzione della variabile, osservabile dal suo istogramma e legata all’asimmetria negativa. È probabile che questo abbia degli effetti sulla normalità dei residui.

Tramite la seguente funzione si indagano le corelazioni lineari tra le variabili.

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 <- 1.4
  text(0.5, 0.5, txt, cex = cex.cor * r)
}

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

Dalla tabella che mette in relazione le variabili, che mostra sia il coefficiente di correlazione di Pearson sia gli scatterplot, si evidenzia una correlazione significativa della variabile risposta Peso con le variabili Lunghezza (0.80), Cranio (0.70) e Gestazione (0.59). Non sembra invece significativa la relazione con le altre variabili quantitative Anni.madre e N.gravidanze.

Inoltre, la variabile Lunghezza sembra correlata in maniera moderata con Gestazione (0.62) e Cranio (0.60).

Per capire meglio l’eventuale correlazione di Peso con le variabili categoriche, si analizzano i relativi boxplot rispetto a Ospedale, Tipo.parto e Fumatrici.

La correlazione con il genere è già stata evidenziata ed è risultata significativa. Quindi Sesso sarà una variabile esplicativa importante da considerare nel modello lineare, quantomeno come variabile di controllo.

4.1 Correlazione tra Peso e Tipo.parto

boxplot(dati$Peso~dati$Tipo.parto, main = "Peso condizionato al tipo di parto", xlab = "Tipo di parto", ylab = "Peso (g)")

t.test(dati$Peso~dati$Tipo.parto)
## 
##  Welch Two Sample t-test
## 
## data:  dati$Peso by dati$Tipo.parto
## t = -0.12968, df = 1493, p-value = 0.8968
## alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
## 95 percent confidence interval:
##  -46.27992  40.54037
## sample estimates:
## mean in group Ces mean in group Nat 
##          3282.047          3284.916

Al contrario, Tipo.parto (cesareo o naturale) non influisce sensibilmente sul Peso del neonato. Il test restituisce un p-value pari a 0.90, quindi non si rifiuta l’ipotesi nulla di uguaglianza tra le medie nei due casi. La differenza tra le medie delle due modalità di Tipo.parto non è quindi significativamente diversa da zero. Questa variabile potrà essere rimossa dal modello lineare.

4.2 Correlazione tra Peso e Ospedali

In questo caso vanno eseguiti confronti multipli, essendo gli ospedali tre.

boxplot(dati$Peso~dati$Ospedale, main = "Peso condizionato all'ospedale", xlab = "Ospedale", ylab = "Peso (g)")

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

I p-value dei confronti multipli sono tutti superiori al livello di significatività del 5%, quindi non emergono differenze statisticamente significative tra i pesi medi dei neonati nei tre livelli della variabile Ospedale.

Anche questa variabile potrà quindi essere ignorata nel modello lineare.

4.3 Correlazione tra Peso e Fumatrici

boxplot(dati$Peso~dati$Fumatrici, main = "Peso condizionato a Fumatrici", xlab = "Madre fumatrice", ylab = "Peso (g)")

t.test(dati$Peso~dati$Fumatrici)
## 
##  Welch Two Sample t-test
## 
## data:  dati$Peso by dati$Fumatrici
## t = 1.034, df = 114.1, p-value = 0.3033
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -45.61354 145.22674
## sample estimates:
## mean in group 0 mean in group 1 
##        3286.153        3236.346

Con un p-value pari a 0.30 non si rifiuta l’ipotesi nulla di uguaglianza tra le medie del peso dei neonati rispetto alla variabile Fumatrici.

In questo campione, quindi, non emerge un effetto statisticamente significativo del fumo materno sul peso alla nascita.

4.4 Correlazione tra Peso e Gestazione

Sono entrambe variaibili quantitative, quindi si procede con la correlazione lineare.

mod_lin <- lm(Peso ~ Gestazione, data=dati)
summary(mod_lin)
## 
## Call:
## lm(formula = Peso ~ Gestazione, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1603.61  -287.34   -11.07   280.46  1897.75 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3197.250    176.851  -18.08   <2e-16 ***
## Gestazione    166.272      4.532   36.69   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 423.3 on 2498 degrees of freedom
## Multiple R-squared:  0.3502, Adjusted R-squared:  0.3499 
## F-statistic:  1346 on 1 and 2498 DF,  p-value: < 2.2e-16
plot(dati$Gestazione, dati$Peso,
     xlab = "Settimane di gestazione",
     ylab = "Peso (g)",
     main = "Peso vs Gestazione",
     pch = 20, col = "steelblue")
abline(mod_lin, col="red", lwd=2)

Lo scatterplot mette in relazione le due variabili Peso e Gestazione. Il valore di \(\beta_1\) pari a 166 rappresenta il coefficiente angolare della retta di regressione: ciò significa che, in media, per ogni settimana aggiuntiva di gestazione il peso medio dei neonati aumenta di circa 166 grammi.

5 Regressione lineare multipla

5.1 Modello 1

Si parte, comunque, inserendo tutte le variabili esplicative nel modello, che verranno eliminate progressivamente nei successivi modelli in base alle analisi fin qui svolte. Si utilizza quindi la procedura stepwise.

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

Si nota un Adjusted \(R^2\) del modello 1 pari a 0.7278 e una forte relazione con le variabili Gestazione, Lunghezza, Cranio e Sesso. In particolare, a parità delle altre variabili, un neonato maschio pesa in media circa 77 grammi in più rispetto a una femmina.

Emerge inoltre una significatività più debole per N.gravidanze, Tipo.partoNat (con circa 29 grammi di differenza rispetto al parto cesareo) e Ospedale3.

Non si osserva invece una relazione statisticamente significativa con Anni.madre e Fumatrici.

5.2 Modello 2

Come anticipato, si decide quindi di eliminare le due variabili qualitative Ospedale e Tipo.parto, ritenute poco informative sulla base dell’analisi precedente.

mod2 <- update(mod1, ~.-Ospedale - Tipo.parto)
summary(mod2)
## 
## Call:
## lm(formula = Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + 
##     Lunghezza + Cranio + Sesso, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1161.56  -181.19   -15.75   163.70  2630.75 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6714.4109   141.1515 -47.569  < 2e-16 ***
## Anni.madre       0.9585     1.1347   0.845   0.3984    
## N.gravidanze    11.2756     4.6690   2.415   0.0158 *  
## Fumatrici      -30.2959    27.5971  -1.098   0.2724    
## Gestazione      32.9331     3.8267   8.606  < 2e-16 ***
## Lunghezza       10.2342     0.3009  34.009  < 2e-16 ***
## Cranio          10.5177     0.4268  24.642  < 2e-16 ***
## SessoM          78.0845    11.2039   6.969 4.06e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2492 degrees of freedom
## Multiple R-squared:  0.7272, Adjusted R-squared:  0.7264 
## F-statistic:   949 on 7 and 2492 DF,  p-value: < 2.2e-16
anova(mod2, mod1)
## Analysis of Variance Table
## 
## Model 1: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso
## Model 2: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Tipo.parto + Ospedale + Sesso
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)   
## 1   2492 187919851                               
## 2   2489 186762521  3   1157330 5.1413 0.00152 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BIC(mod2,mod1)
##      df      BIC
## mod2  9 35233.81
## mod1 12 35241.84

Il modello 2 ha un valore di Adjusted \(R^2\) pari a 0.7264, quindi inferiore di meno di due millesimi rispetto al modello 1, e presenta inoltre un BIC migliore.

Dal confronto tramite ANOVA risulta che l’inclusione delle variabili Tipo.parto e Ospedale migliora significativamente il modello (\(p = 0.0015\)). In particolare Tipo.parto risulta debolmente significativo, così come uno dei livelli della variabile Ospedale (osp3), mentre osp2 non risulta significativo.

Tuttavia il modello completo introduce 3 gradi di libertà aggiuntivi e l’incremento di Adjusted \(R^2\) è trascurabile (0.0014). Per principio di parsimonia si preferisce quindi continuare con il modello ridotto.

5.3 Modello 3

Si procede quindi con l’eliminazione della variabile Anni.madre.

mod3 <- update(mod2, ~.-Anni.madre)
summary(mod3)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso, data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1150.3  -181.3   -15.7   163.0  2636.3 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.6714   135.7178 -49.232  < 2e-16 ***
## N.gravidanze    12.7185     4.3450   2.927  0.00345 ** 
## Fumatrici      -30.4634    27.5948  -1.104  0.26972    
## Gestazione      32.5914     3.8051   8.565  < 2e-16 ***
## Lunghezza       10.2341     0.3009  34.011  < 2e-16 ***
## Cranio          10.5359     0.4262  24.718  < 2e-16 ***
## SessoM          78.1713    11.2028   6.978 3.83e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.6 on 2493 degrees of freedom
## Multiple R-squared:  0.7271, Adjusted R-squared:  0.7265 
## F-statistic:  1107 on 6 and 2493 DF,  p-value: < 2.2e-16
anova(mod3, mod2)
## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Sesso
## Model 2: Peso ~ Anni.madre + N.gravidanze + Fumatrici + Gestazione + Lunghezza + 
##     Cranio + Sesso
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1   2493 187973654                           
## 2   2492 187919851  1     53803 0.7135 0.3984
BIC(mod3, mod2)
##      df      BIC
## mod3  8 35226.70
## mod2  9 35233.81

L’Adjusted \(R^2\) ora vale 0.7265, sostanzialmente identico al modello 2. L’ANOVA mostra che la rimozione della variabile Anni.madre non modifica significativamente il modello lineare, mentre il BIC migliora leggermente. Si preferisce quindi proseguire con il modello 3 e si procede all’eliminazione della variabile Fumatrici.

5.4 Modello 4

mod4 <- update(mod3, ~.-Fumatrici)
summary(mod4)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso, data = dati)
## 
## 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
anova(mod4, mod3)
## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Fumatrici + Gestazione + Lunghezza + Cranio + 
##     Sesso
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1   2494 188065546                           
## 2   2493 187973654  1     91892 1.2187 0.2697
BIC(mod4, mod3)
##      df     BIC
## mod4  7 35220.1
## mod3  8 35226.7

Anche per questo modello, il valore di Adjusted \(R^2\) rimane sostanzialmente costante rispetto al precedente (0.7265). L’ANOVA conferma la non significatività del parametro (p-value = 0.27) e il BIC migliora ulteriormente. La variabile Fumatrici, così come Anni.madre, non risulta statisticamente significativa.

Durante la selezione del modello, alcune variabili poco informative (Anni della madre, Fumatrici, Tipo parto e Ospedale) sono quindi state eliminate seguendo i criteri di parsimonia e BIC. Questa riduzione ha rafforzato l’effetto stimato di N.gravidanze, che ora risulta più significativo.

5.5 Modello 5

Si prova a eliminare la variabile N.gravidanze.

mod5 <- update(mod4, ~.-N.gravidanze)
summary(mod5)
## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso, 
##     data = dati)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1138.2  -184.3   -17.6   163.3  2627.3 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6651.1188   135.5172 -49.080  < 2e-16 ***
## Gestazione     31.2737     3.7856   8.261 2.31e-16 ***
## Lunghezza      10.2054     0.3007  33.939  < 2e-16 ***
## Cranio         10.6704     0.4245  25.139  < 2e-16 ***
## SessoM         79.1049    11.2117   7.056 2.22e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 275 on 2495 degrees of freedom
## Multiple R-squared:  0.7261, Adjusted R-squared:  0.7257 
## F-statistic:  1654 on 4 and 2495 DF,  p-value: < 2.2e-16
anova(mod5, mod4)
## Analysis of Variance Table
## 
## Model 1: Peso ~ Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##   Res.Df       RSS Df Sum of Sq      F   Pr(>F)   
## 1   2495 188688687                                
## 2   2494 188065546  1    623141 8.2637 0.004079 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BIC(mod5, mod4)
##      df      BIC
## mod5  6 35220.54
## mod4  7 35220.10

Eliminando N.gravidanze, il valore di Adjusted \(R^2\) rimane sostanzialmente costante, ma il BIC peggiora leggermente, quindi non è giustificata la rimozione di questa variabile. Si mantiene il modello 4.

5.6 Modello 6

Dagli scatterplot sembrava emergere un effetto non lineare tra Peso e Gestazione. Per verificarlo, si aggiunge una dipendenza quadratica della Gestazione nel modello.

mod6 <- update(mod4, ~.+I(Gestazione^2))
summary(mod6)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + I(Gestazione^2), data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1144.11  -181.17   -13.16   165.73  2662.56 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -4650.2268   898.1706  -5.177 2.43e-07 ***
## N.gravidanze       12.5660     4.3361   2.898  0.00379 ** 
## Gestazione        -81.0486    49.7127  -1.630  0.10316    
## Lunghezza          10.3534     0.3039  34.074  < 2e-16 ***
## Cranio             10.6363     0.4280  24.854  < 2e-16 ***
## SessoM             75.7900    11.2340   6.747 1.88e-11 ***
## I(Gestazione^2)     1.5136     0.6617   2.287  0.02226 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.4 on 2493 degrees of freedom
## Multiple R-squared:  0.7276, Adjusted R-squared:  0.7269 
## F-statistic:  1110 on 6 and 2493 DF,  p-value: < 2.2e-16
anova(mod4, mod6)
## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso + 
##     I(Gestazione^2)
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
## 1   2494 188065546                              
## 2   2493 187671672  1    393875 5.2322 0.02226 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BIC(mod6, mod4)
##      df      BIC
## mod6  8 35222.68
## mod4  7 35220.10

L’Adjusted \(R^2\) vale ora 0.7269, migliorando di poco, e anche il BIC aumenta leggermente, dimostrando che l’effetto curvilineo è debole. Si decide quindi di ignorare l’effetto quadratico per mantenere un modello più semplice e parsimonioso.

5.7 Modello 7

Data la moderata correlazione tra Lunghezza e Cranio, si verifica se esiste un effetto interattivo tra le due variabili esplicative sul Peso. Questo permette di testare se l’effetto della Lunghezza sul risultato dipende dal valore di Cranio (o viceversa), cioè se esiste un effetto moltiplicativo tra le due variabili.

mod7 <- update(mod4, ~.+Lunghezza:Cranio)
summary(mod7)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Lunghezza:Cranio, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1150.74  -180.48   -13.62   165.82  2866.17 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -1.801e+03  1.018e+03  -1.770  0.07685 .  
## N.gravidanze      1.295e+01  4.321e+00   2.996  0.00276 ** 
## Gestazione        3.810e+01  3.964e+00   9.610  < 2e-16 ***
## Lunghezza        -3.047e-01  2.202e+00  -0.138  0.88996    
## Cranio           -4.759e+00  3.191e+00  -1.491  0.13601    
## SessoM            7.327e+01  1.119e+01   6.545 7.20e-11 ***
## Lunghezza:Cranio  3.158e-02  6.528e-03   4.838 1.39e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.4 on 2493 degrees of freedom
## Multiple R-squared:  0.7295, Adjusted R-squared:  0.7289 
## F-statistic:  1121 on 6 and 2493 DF,  p-value: < 2.2e-16
BIC(mod7, mod4)
##      df      BIC
## mod7  8 35204.57
## mod4  7 35220.10
anova(mod7, mod4)
## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso + 
##     Lunghezza:Cranio
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1   2493 186316621                                  
## 2   2494 188065546 -1  -1748926 23.401 1.395e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

L’effetto combinato di Lunghezza e Cranio sul Peso non è semplicemente additivo. Neonati più lunghi con cranio più grande hanno un peso maggiore di quanto ci si aspetterebbe sommando gli effetti singoli. Questo giustifica l’inclusione dell’interazione tra le due variabili nel modello.

5.8 Modello 8

Intuendo però una dipendenza della Lunghezza dalla Gestazione, si verifica se l’effetto combinato di questa interazione sia superiore a quella appena osservata.

mod8 <- update(mod4, ~.+Gestazione:Lunghezza)
summary(mod8)
## 
## Call:
## lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
##     Sesso + Gestazione:Lunghezza, data = dati)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1133.50  -179.45   -11.74   168.77  2653.34 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -1.991e+03  9.202e+02  -2.163 0.030619 *  
## N.gravidanze          1.304e+01  4.319e+00   3.020 0.002551 ** 
## Gestazione           -9.396e+01  2.480e+01  -3.789 0.000155 ***
## Lunghezza            -8.111e-02  2.027e+00  -0.040 0.968082    
## Cranio                1.076e+01  4.262e-01  25.245  < 2e-16 ***
## SessoM                7.228e+01  1.120e+01   6.454 1.31e-10 ***
## Gestazione:Lunghezza  2.729e-01  5.295e-02   5.153 2.76e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 273.2 on 2493 degrees of freedom
## Multiple R-squared:  0.7299, Adjusted R-squared:  0.7292 
## F-statistic:  1123 on 6 and 2493 DF,  p-value: < 2.2e-16
anova(mod8, mod4)
## Analysis of Variance Table
## 
## Model 1: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso + 
##     Gestazione:Lunghezza
## Model 2: Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1   2493 186083565                                  
## 2   2494 188065546 -1  -1981981 26.553 2.765e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

L’Adjusted \(R^2\) vale 0.7292 (il modello spiega quasi il 73% della variabilità) e l’ANOVA conferma la significatività dell’interazione. L’effetto combinato di Gestazione e Lunghezza è più forte e più significativo.

BIC(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8)
##      df      BIC
## mod1 12 35241.84
## mod2  9 35233.81
## mod3  8 35226.70
## mod4  7 35220.10
## mod5  6 35220.54
## mod6  8 35222.68
## mod7  8 35204.57
## mod8  8 35201.44

Il confronto multiplo BIC tra tutti i modelli fa preferire l’ultimo, il modello 8, che ha un migliore Adjusted \(R^2\) rispetto al modello 7 e il valor BIC migliore tra tutti.

Si verifica, infine, se ci sono problemi di multicollinearità nei modelli.

vif(mod4)
## N.gravidanze   Gestazione    Lunghezza       Cranio        Sesso 
##     1.023475     1.669189     2.074689     1.624465     1.040054
vif(mod8)
##         N.gravidanze           Gestazione            Lunghezza 
##             1.024145            71.896088            95.264171 
##               Cranio                Sesso Gestazione:Lunghezza 
##             1.640863             1.050345           262.766336

Nel modello finale (mod8) le variabili Gestazione e Lunghezza (e la loro interazione) mostrano una multicollinearità elevata (VIF molto alti) proprio a causa della presenza dell’inserimento del fattore di interazione. Le altre variabili hanno valori inferiori a 5.

Nel modello senza interazioni (mod4), invece, non si osservano problemi di multicollinearità, con tutti i VIF minori di 5.

Si decide, comunque, di utilizzare mod8 sia per il miglior Adjusted \(R^2\), sia per il BIC più parsimonioso.

6 Previsione

Scelto il modello 8 si procede alla previsione del peso di una neonata noti il sesso (femmina), una madre alla terza gravidanza e una gestazione di 39 settimane. Dato che il modello scelto prevede anche la dipendenza dalle variabili Lunghezza e Cranio si inserisce la media femminile.

cat(round(predict(mod8, newdata = data.frame(
  N.gravidanze = 2,
  Gestazione   = 39,
  Lunghezza    = mean(dati$Lunghezza[dati$Sesso=="F"]),
  Cranio       = mean(dati$Cranio[dati$Sesso=="F"]),
  Sesso        = "F"))))
## 3176

Il modello stima che il peso della neonata sarà di 3176 grammi.

7 Analisi dei residui

I residui devono rispettare criteri di normalità, omoschedasticità, incorrelazione e media di zero.

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

Nel primo grafico i residui dovrebbero essere sparsi casualmente attorno a zero e sembrano esserlo. Nel secondo, i residui confrontati con i quantili di una normale seguono la bisettrice, con una netta eccezione come l’osservazione 1551 e qualche valore iniziale e finale (meno evidenti). Nel terzo grafico la dispersione dei residui permette di valutare se la varianza è costante. Nell’ultimo grafico, alcuni punti presentano valori di leverage elevati: valori sopra 0.5 indicano avvertimento, sopra 1 rappresentano allarme. Sempre l’osservazione 1551 sembra dare problemi.

7.1 Valori di leva

I valori di leva (leverage) misurano quanto un’osservazione è distante dalle altre nello spazio delle variabili esplicative. Un leverage elevato indica che l’osservazione possiede combinazioni di valori delle variabili indipendenti relativamente rare nel campione e potrebbe quindi avere una maggiore capacità di influenzare la stima dei coefficienti del modello.

Utilizzando la soglia teorica è possibile individuare le osservazioni potenzialmente influenti. Dal grafico emerge che solo un numero limitato di punti supera tale soglia.

lev <- hatvalues(mod8)
plot(lev)
p = sum(lev)
soglia = 2*p/n
abline(h=soglia, col=2)

cat(sum(lev > soglia))
## 125

125 osservazioni su 2500 osservazioni totali superano la soglia di leverage. Tuttavia, la presenza di valori di leva elevati non implica necessariamente che l’osservazione sia problematica: è necessario verificare anche il loro impatto effettivo sulle stime del modello attraverso altre misure di influenza, come la distanza di Cook.

7.2 Outliers

Gli outliers nei residui rappresentano osservazioni per le quali il valore osservato della variabile risposta si discosta in modo rilevante da quello previsto dal modello. In genere, valori dei residui studentizzati superiori a 2 in valore assoluto indicano potenziali anomalie.

plot(rstudent(mod8))
abline(h=c(-2,2), col=2)

outlierTest(mod8)
##       rstudent unadjusted p-value Bonferroni p
## 1551 10.160278         8.6136e-24   2.1534e-20
## 155   5.042472         4.9259e-07   1.2315e-03
## 1306  4.751668         2.1321e-06   5.3303e-03

Nel modello analizzato emergono solo poche osservazioni che superano queste soglie, indicando che la maggior parte dei dati è ben spiegata dalla relazione lineare stimata. L’outlierTest conferma la presenza di un numero molto limitato di osservazioni anomale (precisamente 3).

In un dataset di dimensioni relativamente grandi come quello analizzato (2500 osservazioni), la presenza di pochi outliers è fisiologica e non rappresenta necessariamente un problema per la stabilità del modello.

7.3 Distanza di Cook

Per valutare l’influenza delle singole osservazioni sulle stime del modello di regressione si utilizza la distanza di Cook, una misura che combina l’informazione proveniente dai residui e dai valori di leva.

cook<-cooks.distance(mod8)
plot(cook)
abline(h=0.5, col=2)

cat(max(cook))
## 0.7277768

Una regola empirica comunemente utilizzata considera potenzialmente influenti le osservazioni con distanza di Cook superiore a 0.5, mentre valori superiori a 1 indicano casi fortemente influenti. Nel modello analizzato si osserva un solo punto con valore relativamente elevato (0.73), ma comunque inferiore alla soglia critica pari a 1 (sempre l’osservazione 1551).

7.4 Test sui residui

Infine si valutano i residui del modello tramite tre test classici: il Breusch-Pagan per l’omoschedasticità, il Durbin-Watson per l’autocorrelazione e il Shapiro-Wilk per la normalità.

bptest(mod8)
## 
##  studentized Breusch-Pagan test
## 
## data:  mod8
## BP = 90.217, df = 6, p-value < 2.2e-16
dwtest(mod8)
## 
##  Durbin-Watson test
## 
## data:  mod8
## DW = 1.9529, p-value = 0.1192
## alternative hypothesis: true autocorrelation is greater than 0
shapiro.test(residuals(mod8))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod8)
## W = 0.97344, p-value < 2.2e-16
plot(density(residuals(mod8)))

Dai risultati emerge che la varianza dei residui non è costante (caso di eteroschedasticità), i residui non mostrano autocorrelazione, ma la loro distribuzione non è perfettamente normale, con una coda destra più lunga Si rifiutano dunque le ipotesi di normalità e di omoschedasticità sui residui.

8 Grafici aggiuntivi

8.1 Grafico Peso-Gestazione

ggplot(data=dati)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col=Sesso), position="jitter")+
  geom_smooth(aes(x=Gestazione,
                  y=Peso,
                  col=Sesso), se=F, method = "lm")+
  geom_smooth(aes(x=Gestazione,
                  y=Peso), col="black", se=F, method = "lm")+
  theme_classic()

Questo grafico mette in relazione il Peso con la Gestazione. Le rette di regressione, per i maschi e per le femmine, hanno la stessa pendenza ma si deduce che i maschi pesano mediamente di più, avendo l’intercetta più in alto.

8.2 Grafico Peso-Cranio

ggplot(data=dati)+
  geom_point(aes(x=Cranio,
                 y=Peso,
                 col=Sesso), position="jitter")+
  geom_smooth(aes(x=Cranio,
                  y=Peso,
                  col=Sesso), se=F, method = "lm")+
  geom_smooth(aes(x=Cranio,
                  y=Peso), col="black", se=F, method = "lm")+
  theme_classic()

Anche per questo grafico, che mette in relazione Peso e Cranio, le pendenze delle rette per i maschi e per le femmine sembrano simili e si conferma la correlazione positiva le due variabili.

9 Conclusioni

L’analisi condotta sul dataset di 2500 neonati ha permesso di sviluppare un modello di regressione lineare multipla per stimare il peso alla nascita basandosi su variabili cliniche e demografiche. Dalla fase di analisi preliminare emerge che alcune variabili quantitative, come Gestazione, Lunghezza e Cranio, sono fortemente correlate con il peso del neonato, mentre variabili come Anni.madre, Fumatrici, Tipo di parto e Ospedale risultano poco influenti. Tra le variabili qualitative, il Sesso del neonato mostra un impatto significativo sul peso, con i maschi mediamente più pesanti delle femmine.

Il modello finale (mod8) include le variabili principali e l’interazione tra Gestazione e Lunghezza, ottenendo il miglior Adjusted \(R^2\) e il valore BIC più parsimonioso tra tutti i modelli implementati. Nonostante la presenza di multicollinearità tra le variabili principali e l’interazione, il modello consente stime interpretabili e predizioni affidabili.

Dall’analisi dei residui si osserva che non vi sono problemi significativi di autocorrelazione (Durbin-Watson). La distribuzione dei residui non è però perfettamente normale, con una coda destra leggermente più pesante (Shapiro-Wilk). La varianza dei residui non è costante (Breusch-Pagan), indicando eteroschedasticità. I valori di leverage e la distanza di Cook evidenziano pochi punti con impatto moderato, ma nessuno influenza in modo critico le stime del modello.

In termini di previsione, il modello permette di stimare con precisione il peso dei neonati in scenari pratici. Ad esempio, stimando il peso di una neonata di una madre alla terza gravidanza con 39 settimane di gestazione, il modello predice circa 3176 grammi, in linea con le medie di riferimento della popolazione.

Nonostante alcune violazioni delle ipotesi classiche della regressione, il modello sviluppato fornisce uno strumento utile per Neonatal Health Solutions per prevedere il peso dei neonati e identificare casi a rischio, pianificare risorse ospedaliere e interventi clinici in modo più mirato.