Aplicación Modelos Lineales Mixtos
Carlos Mario Calle González
25/3/2020
Base de datos
En este documneto se utilizará la base de datos Theoph del paquete base de R, que contiene información sobre el proceso desde la administración hasta la eliminación (farmacocinética) del medicamento antiasmático teofilina. La base de datos tiene 132 datos y cinco variables, las cuales son:
Subject: un factor ordenado con los niveles del 1, …, 12 que identifica al sujeto sobre quien se hizo la observación. El orden esta dado por el incremento de la concentración máxima de teofilina observada.
Wt: peso del sujeto (kg).
Dose: dosis de teofilina administrada oralmente al sujeto (mg/kg).
Time: tiempo transcurrido desde la administración del fármaco cuando se extrajo la muestra (h).
conc: concentración de teofilina en la muestra (mg / L).
head(Theoph)## Subject Wt Dose Time conc
## 1 1 79.6 4.02 0.00 0.74
## 2 1 79.6 4.02 0.25 2.84
## 3 1 79.6 4.02 0.57 6.57
## 4 1 79.6 4.02 1.12 10.50
## 5 1 79.6 4.02 2.02 9.66
## 6 1 79.6 4.02 3.82 8.58Objetivo
El objetivo del documento será ajustar y comparar algunos modelos multinivel para la variable respuesta concentración de teofilina, así como identificar algunas relaciones entre las diferentes variables y la variable respuesta.
Análsis descriptivo
Una primer relación a identificar será como se comporta la concentración del fármaco a medida que transcurren las horas.
attach(Theoph)
library(ggplot2)
ggplot(data = Theoph, aes(x = Time, y = conc, color = Subject)) +
geom_point() +
theme_bw() +
facet_wrap(~ Subject) +
theme(legend.position = "none")
Se observa que para todos los sujetos, la concentración de teofilina comienza a aumentar a medida que transcurre el tiempo hasta alcanzar un punto máximo, a partir del cual comienza a decrecer de forma lenta.
Ahora se realizará un gráfico para observar que posible relación existe entre la concentración del fármaco y la dosis administrada.
library(ggplot2)
library(RColorBrewer)
ggplot(Theoph, aes(Dose, conc, color= Subject) ) +
geom_point() +
labs(colour="Sujeto")+ scale_color_brewer(palette="Paired") + theme_bw()
De la figura anterior se puede observar que para los sujetos que recibieron dosis de teofilina mayores 4 y menores a 5 mg por kg, la concentración del fármaco es menos variante. Mientras que para aquellos sujetos que se encuentran fuera de ese rango, la concentración tiende a alcanzar máximos más altos y a ser más dispersa.
Modelos a considerar
Con lo planteado anteriormente se pueden proponer modelos en los cuales se introduzca efectos aleatorios por sujeto, además de componentes de segundo orden para la variable tiempo. Se tendra como modelo inicial aquel en el que solo se considera la variable tiempo desde la administracoión de la teofilina, luego se le agregara la variable dosis. Los modelos propuestos se presentan a continuación:
Modelo 1
\(Concentration_{ij} \sim N(\mu_{ij}, \sigma_{Concentration}^2) \\ \mu_{ij} = \beta_0 + \beta_1Tiempo_{ij} + \beta_2Tiempo_{ij}^2+ b_{0i} \\ b_0 \sim N(0, \sigma_{b_0}^2)\)
Con vector de parámetros \(\Theta = (\beta_0, \beta_1, \beta_2, \sigma_{Concentration}, \sigma_{b_0})^\top\)
Modelo 2
\(Concentration_{ij} \sim N(\mu_{ij}, \sigma_{Concentration}^2) \\ \mu_{ij} = \beta_0 + \beta_1Time_{ij} + \beta_2Time_{ij}^2 + \beta_3Dose_i + b_{0i} \\ b_0 \sim N(0, \sigma_{b_0}^2)\)
Con vector de parámetros \(\Theta = (\beta_0, \beta_1, \beta_2, \beta_3, \sigma_{Concentration}, \sigma_{b_0})^\top\)
Construcción de los modelos
A continuación se presenta el código del ajuste de los modelos antes presentados. Para ello se usará el paquete nlme.
library(nlme)
mod1 <- lme(conc ~ Time + I(Time^2), random = ~1| Subject, method = 'REML')
summary(mod1)## Linear mixed-effects model fit by REML
## Data: NULL
## AIC BIC logLik
## 643.2821 657.5811 -316.641
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 0.3757623 2.51871
##
## Fixed effects: conc ~ Time + I(Time^2)
## Value Std.Error DF t-value p-value
## (Intercept) 4.705009 0.3731011 118 12.610546 0.0000
## Time 0.315840 0.0985189 118 3.205884 0.0017
## I(Time^2) -0.019505 0.0041192 118 -4.735236 0.0000
## Correlation:
## (Intr) Time
## Time -0.696
## I(Time^2) 0.564 -0.947
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.93309081 -0.65903989 0.04447066 0.70020071 2.47538169
##
## Number of Observations: 132
## Number of Groups: 12mod2 <- lme(conc ~ Time + I(Time^2)+ Dose, random = ~1| Subject, method = 'REML')
summary(mod2)## Linear mixed-effects model fit by REML
## Data: NULL
## AIC BIC logLik
## 644.7948 661.907 -316.3974
##
## Random effects:
## Formula: ~1 | Subject
## (Intercept) Residual
## StdDev: 0.3949838 2.51871
##
## Fixed effects: conc ~ Time + I(Time^2) + Dose
## Value Std.Error DF t-value p-value
## (Intercept) 3.296271 1.6409586 118 2.008748 0.0468
## Time 0.315700 0.0985191 118 3.204450 0.0017
## I(Time^2) -0.019499 0.0041192 118 -4.733737 0.0000
## Dose 0.304607 0.3454332 10 0.881812 0.3986
## Correlation:
## (Intr) Time I(T^2)
## Time -0.157
## I(Time^2) 0.127 -0.947
## Dose -0.974 -0.001 0.002
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.0560701 -0.5975322 0.0782864 0.6806460 2.3524557
##
## Number of Observations: 132
## Number of Groups: 12Bajo los resultados anteriores se tienen los siguientes vectores de parámetros:
Para el modelo 1 \(\Theta = (\beta_0 = 4.705009, \beta_1 = 0.315840, \beta_2 = -0.019505, \sigma_{Concentration} = 2.51871, \sigma_{b_0} = 0.3757623)^\top\)
Para el modelo 2 \(\Theta = (\beta_0 = 3.296271, \beta_1 = 0.315700, \beta_2 = -0.019499, \beta_3 = 0.304607, \sigma_{Concentration} = 2.51871, \sigma_{b_0} = 0.3949838)^\top\)
Comparación de los modelos
La prueba de hipotesis a contrastar para comparar los modelos es la siguientes:
\(H_0\): La variable Dose no aporta al modelo, versus, \(H_a\): La variable Dose si aporta al modelo. Para obtener un Valor-P más preciso se utilizará simulación con la función simulate.lme
simul <- simulate.lme(object=mod1, m2=mod2, method = 'ML', nsim=10000)
lrts_nlme <- -2 * (simul$null$ML[, 2] - simul$alt$ML[, 2])
acumulada1 <- ecdf(x=lrts_nlme)
lrt <- -2 * (logLik(mod1) - logLik(mod2))
1 - acumulada1(lrt)## [1] 0.4888Del resultado anteriore se concluye que la variable Dose no aporta al modelos, pues con un Valor-P de \(0.4935\) no se puede rechazar la hipotesis nula.
Análisis de residuales
plot(mod1)
En el gráfico de los residuales contra los valores ajustados se evidencia un conjunto de datos atípicos. Sin embargo al no ser pocos se puede elegir el modelo 1.
Modelo Final
\(Concentration_{ij} \sim N(\widehat\mu_{ij}, \widehat\sigma_{Concentration}^2) \\ \widehat\mu_{ij} = 4.705009 + 0.315840*Tiempo_{ij} + -0.019505*Tiempo_{ij}^2+ b_{0i} \\ \widehat\sigma_{Concentration} = 2.51871\)