Modelos Jerárquicos para la Estimación del Peso de Bebes diagnosticados con BPN en la ciudad de Medellín

El bajo peso al nacer (BPN) ha constituido un enigma en la ciencia a través de los tiempos. Múltiples son las investigaciones realizadas acerca de las causas que lo producen y las consecuencias que provoca. Su importancia no solo radica en lo que significa en la morbilidad y la mortalidad infantil, sino que estos niños tienen habitualmente múltiples problemas posteriores. Dentro de los factores de riesgo del BPN se han encontrado con mayor frecuencia en estudios realizados por diferentes autores, los siguientes: el embarazo en la adolescencia, la desnutrición en la madre, el hábito de fumar, la hipertensión arterial durante el embarazo, la sepsis cervicovaginal, la anemia y los embarazos gemelares, entre otros.

Descripción de la Base de Datos

La base de datos utilizada para el analisis corresponde al Registro de pacientes atendidos en las Instituciones Prestadoras de Servicios de Salud en la ciudad de Medellín con diagnóstico confirmado de BAJO PESO AL NACER. Se puede consultar la informacion en el siguiente enlace: http://medata.gov.co/dataset/bajo-peso-al-nacer/resource/51765ef7-f2e6-4795-9ea1-dd6e4e5852f6#{}

Inicialmente el conjuto de datos cuenta con 8174 observaciones y 19 varibles, se realiza un preprocesamiento de la informacion eliminando las observaciones que no tenían registro en la variable respuesta y se hace un renombramiento de las columnas, y se cambia por variables tipo factor y numericas algunas variables que eran de tipo caracter, entre otros ajustre otros quedandando asi con una base de 6253 observaciones y 10 variables para el análisis.

## LECTURA DE LA BASE DE DATOS
data <- read.csv("C:/Users/vreng/Desktop/Mod_Jerárquicos/sivigila_bpn.csv",sep=";")
data <- data[,-c(2,4,6,9,10,11,12,13,18)]

## RENOMBRANDO LAS VARIABLES
names(data)<-c("ID","EdadMadre","Sexo","Comuna", "TipoSalud", 
               "Peso","Talla", "SemanaGest", "EduMadre", "Año")
# SEXO
data$Sexo<-factor(data$Sexo,
                  levels = c("M","F"),labels = c("1","0"))
# COMUNA
data$Comuna<-factor(data$Comuna,
                    levels=c("Belen","Castilla","Aranjuez","Buenos Aires",
                             "Doce de Octubre","El Poblado","Guayabal","La America",
                             "LaCandelaria","Laureles","Manrique","Popular","Robledo",
                             "San Javier","Santa Cruz","Villa Hermosa"))
colSums(is.na(data))
data<-na.omit(data) # Elimina las obs que no tienen información respecto a la comuna.

# TIPOSALUD
data$TipoSalud<- as.factor(data$TipoSalud)

# PESO: hay 1 entrada sin peso registrado
data[data$Peso == "SD",]
data<-data[!(data$ID==734),]
data$Peso <- as.numeric(data$Peso)

# TALLA
data$Talla <- as.numeric(data$Talla)

# SEMANAGEST
data$SemanaGest <- as.numeric(data$SemanaGest)

dim(data)

## [1] 6253   10

	ID	EdadMadre	Comuna	TipoSalud	Peso	Talla	SemanaGest	EduMadre	Año
20	20	19	Doce de Octubre	S	2240	47	38	2	1
31	31	29	Doce de Octubre	C	2361	47	39	2	1
32	32	27	San Javier	S	2440	50	40	1	1
33	33	40	Buenos Aires	S	2140	45	37	2	1
35	35	25	Belen	C	2310	47	38	3	1
36	36	19	Villa Hermosa	S	2420	45	38	1	1
37	37	30	Robledo	C	2340	46	38	2	1
39	39	18	La America	S	2460	48	37	1	1

Descripción de las variables

ID: Número de Registro
EdadMadre: Edad de la madre en años
Sexo: Sexo del bebé (Masculino = 1. Femenino = 0)
Talla: Talla del bebé en centimetros
Peso: Peso del bebé en gramos
SemanaGest: Duración del embarazo
TipoSalud: C= Contributivo, S=Subsidiado, P=Excepción, E=Especial, N= No asegurado, I= Indeterminado/Pendiente
EduMadre: Nivel de educación de la madre(1= Primaria, 2= Secundaria, 3= Técnico o Superior, 4= Ninguna).
Comuna: Comuna en la que reside la madre
Año:Año en el que se tomó el registro

Objetivo.

