Adding reference levels without having continuous covariates go through zero

Hello! I’m quite new to brms so forgive me if any of the following doesn’t make sense. I’m happy to answer clarifying questions.

I’m fitting a survival model with two continuous covariates of weight and temperature, and two fixed ones of bander (1-54) and day (0-14). Basically, for each day of a bird’s life through 14 days after banding, which of the following impacts survival?

In order to interpret results easier, I don’t want to have a reference category. Currently, the way the model is coded means that bander1 and day0 are the reference levels. Therefore, the other banders and days are dependent on them. That is what I don’t want, as I want to be able to say, “this bander had a negative effect keeping everything else constant” instead of “this bander had a more negative effect than bander 1 on day 0”.

I want to additionally make sure that the weight and temperature aren’t going to go through zero and they will have their own effect size as well.

Here is my current model.

surv_model <- brm(survived ~ weight_scaled + 
                    mean_mintemp_prev5_scaled + 
                    bander_index + 
                    day_index + 
                    (1 | year_index) +
                    (1 | state_index),
                  data = data_model,
                  family = bernoulli(),
                  prior = new_priors,
                  iter = 20000,
                  warmup = 10000, 
                  chains = 3,
                  cores = 3,
                  thin = 10,
                  control = list(adapt_delta = 0.99))

My current computer:

  • Operating System: MacOS Ventura 13.0
  • brms Version: 2.22.0

Thank you! I look forward to hearing some ideas.

My only exposure to the type of model you are fitting was while reading this paper. Appendix B (which I am linking to) has a discrete time-to-event model with censored observations, which I believe is similar to your case. With that caveat in mind, here are some suggestions.

It looks like you’ve standardized your continuous predictors, which takes care of what you mention in the title of your post? Or did I misunderstand that?

One way to get bander deviations from the average is with partial pooling, i.e (1|bander_index) so that ranef(surv_model)$bander_index would give you the desired deviations.

If you want to do this without partial pooling, then regardless of the fact that your model is coded with a reference bander, as you say, you can always get bander differences from the average using tidybayes+emmeans:

library(tidyverse)
library(tidybayes)
library(emmeans)

surv_model %>%
  emmeans(~ bander_index) %>%
  contrast(method = "eff") %>%
  gather_emmeans_draws() %>%
  mean_qi()

Again, not really familiar with these models, but based on the tutorial I linked to, shouldn’t you be using a “cloglog” link function? And shouldn’t your binary response be coded 1=died, 0=did not die (yet)?

Finally, I do not understand why you choose to go with 20000 iterations and then thin them rather than do 2000 and keep everything.

To give the best advice possible, I’ll ask a few clarifying questions here.

  1. What effects are of interest here and what effects are supposed to control for nuisance variation? That is, are you primarily interested in how much variation there is across banders, while controlling for the effects of temperature and weight? Or are you more interested in the effect of e.g. temperature while controlling for bander and weight?
  2. Are most of the banding captures at around the same day, or are they spread out? That is, does your day covariate capture variation in days-since-banding that is primarily independent of variation in day-of-year, or closely related?
  3. Do you expect the impacts of day to be smooth and/or monotonic across days?
  4. Do you have reliable daily observations of whether each individual survived, or just detection histories that leave ambiguity about whether an individual has died or simply hasn’t been seen for a few days in a row?