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