Modello Statistico per la Previsione del Peso Neonatale

Analisi descrittiva delle variabili

Il dataset neonati.csv contiene dati relativi a 2500 neonati provenienti da tre ospedali.

Si compone di 10 variabili: 6 quantitative e 4 qualitative.

setwd("C:/Users/guido/Desktop/Python/R/Dataset")
neonati = read.csv("neonati.csv", stringsAsFactors=T)
attach(neonati)

library(knitr)
library(ggplot2)
library(patchwork)
library(TeachingDemos)
library(dplyr)
library(DT)
library(kableExtra)
library(hexbin)
library(car)
library(MASS)
library(lmtest)
library(rgl)
library(plotly)

Per ogni variabile quantitativa verrà calcolato

minimo
1° quartile
mediana
media
3° quartile
massimo

Per le variabili qualitative verranno calcolate invece le frequenze assolute e relative.

Verranno poi fatte delle rappresentazioni visive con distribuzioni di frequenza, istogrammi e boxplot.

Anni.madre

Variabile quantitativa: misura dell’età in anni della madre

# Funzione per visualizzare summary
stat <- function(variable){
summ = summary(variable)
df_summ = data.frame(Statistica = names(summ), Valore = as.vector(summ))
kable(df_summ)
}

# Funzione per visualizzare density e boxplot
graf <- function(variable,label){

density = ggplot(data=neonati) +
  geom_density( aes(x=variable), color="red") +
  labs(x=label, y="Densità") + 
  theme_bw()

boxplot = ggplot(data=neonati) +
  geom_boxplot(aes(y=variable), color="blue") +
  labs(y=label) +
  scale_x_continuous(breaks = NULL) +
  theme_bw()

density + boxplot
}

stat(Anni.madre)

Statistica	Valore
Min.	0.000
1st Qu.	25.000
Median	28.000
Mean	28.164
3rd Qu.	32.000
Max.	46.000

graf(Anni.madre,"N° di anni")

L’età media delle madri è di 28 anni con una distribuzione leggermente assimmetrica verso destra, con una coda più lunga di madri più anziane.

La variabile contiene sicuramente 2 anomalie individuate dai valori 0 ed 1.

N.gravidanze

Variabile quantitativa: indica quante gravidanze ha avuto la madre

stat(N.gravidanze)

Statistica	Valore
Min.	0.0000
1st Qu.	0.0000
Median	1.0000
Mean	0.9812
3rd Qu.	1.0000
Max.	12.0000

# Funzione per visualizzare un grafico ad istogramma
ist <- function(variable,label){

ggplot(data=neonati)+
  geom_bar(aes(x=variable),
           stat="count",
           color = "black",
           fill = "grey") + 
  labs(title="",
       x=label,
       y="Frequenze assolute")+
  scale_x_discrete(breaks = unique(variable), 
                   labels = unique(variable)) +
  theme_light() + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

}

ist(N.gravidanze,"N° di gravidanze (0 --> 12)")

La maggior parte delle donne hanno avuto un numero di gravidanze precedenti a quella esaminata compreso tra 0 e 3. La distribuzione poi si allunga e si schiaccia verso destra fino ad arrivare ad una valore massimo di 12.

Fumatrici

Variabile qualitativa: indica se la madre è fumatrice “1” o non fumatrice “0”

# Funzione per visualizzare il n° di record per categoria
freq <- function(variable){
tab = table(variable)
df_tab = data.frame(Categoria = names(tab), Record = as.vector(tab))
df_tab$Percentuale <- round((df_tab$Record / sum(df_tab$Record)) * 100, 1)
kable(df_tab)
}

freq(Fumatrici)

Categoria	Record	Percentuale
0	2396	95.8
1	104	4.2

ist(Fumatrici,"Non fumatrici / fumatrici")

La maggior parte delle madri non sono fumatrici.

Gestazione

Variabile quantitativa: durata della gravidanza in settimane

stat(Gestazione)

Statistica	Valore
Min.	25.0000
1st Qu.	38.0000
Median	39.0000
Mean	38.9804
3rd Qu.	40.0000
Max.	43.0000

graf(Gestazione,"N° di settimane")

Il tempo di gestazione medio è di circa 39 settimane. Distribuzione con tanti outlier nella parte bassa corrispondenti ai parti prematuri.

Peso

Variabile quantitativa: peso del neonato in grammi

stat(Peso)

Statistica	Valore
Min.	830.000
1st Qu.	2990.000
Median	3300.000
Mean	3284.081
3rd Qu.	3620.000
Max.	4930.000

graf(Peso,"Grammi")

Distribuzione abbastanza simmetrica con la maggior parte di bambini che pesa fra 3 e 3,6 kg circa, tanti outlier nella parte bassa corrispondenti ai parti prematuri.

Da notare 6 record con peso inferiore al Kg.

Il peso medio del campione osservato è di 3284.1 grammi che è significativamente uguale a quello della popolazione che è di 3300 grammi (ipotesi nulla H0). Quanto detto può essere verificato mediante uno Z-test in quanto è nota sia la media (3300 grammi) sia la deviazione standard (500 grammi) dell’intera popolazione. Considero un valore di significatività dell 1%.

mean_pop = 3300 # Media della popolazione
stdev_pop = 500 # Deviazione standard della popolazione

ztest = z.test(Peso,
        mean_pop,
        stdev=stdev_pop,
        alternative="two.sided",
        conf.level=0.99)

conf_int <- ztest$conf.int

if (mean_pop >= conf_int[1] && mean_pop <= conf_int[2]) {
  cat("Non rifiuto l'ipotesi nulla H0")
} else {
  cat("Rifiuto l'ipotesi nulla H0")
}

## Non rifiuto l'ipotesi nulla H0

