I am modeling heroin use during the first year of opioid agonist treatment based on different levels of amphetamine use at start of treatment. I am using brms. Days of amphetamine use at baseline in the 28 days prior to start of treatment (ats_baseline) is a bounded count variable. Time from start of treatment is in years (variable yearsFromStart) modelled using a thin-plate regression spline to allow for departures from linearity. The outcome, days of heroin use in the 28 days prior to each measurement (heroin28), is also a bounded 0-28 count, hence I am using a betabinomial regression model.
Here is the dataset
atsUseInOTP_alloutcomes_noMissing_inf.RData (17.3 KB)
I want to perform a prior sensitivity check using the priorsense package. As @n-kall describes in this vignette for the priorsense package when using splines it is better to powerscale the predictions than the parameters themselves. This requires tagging each prior using the new prior syntax in brms. Like so.
fit_betabinomial_heroin_s2_ps <- brm(formula = heroin28 | trials(set) ~ ats_baseline + s(yearsFromStart) + s(yearsFromStart, ats_baseline) + (1 | encID),
data = d,
family = beta_binomial(),
prior = c(prior(prior = normal(0,4),
class = "b",
coef = "ats_baseline",
tag = "b"),
prior(prior = normal(0,4),
class = "b",
coef = "syearsFromStart_1",
tag = "b"),
prior(prior = normal(0,4),
class = "b",
coef = "syearsFromStartats_baseline_1",
tag = "b"),
prior(prior = normal(0,4),
class = "b",
coef = "syearsFromStartats_baseline_2",
tag = "b"),
prior(prior = normal(0,10),
class = "sds",
coef = "s(yearsFromStart)",
tag = "sds"),
prior(prior = normal(0,10),
class = "sds",
coef = "s(yearsFromStart, ats_baseline)",
tag = "sds")),
save_pars = save_pars(all=TRUE),
seed = 1234,
refresh = 0,
chains = 4,
warmup = 1e3,
iter = 2e3)
The predictions I would like to extract are levels of heroin28 at
each of five candidate levels of ats_baseline: 0, 7, 14, 21, and 28 days use, at each of five timepoints: 0 years, 0.25, 0.5, 0.75, and 1 year. Now we follow the procedure in the vignette. We bind the draws of the predictions and measures to the posterior draws using posterior::bind_draws(). For this we will use the predictions_as_draws helper function provided by priorsense which transforms predictions from brms to draws objects. The idea is that we are attaching predictions at various levels of the outcome to the posterior draws object, and then doing power-scaling on THOSE predictions to assess prior sensitivity.
# We first define a log-score function.
logscore <- function(x) {
as.matrix(rowSums(log_lik(x)))
}
# now let's get the grid of predictions
tibble(encID=nEnc+1, set=28) %>%
data_grid(yearsFromStart = seq(0,1,by=0.25),
ats_baseline=c(0,7,14,21,28),
encID,
set) -> predGrid
# how many cells?
nGrid <- nrow(predGrid)
# and the names of each cell
expand_grid(paste0("yearsFromStart",seq(0,1,by=0.25)),
paste0("ats_baseline",c(0,7,14,21,28))) %>%
unite(col = "cell") %>%
pull -> cellNames
# now to add the predictions to the posterior draws
post_draws_s2_ps <- as_draws_df(fit_betabinomial_heroin_s2_ps) %>%
bind_draws(log_lik_draws(fit_betabinomial_heroin_s2_ps)) %>%
bind_draws(predictions_as_draws(x = fit_betabinomial_heroin_s2_ps,
predict_fn = bayes_R2,
prediction_names = "R2",
summary = FALSE)) %>%
bind_draws(predictions_as_draws(x = fit_betabinomial_heroin_s2_ps,
predict_fn = logscore,
prediction_names = "logscore")) %>%
bind_draws(predictions_as_draws(x = fit_betabinomial_heroin_s2_ps,
predict_fn = posterior_epred,
prediction_names = cellNames,
newdata = predGrid,
allow_new_levels = TRUE))
# number of variables created in the function above
nVarPredDraws <- length(variables(as_draws(post_draws_s2_ps)))
# The first ten variables are the posterior estimates, and the variables we want to power scale, that we created in the post_draws_s2_ps object immediately above, were pasted onto the end of the posterior draws objects AFTER the R squared and logscore variable (i.e. the last 25 columns). So we want R-squared, logscore, plus all the new prediction columns. We'll extract their names first
variables(as_draws(post_draws_s2_ps))[c(1:10,(nVarPredDraws-(nGrid+1)):nVarPredDraws)] -> varsToPS
# finally to get the powerscale assessment for the coefs and the predictions. The function actually power-scales all priors unless we tell it what we are specifically interested in (specified above in varsToPS object), the fixed effects coefficients, R-squared, logscore, and the predictions
powerscale_sensitivity(x = post_draws_s2_ps,
variable=varsToPS) %>%
tt() %>%
format_tt(num_fmt="decimal") -> powerscale_sens_heroin_s2_ps
powerscale_sens_heroin_s2_ps
This is what should print in the viewer. Lots of variables so please excuse the lengthy output. You can see that there are quiet a few prior-data conflicts on the parameters
…but, as pointed out in the vignette, these are not what we are interested in. We are interested in the predictions. In the past I have had a lot of success with powerscaling these predictions using similar betabinomial models and the same priors. But this time, unfortunately there are some prior-data conflicts on the predictions as well.
Can anyone (e.g. @n-kall @avehtari?) give me some advice on how to resolve these conflicts. potentially on what priors to set on the spline parameters?
