Modello Statistico per la Previsione del Peso Neonatale v2

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.

record_anomali <- neonati[Anni.madre %in% c(0, 1), ]
kable(record_anomali)

	Anni.madre	N.gravidanze	Fumatrici	Gestazione	Peso	Lunghezza	Cranio	Tipo.parto	Ospedale	Sesso
1152	1	1	0	41	3250	490	350	Nat	osp2	F
1380	0	0	0	39	3060	490	330	Nat	osp3	M

Elimino dal dataset i 2 valori anomali:

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

N.gravidanze

Variabile quantitativa: indica quante gravidanze ha avuto la madre

stat(N.gravidanze)

Statistica	Valore
Min.	0.0000000
1st Qu.	0.0000000
Median	1.0000000
Mean	0.9815853
3rd Qu.	1.0000000
Max.	12.0000000

# 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, 2)
kable(df_tab)
}

freq(Fumatrici)

Categoria	Record	Percentuale
0	2394	95.84
1	104	4.16

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.00000
1st Qu.	38.00000
Median	39.00000
Mean	38.97958
3rd Qu.	40.00000
Max.	43.00000

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.184
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%.

Nota: tutti i valori di riferimento della popolazione sono stati estrapolati della seguente pagina OMS: https://www.who.int/tools/child-growth-standards

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

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

Osservando l’intervallo di confidenza non rifiuto l’ipotesi nulla H0.

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.0000
1st Qu.	480.0000
Median	500.0000
Mean	494.6958
3rd Qu.	510.0000
Max.	565.0000

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

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

Osservando l’intervallo di confidenza rifiuto l’ipotesi nulla H0.

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.14
Nat	1770	70.86

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. L’ospedale non sembra quindi influenzare in modo significativo la tipologia di parto, questo può essere verificato con un test Chi-quadro:

tab_contingenza <- table(Tipo.parto, Ospedale)
kable(tab_contingenza)

	osp1	osp2	osp3
Ces	242	254	232
Nat	574	594	602

test_chi_quadro <- chisq.test(tab_contingenza)
print(test_chi_quadro)

## 
##  Pearson's Chi-squared test
## 
## data:  tab_contingenza
## X-squared = 1.083, df = 2, p-value = 0.5819

Il valore del p-value è alto (0.5819), non c’è quindi alcuna relazione significativa fra le variabili Tipo.parto e Ospedale.

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.67
osp2	848	33.95
osp3	834	33.39

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	1255	50.24
M	1243	49.76

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),2), 
            Lunghezza=round(mean(Lunghezza),2), 
            Cranio=round(mean(Cranio),2))

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

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)
t_test$p.value

## [1] 7.275684e-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)
t_test$p.value

## [1] 2.232328e-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)
t_test$p.value

## [1] 1.414019e-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, 3)
# 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.383	-0.135
N.gravidanze	0.383	1	-0.102
Gestazione	-0.135	-0.102	1	0.619	0.461	0.592
Lunghezza			0.619	1	0.603	0.796
Cranio			0.461	0.603	1	0.705
Peso			0.592	0.796	0.705	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)
t_test$p.value

## [1] 0.3022785

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

ggplot(neonati, aes(x=Ospedale, y=Peso, fill=Ospedale)) +
  geom_boxplot() +
  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]

Ha senso, secondo me, effettuare un test ANOVA perchè ci potrebbero essere delle cause che potrebbero portare ad un’influenza fra ospedale e peso del neonato, per esempio:

Qualità dell’assistenza prenatale: la qualità del monitoraggio e delle cure prestate alla madre durante la gravidanza può variare da ospedale a ospedale. Un’assistenza prenatale più accurata e personalizzata potrebbe contribuire a un migliore sviluppo del feto e a un peso alla nascita più adeguato.
Fattori socio-economici: gli ospedali situati in aree con diverse condizioni socio-economiche potrebbero assistere a una diversa distribuzione del peso dei neonati.

Modello di regressione multipla

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

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)

# 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(" ")
  }
}

