I’m trying to model the following toy dataset as part of my self-study:

running.csv (4.2 KB)

You can read this data in using this function:

```
get_normal_hierarchical_predictor_data <- function(integer_encode_groups=F) {
df <- read.csv("./running.csv")
if (integer_encode_groups) {
df$runner <- factor(df$runner)
df$runner <- as.numeric(df$runner)
}
df
}
df <- get_normal_hierarchical_predictor_data()
df$age <- (df$age - mean(df$age)) / sd(df$age) # normalize predictor
```

In RStan, I defined the following hierarchical models:

```
modelString = "
data {
// parameters for specifying prior for b0
real m_0;
real<lower=0> s_0;
// parameters for specifying prior for sigma_0
real<lower=0> lambda_0;
// parameters for specifying prior for b1
real m_1;
real<lower=0> s_1;
// parameters for specifying prior for sigma_1
real<lower=0> lambda_1;
// parameters for specifying prior for sigma_y
real<lower=0> lambda_y;
// parameters for data
int<lower=0> n_obs;
int<lower=0> n_group;
vector[n_obs] x;
real y[n_obs];
int group[n_obs];
int prior_only;
}
parameters {
real b0;
real<lower=0> sigma_0;
vector[n_group] b0_j;
real b1;
real<lower=0> sigma_1;
vector[n_group] b1_j;
real<lower=0> sigma_y;
}
model {
b0 ~ normal(m_0, s_0);
sigma_0 ~ exponential(lambda_0);
b0_j ~ normal(b0, sigma_0);
b1 ~ normal(m_1, s_1);
sigma_1 ~ exponential(lambda_1);
b1_j ~ normal(b1, sigma_1);
sigma_y ~ exponential(lambda_y);
if (!prior_only) {
y ~ normal(b0_j[group] + b1_j[group] .* x, sigma_y);
}
}
generated quantities {
real y_rep[n_obs] = normal_rng(b0_j[group] + b1_j[group] .* x, sigma_y);
}
"
```

Compiling and sampling are pretty fast:

```
library(rstan)
stanDso = stan_model( model_code=modelString )
dataList = list(
x = df$age,
y = df$net,
group = df$runner,
n_obs = n_obs,
n_group = n_group,
m_0 = 0, s_0 = 100,
lambda_0 = 0.1,
m_1 = 0, s_1 = 100,
lambda_1 = 0.1,
lambda_y = 0.1,
prior_only = 0
)
fit = sampling(
object=stanDso,
data=dataList,
chains=3,
warmup=2000,
iter=2000+10000
)
```

But the same (almost) model in rstanarm takes way (at least 5x time) longer to sample (for the same number of warmup and iter):

```
library(rstanarm)
fit <- stan_glmer(
net ~ 1 + age + (1 + age | runner),
data = df,
family = gaussian(),
prior_intercept = rstanarm::normal(0, 100),
prior = rstanarm::normal(0, 100),
prior_aux = rstanarm::exponential(0.1),
prior_covariance = rstanarm::decov(reg=1, conc=1, shape=1, scale=1),
chains=3,
warmup=2000,
iter=2000+10000
)
```

Does anyone know the reason behind this? Could it be because that rstanarm defines a joint prior via prior covariance, while my stan code uses independent priors?

Note: I only observe the slowdown of rstanarm for this hierarchical model; for binomial and simple linear regression models, rstanarm is consistently faster than rstan (time measured for compilation + sampling).