Minimo dell’intervallo confidenza	Media della popolazione	Massimo dell’intervallo confidenza
3258.3 grammi	3300 grammi	3309.8 grammi

mean_pop = 3300 # Media della popolazione
stdev_pop = 500 # Deviazione standard della popolazione

peso_popolazione <- data.frame(peso = rnorm(10000, 
                                            mean=mean_pop, 
                                            sd=stdev_pop))

ggplot(peso_popolazione, aes(x = peso)) +
  geom_density(fill = "skyblue", 
               color = "black", 
               alpha = 0.5) + 
  labs(title = "Distribuzione del peso dei neonati",
       x = "Peso (grammi)",
       y = "Densità") +
  theme_minimal() +
  geom_vline(xintercept = conf_int[1], linetype = "dashed", color = "red") +
  geom_vline(xintercept = conf_int[2], linetype = "dashed", color = "red") +
  geom_point(x = mean_pop, y = 0, color = "blue", size = 3)

Si può osservare che il valore medio del peso della popolazione (punto verde) si trova all’interno dell’intervallo di confidenza, non rifiuto quindi l’ipotesi nulla H0.

Lunghezza

Variabile quantitativa: lunghezza del neonato in millimetri

stat(Lunghezza)

Statistica	Valore
Min.	310.000
1st Qu.	480.000
Median	500.000
Mean	494.692
3rd Qu.	510.000
Max.	565.000

graf(Lunghezza,"Millimetri")

Lunghezza media di esattamente 50 cm con la maggior parte dei valori compresi tra 48 e 51 cm.

Distribuzione con tanti outlier nella parte bassa corrispondenti ai parti prematuri.

La lunghezza media del campione osservato è di 494.7 millimetri che non è significativamente uguale a quello della popolazione che è di 500 millimetri (ipotesi nulla H0). Quanto detto può essere verificato mediante uno Z-test in quanto è nota sia la media (500 millimetri) sia la deviazione standard (30 millimetri) dell’intera popolazione. Considero un valore di significatività dell 1%.

mean_pop = 500 # Media della popolazione
stdev_pop = 30 # Deviazione standard della popolazione

ztest = z.test(Lunghezza,
        mean_pop,
        stdev=stdev_pop,
        alternative="two.sided",
        conf.level=0.99)

conf_int <- ztest$conf.int

if (mean_pop >= conf_int[1] && mean_pop <= conf_int[2]) {
  cat("Non rifiuto l'ipotesi nulla H0")
} else {
  cat("Rifiuto l'ipotesi nulla H0")
}

## Rifiuto l'ipotesi nulla H0

Minimo dell’intervallo confidenza	Massimo dell’intervallo confidenza	Media della popolazione
493.1 millimetri	496.2 millimetri	500 millimetri

mean_pop = 500 # Media della popolazione
stdev_pop = 30 # Deviazione standard della popolazione

lunghezza_popolazione <- data.frame(peso = rnorm(10000, 
                                            mean=mean_pop, 
                                            sd=stdev_pop))

ggplot(lunghezza_popolazione, aes(x = peso)) +
  geom_density(fill = "skyblue", 
               color = "black", 
               alpha = 0.5) + 
  labs(title = "Distribuzione della lunghezza dei neonati",
       x = "Lunghezza (millimetri)",
       y = "Densità") +
  theme_minimal() +
  geom_vline(xintercept = conf_int[1], linetype = "dashed", color = "red") +
  geom_vline(xintercept = conf_int[2], linetype = "dashed", color = "red") +
  geom_point(x = mean_pop, y = 0, color = "blue", size = 3)

Cranio

Variabile quantitativa: circonferenza cranica del neonato in millimetri

stat(Cranio)

Statistica	Valore
Min.	235.0000
1st Qu.	330.0000
Median	340.0000
Mean	340.0292
3rd Qu.	350.0000
Max.	390.0000

graf(Cranio,"Millimetri")

La circonferenza media è di 34 cm. Distribuzione con tanti outlier nella parte bassa corrispondenti ai parti prematuri.

Tipo.parto

Variabile qualitativa: indica il tipo di parto. Naturale “Nat” o cesareo “Ces”.

freq(Tipo.parto)

Categoria	Record	Percentuale
Ces	728	29.1
Nat	1772	70.9

ist(Tipo.parto,"Cesareo / Naturale")

Predominanza del parto naturale.

ggplot(data=neonati) +
  geom_bar(aes(x=Tipo.parto, fill=Ospedale),
           position="dodge",
           stat="count",
           col="black") +
  labs(title="",
       x="Cesareo / Naturale",
       y="Frequenze assolute")+
  scale_x_discrete(breaks = unique(Tipo.parto), 
                   labels = unique(Tipo.parto)) +
  theme_light()

I parti cesarei e naturali sono, seppur con qualche piccolo scostamento, equamente distribuiti fra i 3 ospedali.

Ospedale

Variabile qualitativa: indica l’ospedale in cui è avvenuto il parto. Ospedale 1 “osp1”, ospedale 2 “osp2” oppure ospedale 3 “osp3”.

freq(Ospedale)

Categoria	Record	Percentuale
osp1	816	32.6
osp2	849	34.0
osp3	835	33.4

ist(Ospedale,"Ospedale")

I record sono, seppur con qualche piccolo scostamento, equamente distribuiti fra i 3 ospedali esaminati.

Sesso

Variabile qualitativa: sesso del neonato. “F” Femmina e “M” maschio.

freq(Sesso)

Categoria	Record	Percentuale
F	1256	50.2
M	1244	49.8

ist(Sesso,"Sesso")

Record equamente distribuiti fra masci e femmine.

Misure antropometriche fra i due sessi