# Funzione per estrarre tutti i coefficienti
tabella_coefficienti <- function(modello) {
  tab_coef <- summary(modello)$coefficients
  p_values <- tab_coef[, "Pr(>|t|)"]
  estimate <- round(tab_coef[, "Estimate"],5)
  std_error <- round(tab_coef[, "Std. Error"],5)
  t_value <- round(tab_coef[, "t value"],5)

  signif_codes <- sapply(p_values, assign_significance)
  tab_coef <- cbind(estimate, 
                    std_error, 
                    t_value, 
                    p_values, 
                    "Signif." = signif_codes)
  kable(tab_coef, align = "r")
}

tabella_coefficienti(mod_1)

	estimate	std_error	t_value	p_values	Signif.
(Intercept)	-6735.79597	141.47895	-47.60988	0	< 0.1 %
Gestazione	32.57729	3.82075	8.52641	2.57878810763896e-17	< 0.1 %
Lunghezza	10.29218	0.30088	34.20703	1.62284443764329e-210	< 0.1 %
Cranio	10.47221	0.42627	24.56707	1.67101441002222e-119	< 0.1 %
N.gravidanze	11.38117	4.66858	2.43782	0.0148457496951216	< 5 %
Anni.madre	0.80179	1.14671	0.69921	0.484486122210521
SessoM	77.57227	11.18652	6.93444	5.17920945327777e-12	< 0.1 %
Fumatrici	-30.27413	27.54918	-1.09891	0.271912743666334
Tipo.partoNat	29.63351	12.0905	2.45097	0.014315417986905	< 5 %
Ospedaleosp2	-11.09121	13.44708	-0.8248	0.40956171544463
Ospedaleosp3	28.24955	13.50545	2.09171	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)

tabella_coefficienti(mod_2)

	estimate	std_error	t_value	p_values	Signif.
(Intercept)	-6707.92522	136.02571	-49.31366	0	< 0.1 %
Gestazione	32.03856	3.79247	8.44793	4.96777633897981e-17	< 0.1 %
Lunghezza	10.30586	0.30058	34.28636	2.39356143820422e-211	< 0.1 %
Cranio	10.492	0.42567	24.64829	3.26642227161212e-120	< 0.1 %
N.gravidanze	12.33599	4.33444	2.84604	0.0044628604375385	< 1 %
SessoM	77.46565	11.1842	6.92635	5.47721175660367e-12	< 0.1 %
Tipo.partoNat	29.40797	12.08746	2.43293	0.0150471318815552	< 5 %
Ospedaleosp2	-10.89393	13.44469	-0.81028	0.417858322767692
Ospedaleosp3	28.79168	13.49695	2.1332	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)

tabella_coefficienti(mod_3)

	estimate	std_error	t_value	p_values	Signif.
(Intercept)	-6681.72512	135.80361	-49.20138	0	< 0.1 %
Gestazione	32.38273	3.80081	8.51996	2.71939279548013e-17	< 0.1 %
Lunghezza	10.24546	0.30082	34.05898	4.32202191640788e-209	< 0.1 %
Cranio	10.54095	0.42647	24.71668	8.11929059051433e-121	< 0.1 %
N.gravidanze	12.45544	4.34158	2.86887	0.00415407274428675	< 1 %
SessoM	77.98074	11.21108	6.95569	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.

Il modello però risulta decisamente semplificato rispetto a quello iniziale in quanto sono state eliminate ben 4 variabili non significative. Osservando la colonna estimate si possono osservare gli effetti marginali di ciascuna variabile tenendo fisse le altre (coefficienti beta). Per esempio per ogni incremento unitario delle settimane di gestazione un bambino peserà 32.38 grammi in più. Oppure un bambino peserà 78 grammi in più di una bambina. Molto interessante è il contributo positivo del numero di gravidanze: circa 12.5 grammi.

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 dalle 5 varibili significative del modello 3 si testa una possibile interazione fra esse.

Vengono considerate tutte le possibili combinazioni esaminandole con un modello automatico stepwise:

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

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