El objetivo principal de este estudio es estimar el peso de los pacientes diagnosticados con bajo peso al nacer en la ciudad de Medellín, por medio de modelos mixtos.

Análisis descriptivo.

Primero, se estudia el comportamiento de las variables de tipo numérico se evidencia que la media de la variable Peso es de 2340 gramos, la Edad promedio de las madres que participaron en el estudio esta entre 25 y 26 años, por último la Talla promedio de los bebes es de 46.24 centimetros.

##    EdadMadre          Peso          Talla      
##  Min.   :11.00   Min.   :2000   Min.   :38.00  
##  1st Qu.:20.00   1st Qu.:2270   1st Qu.:45.00  
##  Median :24.00   Median :2365   Median :46.00  
##  Mean   :25.37   Mean   :2340   Mean   :46.24  
##  3rd Qu.:30.00   3rd Qu.:2440   3rd Qu.:48.00  
##  Max.   :57.00   Max.   :2499   Max.   :53.00  
##

Analizando las variables categóricas se evidencia que la comuna en la que se presentan un mayor peso promedio de los bebés es en Guayabal mientras que el peso promedio mas bajo se reporta en Buenos Aires y Villa Hermosa.

Modelos a considerar.

Para modelar la variable respuesta Peso se debe tener en cuenta lo siguiente:

Se ajustaran tres modelos lineales generalizados mixtos ulitlizando la función glmer del paquete lme4 con familia Gamma y función de enlace "log" pues esta es una generalización de la distribución exponencial, específicamente una distribución de probabilidad continua.

Modelo 1

Este modelo se ajusta con las variable Talla mas intercepto aleatorio asociado a la Comuna en la que reside la Madre.

\[Peso=\beta_0+\beta1Talla{ij}+b_{0i}\\ b_0 \sim N(0, \sigma^2_{b0})\]

Modelo 2

Modelo Ajustado con las variables Talla y SemanaGest a las cuales de se les aplica una transformación logaritmo para corregir la variación dentro de cada covariable y se mantiene el intercepto aleatorio asociado a Comuna

\[Peso=\beta_0+\beta1Log(Talla{ij})+\beta_2Log(SemanaGest_{ij})+b_{0i}\\ b_0 \sim N(0, \sigma^2_{b0})\]

Modelo 3

En este caso se realiza el mismo ajuste que se realizó en el segundo modelo pero se añade una tercera variable categórica EduMadre que corresponde al grado de escolaridad de las madres y se mantiene el intercepto aleatorio asociado a la Comuna.

\[Peso=\beta_0+\beta1Log(Talla_{ij})+\beta_2Log(SemanaGest_{ij})+ \beta_3EduMadre_{ij}+b_{0i}\\ b_0 \sim N(0, \sigma^2_{b0})\]

Construcción de los modelos.

Modelo 1

