Mixed effect beta regression in RStan

Modeling
, ,
1

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 + beta1x1[i] + beta2week[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);
}
``

2

I have a suggestion and question:

parameters{

real<lower=0> kappa;

}
transformed parameters{ 

vector[N] mu = beta0+beta1*x1; 

}

model{

score ~ beta_proportion(inv_logit(mu),kappa);` 
kappa ~ normal(0, 50); //change as needed
}

Reference:
https://mc-stan.org/docs/functions-reference/continuous_distributions_on_0_1.html#beta-proportion-distribution

  • Can you elaborate on what the covariate “week” represents and why you’re trying to add a quadratic term to the model?
1 Like
3

Thank you so much for your response, fsdias!
My code works perfectly now.

The covariate “week” predicts a patient’s score at each week. I am running a number of simulation scenarios using this model, one being a patient’s score plateauing beyond a certain week and so I have used a quadratic term to account for this.
Many thanks once again