summary(mod.stepwise.inter)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + N.gravidanze + 
##     Sesso + Gestazione:Lunghezza + Gestazione:Cranio + Lunghezza:Cranio + 
##     Cranio:N.gravidanze + Gestazione:Lunghezza:Cranio)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1158.14  -180.69   -13.62   163.34  2596.03 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  3.817e+04  8.340e+03   4.577 4.95e-06 ***
## Gestazione                  -1.238e+03  2.365e+02  -5.234 1.80e-07 ***
## Lunghezza                   -7.291e+01  1.954e+01  -3.731 0.000195 ***
## Cranio                      -1.318e+02  2.754e+01  -4.787 1.79e-06 ***
## N.gravidanze                -2.261e+02  7.829e+01  -2.889 0.003903 ** 
## SessoM                       7.235e+01  1.114e+01   6.496 9.95e-11 ***
## Gestazione:Lunghezza         2.376e+00  5.268e-01   4.510 6.78e-06 ***
## Gestazione:Cranio            3.992e+00  7.575e-01   5.270 1.48e-07 ***
## Lunghezza:Cranio             2.640e-01  6.193e-02   4.263 2.09e-05 ***
## Cranio:N.gravidanze          7.019e-01  2.299e-01   3.053 0.002288 ** 
## Gestazione:Lunghezza:Cranio -7.463e-03  1.642e-03  -4.546 5.73e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 271.5 on 2487 degrees of freedom
## Multiple R-squared:  0.7339, Adjusted R-squared:  0.7328 
## F-statistic: 685.8 on 10 and 2487 DF,  p-value: < 2.2e-16

tabella_coefficienti(mod.stepwise.inter)

	estimate	std_error	t_value	p_values	Signif.
(Intercept)	38172.99717	8340.06762	4.57706	4.94737525878854e-06	< 0.1 %
Gestazione	-1237.83254	236.49096	-5.23416	1.79590279197479e-07	< 0.1 %
Lunghezza	-72.90948	19.53954	-3.73138	0.000194695697552472	< 0.1 %
Cranio	-131.81689	27.53681	-4.78693	1.79276629745083e-06	< 0.1 %
N.gravidanze	-226.1458	78.2897	-2.88858	0.00390335409090194	< 1 %
SessoM	72.35114	11.13836	6.49567	9.94911165916199e-11	< 0.1 %
Gestazione:Lunghezza	2.37591	0.52679	4.5102	6.77716311068207e-06	< 0.1 %
Gestazione:Cranio	3.99211	0.7575	5.27012	1.48047852902535e-07	< 0.1 %
Lunghezza:Cranio	0.26402	0.06194	4.26286	2.09354109167527e-05	< 0.1 %
Cranio:N.gravidanze	0.70192	0.2299	3.0532	0.00228815362549761	< 1 %
Gestazione:Lunghezza:Cranio	-0.00746	0.00164	-4.54605	5.7278903932072e-06	< 0.1 %

r_squared_adjusted <- round(summary(mod.stepwise.inter)$adj.r.squared,3)

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

Modello con relazioni non lineari:

Partendo dall’osservazione dei grafici scatterplot si ipotizzano delle relazioni quadratiche fra la variabile Peso con le variabili Gestazione, Cranio e Lunghezza

Per la selezione del modello ottimale si utilizza ancora la funzione stepAIC

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

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

summary(mod.stepwise.nl)

## 
## Call:
## lm(formula = Peso ~ Gestazione + Cranio + Lunghezza + 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 ***
## Cranio          10.5410     0.4265  24.717  < 2e-16 ***
## Lunghezza       10.2455     0.3008  34.059  < 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

tabella_coefficienti(mod.stepwise.nl)

	estimate	std_error	t_value	p_values	Signif.
(Intercept)	-6681.72512	135.80361	-49.20138	0	< 0.1 %
Gestazione	32.38273	3.80081	8.51996	2.71939279548206e-17	< 0.1 %
Cranio	10.54095	0.42647	24.71668	8.11929059050599e-121	< 0.1 %
Lunghezza	10.24546	0.30082	34.05898	4.32202191641304e-209	< 0.1 %
N.gravidanze	12.45544	4.34158	2.86887	0.00415407274428663	< 1 %
SessoM	77.98074	11.21108	6.95569	4.46626455091833e-12	< 0.1 %

r_squared_adjusted <- round(summary(mod.stepwise.nl)$adj.r.squared,3)

