Hello, @bgoodri and @Bob_Carpenter,

Again, I truly appreciate your insights. I reviewed several articles on non-centered parameterizations, and also the relevant Stan documentation, and made several changes to the code.

I was able, with a two-level hierarchical model to successfully recover mean and variance parameters of simulated data after vectorizing my code and using the non-centered parameterization as follows:

```
parameters {
vector[J] uj;
real<lower=0> sigma_item;
vector[P] family_mu;
vector[K] betai_offset;
}
transformed parameters{
vector[N] eta;
vector[K] betai;
betai = family_mu[parent] + betai_offset*sigma_item;
eta = uj[peeps]-betai[item];
}
model {
betai_offset ~ normal(0,1);
uj ~ normal(0,1);
sigma_item ~ cauchy(0,2.5);
score ~ bernoulli_logit(eta);
}
```

I am still having difficulty with the decomposed means in the more complex model, however. I did define the arrays of indices as suggested by the previous post:

```
mydat<-mydat[with(mydat,order(family)),]
multten<-unique(mydat[mydat$m10==1,c("family")])
threedig<-unique(mydat[mydat$d3==1,c("family")])
```

However, when I run the code to estimate the model with decomposed means at that second level, I get the following error with the arithmetic operations:

```
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
error in 'model3070257e244c_rpri' at line 67, column 24
65:
66: family_mu = model_mu[mm] + fammu_offset*fam_resid;
67: family_mu[multten] += char_m10*ident_x;
^
68: family_mu[threedig] += char_d3*ident_y;
PARSER EXPECTED: <expression>
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model 'rpri' due to the above error.
```

The code that I am using to run the model is here (and I think, @bgoodri, that I have successfully used the non-centered parameterizations and also defined everything I am looking to estimate?) ā any additional suggestions would be greatly appreciated. I believe I have coded the model correctly, based on my reading of the documentation, but I am sure there is something I am overlooking.

Thank you again!

```
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;
real<lower=0> sigma_item;
real<lower=0> fam_resid;
vector[K] betai_offset;
vector[P] fammu_offset;
vector[M] model_mu;
real char_m10; //fixed effects of content characteristics
real char_d3; //fixed effects of content characteristics
}
transformed parameters{
vector[N] eta;
vector[K] betai;
vector[P] family_mu;
// 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;
family_mu[multten] += char_m10;
family_mu[threedig] += char_d3;
// item difficulties parameterized as random, with parent-specific means
betai = family_mu[parent] + betai_offset*sigma_item;
//log odds of a correct probability
eta = uj[peeps]-betai[item];
}
model {
//hyperprior
betai_offset ~ normal(0,1);
fammu_offset ~ normal(0,1);
//prior on random variable theta to scale the item difficulty parameters
uj ~ normal(0,1);
//weakly informative prior on variances
sigma_item ~ cauchy(0,2.5);
fam_resid ~ Cauchy(0,2.5);
//likelihood function
score ~ bernoulli_logit(eta);
}
"
mymodel<-stan(model_code = mymodel, model_name = "rpri",
data = c("N","J","K","P","M","Q",
"A","B","X","Y",
"ident_x","ident_y",
"peeps","item","pid",
"iid","fid","mid","tid",
"parent","score","mm","tt",
"multten","threedig",
"vv","hh"),
control = list(stepsize=0.01, adapt_delta=0.99, max_treedepth=15),
iter = 2000, chains = 3, thin=1, save_warmup=TRUE,
verbose = T)
##### Check recovery
mydat<-read.csv("empiricaldataforestimation_model1_v3_params.csv",header=TRUE)
fammu<-cbind(unique(mydat[order(mydat$fid),c("fid","family_mu")]),summary(mymodel,pars="family_mu"))
fammu[,c(2,3,6,10)]
sigma_item
summary(mymodel,pars="sigma_item")
```