Vectorize or optimizing a BUGS style random effects model

Modeling
1
  1. The model in statistical language:

I am trying to model the probability of having at least one event during a trip, which can be modeled using a Poisson distribution.

P_i = P(\text{at least one event in trip i})\\ = 1 - P(\text{no event in trip i})\\ = 1 - \frac{e^{-t_i\lambda_i}(t_i\lambda_i)^0)}{0!}\\ = 1 - \text{EXP}(-t_i\lambda_i)

The \lambda_i here, as the parameter of a Poisson distribution, can be modeled using both trip and driver level variables. We assume:

\lambda_i = \exp(\beta_{0, d(i)} + \beta_{1,d(i)}x_{1i} + \beta_2x_{2i})

Where i is the trip number (also the index of each row), d(i) is the driver ID, indicating driver level random intercepts \beta_{0, d(i)} and random slopes \beta_{1,d(i)} for x_{1i}, and there is also a fixed effect for x_{2i}.

  1. My question

I regularly used BUGS to write Bayesian codes, and I am trying to change to Stan since it is faster and more efficient. I read several papers and the optimizing section of Stan manual, but I still have difficult in vectorizing or optimizing the described the model.

The primary problem for me is that there are both random and fixed effects in the model. There is no way for me to avoid the for loop since I need to index both driver level \beta_{0, d(i)}, \beta_{1,d(i)} and trip level \beta_2. Therefore, I wrote a quite BUGS style code in the model section, although I vectorized b0, b1 and their hyperpriors.

I wonder if anyone could help vectorize or optimize the Stan code to make run faster and more efficient than JAGS. Thanks a lot!

  1. Here are my Stan code:
data {
  int<lower=0> n; // # of trips/observations
  int<lower=0> k; // # of drivers

  int<lower=0> Y[n]; //binary outcome
  real<lower=0> T[n]; //trip time
  int<lower=0> dID[n]; //driver id
  real<lower=0> x1[n]; //x1 with random effects at driver level
  int<lower=0> x2[n]; //x2 with fixed effect
}
parameters{
  vector[k] b0;
  vector[k] b1;
  real b2;
  real mu0;
  real mu1;
  real<lower=0> sigma0;
  real<lower=0> sigma1;
}

model{
  for(i in 1:n){
  Y[i] ~ bernoulli(1 - exp(-T[i]* exp(beta0[dID[i]] + b1[dID[i]]*x1[i] + b2*x2[i])));
  }
  //PRIORS
  b0 ~ normal(mu0, sigma0);
  b1 ~ normal(mu1, sigma1);
  b2 ~ normal(0, 10);
  //HYPERPRIORS
  mu0 ~ normal(0, 10);
  mu1 ~ normal(0, 10);
  sigma0 ~ gamma(1, 1);
  sigma1 ~ gamma(1, 1);
}
2 
model {
  vector[n] mu;
  for (i in 1:n) {
    mu[i] = 1 - exp(-T[i] * exp(beta0[dID[i]] + b1[dID[i]]*x1[i] + b2*x2[i]));
  }
  Y ~ bernoulli(mu);
  ...
}

But you may ultimately need some non-centered parameterization. Also, if that is a cloglog inverse link function or something, it would be better to figure out (analytically) what its logarithm is and use the bernoulli_logit version of the log-likelihood for numerical stability.

1 Like
3

Thanks Ben. Is there any tutorial on how to make non-centered parameterization or transforming the inverse link function? Thanks a lot.

4
5

Thanks Ben. That help me a lot.