Multi-level and cross-classified model... and fixed effects are being pulled toward zero?

[Shauna Sweet]

Hello all,

I am new to Stan and also very new to Bayesian modeling. I think that I have successfully translated a multi-level model into the Bayesian estimation framework (though I'm totally open to suggestions/comments/helpful tips!) - but when I run the model I notice that all of my fixed effects seem to be "pulled" toward zero, and the residual variance at level one gets bizarrely large?

Might there be something with how I'm specifying the structure of the model in Stan, or how I'm specifying my priors, that could be leading to this result?

My model code is attached below.

Thanks so much!
Shauna

[Shauna Sweet]

I wanted to follow up on my own message and say that the problems that I have been getting -- which is unusually large variance parameters is absolutely consistent with my (likely improper! eek!) hyperpriors.

Looking again through the literature, I think a half-cauchy distribution would be more appropriate, rather than an inverse gamma (and certainly not the one in the code...), because the number of groups that I have at each level are relatively small - 3 templates and 4 "parents" within each template.

If anyone sees anything else in the code, particularly in how I've attempted to specify regressions at each level (see the transformed parameters block), I would really appreciate it.

Thanks!

[Bob Carpenter]

Why do you restrict this parameter to be between 0 and 1?

vector<lower=0,upper=1>[4] kc[P]; //key content characteristics

This isn't right unless mid is a permutation

for (m in 1:M){
mum[mid[m]] <- mut[itemplate[mid[m]]] + bd*display[mid[m]];
}

The same problem's repeated throughout. So this might give mum two
priors:

for (m in 1:M){
mum[mid[m]] ~ normal(mut[itemplate[mid[m]]],sigmam[mid[m]]);
}

And please, please, indent code so it's easy to line up open and
close brackets.

You can also vectorize.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Shauna Sweet]

Bob,

Thank you so much for your reply.

Why do you restrict this parameter to be between 0 and 1?

  vector<lower=0,upper=1>[4] kc[P];                //key content characteristics

kc is a vector of zeros and ones, essentially a set of dummy/design variables which activate effects that are captured by bkc, in the parameters block.

  for (m in 1:M){    
  mum[mid[m]] <- mut[itemplate[mid[m]]] + bd*display[mid[m]];
  }

The same problem's repeated throughout.  So this might give mum two
priors:

for (m in 1:M){
mum[mid[m]] ~ normal(mut[itemplate[mid[m]]],sigmam[mid[m]]);
}

Regarding the two priors, I will look back at my code, because what I'm trying to specify are regressions at multiple levels, and I'm not sure how/where to specify the residuals. So for example, the mean difficulty of item models (mum) is decomposed into the difficulty of the template which is used to generate that model (which has a mean difficulty mut), and a set of effects (bd). sigmam is supposed to be the residual variance, so what is not explained by the fixed effects... Does this make sense?

Thank you again. 

[Shauna Sweet]

Might the solution be (in part) to remove the parameters from the transformed block, declaring all of the parameters in the parameter block, and writing the following in the modeling block:

