# Error message executing powerscale\_sensitivity() function

**URL:** <https://discourse.mc-stan.org/t/error-message-executing-powerscale-sensitivity-function/39885>\
**Category:** General\
**Tags:** brms, techniques\
**Created:** [July 14, 2025, 2:36am UTC](https://discourse.mc-stan.org/t/error-message-executing-powerscale-sensitivity-function/39885 "2025-07-14T02:36:56Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![llewmills](https://avatars.discourse-cdn.com/v4/letter/l/4da419/32.png) [@llewmills](https://discourse.mc-stan.org/u/llewmills)\
**Post date:** [July 14, 2025, 2:36am UTC](https://discourse.mc-stan.org/t/error-message-executing-powerscale-sensitivity-function/39885/1 "2025-07-14T02:36:56Z")

</div>

I am running an analysis based closely on [this notebook](https://users.aalto.fi/~ave/casestudies/Nabiximols/nabiximols.html) by Aki Vehtari. I want to estimate the effect of amphetamine use prior to start of treatment on number of days of cannabis use during treatment among people enrolled in an opioid treatment program.

The data come from routinely collected outcomes data from clients enrolled in opioid agonist treatment at public drug and alcohol clinics. The goal is to examine the effects of amphetamine use prior to starting treatment on clinical outcomes during treatment in this cohort. The outcome variables in this model is `outcome` and is a bounded count indicating days of cannabis use in the previous 28 day period. The predictors are `ats_baseline` indicating days used in the 28 days prior to starting treatment and `yearsFromStart` indicating time in years that the outcome variable was observed. Using a betabinomial model the goal is to estimate differences in days of cannabis use at various points throughout the first year of treatment based on candidate values of `ats_baseline` at baseline and then three-monthly intervals during the first year of treatment. `yearsFromStart` and the interaction between `yearsFromStart` and `ats_baseline` are represented by a spline. `set` is a variable indicating the number of trials for each measure (always 28).

When I originally ran this model I used the default priors in brms. The model converged ok but when I did a prior sensitivity check using the priorsense package I got several prior-data conflicts. Here is the output

 ![priorSensitivity_spline2](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/1/5/15f4bd5ea809a363204601cd4505b36c22d16630.jpeg)

Professor Vehtari advised than in the case of a spline model is was better to perform the prior-scaling sensitivity analysis on the predictions rather than the posterior estimates, and referred me to [this vignette](https://n-kall.github.io/priorsense/articles/selecting_priors_and_quantities.html) detailing how to power-scale priors in a model with splines. This involves specifying priors on betas and the sds parameters for the splines and inclusion of the `tag = ` argument in each of the prior specifications.

Here is the data.  
[amphUseInOTP\_cann.RData](https://discourse.mc-stan.org/uploads/short-url/bkKk9uTxc9gXkOdtqUVUvT9Y2cQ.RData) (22.7 KB)

Code to run the model is below

```no-highlight
load(file = "amphUseInOTP_cann.RData") # object is amphUseInOTP_cann

# give the data object a generic name. Note the days' cannabis use variable is called `outcome` abd is a bounded count between 0 and 28

workingDF <- amphUseInOTP_cann

workingDF %>%
  distinct(encID) %>%
    nrow -> nEnc

# run betabinomial spline model
fit_betabinomial_cann_s2_ps <- brm(formula = outcome | trials(set) ~ ats_baseline + s(yearsFromStart) + s(yearsFromStart, ats_baseline) + (1 | encID),
                                   data = workingDF,
                                   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)

# define a logscore function
logscore <- function(x) {
  as.matrix(rowSums(log_lik(x)))
}

```

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.

First we get the grid of predictions (using the data\_grid function from modelr), for each of five candidate levels of amphetamine use at baseline - 0 days amphetamine use in the 28 days prior to baseline, 7 days, 14 days, 21 days, and 28 dyas - at each of five timepoints in the first year of treatment: 0 years from start of treatment, 0.25 years, 0.5 years, 0.75 years, 1 year.

```no-highlight
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 combinations of amphetamine use at baseline and yearsFromStart?

```no-highlight
nGrid <- nrow(predGrid)

```

We also need to create names for each of the 25 `ats_baseline` by `yearsFromStart` combinations .

```no-highlight
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 we add the predictions to the posterior draws using bind\_draws and predictions\_as\_draws

```no-highlight
post_draws_s2_ps <- as_draws_df(fit_betabinomial_cann_s2_ps) %>%
  bind_draws(log_lik_draws(fit_betabinomial_cann_s2_ps)) %>%
  bind_draws(predictions_as_draws(x = fit_betabinomial_cann_s2_ps,
                                  predict_fn = bayes_R2,
                                  prediction_names = "R2", 
                                  summary = FALSE)) %>%
  bind_draws(predictions_as_draws(x = fit_betabinomial_cann_s2_ps,
                                  predict_fn = logscore,
                                  prediction_names = "logscore")) %>%
  bind_draws(predictions_as_draws(x = fit_betabinomial_cann_s2_ps,
                                  predict_fn = posterior_epred,
                                  prediction_names = cellNames,
                                  newdata = predGrid,
                                  allow_new_levels = TRUE))

```

How many variables created in the function above?

```no-highlight
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

```no-highlight
variables(as_draws(post_draws_s2_ps))[c(1:10,(nVarPredDraws-(nGrid+1)):nVarPredDraws)] -> varsToPS

```

Now 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

```no-highlight
powerscale_sensitivity(x = post_draws_s2_ps,
                       variable=varsToPS) %>%
  tt() %>%
  format_tt(num_fmt="decimal") -> powerscale_sens_cann_s2_ps

powerscale_sens_cann_s2_ps

```

This is great! We see the sensitivity in the earlier model on the posterior estimates

 ![postParEstCoefs](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/5/b/5b3546883a9c71ab7b83ab78c994699f1c81c8f8.png)

but NOT in R-squared, logscore, or any of the predictions.

 ![preds](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/f/2/f263b83e78c81947f98c0ea4a1eb8a873c0dc1b8.png)

I’m new to all of this but I gather this means that the predictions generated by the model are not susceptible to the prior sensitivity evident in the posterior parameter estimates.

Now for the next step in the workflow outlined in the vignette: selectively power-scaling different priors by specifying the corresponding tag in prior\_selection. This is where I get stuck

We introduced a prior on the sds term, so let’s try that

```no-highlight
powerscale_plot_quantities(
  x = post_draws_s2_ps,
  div_measure = NULL,
  variable = c("R2",
               "logscore",                        
               "yearsFromStart0_ats_baseline0",   
               "yearsFromStart0_ats_baseline7"),
  prior_selection = "sds"
)

```

This throws the error

```no-highlight
Error in powerscale_sequence.priorsense_data(psd, lower_alpha = lower_alpha, : 
  Assertion on 'prior_selection' failed: Must be of type 'numeric' (or 'NULL'), not 'character'.

```

So it seems like the `prior_selection = "sds"` argument is the problem and the problem appears to be that it does not like strings. This cannot be though because when I run the code in the vignette, using a string, I get no such error.

Clearly I am doing something wrong.

Any help much appreciated

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [July 15, 2025, 3:01pm UTC](https://discourse.mc-stan.org/t/error-message-executing-powerscale-sensitivity-function/39885/2 "2025-07-15T15:01:52Z")

</div>

> [@llewmills](#):
>
> This throws the error
> 
> ```no-highlight
> Error in powerscale_sequence.priorsense_data(psd, lower_alpha = lower_alpha, : 
> Assertion on 'prior_selection' failed: Must be of type 'numeric' (or 'NULL'), not 'character'.
> 
> ```
> 
> So it seems like the `prior_selection = "sds"` argument is the problem and the problem appears to be that it does not like strings. This cannot be though because when I run the code in the vignette, using a string, I get no such error.

Tagging the priorsense maintainer @n-kall. Hopefully he can explain the error.

---

<div class="post-metadata">

**Author:** ![n-kall](https://avatars.discourse-cdn.com/v4/letter/n/e5b9ba/32.png) [@n-kall](https://discourse.mc-stan.org/u/n-kall)\
**Post date:** [July 15, 2025, 6:36pm UTC](https://discourse.mc-stan.org/t/error-message-executing-powerscale-sensitivity-function/39885/3 "2025-07-15T18:36:05Z")

</div>

Thanks @llewmills for reporting this and for a nice write-up of your workflow. I’m not able to reproduce the error though. It looks like an error that would occur in an older version of priorsense which was checking that `prior_selection` was numeric, but that changed in 1.1.0. Can you post the package version as shown in `sessionInfo()`

Below is the plot I get when running your code with the latest priorsense:

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/8/7/87ccbd19e868d18716a293f56092bccdeeea98d6.png)

---

<div class="post-metadata">

**Author:** ![llewmills](https://avatars.discourse-cdn.com/v4/letter/l/4da419/32.png) [@llewmills](https://discourse.mc-stan.org/u/llewmills)\
**Post date:** [July 17, 2025, 1:19am UTC](https://discourse.mc-stan.org/t/error-message-executing-powerscale-sensitivity-function/39885/4 "2025-07-17T01:19:12Z")

</div>

Dammit I wondered if it was a old version issue. And it was. After updating the graphs work. Now I just need to figure out what it means. A subject for another post I think