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.