En el siguiente post explicaremos como funciona el paquete bmrs usando modelos multinivel, para este caso utilizaremos el ejemplo sobre el tiempo de recurrencia de una infección en pacientes renales publicado por McGilchrist and Aisbett (1991)
La base de datos kidney contiene la siguientes variables:
time: representa el tiempo de recurrrencia de la infección
censored: indica si el tiempo está censurado a la derecha (1) o no censurado (0)
patient: es el id del paciente
recur : indica si es la primera o segunda recurrencia en ese paciente
age: edad del paciente
sex: sexo del paciente
disease: tipo de enfermedad
library("Rcpp")
library("brms")
data("kidney")
head(kidney, n=5) time censored patient recur age sex disease
1 8 0 1 1 28 male other
2 23 0 2 1 48 female GN
3 22 0 3 1 32 male other
4 447 0 4 1 31 female other
5 30 0 5 1 10 male otherbmrsEl núcleo para realizar el ajuste de los modelos con este paquete es la función brm. Para la base de datos kidney supongamos que se quiere predecir el tiempo de recurrencia (posiblemente censurado) utilizando un modelo log-normal en el que intercepto así como el efecto de la edad están anidados dentro de los pacientes. Entonces escribimos el siguiente código para representar el ajuste del modelo.
results=‘hide’
fit1 <- brm(formula = time | cens(censored) ~ age * sex + disease
+ (1 + age|patient),
data = kidney, family = lognormal(),
prior = c(set_prior("normal(0,5)", class = "b"),
set_prior("cauchy(0,2)", class = "sd"),
set_prior("lkj(2)", class = "cor")),
warmup = 1000, iter = 2000, chains = 4,
control = list(adapt_delta = 0.95))Todo antes del símbolo \(\sim\) se refiere a la parte de la variable de respuesta.
time | cens(censored) indica lo siguiente cens constituye la función interna que maneja los datos censurados, y censored es la variable que contiene información sobre la censura.Otras funciones disponibles en este contexto son weights y disp para permitir diferentes tipos de ponderación.
Todo lo que hay a la derecha de \(\sim\) especifica los predictores.
Los términos a nivel de grupo son de la forma (coefs | grupo), donde coefs contiene una o más variables cuyos efectos se supone que varían con los niveles del factor de agrupación dado en grupo Bürkner (2017).
Para este caso solo se especifica un término a nivel de grupo en el que 1+age son los coeficientes que varían con el factor de agrupación patient. Asimismo,lo que no se reconoce como parte de un término a nivel de grupo se trata como un efecto a nivel de población. En este ejemplo, los efectos a nivel de población son la age, el sex y la disease.
Aquí brms tiene una gran variedad de familias por ejemplo para la regresión lineal se puede utilizar la familia gaussiana o student combinada con el enlace identidad, los datos dicótomicos o categóricos las familias son bernoulli, binomial y categórica combinadas con el enlace logit, cuando tenemos datos de recuento se utiliza la poisson, negbinomial y geométrica. Para nuestro ejemplo utilizamos la family=lognormal() que implica una “supervivencia” log-normal para la variable respuesta time.
Cada efecto a nivel de población tiene su correspondiente parámetro de regresión, estos parámetros se denominan como b_<coef>, where <coef> representa el nombre del correspondiente a nivel de población. Supongamos, por ejemplo, que queremos establecer una prior normal con media \(0\) y desviación estándar \(10\) sobre el efecto de la age y una prior Cauchy con localización \(1\) y escala \(2\) sobre sexfemale. Entonces, podemos escribir
prior <- c(set_prior("normal(0,10)", class = "b", coef = "age"),
set_prior("cauchy(1,2)", class = "b", coef = "sexfemale"))Ahora bien, en los efectos a nivel de grupo de cada factor de agrupamiento se tiene un parámetro de deviacion estándar , en brms estos parámetros se definen como sd<group>_<coef> por lo cual sd_patient_Intercept y sd_patient_age son los nombres para nuestro ejemplo, si desea establecer prior diferentes para cada uno de los parametros utilizando sentencias como
set_prior("student_t(3,0,5)", class = "sd", group = "patient")
set_prior("student_t(3,0,5)", class = "sd")Si hay más de un efecto a nivel de grupo por factor de agrupación, entonces se realizan las estimaciones de las correlaciones entre los efectos a nivel de grupo y aquí utilizamos la prior de correlación \(LKJ\) con parámetro \(\zeta>0\) en brms esto se abrevia como \(lkj(zeta)\) y los parámetros de la matriz de correlación se denominan cor_<grupo>, (por ejemplo, cor_paciente), de modo que
set_prior("lkj(2)", class = "cor", group = "patient")StanAquí lo importante, es que ademas de elegir el numero de iteraciones, muestras de calentamiento y cadenas el usuario puede controlar el comportamiento del muestreador NUTS utilizando el argumento control. Con esto podemos reducir (o eliminar en la medida de lo posible) la cantidad de transiciones divergentes que pueden introducir sesgos en las muestras posteriores obtenidas.
Para ello utilizamos la siguiente sintaxis control = list(adapt_delta = <x>) , donde <x> normalmente debe tomar un valor por defecto entre \(0.8\) y \(1\).
El modelo fit1 se ajusta utilizando \(4\) cadenas , cada una con \(2000\) iteraciones las cuales las primeras \(1000\) son calentamiento para calibrar el muestreador , lo que lleva a un total de \(4000\) muestra posteriores.
summary(fit1, waic=TRUE) Family: lognormal
Links: mu = identity; sigma = identity
Formula: time | cens(censored) ~ age * sex + disease + (1 + age | patient)
Data: kidney (Number of observations: 76)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~patient (Number of levels: 38)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.43 0.35 0.02 1.23 1.01 718 266
sd(age) 0.01 0.01 0.00 0.03 1.01 571 247
cor(Intercept,age) -0.16 0.46 -0.92 0.74 1.01 623 230
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.78 0.95 0.95 4.62 1.00 1063 1381
age 0.01 0.02 -0.03 0.06 1.00 1303 2323
sexfemale 2.39 1.10 0.26 4.60 1.00 1609 2380
diseaseGN -0.40 0.51 -1.40 0.62 1.00 2432 2720
diseaseAN -0.51 0.51 -1.52 0.48 1.00 2460 2441
diseasePKD 0.62 0.72 -0.79 2.00 1.00 2435 2488
age:sexfemale -0.02 0.03 -0.07 0.03 1.00 1568 2321
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.14 0.13 0.91 1.43 1.00 2161 2354
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Observando los efectos a nivel de grupo , el parámetro de desviación estándar de age es sospechosamente pequeño, para ello se realizamos una prueba de hipótesis
hypothesis(fit1, "Intercept - age > 0", class = "sd", group = "patient")Hypothesis Tests for class sd_patient:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
1 (Intercept-age) > 0 0.42 0.35 0.03 0.99 104.26 0.99
Star
1 *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.Ahora veamos los resultados de manera gráfica
Gráficos de trazas y densidades de todos los parámetros relevantes del modelo kidney analizado anteriormente.
pp_check(fit1,type = "intervals")