Hi Stan forum,
If possible, I would love to have someone review my brms model to make sure it is doing what I want it to do! The model is fairly simple (it’s just 2 regressions), but it looks kind of weird the way I’ve implemented it, so I thought it would be worth double-checking my work.
I am attempting to test whether the slope of a response variable resp by it’s explanatory variable expl changes between two time periods. I’m looking at multiple species at once, and I’m hoping to see whether the overall slope changes between these two time periods. Since there are only two time periods, I am simply doing the regression twice, and taking the difference between the slope posteriors for each period.
The regression I’m doing is log(resp) ~ log(expl), so in order to put all of the observations on the same scale, I’ve opted to pre-translate my data, so that I take the log then rescale (to z-scores) both resp and expl to produce log_z.resp and log_z.expl (despite my naming convention, I log first, then take z-scores).
I have rbind’d together the dataset for the two time periods (before and after) with different column names and a boolean for when to use each column. I add together the two formulas using bf, and then add a set_rescor term—I couldn’t set it to be TRUE (I don’t actually want any correlation between the two epochs), so I set it to FALSE and added || bars to avoid any cross-epoch within-species correlations. (I am not sure if that was correct or not.)
Here is the model:
bf(
log_z.resp_before | subset(before) ~ 1 + log_z.expl_before + (1 + log_z.expl_before || sp)
)+
bf(
log_z.resp_after | subset(after) ~ 1 + log_z.expl_after + (1 + log_z.expl_after || sp)
)+
set_rescor(FALSE)
Questions
- I am hoping to get a species-structured regression of two log and rescaled variables across two time periods independently. Does my model currently do this?
- Is my
z-score(log(trait))transformation reasonable? - Will my use of
set_rescormess up the posterior? - I don’t have any manually set priors—weakly or non-informative priors work for what I’m doing anyway. For the sake of completeness, should I do this?
- I am going to replace the
|| spterm with|| gr(species, cov = phy)to make this analysis phylogenetically controlled. Is this the appropriate thing to do? - Is there anything else out of place here?
Full R code
# Get data
dat_raw <- read.csv("path/to/file") # Filled in with path
keep_sp <- dat_raw %>% # Make list of species which have at least 30 observations during both time periods
group_by(epoch,sp) %>%
tally %>%
pivot_wider(names_from="epoch",values_from="n") %>%
filter(before>=30,after>=30) %>%
pull(sp)
dat <- dat_raw %>%
filter(sp%in%keep_sp) %>%
filter(!is.na(resp),!is.na(expl),!is.na(sp)) %>%
group_by(sp) %>%
mutate(
log_z.resp=as.vector(scale(log(resp))),
log_z.expl=as.vector(scale(log(expl)))
) %>%
ungroup %>%
mutate(
sp=as.factor(sp)
) %>%
select(sp,log_z.resp,log_z.expl,epoch)
# Format data
dat_rbind <- rbind( # Combine the two time periods into a single data frame
dat %>%
filter(epoch=="before") %>%
mutate(before=T,after=F) %>%
rename(
log_z.resp_before=log_z.resp,
log_z.expl_before=log_z.expl
) %>%
mutate(
log_z.resp_after=NA,
log_z.expl_after=NA
),
dat %>%
filter(epoch=="after") %>%
mutate(before=F,after=T) %>%
rename(
log_z.resp_after=log_z.resp,
log_z.expl_after=log_z.expl
) %>%
mutate(
log_z.resp_before=NA,
log_z.expl_before=NA
)
)
# Load model
model <- brm(
bf(
log_z.resp_before | subset(before) ~ 1 + log_z.expl_before + (1 + log_z.expl_before || sp)
)+
bf(
log_z.resp_after | subset(after) ~ 1 + log_z.expl_after + (1 + log_z.expl_after || sp)
)+
set_rescor(FALSE),
data=dat_rbind,
chains=0
)
model_copy <- model # To avoid updating `model` itself in the next step
# Fit model
fit <- update(
model_copy,recompile=FALSE,
chains=4,
cores=4,
iter=6000,
warmup=2000
)
# View model
fit %>%
as.data.frame(variable=c("b_logzrespbefore_log_z.expl_before","b_logzrespafter_log_z.expl_after")) %>%
mutate(delta=b_logzrespafter_log_z.expl_after-b_logzrespbefore_log_z.expl_before) %>%
pull(delta) %>%
quantile(., probs=c(0.055, 0.945)) # Get 89% CI