Random intercept and slopes with multilevel generalized linear model

Modeling
1

Hi all,

I am trying to construct a multilevel model with random intercepts and slopes, but just can’t figure out how to add the slopes. I know this has been covered before, but it doesn’t quite sink in. Below is a toy example with corresponding STAN code for a random intercepts model. The code is adopted from here. I can recover the parameters, but not sure how to add in varying slopes. Any help would be appreciated.

library(lme4)
library(rstan)
library(rstanarm)

### get data
data(jsp, package="faraway")
jspr <- jsp[jsp$year==2,]
jspr$school = droplevels(jspr$school)
jspr$craven <- jspr$raven-mean(jspr$raven)

jspr$math_class = 0
jspr$math_class[jspr$math > mean(jspr$math,na.rm=T)] = 1

### parameters for stan
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
set.seed(123)

## MCMC settings
ni <- 2000
nt <- 1
nb <- 1000
nc <- 4

mod_randInt.glmm = stan_glmer(math_class ~ craven + (1 | school/class),family=binomial,data=jspr)
print(mod_randInt.glmm,pars=c("beta","sigmalev1","sigmalev2","sigmaeps"))

### now with Stan
Xmat = model.matrix(glmer(math_class ~ craven + (1 | school/class),family=binomial,data=jspr))
Nobs = nrow(Xmat)
Npreds=ncol(Xmat)
Nlev1=length(unique(jspr$school))
Nlev2=length(unique(jspr$class))
y=jspr$math_class
levind1=as.numeric(jspr$school)
levind2=as.numeric(jspr$class)

data.list = list(Nobs = Nobs, Npreds = Npreds, Nlev1 = Nlev1, Nlev2 = Nlev2, y = y, x = Xmat,levind1 = levind1,
                 levind2 = levind2)

multilevel_glm.stan = '
data {
  int<lower=0> Nobs;
  int<lower=0> Npreds;
  int<lower=0> Nlev1;
  int<lower=0> Nlev2;
  int y[Nobs];
  matrix[Nobs,Npreds] x;
  int<lower=1,upper=Nlev1> levind1[Nobs];
  int<lower=1,upper=Nlev2> levind2[Nobs];
}

parameters {
  vector[Npreds] beta;
  real<lower=0> sigmalev1;
  real<lower=0> sigmalev2;
  real<lower=0> sigmaeps;
  
  vector[Nlev1] eta1;
  vector[Nlev2] eta2;
}
transformed parameters {
  vector[Nlev1] ran1;
  vector[Nlev2] ran2;
  vector[Nobs] yhat;
  
  ran1  = sigmalev1 * eta1;
  ran2  = sigmalev2 * eta2;
  
  for (i in 1:Nobs)
    yhat[i] = x[i]*beta+ran1[levind1[i]]+ran2[levind2[i]];
  
}
model {
  eta1 ~ normal(0, 2.5);
  eta2 ~ normal(0, 2.5);
  sigmalev1 ~ cauchy(0, 2.5);
  sigmalev2 ~ cauchy(0, 2.5);
  sigmaeps ~ cauchy(0, 2.5);
  y ~ bernoulli_logit(yhat);
}
'

mod_randInt2.glmm = stan(model_code = multilevel_glm.stan,data=data.list,chains = nc,thin = nt,iter = ni, warmup = nb,
                         seed = 1)

print(mod_randInt3.glmm,pars=c("beta","sigmalev1","sigmalev2","sigmaeps"))

### Now with random intercept and slopes
mod_randslope2.glmm = stan_glmer(math_class ~ craven + (1 + craven  | school/class),family=binomial,data=jspr)
print(mod_randslope2.glmm,pars=c("beta","sigmalev1","sigmalev2","sigmaeps"))

At this point, I am just not sure how to code the random slopes. I have looked at examples, but not sure how to do it.

2 Likes
2

Okay. I have cobbled together some code based on this paper, and I at least can get it to run. I receive error messages about divergent chains etc, but it returns parameters estimates similar to stan_glmer() results.

Here is the R code with embedded STAN code:

mod_randslope2.glmm = stan_glmer(math_class ~ craven + (1 + craven  | school/class),family=binomial,data=jspr)
print(mod_randslope2.glmm)