Nella seguente tabella vengono riportate le medie delle misure antropometriche suddivise per i due sessi.

neonati_mf = neonati %>%
  group_by(Sesso) %>%
  summarise(Peso=round(mean(Peso),1), 
            Lunghezza=round(mean(Lunghezza),1), 
            Cranio=round(mean(Cranio),1))

delta <- neonati_mf %>%
  summarise(
    Peso = round(diff(Peso),1),
    Lunghezza = round(diff(Lunghezza),1),
    Cranio = round(diff(Cranio),1)
  )

neonati_mf <- neonati_mf %>%
  add_row(Sesso = "Delta", 
          Peso = delta$Peso, 
          Lunghezza = delta$Lunghezza, 
          Cranio = delta$Cranio)

datatable(neonati_mf)

Peso:

Dalla letteratura si desume che la differenza di peso fra maschi e femmine alla nascita può essere approssimata ad una distribuzione normale con media 200 grammi e deviazione standard 150 grammi.

Per quanto riguarda il campione esaminato è possibile dire se questa differenza è significativamente rilevante effettuando un T-test.

Ipotesi nulla (H0): Non c’è differenza significativa nel peso alla nascita tra maschi e femmine

Ipotesi alternativa (H1): C’è una differenza significativa nel peso alla nascita tra maschi e femmine

Considero un valore di significatività dell’ 1%.

t_test <- t.test(Peso~Sesso, data=neonati, conf.level=0.99)

p_value <- t_test$p.value

