Hi all,
I’m looking to get a deeper understanding of the hierarchical linear model I’ve built in brms
to see if there’s any way I can improve it / to see if it’s doing what I think it’s doing. To describe the model and data, I have a variable, CHG, which I assume changes as a function of plasma concentration, PKCONC, and visit, AVISIT, for a series of subjects, who all have unique SUBJIDs.
For anyone trying to get a better sense of the data, PKCONC and CHG are both continuous variables, while AVISIT is categorical and can take on one of two values, “Visit 1” or “Visit 2”.
I started with a simple model bf(CHG ~ 1 + PKCONC + (1 + PKCONC | SUBJID)
(full code below) which to me means that I’m estimating a population level intercept and slope (as a function of PKCONC) and that I’m also estimating individual level intercepts and slopes (as a function of PKCONC), for each individual (each unique SUBJID). I then tried a bunch of other variations on that theme, and selected the best model using LOO-CV and WAIC.
# Define the model formula
base_formula <- bf(CHG ~ 1 + PKCONC + (1 + PKCONC | SUBJID))
# Fit the model using the brm function
base_fit <- brm(base_formula,
data = qtc2,
family = "gaussian",
cores = 4,
chains = 4,
file = "base",
save_all_pars = TRUE)
The final model was bf(CHG ~ 1 + PKCONC * AVISIT + (1 + PKCONC * AVISIT | SUBJID)
with a standard deviation estimated for each visit level (below) but I’m wondering if I should use a different syntax to specify the IDs following here. Would it make sense to instead use CHG ~ 1 + PKCONC * AVISIT + (1 + PKCONC * AVISIT | gr(SUBJID, id = "i"))
, and how does this structurally change the model (I’m assuming that it now treats data for each SUBJID as correlated)?
final_fit <- brm(bf(CHG ~ 1 + PKCONC * AVISIT + (1 + PKCONC * AVISIT | SUBJID),
sigma ~ AVISIT),
data = qtc2,
family = gaussian(),
# prior = NULL,
prior = c(set_prior("normal(0, 1)", class = "b"), # Fixed effects, usually small, around 0
set_prior("exponential(1)", class = "sd"), # Random effects standard deviations, usually not very uncertain
set_prior("lkj(1)", class = "cor")), # Random effects correlations, not a ton of correlation to be expected
warmup = 4000,
iter = 8000,
chains = 4,
control = list(adapt_delta = 0.99),
cores = 4) #,
# file = "final_fit")
Would it be syntactically better to instead group the individuals by ~ x + (1 | gr(subject, by = x))
or am I totally off the mark on that one?
Thank you very much and I really appreciate the package.
- Operating System: Windows 10
- brms Version: 2.14.4