Hello all,
I am struggling with the expansion of my code designed to estimate a cross-classified hierarchical model. When estimating a model with homogeneous variances and fixed effects that are constant across groups (when the assumptions are too stringent to be realistic), I have code that runs in approximately 12-16 hours per replication, with four chains of 50000 (on a virtual linux machine with 6 cores). Parameter recovery is good.
I am trying to modify that code so that the variances are heterogeneous across groups, and to also allow covariates to have differential effects across groups (a model that might be more realistic).
Unfortunately, the code is exceptionally slow. After some trial and error, it appears as though the fixed effects are not an issue, but the variances are an issue. It seems as though the multiplication of an offset by a vector of values at each level (as opposed to situations where, with homogeneous variances, the residual was a single value), i.e. stan modmu_offset*mod_resid[tt] Any suggestions for making the code more efficient?
I would greatly appreciate any suggestions.
The code is below, and a trial data set is attached that I am using for parameter recoverypr_hetfe_hetre_2_respfile.csv (604.8 KB)
.
mymodel<- "
data {
int<lower=1> N; // number of observations
int<lower=1> J; // number of students
int<lower=1> K; // number of items
int<lower=1> P; // number of families
int<lower=1> M; // number of item models
int<lower=1> Q; // number of templates
int<lower=1> X; // length of multten array
int<lower=1> Y; // length of threedig array
int<lower=1> A; // length of vv array
int<lower=1> B; // length of hh array
int peeps[N]; // student giving response n
int<lower=1,upper=K> item[N]; // item for response n
int<lower=0,upper=1> score[N]; // correctness for response n
int<lower=0,upper=J> pid[J]; // Person ID number
int<lower=0,upper=K> iid[K]; // Item ID number
int<lower=0,upper=P> fid[P]; // family ID number
int<lower=0,upper=M> mid[M]; // item model ID number
int<lower=0,upper=Q> tid[Q]; // template ID number
int<lower=1,upper=P> parent[K]; //indexes items to families
int<lower=1,upper=M> mm[P]; //indexes families to item models
int<lower=1,upper=Q> tt[M]; //indexes item models to templates
int multten[X]; // Array of indices for families - numbers are some multiple of ten
int threedig[Y]; // Array of indices for families - numbers are maximum of three digits
int vv[A]; //Array of indices for imodels - display format is verbal
int hh[B]; //Array of indices for imodels - display format is horizontal
}
parameters {
vector[J] uj;
vector <lower=0> [P] sigma_item;
vector <lower=0> [M] fam_resid;
vector <lower=0> [Q] mod_resid;
real betai_offset;
real fammu_offset;
real modmu_offset;
vector[Q] template_mu;
vector[Q] disp_horiz;
vector[Q] disp_verb;
vector[M] char_m10; //fixed effects of content characteristics
vector[M] char_d3; //fixed effects of content characteristics
}
transformed parameters{
vector[N] eta;
vector[K] betai;
vector[P] family_mu;
vector[M] model_mu;
// varying item family difficulty across item families within models
// decomposition of family means into a model mean and fixed effects
model_mu = template_mu[tt] + modmu_offset*mod_resid[tt];
model_mu[vv] = model_mu[vv] + disp_verb[tt[hh]];
model_mu[hh] = model_mu[hh] + disp_horiz[tt[vv]];
// varying item family diffculty across item families within models
// decomposition of family means into a model mean and fixed effects
family_mu = model_mu[mm] + fammu_offset*fam_resid[mm];
family_mu[multten] = family_mu[multten] + char_m10[mm[multten]];
family_mu[threedig] = family_mu[threedig] + char_d3[mm[threedig]];
// item difficulties parameterized as random, with parent-specific means
betai = family_mu[parent] + betai_offset*sigma_item[parent];
//log odds of a correct probability
eta = uj[peeps]-betai[item];
}
model {
//hyperprior
betai_offset ~ normal(0,1);
fammu_offset ~ normal(0,1);
modmu_offset ~ normal(0,1);
//prior on random variable theta to scale the item difficulty parameters
uj ~ normal(0,1);
//weakly informative prior on item variance
sigma_item ~ cauchy(0,2.5);
fam_resid ~ cauchy(0,2.5);
mod_resid ~ cauchy(0,2.5);
//likelihood function
score ~ bernoulli_logit(eta);
}
"
my_stan_code <- stanc(model_code=mymodel)
my_compiled_model <- stan_model(stanc_ret = my_stan_code,verbose=FALSE)
## READ IN DATA FILE AS DD
N <- length(unique(dd$obs))
K <- length(unique(dd$item))
J <- length(unique(dd$peeps))
P <- length(unique(dd$family))
M <- length(unique(dd$imodel))
Q <- length(unique(dd$template))
peeps <- dd$peeps
item <- dd$item
score <- dd$score
### People ID - Length J
pid <- sort(unique(dd[c("peeps")])[,"peeps"])
### Item ID - Length K
iid <- sort(unique(dd[c("item")])[,"item"])
### Family ID - Length P
fid <- sort(unique(dd[c("family")])[,"family"])
### Item Model ID - Length M
mid <- sort(unique(dd[c("imodel")])[,"imodel"])
### Template ID - Length T
tid <- sort(unique(dd[c("template")])[,"template"])
### Item - Family Index - Length K, P Unique Values
### For each row in parent[], it returns the family ID number
dd<-dd[with(dd,order(item)),]
parent <- unique(dd[,c("item","family")])[,"family"]
### Family - Item Model Index - Length P, M Unique Values
### For each row in mm[], it returns the model ID number
dd<-dd[with(dd,order(family)),]
mm <- unique(dd[c("family","imodel")])[,"imodel"]
### Item Model - Template Index - Length M, T Unique Values
dd<-dd[with(dd,order(imodel)),]
tt <- unique(dd[c("imodel","template")])[,"template"]
### Content Characteristics
### For each row denoting a family, are the numbers multiples of 10 or are they three digits
### Array of indices
dd<-dd[with(dd,order(family)),]
multten<-unique(dd[dd$m10==1,c("family")])
threedig<-unique(dd[dd$d3==1,c("family")])
X<-length(multten)
Y<-length(threedig)
### Form Characteristics
### For each row denoting a template, is the form a verbal representation or is it arranged horizontally
### Base category is comprised of items displayed in numeric format, numbers arranged
### array of indices
dd<-dd[with(dd,order(imodel)),]
vv<-unique(dd[dd$verbal==1,c("imodel")])
hh<-unique(dd[dd$horiz==1,c("imodel")])
A<-length(vv)
B<-length(hh)
#########################
### Estimate Model
#########################
myfit <-sampling(my_compiled_model,
data = c("N","J","K","P","M","Q",
"A","B","X","Y",
"peeps","item","pid",
"iid","fid","mid","tid",
"parent","score","mm","tt",
"multten","threedig",
"vv","hh"),
pars = c("uj","sigma_item","fam_resid","mod_resid","template_mu","char_m10","char_d3","disp_verb","disp_horiz"),
control = list(stepsize=0.001, adapt_delta=0.999, max_treedepth=15),
iter = 50, warmup=25, chains = 1, thin=1, save_warmup=TRUE,
verbose = T)