Si può notare come l’ipotesi di relazioni quadratiche non venga confermata e la funzione stepAIC individui nuovamente il modello 3 come ottimale.

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)),3)
  
  term_names <- names(coef(model))
  interaction_terms <- grep(":", term_names, value = TRUE)
  categorical_terms <- grep("Ospedale", term_names, value = TRUE)
  
  correction = length(categorical_terms) / 2
  num_indep_vars <- length(coef(model)) - correction - 1
  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.stepwise, 
               mod.stepwise.inter, 
               mod.stepwise.nl)

model_descriptions <- c(
  "Modello 1",
  "Modello 2",
  "Modello 3",
  "Modello stepwise",
  "Modello stepwise con interazione",
  "Modello stepwise 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.417	9	0
Modello 2	35143.29	35201.52	0.728	273.511	7	0
Modello 3	35152.89	35193.65	0.726	274.367	5	0
Modello stepwise	35152.89	35193.65	0.726	274.367	5	0
Modello stepwise con interazione	35099.25	35169.13	0.733	270.894	10	5
Modello stepwise non lineare	35152.89	35193.65	0.726	274.367	5	0

Il miglior modello sembrerebbe essere il modello stepwise con interazione

per la sua semplicità: utilizza solo 5 delle 10 variabili del dataset

(più altre 5 variabili di interazione)
ottiene AIC, BIC, RMSE più bassi rispetto a tutti gli altri modelli
R^2 adjusted più elevata rispetto a tutti gli altri modelli

Non può però essere scelto come migliore perchè presenta problemi di multicollinearità (valori VIF superiori a 5).

vif_mod.stepwise.inter = vif(mod.stepwise.inter, type = "predictor")

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

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

Variable	VIF
Gestazione	1.734655
Lunghezza	1.734655
Cranio	1.051114
N.gravidanze	7509.377144
Sesso	1.051114

Cosa che non succede con il modello 3 dove tutti i valori sono decisamente minori di 5:

vif(mod_3, type = "predictor")

## [1] 1.669779 2.075747 1.624568 1.023462 1.040184

Si sceglie quindi il modello 3 come migliore, perde 0.007 in R^2 adjusted, ma è più robusto e semplice.

Analisi dei residui:

Residui rispetto ai valori predetti:

plot(mod_3, which = 1)

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

Si possono osservare 3 punti outlier:

outlierTest(mod_3)

##       rstudent unadjusted p-value Bonferroni p
## 1549 10.046230         2.6345e-23   6.5810e-20
## 155   5.025345         5.3818e-07   1.3444e-03
## 1305  4.824963         1.4848e-06   3.7092e-03

outlier_data <- neonati[c(155,1305,1549), ] 
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
1306	23	0	Non fumatrici	41	4900	510	352	Nat	osp2	F
1551	35	1	Non fumatrici	38	4370	315	374	Nat	osp3	F

Test Durbin-Watson per testare la non correlazione:

dwtest(mod_3)

## 
##  Durbin-Watson test
## 
## data:  mod_3
## DW = 1.9532, p-value = 0.1209
## 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, 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))

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

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

Da un’analisi grafica si può però dire che la distribuzione dei residui non si discosta molto da una normale con media 0 e deviazione standard 255 (linea rossa).

residui <- residuals(mod_3)
df_residui <- data.frame(residui = residui)
dati_normali <- data.frame(normali = rnorm(100000, 0, 260))

ggplot() +
  geom_density(data = df_residui, 
               aes(x = residui), fill = "lightblue", alpha = 0.5, color = NA) +
  geom_density(data = dati_normali, 
               aes(x = normali), color = "red") +
  labs(title = "Densità dei residui vs. Normale") +
  theme_minimal()

Omoschedasticità (varianza costante dei residui):

plot(mod_3, which = 3)

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

bptest(mod_3)

## 
##  studentized Breusch-Pagan test
## 
## data:  mod_3
## BP = 90.297, df = 5, 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, which = 4)

Solo il campione 1549 prova ad avvicinarsi alla 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
N.gravidanze madre: 3
Sesso neonato: F

Per le variabili Lunghezza e Cranio di cui non è specificato un valore considero quello medio:

Lunghezza neonato : 490 millimetri
Cranio neonato: 350 millimetri

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

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

Il peso previsto dal modello è 3328 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 v2

Guido Miglio

2025-03-05

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 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: