Hi all,

I would be grateful for your assistance. I am trying to predict a patient’s score (on interval [0,1]) using a mixed effect beta regression model.

\text{logit}(\mu_{ij}) = \beta_0 + \beta_1 (x_1)_i + \beta_2 (x_2)_{ij}+ \beta_3 (x_2)_{ij}^2 +u_i

Y_{ij} \sim \text{Beta}(\mu_{ij}\tau, (1-\mu_{ij})\tau)

with Y_{ij} the score for patient i at time point j (in interval [0,1]), (x_1)_i is a covariate, (x_2)_{ij} is time point the j-th data point is collected for patient i, and u_i the random effect for patient i.

I have done this successfully in rjags but it’s slow - I’m hoping it can be made faster in RStan.

Currently when I run my Rstan model fitting, I find that it is taking a very long to compile and my output is incorrect with unsatisfactory (ESS). I believe all these problems are occurring because I’ve made a mistake in my code.

I’ve built the code up from a simple linear mixed effect model and this simpler case works well. Therefore I am thinking that there must be a problem with the specific link functions and “transformed parameters” section I have written.

I have attached my .stan script below. I would be very grateful for your thoughts. Thank you in advance!

``

data {

int<lower=1> N; //number of datapoints

array[N] int subj; //subject id

vector[N] week; //week id

vector[N] week2; //week^2 id

vector[N] x1; //covariate 1 id

vector[N] score; //score at datapoint N

}

parameters {

real beta0;

real beta1;

real beta2;

real beta3;

real<lower=0> tau;

vector[I] u;

real<lower=0> sigma_u;

}

transformed parameters{

vector<lower=0,upper=1>[N] mu; // transformed linear predictor for mean of beta distribution

vector<lower=0>[N] A; // parameter for beta distn

vector<lower=0>[N] B; // parameter for beta distn

for (i in 1:N) {

mu[i] = inv_logit(beta0 + beta1*x1[i] + beta2*week[i] + beta3*(week2[i])^2 + u[subj[i]]);

}

A = mu .* tau;

B = (1.0 - mu) .* tau;

}

model {

// priors

beta0 ~ normal(0,1);

beta1 ~ normal(0,1);

beta2 ~ normal(0,1);

beta3 ~ normal(0,1);

sigma_u ~ gamma(1,1);

u ~ normal(0,sigma_u);

tau ~ gamma(1,1);

// likelihood

score ~ beta(A, B);

}

``