# Help with understanding my results

**URL:** <https://discourse.mc-stan.org/t/help-with-understanding-my-results/7549>\
**Category:** Modeling\
**Tags:** performance, specification, fitting-issues\
**Created:** [February 1, 2019, 1:16am UTC](https://discourse.mc-stan.org/t/help-with-understanding-my-results/7549 "2019-02-01T01:16:56Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![ffwang2011](https://avatars.discourse-cdn.com/v4/letter/f/8c91f0/32.png) [@ffwang2011](https://discourse.mc-stan.org/u/ffwang2011)\
**Post date:** [February 1, 2019, 1:16am UTC](https://discourse.mc-stan.org/t/help-with-understanding-my-results/7549/1 "2019-02-01T01:16:56Z")

</div>

Please share your Stan program and accompanying data if possible.

I am fitting a multilevel 2PL IRT model to some simulated data with a multilevel 2PL model with rstan. I have tweaked the parameters for a week and find that I cannot really understand what’s going on with rstan and the results.  
My data looks like this: [20\_10\_0.3\_TwoPL\_x1.txt](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/4/45b149afc4f98a1da21fe33c30b4af5656438512.txt) (797.1 KB)

my rstan program looks like this:

* * *

data {  
int\<lower=1\> N; // number of students  
int\<lower=1\> NFam; // number of item families  
int\<lower=1\> NObs; // number of observations  
int\<lower=1,upper=N\> ID\_Stu[NObs]; // student id  
int\<lower=1,upper=NFam\> ID\_Fam[NObs]; // item family id  
int\<lower=0,upper=1\> Resp[NObs];// resp  
//real x[NObs]; // design matrix  
int NIC; //number of item clones, fixed in this syntax  
int\<lower=1,upper=NIC\> ID\_Clone[NObs];//item clone id  
}  
parameters {  
real theta[N]; // ability  
real mu[NFam]; // difficulty  
real\<lower=0, upper=5\> alpha[NFam]; // discrimination  
//real lambda[I]; // covariate  
real\<lower=0\> sigma[NFam]; // random variation  
vector[NIC] beta[NFam]; // item clone difficulties  
}  
model {  
row\_vector[NObs] alpha\_local; // temporary parameters  
row\_vector[NObs] diff\_local;  
int fam\_n;  
theta~normal(0,1); // priors  
mu~normal(0,10);  
//lambda ~normal(0,10);  
sigma~cauchy(0,.3);  
for(i in 1:NFam){  
beta[i] ~normal(mu[i], sigma[i]);  
}  
for (n in 1:NObs) {  
fam\_n \<- ID\_Fam[n];  
alpha\_local[n] \<- alpha[fam\_n];  
diff\_local[n] \<- theta[ID\_Stu[n]] - beta[fam\_n, ID\_Clone[n]];  
}  
Resp~bernoulli\_logit(alpha\_local .\* diff\_local);  
}

* * *

The rstan model fitting looks like this:

data \<- list(  
N = N, #number of students  
NFam = NFam, #number of item families  
NObs = N \* NFam \* NItemClone, #number of observations  
Resp = as.vector(t(irv)),#person by person response vector

```
ID_Stu = seq(0,length(Resp)-1)%/% (NFam*NItemClone) + 1, #student ID 

ID_Fam = rep(seq(0,(NFam*NItemClone-1)) %/% NItemClone + 1,N), #item family ID 

NIC = NItemClone, #number of item clones 
ID_Clone = rep(rep(1:NItemClone),length(Resp)/NItemClone) #item clone ID

```

)

fit1 \<- stan(  
file = paste0(path,“TwoPL.stan”), # Stan program  
data = data, # named list of data  
chains = 2, # number of Markov chains  
warmup = 200, # number of warmup iterations per chain  
iter = 1200, # total number of iterations per chain  
thin = 10,  
cores = 4, # number of cores (using 2 just for the vignette)  
refresh = 1, # show progress every ‘refresh’ iterations  
control = list(max\_treedepth = 20,adapt\_delta=0.90)  
)

sometimes I have divergent chains, so the warning message recommends increasing acceptance rate, which is already 0.90. I thought if the acceptance rate is too high, that would introduce more noises to the estimates because we are pretty much accepting whatever in the posterior. If not divergence problem, I may have chains with really bad mixing, such as the followings: [10\_10\_0.3\_TwoPL\_Mu\_Traceplot\_1.pdf](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/b/b1599d2640c5869b278c66265ada6d4c376fc728.pdf) (18.0 KB), like mu[2] and mu[7], they band is really wide, while other parameter estimates look ok. Why would that happen?

For the AC plot, some looks really bad too. Not sure why because I picked every 10th value so I thought the AC should be low for all the estimates. [20\_10\_0.3\_TwoPL\_Sigma\_ACplot\_1.pdf](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/2/209a4214e727ba9ee1f5effe1bc660a8f60979c5.pdf) (9.8 KB)

I also spent a couple days just to play with the iterations and number of chains. I did read in the literature that 10000 effective sample size is desirable for social science research, but my effective sample sizes for the 2 chains for each parameter is around 100 or something. When there are divergent chains, some may end up with 20 to 30 effective sample size, which also reflected on very high se of the estimate. This is one output I had. If my chain is 1000 steps long after burned in 200, and I thinned every 10th, with 2 chains, should I have about 200 sample sizes ideally if no AC is presented?

So I would appreciate if anybody would help me take a look at my code and give me some suggestions in terms of how to specify the number of chains, number of iterations for IRT models, step size(reflected in the max\_treedepth), and acceptance rate? Would .95 be too high? How big and small can be max\_treedepth be? My actual step is less than 10 (for example: [10\_10\_0.3\_TwoPL\_Diagnostics\_1.txt](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/c/c70fb81e576092a7b21ac465dd7da12eb2e3ec7d.txt) (1.9 KB)  
) so I specified 20. If I make stan to walk big steps, and accept the posterior with high rate like 0.95, would it mean it actually converges fast? But actually it does not. For the calibration I attached here, it took about 4 hours. :-( just learned a bit about Bayestian computation, so I have so many dump questions. Any feedback would be highly appreciated! I

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [February 2, 2019, 3:51am UTC](https://discourse.mc-stan.org/t/help-with-understanding-my-results/7549/2 "2019-02-02T03:51:26Z")

</div>

> [@ffwang2011](#):
>
> for(i in 1:NFam){  
> beta[i] ~normal(mu[i], sigma[i]);  
> }

Maybe non-center this. So do:

```no-highlight
for(i in 1:NFam) {
  beta[i] ~ mu[i] + beta_z[i] * sigma[i];
}

```

Where you’ve defined beta as a local variable (you gotta put these at the top of the model or transformed parameters block), beta\_z is now a parameter:

```no-highlight
parameters {
  vector[NIC] beta_z[NFam];
  ...
}

```

With prior:

```no-highlight
model {
  for(i in 1:NFam) {
    beta_z[i] ~ normal(0, 1);
  }
}

```

This is usually something to try when you have divergences. Check here for a description of why (and another example of non-centering cause what I wrote above probably isn’t clear) [Diagnosing Biased Inference with Divergences](https://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html) .

> [@ffwang2011](#):
>
> I thought if the acceptance rate is too high, that would introduce more noises to the estimates because we are pretty much accepting whatever in the posterior

That is for a Metropolis sampler. HMC samplers (which Stan uses a variant of) are not like this. It’s perfectly normal to operate them near 100% accept rates. I’m not sure the current Stan sampler actually does rejections haha. The short story is adapt delta=0.9999 is fine if that gets rid of your divergences. Divergences are the thing to worry about with HMC.

> [@ffwang2011](#):
>
> and I thinned every 10th

You shouldn’t need to thin your Stan samples. Use the Neff estimates to gauge the sampler efficiency and for any standard error calculations you need to do.

> [@ffwang2011](#):
>
> How big and small can be max\_treedepth be?

If you need to make it larger than the default 10, there’s something interesting going on. Hopefully you don’t have to. That max\_treedepth is an adaptation limit. You shouldn’t need to make it lower. Ideally if the sampler doesn’t need to use the whole tree, it automatically won’t.

> [@ffwang2011](#):
>
> For the calibration I attached here, it took about 4 hours.

After non-centering, I’d replace the Cauchy prior on sigma and with a normal and see if that helps. Cauchy priors are super non-informative and that can make the sampling really slow.

---

<div class="post-metadata">

**Author:** ![ffwang2011](https://avatars.discourse-cdn.com/v4/letter/f/8c91f0/32.png) [@ffwang2011](https://discourse.mc-stan.org/u/ffwang2011)\
**Post date:** [February 4, 2019, 2:03am UTC](https://discourse.mc-stan.org/t/help-with-understanding-my-results/7549/3 "2019-02-04T02:03:36Z")

</div>

> [@bbbales2](#):
>
> p

Thank you for your response. I will try to get convergence first and then play with the priors.

---

<div class="post-metadata">

**Author:** ![ffwang2011](https://avatars.discourse-cdn.com/v4/letter/f/8c91f0/32.png) [@ffwang2011](https://discourse.mc-stan.org/u/ffwang2011)\
**Post date:** [February 8, 2019, 2:24am UTC](https://discourse.mc-stan.org/t/help-with-understanding-my-results/7549/4 "2019-02-08T02:24:29Z")

</div>

Hi bbbales2, after delving into the non-convergence issue I had, I learned that stan is not efficient in estimating hierarchical models. So I tried non-centering of one of the variable I am estimating, just like what you suggested. It converged and relatively efficient. I wonder if you could help me take a look at my code and my results because I do see some catepillar like traceplot, which I do not understand what that means.

My stan model goes as follows:  
data {  
int\<lower=1\> N; // number of students  
int\<lower=1\> NFam; // number of item families  
int\<lower=1\> NObs; // number of observations  
int\<lower=1,upper=N\> ID\_Stu[NObs]; // student id  
int\<lower=1,upper=NFam\> ID\_Fam[NObs]; // item family id  
int\<lower=0,upper=1\> Resp[NObs];// resp  
int NIC; //number of item clones, fixed in this syntax  
int\<lower=1,upper=NIC\> ID\_Clone[NObs];//item clone id  
}  
parameters {  
real theta[N]; // ability  
real\<lower=0, upper=5\> alpha[NFam]; // discrimination  
real mu[NFam]; // difficulty  
real\<lower=0\> sigma[NFam]; // random variation  
vector[NIC] beta\_z[NFam]; //\*beta reparameterization;  
}

transformed parameters {  
vector[NIC] beta[NFam]; // item clone difficulties  
for(i in 1:NFam) {  
beta[i]= mu[i] + beta\_z[i] \* sigma[i];  
}  
}

model {  
row\_vector[NObs] alpha\_local; // temporary parameters  
row\_vector[NObs] diff\_local;  
int fam\_n;  
theta~normal(0,1); // priors  
mu~normal(0,10);  
sigma~normal(0,1);  
for(i in 1:NFam){  
beta\_z[i] ~ normal(0, 1);  
//beta[i] ~normal(mu[i], sigma[i]);  
}  
for (n in 1:NObs) {  
fam\_n \<- ID\_Fam[n];  
alpha\_local[n] \<- alpha[fam\_n];  
diff\_local[n] \<- theta[ID\_Stu[n]] - beta[fam\_n, ID\_Clone[n]];  
}  
Resp~bernoulli\_logit(alpha\_local .\* diff\_local);  
}

The model I fitted in rstan goes as follows:  
fit1 \<- stan(  
file = paste0(path,“TwoPL1.stan”), # Stan program  
data = data, # named list of data  
seed=4562627,  
chains = 4, # number of Markov chains  
warmup = 2500, # number of warmup iterations per chain  
iter = 5000, # total number of iterations per chain  
#thin = 10,  
cores = 4, # number of cores (using 2 just for the vignette)  
refresh = 100, # show progress every ‘refresh’ iterations  
control = list(max\_treedepth = 10,adapt\_delta=0.80)  
)

The concerning trace plot looks like as follows: [10\_10\_0.15\_TwoPL\_Mu\_Traceplot\_1.pdf](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/a/ae54051917d7ff54a7448aa80d5fed3d91fd4122.pdf) (820.1 KB)  
[10\_10\_0.15\_TwoPL\_Sigma\_Traceplot\_1.pdf](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/3/32f9f1fa369015ec4246b6bde83d7a8dd9c59abf.pdf) (868.0 KB)

The same problem of mu[4] is reflected in the Neff too as it is the lowest (partial results: [10\_10\_0.15\_TwoPL\_Par\_Ests\_1.pdf](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/a/a4c4f1de0989d06f279675c20f08073a45cf03cc.pdf) (85.4 KB))  
I ran 4 chains of 2500 iterations, so it should be around 10,000, but it is only around 1000 for this parameter estimate. Why some parameter estimates are so unstable? I can understand partly because the distribution of items and people are not exactly overlapped so when items are extremely low in terms of alpha, there are not many information(reflected in people) down there to estimate these item parameters. What do you think? Should I be concerned about the unstable estimates after I get a converged solution?

I also have a question about the number of iterations ( relately, Neff). How large the Neff should be to get stable parameter estimates? I did read that it differs from field to field. Is there a rough rule of thumb?

I am running a simulation study. In my field, we always run multiple repetitions and then aggregate the results. How many repetitions would you recommend based on the rstan model I fit? As I already did 4 chains of 2500 iterations, I think this is also repetitions as well. I guess I am all thinking about saving time but getting the most reliable results. :-) Thank you for your attention.

---

<div class="post-metadata">

**Author:** ![bbbales2](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bbbales2/32/77_2.png) [@bbbales2](https://discourse.mc-stan.org/u/bbbales2)\
**Post date:** [February 8, 2019, 2:55am UTC](https://discourse.mc-stan.org/t/help-with-understanding-my-results/7549/5 "2019-02-08T02:55:40Z")

</div>

> [@ffwang2011](#):
>
> but it is only around 1000

Neff of 1000 should be pretty good.

Check the MC\_error column there. That’s an estimate of the Monte Carlo error: [Monte Carlo error analysis](https://web.northeastern.edu/afeiguin/phys5870/phys5870/node71.html).

The posterior mean of mu[4] is -2.8. The MC\_error estimate is like 0.05. So do you need to be able to identify the mean of mu[4] more accurately than within a range of like -2.9 to -2.7? If so, you’ll need more samples. If not, you’re done.

We’re usually interested in more complicated things than the mean of one parameter. If you want to avoid going crazy checking Monte Carlo diagnostics on your posterior samples, work on checking your model calibration with posterior predictives and such once the Rhat/Neff diagnostics seem okay.

> [@ffwang2011](#):
>
> I am running a simulation study. In my field, we always run multiple repetitions and then aggregate the results. How many repetitions would you recommend based on the rstan model I fit?

I’ll worm out of this and say it is up to you to convince yourself and your audience :D.

The generic ways to check a model are:

1. Make sure your prior predictions are sane
2. Convince yourself that you can recover parameters with simulated data
3. Make sure your diagnostics are good (4 chains, Rhats = 1.0, no divergences, enough Neffs)
4. Check the calibration of the model with posterior predictions

So if you’ve done all those things and you’re happy with your model, awesome, make a nice presentation and convince other folks of your stuff.

---

<div class="post-metadata">

**Author:** ![ffwang2011](https://avatars.discourse-cdn.com/v4/letter/f/8c91f0/32.png) [@ffwang2011](https://discourse.mc-stan.org/u/ffwang2011)\
**Post date:** [February 8, 2019, 4:04am UTC](https://discourse.mc-stan.org/t/help-with-understanding-my-results/7549/6 "2019-02-08T04:04:32Z")

</div>

Thank you very much, Ben. All the comments are very helpful. :-)