if (p_value < 0.01) {
  cat(paste("Rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
} else {
  cat(paste("Non rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
}

## Rifiuto l'ipotesi nulla H0 (p-value = 8.02974271144157e-33 )

Il P-value è estremamente piccolo rifiuto quindi l’ipotesi nulla H0, c’è quindi una differenza significativa di peso alla nascita tra maschi e femmine.

Lunghezza:

Secondo la letteratura scientifica la distribuzione delle differenze di lunghezza dei neonati maschie e femmine alla nascita può essere approssimata ad una curva normale. Da notare che le differenze hanno valori molto piccoli.

Per quanto riguarda il campione esaminato è possibile dire se questa differenza è significativamente rilevante effettuando un T-test.

Ipotesi nulla (H0): Non c’è differenza significativa nella lunghezza alla nascita tra maschi e femmine

Ipotesi alternativa (H1): C’è una differenza significativa nella lunghezza alla nascita tra maschi e femmine

Considero un valore di significatività dell’ 1%.

t_test <- t.test(Lunghezza~Sesso, data=neonati, conf.level=0.99)

p_value <- t_test$p.value

if (p_value < 0.01) {
  cat(paste("Rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
} else {
  cat(paste("Non rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
}

## Rifiuto l'ipotesi nulla H0 (p-value = 2.23724251765558e-21 )

Il P-value è estremamente piccolo rifiuto quindi l’ipotesi nulla H0, c’è quindi una differenza significativa della lunghezza alla nascita tra maschi e femmine.

Cranio:

Secondo la letteratura scientifica la distribuzione delle differenze di diametro del cranio dei neonati maschie e femmine alla nascita può essere approssimata ad una curva normale. Da notare che le differenze hanno valori molto piccoli e spesso sono considerate trascurabili dal punto di vista clinico.

Per quanto riguarda il campione esaminato è possibile dire se questa differenza è significativamente rilevante effettuando un T-test.

Ipotesi nulla (H0): Non c’è differenza significativa nel diametro del cranio alla nascita tra maschi e femmine

Ipotesi alternativa (H1): C’è una differenza significativa nel diametro del cranio alla nascita tra maschi e femmine

Considero un valore di significatività dell’ 1%.

t_test <- t.test(Cranio~Sesso, data=neonati, conf.level=0.99)

p_value <- t_test$p.value

if (p_value < 0.01) {
  cat(paste("Rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
} else {
  cat(paste("Non rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
}

## Rifiuto l'ipotesi nulla H0 (p-value = 1.71769469990384e-13 )

Il P-value è estremamente piccolo rifiuto quindi l’ipotesi nulla H0, c’è quindi una differenza significativa del diametro del cranio alla nascita tra maschi e femmine.

Correlazioni fra le variabili quantitative

Nella tabella sottostante sono stati calcolati tutti i coefficienti di correlazione lineare fra le variabili quantitative, rappresentando per chiarezza solo i valori minimamente significativi, con valore assoluto >= 0.1

# Variabili quantitative
neonati_quant <- neonati[, c("Anni.madre", 
                             "N.gravidanze", 
                             "Gestazione", 
                             "Lunghezza", 
                             "Cranio",
                             "Peso")]

corr <- cor(neonati_quant)
corr_round <- round(corr, 2)
# Correlazione con blank i valori bassi per una migliore leggibilità
corr_blank <- ifelse(abs(corr_round) < 0.1, "", corr_round)
kable(corr_blank)

	Anni.madre	N.gravidanze	Gestazione	Lunghezza	Cranio	Peso
Anni.madre	1	0.38	-0.14
N.gravidanze	0.38	1	-0.1
Gestazione	-0.14	-0.1	1	0.62	0.46	0.59
Lunghezza			0.62	1	0.6	0.8
Cranio			0.46	0.6	1	0.7
Peso			0.59	0.8	0.7	1

Nei seguanti grafici scatterplot vengono rappresentate tutte le combinazioni fra le variabili quantitative.

var <- names(neonati_quant)
n_var <- length(var)

# Asse x
for (i in 1:(n_var - 1)) {
  # Asse y
  for (j in (i + 1):n_var) {

    grafico <- ggplot(neonati_quant, 
                      aes_string(x=var[i], y=var[j])) +
      geom_hex() +
      scale_fill_gradient(low = "blue4", high = "yellow") +
      labs(title = paste(var[i], " - ", var[j]),
           x = var[i],
           y = var[j]) +
      theme_light()

    print(grafico)
  }
}

Analisi dei risultati grafici:

Anni.madre - N.gravidanze: evidenza la correlazione positiva con coefficiente di correlazione di 0.38
Anni.madre - Gestazione: non si evidenzia alcuna correlazione lineare
Anni.madre - Lunghezza: non si evidenzia alcuna correlazione lineare
Anni.madre - Cranio: non si evidenzia alcuna correlazione lineare
Anni.madre - Peso: non si evidenzia alcuna correlazione lineare
N.gravidanze - Gestazione: all’aumentare del n° di gravidanze i valori minimi delle settimane di gestazione aumentano. Può essere una variabile da tenere in considerazione nel modello di regressione lineare multipla pur avendo un coefficiente di - 0.1
N.gravidanze - Lunghezza: andamento simile a “N.gravidanze - Gestazione”. Le settimane di gestazione sono positivamente correlate alle misure antropometriche
N.gravidanze - Cranio: andamento simile a “N.gravidanze - Gestazione”. Le settimane di gestazione sono positivamente correlate alle misure antropometriche
N.gravidanze - Peso: andamento simile a “N.gravidanze - Gestazione”. Le settimane di gestazione sono positivamente correlate alle misure antropometriche
Gestazione - Lunghezza: evidenza la correlazione positiva con coefficiente di correlazione di 0.62
Gestazione - Cranio: evidenza la correlazione positiva con coefficiciente di correlazione di 0.46
Gestazione - Peso: evidenza la correlazione positiva con coefficiente di correlazione di 0.59
Lunghezza - Cranio: evidenza la correlazione positiva con coefficiente di correlazione di 0.6
Lunghezza - Peso: evidenza la correlazione positiva con coefficiente di correlazione di 0.8
Cranio - Peso: evidenza la correlazione positiva con coefficiente di correlazione di 0.7

Variabile peso in funzione delle variabili qualitative

Mediante grafici boxplot si evidenza la relazione delle variabili qualitative alla variabile peso di cui si vuole fare un modello di regressione lineare.

ggplot(neonati, aes(x=Sesso, y=Peso, fill=Sesso)) +
  geom_boxplot() +
  scale_fill_manual(values = c("F" = "pink", "M" = "skyblue")) +
  labs(title = "Boxplot del peso per sesso",
       x = "Sesso",
       y = "Peso") +
  theme_light()

Il peso del neonato, come già visto precedentemente, è significativamente influenzato dal sesso.

neonati$Fumatrici <- factor(neonati$Fumatrici, 
                            levels = c(0, 1), 
                            labels = c("Non fumatrici", "Fumatrici"))

ggplot(neonati, aes(x=Fumatrici, y=Peso, fill=Fumatrici)) +
  geom_boxplot() +
  scale_fill_manual(values = c("Non fumatrici" = "lightgreen", 
                               "Fumatrici" = "grey")) +
  labs(title = "Boxplot del peso per fumatrici / non fumatrici",
       x = "Non fumatrici / fumatrici",
       y = "Peso") +
  theme_light()

La variabile “Fumatrici” non influenza significativamente il peso del neonato, questo può essere verificato con un test T-test:

t_test <- t.test(Peso~Fumatrici, data=neonati, conf.level=0.95)

p_value <- t_test$p.value

if (p_value < 0.05) {
  cat(paste("Rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
} else {
  cat(paste("Non rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
}

## Non rifiuto l'ipotesi nulla H0 (p-value = 0.303317226271264 )

ggplot(neonati, aes(x=Tipo.parto, y=Peso, fill=Tipo.parto)) +
  geom_boxplot() +
  scale_fill_manual(values = c("Nat" = "lightgreen", "Ces" = "red")) +
  labs(title = "Boxplot del peso per tipo parto",
       x = "Cesareo / naturale",
       y = "Peso") +
  theme_light()

Il tipo di parto non influenza significativamente il peso del neonato, questo può essere verificato con un test T-test:

t_test <- t.test(Peso~Tipo.parto, data=neonati, conf.level=0.95)

p_value <- t_test$p.value

if (p_value < 0.05) {
  cat(paste("Rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
} else {
  cat(paste("Non rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
}

## Non rifiuto l'ipotesi nulla H0 (p-value = 0.896840872866718 )

ggplot(neonati, aes(x=Ospedale, y=Peso, fill=Ospedale)) +
  geom_boxplot() +
  #scale_fill_manual(values = c("Nat" = "lightgreen", "Ces" = "red")) +
  labs(title = "Boxplot del peso per ospedale",
       x = "Ospedale",
       y = "Peso") +
  theme_light()

L’ospedale non influenza significativamente il peso del neonato, questo può essere verificato con un test ANOVA:

modello_anova <- aov(Peso ~ Ospedale, data = neonati)

p_value <- summary(modello_anova)[[1]]$`Pr(>F)`[1]

if (p_value < 0.05) {
  cat(paste("Rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
} else {
  cat(paste("Non rifiuto l'ipotesi nulla H0 (p-value =", p_value,")"))
}

## Non rifiuto l'ipotesi nulla H0 (p-value = 0.183061493975835 )

Modello di regressione multipla

Si vuole creare un modello di regressione multipla per determianare il peso del neonato.

La variabile Peso non ha una distribuzione normale, questo può essere affermato effettuando un Shapiro-Wilk test ed verificando che il p-value è estremamente piccolo:

shapiro.test(Peso)

## 
##  Shapiro-Wilk normality test
## 
## data:  Peso
## W = 0.97066, p-value < 2.2e-16

Nel tentativo di normalizzare la variabile peso elimino dal dataset i due record palesemente anomali:

neonati <- subset(neonati, Anni.madre >= 2)
attach(neonati)

Valuto se eliminare anche i record corrispondenti ad outlier, considerando dei valori di soglia molto ampi (1.6 * IQR) in modo da non distorcere il dataset originale.

ggplot(neonati, aes(sample = Peso)) +
  stat_qq() +
  stat_qq_line() +
  labs(title = "Grafico Q-Q di Peso considerando gli outlier",
       x = "Quantili Teorici (Normale)",
       y = "Quantili Campione (Peso)") +
  theme_minimal()

maschi <- neonati[neonati$Sesso == "M", ]
femmine <- neonati[neonati$Sesso == "F", ]

Q1_maschi <- quantile(maschi$Peso, 0.25)
Q3_maschi <- quantile(maschi$Peso, 0.75)
IQR_maschi <- Q3_maschi - Q1_maschi
limite_inferiore_maschi <- Q1_maschi - 1.6 * IQR_maschi
limite_superiore_maschi <- Q3_maschi + 1.6 * IQR_maschi
outliers_maschi <- maschi[maschi$Peso < limite_inferiore_maschi | 
                          maschi$Peso > limite_superiore_maschi, ]

Q1_femmine <- quantile(femmine$Peso, 0.25)
Q3_femmine <- quantile(femmine$Peso, 0.75)
IQR_femmine <- Q3_femmine - Q1_femmine
limite_inferiore_femmine <- Q1_femmine - 1.6 * IQR_femmine
limite_superiore_femmine <- Q3_femmine + 1.6 * IQR_femmine
outliers_femmine <- femmine[femmine$Peso < limite_inferiore_femmine | 
                            femmine$Peso > limite_superiore_femmine, ]

outliers <- c(rownames(outliers_maschi), rownames(outliers_femmine))
neonati_mod <- neonati[!(rownames(neonati) %in% outliers), ]
row_del = 2498 - nrow(neonati_mod)
row_del

## [1] 76

ggplot(neonati_mod, aes(x=Sesso, y=Peso, fill=Sesso)) +
  geom_boxplot() +
  scale_fill_manual(values = c("F" = "pink", "M" = "skyblue")) +
  labs(title = "Boxplot del peso per sesso senza outlier",
       x = "Sesso",
       y = "Peso") +
  theme_light()

ggplot(neonati_mod, aes(sample = Peso)) +
  stat_qq() +
  stat_qq_line() +
  labs(title = "Grafico Q-Q di Peso senza outlier",
       x = "Quantili Teorici (Normale)",
       y = "Quantili Campione (Peso)") +
  theme_minimal()

Rieffettuo quindi uno Shapiro-Wilk test:

shapiro.test(neonati_mod$Peso)

## 
##  Shapiro-Wilk normality test
## 
## data:  neonati_mod$Peso
## W = 0.99869, p-value = 0.05793

Ottenedo un p-value di 0.05793 posso considerare che Peso abbia una distribuzione normale.

Considerando però che i record che andrei ad eliminare sono ben 76 senza avere alcuna regionevole certezza che siano dati errati, preferisco mantenere il dataset originale senza effettuare distorsioni forzate.

Nei capitoli precedenti sono già state analizzate le correlazioni tra le variabili mediante T-test e scatterplot. Per realizzare un modello di regressione multipla terrei in considerazione questi punti:

La variabile Peso ha una forte correlazione con Lunghezza, Cranio e Gestazione.
Le variabili Lunghezza, Cranio e Gestazione sono fortemente correlate fra di loro e che quindi potrebbero generare multicollinearità, nel caso di utilizzo simultaneo in un modello di regressione.
Le variabili Anni.madre e N.gravidanze non sembrano correlate con la variabile Peso. Da notare però quanto osservato nell’analisi degli scatterplot: N.gravidanze sembra influenzare positivamente i valori minimi delle settimane di gestazione.
Per quanto riguarda le variabili qualitative solo Sesso influenza significativamente il peso.

Modello 1:

Nel primo modello di regressione vengono considerate tutte le variabili esplicative possibili:

mod_1 = lm(Peso ~ Gestazione + 
                  Lunghezza + 
                  Cranio +
                  N.gravidanze +
                  Anni.madre +
                  Sesso +
                  Fumatrici + 
                  Tipo.parto + 
                  Ospedale)

tab_coef <- summary(mod_1)$coefficients
p_values <- tab_coef[, "Pr(>|t|)"]

# Funzione per assegnare le % di significatività
assign_significance <- function(p) {
  if (p < 0.001) {
    return("< 0.1 %")
  } else if (p < 0.01) {
    return("< 1 %")
  } else if (p < 0.05) {
    return("< 5 %")
  } else if (p < 0.1) {
    return("< 10 %")
  } else {
    return(" ")
  }
}

signif_codes <- sapply(p_values, assign_significance)
tab_coef <- cbind(p_values, "Signif." = signif_codes)
kable(tab_coef)

	p_values	Signif.
(Intercept)	0	< 0.1 %
Gestazione	2.57878810763896e-17	< 0.1 %
Lunghezza	1.62284443764329e-210	< 0.1 %
Cranio	1.67101441002222e-119	< 0.1 %
N.gravidanze	0.0148457496951216	< 5 %
Anni.madre	0.484486122210521
SessoM	5.17920945327777e-12	< 0.1 %
FumatriciFumatrici	0.271912743666334
Tipo.partoNat	0.014315417986905	< 5 %
Ospedaleosp2	0.40956171544463
Ospedaleosp3	0.0365652402589988	< 5 %

r_squared_adjusted <- round(summary(mod_1)$adj.r.squared,3)

Il valore di R^2 adjusted è: 0.728

Si conferma quanti ipotizzato nei capitoli precedenti:

La forte significatività delle variabili Gestazione, Lunghezza, Cranio, Sesso
La bassa significatività di N.gravidanze, Tipo.parto, Ospedale
La non significatività di Anni.madre e Fumatrici

Modello 2:

Elimino dal modello le variabili non significative Anni.madre e Fumatrici:

mod_2 = lm(Peso ~ Gestazione + 
                  Lunghezza + 
                  Cranio +
                  N.gravidanze +
                  Sesso +
                  Tipo.parto + 
                  Ospedale)

tab_coef <- summary(mod_2)$coefficients
p_values <- tab_coef[, "Pr(>|t|)"]

signif_codes <- sapply(p_values, assign_significance)
tab_coef <- cbind(p_values, "Signif." = signif_codes)
kable(tab_coef)

	p_values	Signif.
(Intercept)	0	< 0.1 %
Gestazione	4.96777633897981e-17	< 0.1 %
Lunghezza	2.39356143820422e-211	< 0.1 %
Cranio	3.26642227161212e-120	< 0.1 %
N.gravidanze	0.0044628604375385	< 1 %
SessoM	5.47721175660367e-12	< 0.1 %
Tipo.partoNat	0.0150471318815552	< 5 %
Ospedaleosp2	0.417858322767692
Ospedaleosp3	0.0330059543573244	< 5 %

r_squared_adjusted <- round(summary(mod_2)$adj.r.squared,3)

Il valore di R^2 adjusted è: 0.728 identico a quello del modello 1.

Modello 3:

Provo a semplificare ulteriormente il modello eliminado le variabili meno significative Ospedale e Tipo.parto

mod_3 = lm(Peso ~ Gestazione + 
                  Lunghezza + 
                  Cranio +
                  N.gravidanze +
                  Sesso)

tab_coef <- summary(mod_3)$coefficients
p_values <- tab_coef[, "Pr(>|t|)"]

signif_codes <- sapply(p_values, assign_significance)
tab_coef <- cbind(p_values, "Signif." = signif_codes)
kable(tab_coef)

	p_values	Signif.
(Intercept)	0	< 0.1 %
Gestazione	2.71939279548013e-17	< 0.1 %
Lunghezza	4.32202191640788e-209	< 0.1 %
Cranio	8.11929059051433e-121	< 0.1 %
N.gravidanze	0.00415407274428675	< 1 %
SessoM	4.46626455091865e-12	< 0.1 %

r_squared_adjusted <- round(summary(mod_3)$adj.r.squared,3)

Il valore di R^2 adjusted è: 0.726, quindi sempre pressochè invariato.

Modello 4:

Provo a semplificare ulteriormente il modello eliminado la variabile Gestazione, fortemente correlata a Lunghezza e Cranio

mod_4 = lm(Peso ~ Lunghezza +
                  Cranio + 
                  N.gravidanze +
                  Sesso )

tab_coef <- summary(mod_4)$coefficients
p_values <- tab_coef[, "Pr(>|t|)"]

signif_codes <- sapply(p_values, assign_significance)
tab_coef <- cbind(p_values, "Signif." = signif_codes)
kable(tab_coef)

	p_values	Signif.
(Intercept)	0	< 0.1 %
Lunghezza	1.53594098342368e-297	< 0.1 %
Cranio	4.46503182660008e-131	< 0.1 %
N.gravidanze	0.0430654802936948	< 5 %
SessoM	2.81104866732156e-12	< 0.1 %

r_squared_adjusted <- round(summary(mod_4)$adj.r.squared,3)

Il valore di R^2 adjusted è: 0.719, si è ridotto anche se non di molto.

Modello stepwise:

Viene effettuata una selezione delle variabili utilizzando una procedura automatica.

n = nrow(neonati)

mod.stepwise <- stepAIC(mod_1,
                        direction="both",
                        k=log(n), # criterio BIC 
                        trace = FALSE)

summary(mod.stepwise)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze + 
##     Sesso)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1149.37  -180.98   -15.57   163.69  2639.09 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -6681.7251   135.8036 -49.201  < 2e-16 ***
## Gestazione      32.3827     3.8008   8.520  < 2e-16 ***
## Lunghezza       10.2455     0.3008  34.059  < 2e-16 ***
## Cranio          10.5410     0.4265  24.717  < 2e-16 ***
## N.gravidanze    12.4554     4.3416   2.869  0.00415 ** 
## SessoM          77.9807    11.2111   6.956 4.47e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 274.7 on 2492 degrees of freedom
## Multiple R-squared:  0.727,  Adjusted R-squared:  0.7265 
## F-statistic:  1327 on 5 and 2492 DF,  p-value: < 2.2e-16

Il risultato conferma che il modello con il miglior compromesso tra BIC e semplicità è, per il momento, il modello 3 che contiene le seguenti variabili:

Gestazione
Lunghezza
Cranio
N.gravidanze
Sesso

Modello con interazione fra variabili:

Partendo dal modello 3 si testa un’interazione fra le variabili Gestazione, Lunghezza e Cranio:

mod_3_interazione = lm(Peso ~ Gestazione * Lunghezza * Cranio +
                       N.gravidanze +
                       Sesso)

tab_coef <- summary(mod_3_interazione)$coefficients
p_values <- tab_coef[, "Pr(>|t|)"]

signif_codes <- sapply(p_values, assign_significance)
tab_coef <- cbind(p_values, "Signif." = signif_codes)
kable(tab_coef)

	p_values	Signif.
(Intercept)	1.96909889765174e-05	< 0.1 %
Gestazione	8.04261551521498e-07	< 0.1 %
Lunghezza	0.00028688140023518	< 0.1 %
Cranio	9.34151597640356e-06	< 0.1 %
N.gravidanze	0.00362190970786506	< 1 %
SessoM	8.88732617858961e-11	< 0.1 %
Gestazione:Lunghezza	1.24614577405449e-05	< 0.1 %
Gestazione:Cranio	8.5431458075734e-07	< 0.1 %
Lunghezza:Cranio	4.36500839945457e-05	< 0.1 %
Gestazione:Lunghezza:Cranio	1.41459017728832e-05	< 0.1 %

r_squared_adjusted <- round(summary(mod_3_interazione)$adj.r.squared,3)

Il valore di R^2 adjusted è: 0.732, si è quindi migliorato leggermente il modello 3

Modello con relazioni non lineari:

Partendo dall’osservazione dei grafici scatterplot si prova a modificare il modello 3 con interazione, poco fa calcolato, introducendo una relazione quadratica per la variabili Gestazione, Cranio e Lunghezza.

mod_3_nl = lm(Peso ~ I(Gestazione)^2 +
                     I(Cranio)^2 + 
                     I(Lunghezza)^2 + 
                     N.gravidanze +
                     Sesso)

tab_coef <- summary(mod_3_nl)$coefficients
p_values <- tab_coef[, "Pr(>|t|)"]

signif_codes <- sapply(p_values, assign_significance)
tab_coef <- cbind(p_values, "Signif." = signif_codes)
kable(tab_coef)

	p_values	Signif.
(Intercept)	0	< 0.1 %
I(Gestazione)	2.71939279548206e-17	< 0.1 %
I(Cranio)	8.11929059050599e-121	< 0.1 %
I(Lunghezza)	4.32202191641304e-209	< 0.1 %
N.gravidanze	0.00415407274428663	< 1 %
SessoM	4.46626455091833e-12	< 0.1 %

r_squared_adjusted <- round(summary(mod_3_nl)$adj.r.squared,3)

Il valore di R^2 adjusted è: 0.726, non si ottiene quindi alcun miglioramento.

Raffronto fra modelli:

Per tutti i modelli viene fatto un raffronto utilizzando le seguenti metriche:

AIC : Akaike Information Criterion
BIC : Bayesian Information Criterion
R^2 Adj. : R quadro aggiustato
RMSE : Root Mean Square Error
N° var. : numero di variabili utilizzate
N° inter. var. : numero di variabili di interazione

extract_metrics <- function(model) {
  aic <- AIC(model)
  bic <- BIC(model)
  r_squared_adj <- round(summary(model)$adj.r.squared,3)
  rmse <- round(sqrt(mean(residuals(model)^2)))
  num_indep_vars <- length(coef(model)) - 1

  term_names <- names(coef(model))
  interaction_terms <- grep(":", term_names, value = TRUE)
  num_correlated_vars <- length(interaction_terms)

  return(c(aic = aic, 
           bic = bic, 
           r_squared_adj = r_squared_adj, 
           rmse = rmse,
           num_indep_vars = num_indep_vars,
           num_correlated_vars = num_correlated_vars))
}

models <- list(mod_1,
               mod_2, 
               mod_3, 
               mod_4, 
               mod.stepwise, 
               mod_3_interazione, 
               mod_3_nl)

model_descriptions <- c(
  "Modello 1",
  "Modello 2",
  "Modello 3",
  "Modello 4",
  "Modello stepwise",
  "Modello 3 con interazione",
  "Modello 3 non lineare"
)

metrics_list <- lapply(models, extract_metrics)

metrics_table <- do.call(rbind, metrics_list)
rownames(metrics_table) <-model_descriptions

kable(metrics_table, 
      booktabs = T,
      col.names = c("AIC", 
                    "BIC", 
                    "R^2 Adj.", 
                    "RMSE",
                    "N° var.",
                    "N° inter. var.")) %>%
      kable_styling(bootstrap_options = "striped", full_width = T)

	AIC	BIC	R^2 Adj.	RMSE	N° var.	N° inter. var.
Modello 1	35145.57	35215.45	0.728	273	10	0
Modello 2	35143.29	35201.52	0.728	274	8	0
Modello 3	35152.89	35193.65	0.726	274	5	0
Modello 4	35222.61	35257.55	0.719	278	4	0
Modello stepwise	35152.89	35193.65	0.726	274	5	0
Modello 3 con interazione	35106.59	35170.65	0.732	271	9	4
Modello 3 non lineare	35152.89	35193.65	0.726	274	5	0

Il miglior modello è il 3 con interazione

per la sua semplicità: utilizza solo 5 delle 10 variabili del dataset
ottiene AIC, BIC, RMSE più bassi rispetto a tutti gli altri modelli
R^2 adjusted più elevata rispetto a tutti gli altri modelli

Si verifica che il modello 3 con interazione non presenti problemi di multicollinearità.

vif_mod_3_interazione = vif(mod_3_interazione, type = "predictor")

vif_table <- data.frame(Variable = rownames(vif_mod_3_interazione),
                        VIF = vif_mod_3_interazione[, 1])

kable(vif_table) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), 
                full_width = FALSE)

Variable	VIF
Gestazione	1.075440
Lunghezza	1.075440
Cranio	1.075440
N.gravidanze	1.025337
Sesso	1.051026

Tutti i valori sono decisamente minori di 5.

Analisi dei residui:

Residui rispetto ai valori predetti:

plot(mod_3_interazione, which = 1)

I punti si ditribiscono attorno alla media 0 non in modo casuale, ma assumendo una forma tondeggiante ed allungata.

Si possono osservare 4 punti outlier:

outlierTest(mod_3_interazione)

##       rstudent unadjusted p-value Bonferroni p
## 1549 10.574247         1.3587e-25   3.3941e-22
## 155   5.393128         7.5771e-08   1.8928e-04
## 1305  4.790619         1.7605e-06   4.3976e-03
## 1692  4.288807         1.8649e-05   4.6586e-02

outlier_data <- neonati[c(155,1306,1551,1694), ] 
kable(outlier_data) %>%       
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),                      full_width = T)

	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
155	30	0	Non fumatrici	36	3610	410	330	Nat	osp1	M
1307	29	0	Non fumatrici	42	3560	510	355	Nat	osp1	F
1553	30	4	Non fumatrici	35	4520	520	360	Nat	osp2	F
1696	26	1	Non fumatrici	38	2950	500	340	Nat	osp3	F

Test Durbin-Watson per testare la non correlazione:

dwtest(mod_3_interazione)

## 
##  Durbin-Watson test
## 
## data:  mod_3_interazione
## DW = 1.9599, p-value = 0.158
## alternative hypothesis: true autocorrelation is greater than 0

Non rifiuto l’ipotesi nulla quindi i residui non sono autocorrelati.

Normalità dei residui:

plot(mod_3_interazione, which = 2)

I punti siano posizionati sulla retta che rappresenta la distribuzione normale solo nella parte centrale, mentre agli estremi ci sono degli scostamenti.

shapiro.test(residuals(mod_3_interazione))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(mod_3_interazione)
## W = 0.97407, p-value < 2.2e-16

Test di Shapiro-Wilk: si rifiuta l’ipotesi nulla di normalità.

Omoschedasticità (varianza costante dei residui):

plot(mod_3_interazione, which = 3)

Dall’osservazione del grafico si può osservare una varianza abbastanza costante dei residui.

bptest(mod_3_interazione)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod_3_interazione
## BP = 222.84, df = 9, p-value < 2.2e-16

Test di Breusch-Pagan: rifiuto dell’ipotesi nulla. Quindi il test numerico non conferma l’ipotesi di varianza costante fatta dall’analisi grafica.

Punti leverage:

plot(mod_3_interazione, which = 4)

Solo il campione 1551 supera la soglia di Cook di 1.

Previsione dei risultati:

La previsione dei risultati può essere effettuata con la funzione predict.

Per esempio si vuole stimare il peso di una neonata utilizzando i seguanti valori:

Gestazione madre: 39 settimane
Lunghezza neonato : 499 millimetri
Cranio neonato: 335 millimetri
N.gravidanze madre: 3
Sesso neonato: F

var_previsione <- data.frame(Gestazione = 39,
                         Lunghezza = 499,
                         Cranio = 335,
                         N.gravidanze = 3,
                         Sesso = "F")

peso_previsione <- round(predict(mod_3_interazione, newdata = var_previsione))

Il peso previsto dal modello è 3265 grammi.

Rappresentazioni grafiche:

Rappresento la variabile Peso in funzione delle variabili Gestazione, Lunghezza, Cranio, N.gravidanze e Sesso utilizzate nel modello selezionato come migliore. Mediante una geometria smooth viene aggiunta una rappresentazione, suddivisa per sesso, del modello di regressione lineare.

ggplot(data=neonati)+
  geom_point(aes(x=Gestazione,
                 y=Peso,
                 col=Sesso), position = "jitter", size=1) +
  geom_smooth(aes(x=Gestazione,
                  y=Peso,
                  col=Sesso), se=F, method="lm")+
  labs(title="Peso in funzione dei mesi di gestazione",
       x="Mesi gestazione",
       y="Peso") +
  theme_light()

ggplot(data=neonati)+
  geom_point(aes(x=Lunghezza,
                 y=Peso,
                 col=Sesso), position = "jitter", size=1)+
  geom_smooth(aes(x=Lunghezza,
                  y=Peso,
                  col=Sesso), se=F, method="lm") +
  labs(title="Peso in funzione della lunghezza",
       x="Lunghezza",
       y="Peso") +
  theme_light()

ggplot(data=neonati)+
  geom_point(aes(x=Cranio,
                 y=Peso,
                 col=Sesso), position = "jitter", size=1) +
  geom_smooth(aes(x=Cranio,
                  y=Peso,
                  col=Sesso), se=F, method="lm")+
  labs(title="Peso in funzione del diametro del cranio",
       x="Cranio",
       y="Peso") +
  theme_light()

ggplot(data=neonati)+
  geom_point(aes(x=N.gravidanze,
                 y=Peso,
                 col=Sesso), position = "jitter", size=1) +
  geom_smooth(aes(x=N.gravidanze,
                  y=Peso,
                  col=Sesso), se=F, method="lm")+
  labs(title="Peso in funzione del numero di gravidanze",
       x="N° di gravidanze",
       y="Peso") +
  theme_light()

Modello Statistico per la Previsione del Peso Neonatale

Guido Miglio

2025-02-09

Analisi descrittiva delle variabili

Anni.madre

N.gravidanze

Fumatrici

Gestazione

Peso

Lunghezza

Cranio

Tipo.parto

Ospedale

Sesso

Misure antropometriche fra i due sessi

Peso:

Lunghezza:

Cranio:

Correlazioni fra le variabili quantitative

Variabile peso in funzione delle variabili qualitative

Modello di regressione multipla

Modello 1:

Modello 2:

Modello 3:

Modello 4:

Modello stepwise:

Modello con interazione fra variabili:

Modello con relazioni non lineari:

Raffronto fra modelli:

Analisi dei residui:

Residui rispetto ai valori predetti:

Normalità dei residui:

Omoschedasticità (varianza costante dei residui):

Punti leverage:

Previsione dei risultati:

Rappresentazioni grafiche: