I have a general question regarding a varying intercept / varying slope model in jags/stan:

I have data from a psychophysics experiment, with one covariate, one within-subjects factor and several subjects:

The response variable y is binary, and I want to model the probability of giving a response as a function of the (centred) covariate x in all conditions. I think the lme4 model should be:

glm_fit <- glmer(formula = y ~ x + (1 | cond/subject),
               data = df,
               family = binomial(probit))

where the slope and intercept vary among conditions and among subjects within conditions.

The probit link function could be replaced by a logit link.

My question is: how do I correctly model the correlations between subjects' interecepts and slopes in all conditions in jags or stan?

The data are in long format, and the jags model uses nested indexing for the condition factor. The jags model is:

model {
  # likelihood
  for (n in 1:N) { # loop over N observations
    # varying intercept/slope for every condition*subject combination
    probit(theta[n]) <- alpha[cond[n], subject[n]] + beta[cond[n], subject[n]] * x[n]
    y[n] ~ dbern(theta[n])

  # priors
  for (j in 1:J) { # loop over J conditions
    for (s in 1:S) { # loop over S subjects
      # each subjects intercept/slope in each condition comes from a 
      # group-level prior
      alpha[s, j] ~ dnorm(mu_a[j], tau_a[j])
      beta[s, j] ~ dnorm(mu_b[j], tau_b[j])

   # non-informative group level priors
   mu_a[j] ~ dnorm (0, 1e-03)
   mu_b[j] ~ dnorm (0, 1e-03)
   tau_a[j] <- pow(sigma_a[j], -2)
   tau_b[j] <- pow(sigma_b[j], -2)
   sigma_a[j] ~ dunif (0, 100)
   sigma_b[j] ~ dunif (0, 100)

I have left out all correlations between the parameters on purpose. The problem I'm having is that the intercepts and slope are correlated for each subject in each condition, and the intercept/slope pairs are correlated across conditions.

Does anyone have any ideas? What is the best way to implement this?

  • $\begingroup$ How many conditions do you have? Unless you have like 10+, it probably shouldn't be a random effect (it should be a fixed effect). Are they a sample of conditions among a population of possible conditions? (Stated differently, do you want to generalize to these specific conditions, or to a larger population of conditions) $\endgroup$ – Patrick Coulombe Apr 1 '14 at 0:43
  • $\begingroup$ I wanted a general solution to the problem, but in my dataset I have either 2 or 4 conditions. I see your point when using lme4, but I'm not sure what the difference would be in a Bayesian setting, where everything is basically a random effect. I am specifically asking about how to implement this using jags or stan. $\endgroup$ – Andrew Ellis Apr 1 '14 at 7:20
  • $\begingroup$ I would agree that with a small, fixed number of conditions and a large, conceptually unlimited number of subjects, I would tend to just estimate a flat model. But I think in other circumstances you could do the hierarchical model you want with, for example, a multivariate normal prior over the alphas and betas jointly. $\endgroup$ – Ben Goodrich Apr 4 '14 at 23:20
  • $\begingroup$ For later readers: tau_a[j] <- pow(sigma_a[j], -2) is too much BUGSy. Stan has dnorm parmeterized as standard deviation, not a precison. $\endgroup$ – Dieter Menne May 28 '15 at 9:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.