Fitting latent class models in two-way contingency table

[Charles]

I'm starting to get involved with stan to fit various simple models. I've done few replications for those WinBUGS examples and now I'm having some problem fitting a latent class model for two-way contingency table.

Here is what I've got so far:

############################################

 

data {

 int<lower=0> N; //132,number of patients

 int<lower=0> R; //9,number of rows(patterns from cross-classification table)

 int<lower=0> I; //2,number of items

 int<lower=0> culm[R]; //number of patients at each pattern

 vector<lower=0>[2] a1; //prior

 vector<lower=0>[3] a2; //prior

 }

 parameters {

  simplex[2] theta;

  real<lower=0> alpha1[3]; //probability of item1 responses given class1

  real<lower=0> alpha2[3]; //probability of item1 responses given class2

  real<lower=0> beta1[3];  //probability of item2 responses given class1

  real<lower=0> beta2[3];  //probability of item2 responses given class2

 

 }

 model {

 

 theta ~ dirichlet(a1);

 alpha1 ~ dirichlet(a2);

 alpha2 ~ dirichlet(a2);

 beta1 ~ dirichlet(a2);

 beta2 ~ dirichlet(a2);

 }

 for (j in 1:3) {

  for (k in 1:3) {

   culm ~ multinomial(theta[1]*alpha1[j]*beta1[k]+theta[2]*alpha2[j]*beta2[k]);

  }

 }

LC_dat <- list(N=132,R=9,I=2,a1=c(1,1),a2=c(1,1,1),culm=c(43,16,3,6,11,10,9,18,16))

----------------------------------------------------------------------------------------------------------------------------

When I tried to run the model, I received this message:

----------------------------------------------------------------------------------------------------------------------------

no matches for function name="dirichlet_log"

    arg 0 type=real[1]

    arg 1 type=vector

available function signatures for dirichlet_log:

0.  dirichlet_log(vector, vector) : real

unknown distribution=dirichlet

LOCATION OF PARSING ERROR (line = 26, position = 25):

PARSED:

model {

 

 theta ~ dirichlet(a1);

 alpha1 ~ dirichlet(a2);

EXPECTED: <eps> BUT FOUND: 

 alpha2 ~ dirichlet(a2);

 beta1 ~ dirichlet(a2);

 beta2 ~ dirichlet(a2);

 }

-------------------------------------------------------------------

I quickly went to check Dirichlet distribution from stan reference but couldn't figure out why the model failed to compile.

Any idea of what I'm missing?

Thanks,

Charles

[Daniel Lee]

Hi Charles,

You'll need to declare your parameters as vectors, not arrays. There's a chapter in the manual describing the difference. (I'm on my phone and it's hard to look it up.)

Instead of
real<lower=0> alpha1[3];
use
vector<lower=0>[3] alpha1;

On Sep 30, 2014 6:47 PM, "Charles" <chenq...@gmail.com> wrote:

I'm starting to get involved with stan to fit various simple models. I've done few replications for those WinBUGS examples and now I'm having some problem fitting a latent class model for two-way contingency table.

Here is what I've come up so far:

--
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.
For more options, visit https://groups.google.com/d/optout.

[Charles]

Thanks Daniel,

Does that mean I also have to declare counts as a vector instead of arrays? Now I'm receiving the message which says that there is some problem with multinomial function.

[Bob Carpenter]

It's helpful for us to help diagnose your problem if you
cut and paste the exact message you're getting into your
mail. Otherwise, we have to guess what the "some problem" is.

Your model has a number of bugs.

1. alpha1, alpha2, beta1, and beta2 all need to be declared
as simplexes because they're variants drawn from Dirichlet priors.
Daniel noticed they had to be vectors, but the requirement's
even stricter than that. In Stan, the constraints on the variables
need to match the support.

So that's

simplex[3] alpha1;
...

2. I have no idea what you want to do with this:

culm ~ multinomial(theta[1] * alpha1[j] * beta1[k]
+ theta[2] * alpha2[j] * beta2[k]);

but the argument to multinomial() needs to be a 9-simplex (a
size 9 vector with non-negative entries summing to 1).
As is, you get a single scalar, because all those indexed
variables are simplexes.

I also don't understand why it's a doubly-nested loop. You
wind up sampling the same data, culm, 9 times.

If you're trying to create some kind of mixture model, there's a chapter in
the manual on how to do that.

Rather than using a 2-simplex, it's easier to understand and more efficient
to just use a single variable constrained to <lower=0,upper=1>. Then you
can put a beta prior on it instead of a Dirichlet and get exactly the
same posterior for theta with:

parameters {
real<lower=0,upper=1> theta;
...
model {
theta ~ beta(a[1],a[2]);
...
culm ~ multinomial(theta*alpha1[j]*beta1[k] + (1 - theta)*alpha2[j]*beta2[k]);

The right side of this still won't compile because it's just a scalar
and you need a simplex.

- Bob

[Luc Coffeng]

I while back I posted a working toy example for contingency table sampling. May it's useful for you: https://groups.google.com/d/msg/stan-users/ec7NVGZD46A/PWLPCGvUL3EJ

[Charles]

Hi Bob,

Thanks for the help. Below is the message I received in the beginning:

-----------------------------------------------------------------------------------------------------------------------

TRANSLATING MODEL 'LC_code' FROM Stan CODE TO C++ CODE NOW.

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 

  failed to parse Stan model 'LC_code' with error message:

EXPECTATION FAILURE - DIAGNOSTIC(S) FROM PARSER:

no matches for function name="multinomial_log"

    arg 0 type=real

    arg 1 type=real

available function signatures for multinomial_log:

0.  multinomial_log(int[1], vector) : real

unknown distribution=multinomial

-----------------------------------------------------------------------------------

I use " culm ~ multinomial(theta[1] * alpha1[j] * beta1[k] + theta[2] * alpha2[j] * beta2[k])" 

to represent joint probability Pjk based on latent conditional independence model. I'm working on someone else's research

article so I just use his formula to replicate the results.

[Charles]

Thanks for sharing this. I will check it out and see if I can find any helpful information.

[Bob Carpenter]

There are errors in both your arguments. The error message is
saying you have two real arguments, but what multinomial expects
is an integer array outcome of counts and a vector parameter
whose content is a simplex (non-negative entries summing to 1).

As I said last time, what you wrote down,

theta[1] * alpha1[j] * beta1[k] + theta[2] * alpha2[j] * beta2[k],

is a scalar. And if you have a single outcome, you want categorical(),
not multinomial(), and that outcome needs to be an integer.

Check out the description of multinomial in the manual. I'm afraid
we don't have time to do a tutorial one e-mail post at a time.

If you have a working BUGS or JAGS implementation of the model, that
might help isolate your programming errors.

- Bob

[Charles]

Hi Bob,

Thanks very much for pointing out the problem. I have fixed the issue and the model compiles. 

-Charles