model {

   for (m in 1:M){

      mum[mid[m]] ~ normal([itemplate[mid[m]]] + bd*display[mid[m]],sigmam[mid[m]]);
   }

   for (p in 1:P){
      mup[fid[p]] ~ normal(mum[imodel[fid[p]]] + bkc*kc[fid[p]],sigmap[fid[p]]);
   }

   for (k in 1:K){
      betai[iid[k]] ~ normal(mup[parent[iid[k]]],sigmai);
   }

   for (j in 1:J){
      uj[j] ~ normal(0,1);
   }

// response model

   for (n in 1:N){                                           
      score[n] ~ bernoulli_logit(uj[pid[peeps[n]]]-betai[iid[item[n]]]);
   }


Also, I realize I didn't answer your question: I have a data design where I need to numerically index items and the components used to generate those items, and in some cases the numbers are not sequential, with not all people answering all items, etc. The multiple indexes were my solution (albeit not a very elegant one!) to be able to estimate the parameters, linking between the levels.

Thanks again!
Shauna

[Bob Carpenter]

Sorry

> On Apr 26, 2016, at 9:43 PM, Shauna Sweet <shauna...@gmail.com> wrote:
>
> Bob,
>
> Thank you so much for your reply.
>
> Why do you restrict this parameter to be between 0 and 1?
>
> vector<lower=0,upper=1>[4] kc[P]; //key content characteristics
>
> kc is a vector of zeros and ones, essentially a set of dummy/design variables which activate effects that are captured by bkc, in the parameters block.

You don't want to do this this way. It's much more efficient
to declare it as an int and use a conditional rather than multiplying by 0 (drop) or
1 (include). But I'd misparsed it as a parameter, not data.

So instead of kc[p,i] * x, use if_else(kc[p,i], x, 0). But even
better, just drop the term if kc[p,i] == 0.

> for (m in 1:M){
> mum[mid[m]] <- mut[itemplate[mid[m]]] + bd*display[mid[m]];
> }
>
> The same problem's repeated throughout. So this might give mum two
> priors:
>
> for (m in 1:M){
> mum[mid[m]] ~ normal(mut[itemplate[mid[m]]],sigmam[mid[m]]);
> }
>
> Regarding the two priors, I will look back at my code, because what I'm trying to specify are regressions at multiple levels, and I'm not sure how/where to specify the residuals. So for example, the mean difficulty of item models (mum) is decomposed into the difficulty of the template which is used to generate that model (which has a mean difficulty mut), and a set of effects (bd). sigmam is supposed to be the residual variance, so what is not explained by the fixed effects... Does this make sense?

I get confused with fixed vs. random effects terminology.

The point is that you almost always want to give every
parameter a distribution (that is, all mum[k]), and you
almost never want to give a parameter two distributions.

- Bob

[Bob Carpenter]

Is mid just a permutation of 1:M? If not, the below is
almost certainly wrong. And if it is, you can just drop mid[]
from the calcualtions and just use m.

- Bob

[Shauna Sweet]

Bob,

Thank you so much!

All of these changes seem to both improve running time, and the model does run.

I think that now my question has to do with correctly parameterizing a draw of parameters when those parameters are essentially drawn from independent distributions. In this case, I have a hierarchical 2 PL IRT model, with a regression at levels two and three with factors that contribute to both alpha (or log gamma) ad beta, but I need to specify the residuals for those regressions of the fixed effects on alpha and beta as independent. I was attempting to specify a covariance matrix with 0's on the off-diagonals, but I frequently got warnings about the factorization of that covariance matrix being negative. For example, I attempted to specify a covariance matrix Sigma, and then draw from a multi normal distribution as follows:

  transformed parameters {
 
  cov_matrix[2] Sigmap[P];

 
  for (p in 1:P){   

     Sigmap[p,1,1]<-taup[p,1];
     Sigmap[p,1,2]<-nocov;
     Sigmap[p,2,1]<-nocov;
     Sigmap[p,2,2]<-taup[p,2];
    }

and then in the model block I had:

  for (p in 1:P){

  mup[p] ~ multi_normal(mum[p] + if_else(X1[p], kc1, 0) + if_else(X2[p], kc2, 0), taup);
  }

I read several posts about restricting tau to be less than 1, and this fixed the covariance matrix warnings, but I wondered if it wasn't causing other problems, namely it was inappropriately bounding tau, causing problems elsewhere... So then I attempted the following in this example, and the chains are not converging for taup. You'll see your suggested changes in the code implemented below! (I am sure it could be more elegant... :)

  for (p in 1:P){
 

    mup[p,1] ~ normal(mum[p,1] + if_else(X1[p], kc1[1], 0) + if_else(X2[p], kc2[1], 0), taup[1]);
    mup[p,2] ~ normal(mum[p,2] + if_else(X1[p], kc1[1], 0) + if_else(X2[p], kc2[2], 0), taup[2]);

  }
  
Should I return to bounding tau between 0 and 1, and returning to the multi_normal draw?

Again, I greatly appreciate your help!


Best,
Shauna