### now with Stan
randIntRandSl.stan = '
  data {
    int<lower=0> N;               //no trials
  int<lower=1> Npreds;               //no fixefs
  int<lower=0> Nlev1;               //no items, level 1
  int<lower=1> n_u;             //no level1 random effects
  int<lower=0> Nlev2;               // no items, level 2
  int<lower=1> n_w;             //no level2 random effects
  int<lower=1,upper=Nlev1> lev1[N]; //subject indicator
  int<lower=1,upper=Nlev2> lev2[N]; //item indicator
  row_vector[Npreds] Xmat[N];           //fixef design matrix
  row_vector[n_u] Z_u[N];       //subj ranef design matrix
  row_vector[n_w] Z_w[N];       //item ranef design matrix
  int Y[N];                 // response
  }
  
  parameters {
  vector[Npreds] beta;               //fixef coefs
  cholesky_factor_corr[n_u] L_u;  //cholesky factor of subj ranef corr matrix
  cholesky_factor_corr[n_w] L_w;  //cholesky factor of item ranef corr matrix
  vector<lower=0>[n_u] sigma_u; //subj ranef std
  vector<lower=0>[n_w] sigma_w; //item ranef std
  real<lower=0> sigma_e;        //residual std
  vector[n_u] z_u[Nlev1];           //spherical subj ranef
  vector[n_w] z_w[Nlev2];           //spherical item ranef
  
  }
  
  transformed parameters {
  vector[n_u] u[Nlev1];             //level1 ranefs
  vector[n_w] w[Nlev2];             //level2 ranefs
  vector[N] yhat;
  {
    matrix[n_u,n_u] Sigma_u;    //subj ranef cov matrix
    matrix[n_w,n_w] Sigma_w;    //item ranef cov matrix
    Sigma_u = diag_pre_multiply(sigma_u,L_u);
    Sigma_w = diag_pre_multiply(sigma_w,L_w);
    for(j in 1:Nlev1)
    u[j] = Sigma_u * z_u[j];
    for(k in 1:Nlev2)
    w[k] = Sigma_w * z_w[k];
  }
  
  for (l in 1:N)
  yhat[l] = Xmat[l] * beta + Z_u[l] * u[lev1[l]] + Z_w[l] * w[lev2[l]];
  }
  
  model {
  //priors
  L_u ~ lkj_corr_cholesky(2.0);
  L_w ~ lkj_corr_cholesky(2.0);
  for (j in 1:Nlev1)
  z_u[j] ~ normal(0,1); // need to make sure these priors are correct
  for (k in 1:Nlev2)
  z_w[k] ~ normal(0,1); // need to make sure these priors are correct
  //likelihood
  Y ~ bernoulli_logit(yhat);
  }

'

Xmat = model.matrix(glmer(math_class ~ craven + (1+ craven | school/class),family=binomial,data=jspr))
N = nrow(Xmat)
Npreds=ncol(Xmat)
Nlev1=length(unique(jspr$school))
n_u = ncol(Xmat)
Nlev2=length(unique(jspr$class))
n_w = ncol(Xmat)
Z_u = Xmat
Z_w = Xmat
Y=jspr$math_class
lev1=as.numeric(jspr$school)
lev2=as.numeric(jspr$class)

data.list = list(N = N, Npreds = Npreds, Nlev1 = Nlev1, n_u = n_u,
  Nlev2 = Nlev2,n_w = n_w,Xmat = Xmat, Z_u = Z_u, Z_w = Z_w, Y = Y, lev1=lev1,lev2=lev2)

mod_randslope3.glmm = stan("randIntRandSlope.stan",data=data.list,chains = nc,thin = nt,iter = ni, warmup = nb,
                           seed = 1, control = list(adapt_delta = 0.80))

print(mod_randslope3.glmm,pars=c("beta"))

Just not sure how to solve these problems. It probably has to do with the priors, which I just took a guess because I wasn’t familiar with this way of parameterizing the model. Any help would be appreciated.

1 Like
3

I’m working on a similar problem and just looked at your code @Wade. Did you solve this already? It looks like you are assuming that there is a correlation among the subject and item effects. Is this something that is necessary for your model? It makes it a lot more complex and could be the reason that you are getting divergent chains. Listing 6 in the paper your referenced might work better for you.