fit1<-glmer(Peso~(Talla)+(1|Comuna), data=data, family=Gamma(link = "log"))
summary(fit1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: Peso ~ (Talla) + (1 | Comuna)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##  77225.9  77252.8 -38608.9  77217.9     6249 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5695 -0.5707  0.2141  0.7462  2.8453 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Comuna   (Intercept) 2.754e-06 0.00166 
##  Residual             2.398e-03 0.04897 
## Number of obs: 6253, groups:  Comuna, 16
## 
## Fixed effects:
##              Estimate Std. Error t value Pr(>|z|)    
## (Intercept) 7.3705288  0.0095635  770.69   <2e-16 ***
## Talla       0.0083799  0.0002066   40.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr)
## Talla -0.995
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0139176 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

El modelo que se obtiene con las estimaciones es:

\[Log(\widehat{Peso})=7.3705288+0.0083799(Talla_{ij})+b_{0i}\\ b_0 \sim N(-3.285431e-05,~ 6.675801e-07)\]

Modelo 2

# Modelo 2
fit2<-glmer(Peso~log(Talla)+log(SemanaGest)+(1| Comuna), data=data, family=Gamma(link="log"))
summary(fit2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: Peso ~ log(Talla) + log(SemanaGest) + (1 | Comuna)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##  77136.5  77170.2 -38563.2  77126.5     6248 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5295 -0.5774  0.2138  0.7436  2.9913 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Comuna   (Intercept) 3.156e-06 0.001777
##  Residual             2.365e-03 0.048627
## Number of obs: 6253, groups:  Comuna, 16
## 
## Fixed effects:
##                 Estimate Std. Error t value Pr(>|z|)    
## (Intercept)      5.41887    0.09416   57.55   <2e-16 ***
## log(Talla)       0.35514    0.01421   24.98   <2e-16 ***
## log(SemanaGest)  0.26951    0.02379   11.33   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg(Tl)
## log(Talla)  -0.427       
## log(SmnGst) -0.821 -0.165
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.034672 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

El modelo que se obtiene es el siguiente:

\[log(\widehat{Peso})=5.41887+0.35514Log(Talla_{ij})+0.26951Log(SemanaGest_{ij})+b_{0i}\\ b_0 \sim N(-5.00054e-05,~8.849179e-07)\]

Modelo 3

# Modelo 3
fit3<-glmer(Peso~log(Talla)+log(SemanaGest)+(EduMadre)+(1|Comuna),data=data, family= Gamma(link="log"))
summary(fit3)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: Peso ~ log(Talla) + log(SemanaGest) + (EduMadre) + (1 | Comuna)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##  77134.9  77195.6 -38558.5  77116.9     6244 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5567 -0.5753  0.2120  0.7428  3.0658 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Comuna   (Intercept) 2.709e-06 0.001646
##  Residual             2.361e-03 0.048594
## Number of obs: 6253, groups:  Comuna, 16
## 
## Fixed effects:
##                 Estimate Std. Error t value Pr(>|z|)    
## (Intercept)     5.397281   0.092075  58.618   <2e-16 ***
## log(Talla)      0.353501   0.014050  25.160   <2e-16 ***
## log(SemanaGest) 0.276228   0.023151  11.931   <2e-16 ***
## EduMadre2       0.002889   0.002182   1.324   0.1854    
## EduMadre3       0.005123   0.002344   2.186   0.0288 *  
## EduMadre4       0.006369   0.005977   1.066   0.2866    
## EduMadre5       0.033022   0.015016   2.199   0.0279 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg(Tl) lg(SG) EdMdr2 EdMdr3 EdMdr4
## log(Talla)  -0.433                                   
## log(SmnGst) -0.817 -0.165                            
## EduMadre2   -0.027 -0.063  0.047                     
## EduMadre3   -0.057 -0.061  0.080  0.806              
## EduMadre4   -0.022  0.019  0.004  0.311  0.293       
## EduMadre5   -0.013 -0.007  0.016  0.126  0.122  0.046
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.127231 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?

\[Log(\widehat{Peso})=5.397281+0.353501Log(Talla_{ij})+\\0.276228Log(SemanaGest_{ij})+ 0.002889EduMadre_{ij}[2]+0.005123EduMadre_{ij}[3]\\+0.006369EduMadre_{ij}[4]+ 0.033022EduMadre_{ij}[5]+ b_{0i}\] \[ b_0 \sim N(-2.889884e-05, 6.582704e-07)\]

Comparación de los modelos

Se evidencia que el modelo 3 es el que tiene mejor desempeño, si se fija la atención en el criterio de información de Akaike (AIC), vemos el menor valor corresponde tambien al tercer modelo.

Análisis de residuales.

A continuación se presenta el análisis de residuales para cada uno de los modelos ajustados, el análisis se compone de gráficos de normalidad para los errores de cada modelo, grafico de residuales para evaluar el supuesto de homocedasticidad y gráficos de valores ajustados vs. valores reales.

Como en el apartado anterior a partir de la prueba de razón de verosimilitud y el AIC se seleccionó el modelo numero tres como el modelo con mejor desempeño se procede a evaluar la distribución del intercepto aleatorio para el mismo.

Se concluye que no hay evidencia muestral suficiente para rechazar el supuesto de normalidad de \(b_{0i}\)

Para corregir el supuesto de normalidad de los errores de los modelos y de hocedasticidad se aplicaron algunas transformaciones que no generaron muchos cambios en el ajuste, finalmente el mejor ajuste se da sólo con la transformación logaritmo para algunas covariables, dicho esto el modelo con mejor desempeño resulta ser el modelo número 3.

Presentación del Modelo Seleccionado.

\[Log(\widehat{Peso})=5.397281+0.353501Log(Talla_{ij})+\\0.276228Log(SemanaGest_{ij})+ 0.002889EduMadre_{ij}[2]+0.005123EduMadre_{ij}[3]\\+0.006369EduMadre_{ij}[4]+ 0.033022EduMadre_{ij}[5]+ b_{0i}\]

Referencias.

Barajas, F. H., & Martínez, J. L. L. (2022, January 26). Modelos Mixtos Con R. index.knit. Retrieved January 31, 2022, from https://fhernanb.github.io/libro_modelos_mixtos/

Gelman, A., & Hill, J. (2009). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press.

R: A Language and Environment for Statistical Computing,author = R Core Team, organization = R Foundation for Statistical Computing,address = Vienna, Austria, year = 2020, url = https://www.R-project.org/