Graded Response Model

[Jackson King]

Hello, 

I am estimating a simple graded response model (4 items, each with four categories). I have slightly modified a model that runs in JAGS to run in STAN, but am encountering some problems. I am wondering if anybody can provide some pointers.

Below is everything, including the data from the ltm package.

I would be greatly appreciative for some help. 

Thanks

require(ltm)

data<-Science[c(1,3,4,7)]

data<-cbind(as.numeric(data[,1]),as.numeric(data[,2]), as.numeric(data[,3]), as.numeric(data[,4]) )

n<-length(data[,1])

p<-length(data[1,])

K<-apply(data,2,max, na.rm=T)

data.1 <- list(y=data,n=n,p=p,K=K)

grm_code<-"

data {

       int<lower=0> n; //Number of participants

       int<lower=0> p; //Number of items

       int<lower=0> K; //Number of categories per item

       int<lower=1,upper=4> y[n]; //Data

}

parameters {

  real theta[n]; // ability for n

  ordered[K-1] tau[p]; // difficulty for the item

  real<lower=0> alpha[p];// discrimination parameter for the item

  real<lower=0> sigma_tau; // sd of abilities

  real<lower=0> sigma_alpha; // sd of difficulties

  }

model{

real Q[n, p, 3];

real prob[n, p, 4];

theta~normal(0, 1); 

tau ~ cauchy(0,sigma_tau);

alpha~cauchy(0, sigma_alpha);

sigma_tau~cauchy(0,5);

sigma_alpha~cauchy(0,5);

 for(i in 1:n){

    for (j in 1:p){ 

        for (k in 1:(K[j]-1)){ 

        Q[i,j,k]<-inv_logit(alpha[j]*(theta[i]-tau[j,k]));

      }

  prob[i,j,1]<-1-Q[i,j,1];

         for (k in 2:(K[j]-1))

  prob[i,j,k]<-prob[i,j,k-1]-prob[i,j,k];

  prob[i,j,K[j]]<-Q[i,j,K[j]-1];

  y[i,j]~categorical(prob[i,j,1:K[j]]);

}

  }

 }

"

fit <- stan(model_code = grm_code,

            data = data.1, iter = 1000, chains = 2) 

 

[Bob Carpenter]

On 6/23/13 3:23 PM, Jackson King wrote:
> Hello,
>
>
> I am estimating a simple graded response model (4 items, each with four categories). I have slightly modified a model
> that runs in JAGS to run in STAN, but am encountering some problems. I am wondering if anybody can provide some pointers.

It usually helps if you provide an error message indicating the
kind of trouble you're encountering.

In your case, it looks like your model won't compile, because
you declare y as one-dimensional:

int<lower=1,upper=4> y[n]; //Data

then use it as two-dimensional:

y[i,j]~categorical(prob[i,j,1:K[j]]);

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

[Jackson King]

Thanks. That is obviously incorrect and is something I should have caught. However, I still get an error message:

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

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

EXPECTATION FAILURE LOCATION: file=input; line=19, column=5

tau ~ cauchy(0,sigma_tau);

    ^-- here

DIAGNOSTIC(S) FROM PARSER:

no matches for function name="cauchy_log"

    arg 0 type=vector[1]

    arg 1 type=int

    arg 2 type=real

available function signatures for cauchy_log:

0.  cauchy_log(real, real, real) : real

1.  cauchy_log(real, real, real[1]) : real

2.  cauchy_log(real, real, vector) : real

3.  cauchy_log(real, real, row vector) : real

4.  cauchy_log(real, real[1], real) : real

5.  cauchy_log(real, real[1], real[1]) : real

6.  cauchy_log(real, real[1], vector) : real

7.  cauchy_log(real, real[1], row vector) : real

8.  cauchy_log(real, vector, real) : real

9.  cauchy_log(real, vector, real[1]) : real

10.  cauchy_log(real, vector, vector) : real

11.  cauchy_log(real, vector, row vector) : real

12.  cauchy_log(real, row vector, real) : real

13.  cauchy_log(real, row vector, real[1]) : real

14.  cauchy_log(r

Below I have attached the modified code.

require(ltm)

data<-Science[c(1,3,4,7)]

data<-cbind(as.numeric(data[,1]),as.numeric(data[,2]), as.numeric(data[,3]), as.numeric(data[,4]) )

n<-length(data[,1])

p<-length(data[1,])

K<-apply(data,2,max, na.rm=T)

data.1 <- list(y=data,n=n,p=p,K=K)

grm_code<-"

data {

       int<lower=0> n; //Number of participants

       int<lower=0> p; //Number of items

       int<lower=0> K; //Number of categories per item

       int<lower=1,upper=4> y[n,p]; //Data

[Bob Carpenter]

Answer inline below.

What vector[1] is telling you is that tau is declared as an
array of vectors. The cauchy wants either vectors, row
vectors, or 1D arrays.

The following will do what you intended:

for (i in 1:p)
tau[i] ~ cauchy(0,sigma_tau);

Given that tau[i] is an ordered vector, I'm not exactly
sure what that implies for the distribution of tau because
of the order constraints.

You could provide a proper model with a prior on tau[1]
and (tau[p] - tau[1]), with the interior points uniform,
or you can put a prior on tau[1] and then all the differences.
Either way, you'd need to account for the log Jacobian implicated
by the variable transform.

I was asking Andrew about what we should be using for priors
here, but don't recall his response.

- Bob

> stan-users+...@googlegroups.com <javascript:>.
> > For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

[Jackson King]

Thanks, Bob Unfortunately, this still isn't working (setting aside tau's distribution). Now I get:

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

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

EXPECTATION FAILURE LOCATION: file=input; line=33, column=42

  y[i,j]~categorical(prob[i,j,1:K[j]])

                                         ^-- here

DIAGNOSTIC(S) FROM PARSER:

Parser expecting: "]"

Given how the Q and prob matrices are specified, I'm not sure why I am receiving the error. Any advice would be greatly appreciated. Thanks so much.

require(ltm)

data<-Science[c(1,3,4,7)]

data<-cbind(as.numeric(data[,1]),as.numeric(data[,2]), as.numeric(data[,3]), as.numeric(data[,4]) )

n<-length(data[,1])

p<-length(data[1,])

K<-apply(data,2,max, na.rm=T)

data.1 <- list(y=data,n=n,p=p,K=K)

grm_code<-"

data {

       int<lower=0> n; //Number of participants

       int<lower=0> p; //Number of items

       int<lower=0> K[p]; //Number of categories per item

       int<lower=1,upper=4> y[n,p]; //Data

}

parameters {

  real theta[n]; // ability for n

  ordered[K[p]-1] tau[p]; // difficulty for the item

  real<lower=0> alpha[p];// discrimination parameter for the item

  real<lower=0> sigma_tau; // sd of abilities

  real<lower=0> sigma_alpha; // sd of difficulties

  }

model{

real Q[n, p, 3];

real prob[n, p, 4];

theta~normal(0, 1); 

for (i in 1:p)

tau[i] ~ cauchy(0,sigma_tau);

alpha~cauchy(0, sigma_alpha);

sigma_tau~cauchy(0,5);

sigma_alpha~cauchy(0,5);

 for(i in 1:n){

    for (j in 1:p){ 

        for (k in 1: (K[j]-1)){ 

         Q[i, j, k]<-inv_logit(alpha[j]*(theta[i]-tau[j,k]));

      }

  prob[i,j,1]<-1-Q[i,j,1];

         for (k in 2: (K[j]-1))

  prob[i,j,k]<-Q[i,j,k-1]-Q[i,j,k];

  prob[i,j,K[j]]<-Q[i,j,(K[j]-1)];

[Bob Carpenter]

I'm not sure how your code got this far.

Stan isn't BUGS and uses a different syntax for
arrays and vectors and matrices. There are many examples
in the manual along with full definitions.

Are you sure the code attached led to this error? The
Stan parser should gripe about this illegal declaration:

real prob[n, p, 1:4];

Unlike R and BUGS, Stan doesn't have array slicing with
the ":" operator.

The line you cite is also illegal for the same reason:

y[i,j]~categorical(prob[i,j,1:K[j]])

You either need to use a loop or use the entire slice.

Stan has local variables, and also can have expressions
inside distributions, so it's otherwise not as restrictive
as BUGS (more like JAGS).

In general, you can debug these models with these error
messages by looking at where the error message is pointing
and checking your syntax against the definition of what's
allowed in the manual.

- Bob

On 6/24/13 5:28 PM, Jackson King wrote:
> Thanks, Bob Unfortunately, this still isn't working (setting aside tau's distribution). Now I get:
>

> Error in stanc(file = file, model_code = model_code, model_name = model_name, :
>
> failed to parse Stan model 'grm_code' with error message:
>

> EXPECTATION FAILURE LOCATION: file=input; line=33, column=42

>
>
> y[i,j]~categorical(prob[i,j,1:K[j]])
>

> ^-- here
>
>
>
> DIAGNOSTIC(S) FROM PARSER:
>

> Parser expecting: "]"
>
>
>
> Given how the Q and prob matrices are specified, I'm not sure why I am receiving the error. Any advice would be greatly
> appreciated. Thanks so much.
>

> require(ltm)
>
> data<-Science[c(1,3,4,7)]
>
> data<-cbind(as.numeric(data[,1]),as.numeric(data[,2]), as.numeric(data[,3]), as.numeric(data[,4]) )
>
> n<-length(data[,1])
>
> p<-length(data[1,])
>
> K<-apply(data,2,max, na.rm=T)
>
>
>
> data.1 <- list(y=data,n=n,p=p,K=K)
>
>
> grm_code<-"
>
> data {
>
> int<lower=0> n; //Number of participants
>
> int<lower=0> p; //Number of items
>

> int<lower=0> K[p]; //Number of categories per item

>
> int<lower=1,upper=4> y[n,p]; //Data
>
> }
>
> parameters {
>
> real theta[n]; // ability for n
>

> ordered[K[p]-1] tau[p]; // difficulty for the item

>
> real<lower=0> alpha[p];// discrimination parameter for the item
>
> real<lower=0> sigma_tau; // sd of abilities
>
> real<lower=0> sigma_alpha; // sd of difficulties
>
> }
>
> model{
>
> real Q[n, p, 3];
>

> real prob[n, p, 1:4];
>
> theta~normal(0, 1);

>
> for (i in 1:p)
>
> tau[i] ~ cauchy(0,sigma_tau);
>

> alpha~cauchy(0, sigma_alpha);
>
> sigma_tau~cauchy(0,5);
>
> sigma_alpha~cauchy(0,5);
>
> for(i in 1:n){
>
> for (j in 1:p){
>
> for (k in 1: (K[j]-1)){
>

> Q[i, j, k]<-inv_logit(alpha[j]*(theta[i]-tau[j,k]));

>
> }
>
> prob[i,j,1]<-1-Q[i,j,1];
>
> for (k in 2: (K[j]-1))
>

> prob[i,j,k]<-Q[i,j,k-1]-Q[i,j,k];
>
> prob[i,j,K[j]]<-Q[i,j,(K[j]-1)];

>
> y[i,j]~categorical(prob[i,j,1:K[j]])
>
> }
>
> }
>
> }
>
> "
>
>
> fit <- stan(model_code = grm_code,
>
> data = data.1, iter = 1000, chains = 2)
>
>
>

> On Monday, June 24, 2013 9:52:28 AM UTC-7, Bob Carpenter wrote:
>

> <https://groups.google.com/groups/opt_out> <https://groups.google.com/groups/opt_out

[Jackson King]

Thanks, Bob. I was able to get the model to work (after several modifications).  It appears to produce reasonable results that are comparable to a GRM in the "ltm" package. That said, I receive an "informational message" now that is not clear to me. I have provided start values -- but with no luck. The message remains. I assume it has something to do with how I've defined the priors for tau, but when I use the default priors, I get the same thing.

Below is the code and the message. Many thanks in advance. I very much appreciate it.

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

COMPILING THE C++ CODE FOR MODEL 'grm_code' NOW.

Iteration:  1 / 100 [  1%]  (Adapting)

Informational Message: The parameter state is about to be Metropolis rejected due to the following underlying, non-fatal (really) issue (and please ignore that what comes next might say 'error'): Error in function stan::prob::uniform_log(N4stan5agrad3varE): Upper bound parameter is -106.83560414904788:0, but must be greater than -100

If the problem persists across multiple draws, you might have a problem with an initial state or a gradient somewhere.

 If the problem does not persist, the resulting samples will still be drawn from the posterior.

Informational Message: The parameter state is about to be Metropolis rejected due to the following underlying, non-fatal (really) issue (and please ignore that what comes next might say 'error'): Error in function stan::prob::ordered_logistic(N4stan5agrad3varE): Cut points parameter is -24.47121125268054:0, but must be greater than -24.4712:0

If the problem persists across multiple draws, you might have a problem with an initial state or a gradient somewhere.

 If the problem does not persist, the resulting samples will still be drawn from the posterior.

*******

require(rstan)

require(ltm)

data<-Science[c(1,3,4,7)]

data<-cbind(as.numeric(data[,1]),as.numeric(data[,2]), as.numeric(data[,3]), as.numeric(data[,4]) )

c<-4 #Number of categories per item (here, 4)

J<-length(data[,1])

K<-length(data[1,])

N<-J*K

jj<-rep(c(1:J), times=K) #Index for persons

kk<-rep(c(1:K), each=J) #Index for items

y<-c(data[,1], data[,2], data[,3], data[,4])

dat<-list(N=N,

  c=c,

  J=J,

  K=K,

  jj=jj,

  kk=kk,

  y=y

  )

grm_code<-"

data {

       int<lower=0> N; //Number of observations

       int<lower=1> c; //Number of categories

       int<lower=1> J; //Number of ids

       int<lower=1> K; //Number of items

       int<lower=1, upper=J> jj[N]; //Identifier

       int<lower=1, upper=J> kk[N]; //question

       int y[N];  //Data

      }

parameters {

  real theta[N];

  ordered[c-1] tau[K]; //Difficulty

 real<lower=0> alpha[K];// discrimination parameter for each item

  }

model{

 alpha~cauchy(0, 5);

 theta~normal(0,1);

 for (p in 1:4){

tau[p, 1]~uniform(-100, tau[p, 2]);

tau[p, 2]~uniform(tau[p, 1], tau[p, 3]);

tau[p, 3]~uniform(tau[p, 2], 100);

}

for (n in 1:N)

y[n]~ordered_logistic(alpha[kk[n]]*theta[jj[n]], tau[kk[n]]);

}

"

fit <- stan(model_code = grm_code,

            data = dat, iter = 100, chains =3)

[Bob Carpenter]

We put that message in for diagnostic reasons and it's been confusing
users ever since.

What is happening is that the leapfrog integrator used for simulating the
Hamiltonian dynamics (path of parameters over time under random initial momentum and
negative log probability potential energy) is based on a discrete time interval
and low-order approximation to the true curvature of the log prob function; when the time
interval (i.e., step size) is too large, a step in the direction of the
momentum can take you out of the support for a function.

The message usually arises because of overflow or underflow, which can
help diagnosing a poorly specified model (not an error per se, but perhaps
a poor choice of priors or variable constraints).

In your case, the first message may indicate that the parameter in question
doesn't have the appropriate constraints defined, but not necessarily. With
an ordered variable, a big enough random momentum step can destroy the
ordering.

What should happen is that adaptation should find a smaller time interval (step
size) and the messages should go away, or at least become relatively infrequent.
It can actually be optimal for them to keep showing up because we want to
balance larger time intervals for computational efficiency with Metropolis
rejection for sampling efficiency.

- Bob

[Ben Goodrich]

On Wednesday, July 10, 2013 11:48:20 AM UTC-4, Bob Carpenter wrote:

We put that message in for diagnostic reasons and it's been confusing
users ever since.

Maybe it should be off by default during warmup? Always on after warmup is reasonable, but the message isn't even that meaningful on iteration 1.

Ben

[Bob Carpenter]

It will diagnose the problem of too few constraints
being declared, which is one of the main reasons it's
there.

Can we try again to make the message clearer? It's so
tied up with HMC I haven't been able to formulate it
clearly. But it should probably say something about

(a) checking that a parameter's constraint match its sampling
distribution's support, and

(b) strength of priors if there is underflow or overflow at boundaries.

Or we can just get rid of it altogether, which may be easiest.
It's done nothing but confuse users from the get-go.

- Bob

[Ben Goodrich]

On Wednesday, July 10, 2013 12:45:40 PM UTC-4, Bob Carpenter wrote:

On 7/10/13 12:19 PM, Ben Goodrich wrote:
> On Wednesday, July 10, 2013 11:48:20 AM UTC-4, Bob Carpenter wrote:
>
>     We put that message in for diagnostic reasons and it's been confusing
>     users ever since.
>
>
> Maybe it should be off by default during warmup? Always on after warmup is reasonable, but the message isn't even that
> meaningful on iteration 1.

It will diagnose the problem of too few constraints
being declared, which is one of the main reasons it's
there.

If there is a problem with lack of constraints, then the message will presumably appear when it ends the warmup. I guess I am not too concerned about how quickly the message gets to the user when they are not specifying the model correctly. I really don't know how to make the message any clearer. Maybe it would be better to say "A parameter proposal is about to get rejected ..." rather than "The parameter state is about the get rejected ...".

Ben

[Jackson King]

Thanks for the feedback. This makes sense. The warning message consistently appears on the first iteration, and only then, despite the priors or starting values I specify. Since it is a relatively simple model which appears to be correctly specified -- though correct me if I'm wrong -- it sounds like I can ignore the message.

[Bob Carpenter]

On 7/10/13 10:24 PM, Jackson King wrote:
> Thanks for the feedback. This makes sense. The warning message consistently appears on the first iteration, and only
> then, despite the priors or starting values I specify. Since it is a relatively simple model which appears to be
> correctly specified -- though correct me if I'm wrong -- it sounds like I can ignore the message.

Yes. It's just reporting on behavior before
adaptation has found a small enough step size.

- Bob

[louis...@gmail.com]

Hi List,

I was working with examples I found here and I believe that the code posted by Jackson King here contains a mistake. In particular 

he declares  real theta[N]; which is incorrect and he intends to have real theta[J];. 

Not a biggie but if someone else reads this thread and wants to build from what can be found here it might convenient to have at least this sorted out.

I would appreciate if anyone could comment whether this code is otherwise correct, but as far as I understand I believe this is the case. I compared with the results from grm in the ltm package and the comparison looks great.

Kind regards, Louis

library(ltm)

library(rstan)

real theta[J];

ordered[c-1] tau[K]; //Difficulty

real<lower=0> alpha[K];// discrimination parameter for each item

}

model{

alpha~cauchy(0, 5);

theta~normal(0,1);

for (p in 1:4){

tau[p, 1]~uniform(-100, tau[p, 2]);

tau[p, 2]~uniform(tau[p, 1], tau[p, 3]);

tau[p, 3]~uniform(tau[p, 2], 100);

}

for (n in 1:N)

y[n]~ordered_logistic(alpha[kk[n]]*theta[jj[n]], tau[kk[n]]);

}

"

fit <- stan(model_code = grm_code,

            data = dat, iter = 300, warmup=100, chains =3)

####Some Code to extract abilities from Stan###

s <- extract(fit, inc_warmup=FALSE, permuted=FALSE)

s1 <- s[,1,] #order of parameters as in parameter declaration; ends with _lp

s2 <- s[,2,]

s3 <- s[,3,]

stot <- rbind(s1,s2,s3)

summaryFunc <- function(x, q = c(0.025, 0.25, 0.5, 0.75, 0.975)){

  quantile(x,q) 

}

stot.summary.standard <- t(apply(stot,2,summaryFunc))

abilities.stan <- stot.summary.standard[1:392,3]

####Use the ltm package####

grm1 <- grm(data, IRT.param=FALSE)

scoresgrm1 <- factor.scores(grm1, resp.patterns=data)[[1]]

abilities <- scoresgrm1[,7]

####Compare#####

plot(abilities, abilities.stan)

[Bob Carpenter]

No, that code's still problematic.

With

parameters {

ordered[c-1] tau[K]; //Difficulty

it's not consistent to have:

model {

for (p in 1:4) {

tau[p, 1] ~ uniform(-100, tau[p, 2]);
tau[p, 2] ~ uniform(tau[p, 1], tau[p, 3]);

tau[p, 3] ~ uniform(tau[p, 2], 100);

because the support in the model is less than that declared
on the parameters. The declaration allows tau[k,1] to take on
any real value, so it can be less than -100; similarly
tau[k,3] may be greater than 100.

You can use improper priors or provide proper priors, but
the support of the prior should be equal to that of the declared
constraints. The prior support in the model can be greater than
the declared constraint if the truncation is included to match
the declared constraint; the truncation may be dropped if
the distribution parameters and the truncation point are constant.

The model clearly has a multiplicative non-identifiability in
alpha and theta, with linear predictor alpha[kk[n]] * theta[jj[n]].

Taking theta ~ normal(0,1) should fix the scale for theta and
identify the model. There are other ways to identify such models,
such as fixing components --- I described some of them in the manual
chapter on problematic posteriors and mitigating non-identifiabilities.

- 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.

> For more options, visit https://groups.google.com/d/optout.

[louis...@gmail.com]

Dear Bob,

Thank you for the feedback. I agree on the identifiability issue. I knew that but must have missed this and somehow my quick check did not reveal this issue.

I am however a bit puzzled by your comment on the support. Could you clarify the following?

-I thought it is ok to put a restriction on the declaration and have an unrestricted prior

e.g.

declare real<lower=0> alpha[K]

and later on 

alpha~cauchy(0, 2);

but the opposite is not ok? 

No restriction in the declaration but specifying a prior with some restriction e.g. uniform(a,b)  ?

-What would you suggest in the example of above? 

(how would you restrict the ordered declaration if that is the suggestion)

Kind regards, Louis (meanwhile I'll turn to manual once again)

[Daniel Lee]

On Thursday, November 27, 2014, <louis...@gmail.com> wrote:

Dear Bob,

Thank you for the feedback. I agree on the identifiability issue. I knew that but must have missed this and somehow my quick check did not reveal this issue.

I am however a bit puzzled by your comment on the support. Could you clarify the following?

-I thought it is ok to put a restriction on the declaration and have an unrestricted prior

e.g.

declare real<lower=0> alpha[K]

and later on 

alpha~cauchy(0, 2);

but the opposite is not ok? 

No restriction in the declaration but specifying a prior with some restriction e.g. uniform(a,b)  ?

That is correct. If the parameters end up outside of the support of the uniform distribution, the log probability is -inf and there is no gradient. So, NUTS / HMC stops working efficiently. There is a chance it stumbles back into an area of support, but I think our implementation won't.

By default, when Stan initializes its parameters, it chooses random values on the unconstrained scale. If the initial parameters evaluate to a -inf log probability calculation, Stan will choose a different random starting point. It does this a max of 100 times then quits if it can't find a reasonable starting value.

If you're insistent on using that sort of construction, you could start sampling by providing initial values, but it would still be very inefficient. Stan does a change of variables under the hood, which make sampling much easier over bounded variables. There's no possible chance of stepping into a place with no support.

 

-What would you suggest in the example of above? 

(how would you restrict the ordered declaration if that is the suggestion)

You might be able to do a change of variables yourself and add the right Jacobian term, maybe use a positive ordered, then subtract 100?

[Daniel Lee]

On Thursday, November 27, 2014, Daniel Lee <bea...@alum.mit.edu> wrote:

On Thursday, November 27, 2014, <louis...@gmail.com> wrote:

Dear Bob,

Thank you for the feedback. I agree on the identifiability issue. I knew that but must have missed this and somehow my quick check did not reveal this issue.

I am however a bit puzzled by your comment on the support. Could you clarify the following?

-I thought it is ok to put a restriction on the declaration and have an unrestricted prior

e.g.

declare real<lower=0> alpha[K]

and later on 

alpha~cauchy(0, 2);

but the opposite is not ok? 

No restriction in the declaration but specifying a prior with some restriction e.g. uniform(a,b)  ?

That is correct. If the parameters end up outside of the support of the uniform distribution, the log probability is -inf and there is no gradient. So, NUTS / HMC stops working efficiently. There is a chance it stumbles back into an area of support, but I think our implementation won't.

By default, when Stan initializes its parameters, it chooses random values on the unconstrained scale. If the initial parameters evaluate to a -inf log probability calculation, Stan will choose a different random starting point. It does this a max of 100 times then quits if it can't find a reasonable starting value.

If you're insistent on using that sort of construction, you could start sampling by providing initial values, but it would still be very inefficient. Stan does a change of variables under the hood, which make sampling much easier over bounded variables. There's no possible chance of stepping into a place with no support.

 

-What would you suggest in the example of above? 

(how would you restrict the ordered declaration if that is the suggestion)

You might be able to do a change of variables yourself and add the right Jacobian term, maybe use a positive ordered, then subtract 100?

Better solution: scale a simplex. I think that'll 

[Daniel Lee]

On Thursday, November 27, 2014, Daniel Lee <bea...@alum.mit.edu> wrote:

On Thursday, November 27, 2014, Daniel Lee <bea...@alum.mit.edu> wrote:

On Thursday, November 27, 2014, <louis...@gmail.com> wrote:

Dear Bob,

Thank you for the feedback. I agree on the identifiability issue. I knew that but must have missed this and somehow my quick check did not reveal this issue.

I am however a bit puzzled by your comment on the support. Could you clarify the following?

-I thought it is ok to put a restriction on the declaration and have an unrestricted prior

e.g.

declare real<lower=0> alpha[K]

and later on 

alpha~cauchy(0, 2);

but the opposite is not ok? 

No restriction in the declaration but specifying a prior with some restriction e.g. uniform(a,b)  ?

That is correct. If the parameters end up outside of the support of the uniform distribution, the log probability is -inf and there is no gradient. So, NUTS / HMC stops working efficiently. There is a chance it stumbles back into an area of support, but I think our implementation won't.

By default, when Stan initializes its parameters, it chooses random values on the unconstrained scale. If the initial parameters evaluate to a -inf log probability calculation, Stan will choose a different random starting point. It does this a max of 100 times then quits if it can't find a reasonable starting value.

If you're insistent on using that sort of construction, you could start sampling by providing initial values, but it would still be very inefficient. Stan does a change of variables under the hood, which make sampling much easier over bounded variables. There's no possible chance of stepping into a place with no support.

 

-What would you suggest in the example of above? 

(how would you restrict the ordered declaration if that is the suggestion)

You might be able to do a change of variables yourself and add the right Jacobian term, maybe use a positive ordered, then subtract 100?

Better solution: scale a simplex. I think that'll 

*sorry, accidentally hit send.

I think that'll work. You should be able to use a stick breaking construction over the -100 to 100 range.

[Jackson King]

After setting this aside for some time, I have decided to come back to the stan implementation of the GRM and am

still really struggling. If there are to be C-1 ordered thresholds for k items that are used in the model, then

how might I implement this using the ordinal logistic density? 

This, for instance: 

grm_code<-"

data {

int<lower=1> J;                // number of students

int<lower=1> K;                // number of questions

int<lower=1> N;                // number of observations

int<lower=1> C;                // number of categories.

int<lower=1,upper=J> jj[N];    // student for observation n

int<lower=1,upper=K> kk[N];    // question for observation n

int<lower=1,upper=4> y[N];     // correctness of observation n

}

parameters { 

for(k in 1:K)

ordered[C-1] tau[k]; 

real delta[K];

real theta[J];

real<lower=0> sigma;

} 

model{

tau ~ normal(0, 10);

theta ~ normal(0,1); 

delta ~ normal(0,sigma);   

for (n in 1:N){

for(k in 1:K){

y[n] ~ ordered_logistic(delta[kk[n]]*theta[jj[n]], tau[k]); 

}

}

}

"

gives this error:

  

  TRANSLATING MODEL 'grm_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 'grm_code' with error message:

                 SYNTAX ERROR, MESSAGE(S) FROM PARSER:

                 

                 ERROR at line 13

               12:    parameters { 

                 13:    for(k in 1:K)

                   ^

                   14:    ordered[C-1] tau[k]; 

If I were to build this manually using a simplex, any pointers on where I would start. This, for instance, doesn't seem to work:

parameters { 

  simplex[C-1] tau_raw; 

} 

transformed parameters { 

  ordered[C-1] tau; 

  tau <- 100 * tau_raw; 

} 

I apologize if these are elementary questions -- I've estimated these models in JAGS, but am struggling to translate 

the code to stan

Thanks.

[Daniel Lee]

What version of RStan are you running? Mind upgrading to the latest? It has better compiler messages.

[Bob Carpenter]

The next version of Stan's going to provide better
error messages.

You can't have loops in the parameters block. You want
to declare tau with

ordered[C - 1] tau[K];

if you want tau to be a size-K array of ordered (C - 1)-vectors.

> parameters {
> simplex[C-1] tau_raw;
> }
> transformed parameters {
> ordered[C-1] tau;
> tau <- 100 * tau_raw;
> }

This won't work because tau isn't guaranteed to be ordered
if tau_raw is a (C - 1)-simplex. But I'm not sure what you're
trying to do.

If you want to create a parameter that's a centered vector,
the way I was previously recommending is a bad idea compared to
just doing this:

parameters {
vector[N-1] beta_raw;
...
}
transformed parameters {
vector[N] beta;
for (n in 1:(N-1))
beta[n] <- beta_raw[n];
beta[N] <- -sum(beta_raw);
...
}

If you (re-)send the JAGS code, we can help translate.

- Bob

[Jackson King]

Thanks. Here is the JAGS code.  

I'm still struggling with how to easily translate it to stan. Thanks again.

model{

  for(i in 1:n){  # n is respondent

    for (j in 1:p){ #Loop over item, p is item 

      y[i,j]~dcat(prob[i,j,1:K[j]]) 

    }

    theta[i]~dnorm(0, 1)

    for (j in 1:p){ #Loop over item, p is item

      for (k in 1:(K[j]-1)){ #Loop over K-1 response categories

        logit(P[i,j,k])<-alpha[j]*theta[i] -kappa[j,k] 

      }

      P[i,j,K[j]]<-1

    }

    for (j in 1:p){

      prob[i,j,1]<-P[i,j,1]

      for (k in 2:K[j])

        prob[i,j,k]<-P[i,j,k]-P[i,j,k-1]

      }

    }

  }

  #Specify priors/hyperpriors for alpha

  for (j in 1:p){

    alpha[j]~dnorm(mu,prec) 

  }

  

  mu~dnorm(1,0.01)I(0,) #alpha should be positive -- this is the discrimination parameter.

  prec~dgamma(1,.1)

  

  for (j in 1:p){

    kappa0[j,1]~dunif(-3,3)

    for (k in 2:K[j]-1){

      kappa0[j,k]~dunif(kappa0[j,k-1],3)

    }

    kappa[j,1:K[j]-1]<-sort(kappa0[j,1:K[j]-1])

  }

}

On Sunday, June 23, 2013 at 12:23:28 PM UTC-7, Jackson King wrote:

[Jackson King]

Sorry -- I should also note I am using the most recent version of stan. Thanks

On Sunday, June 23, 2013 at 12:23:28 PM UTC-7, Jackson King wrote:

[Bob Carpenter]

It would help if you could let us know which variables are data
and which are parameters and what the sizes are. Both are implicit
in how a JAGS program is called.

I find these JAGS models with sorting and explicit calculations
of cumulative distributions (I think that's what's going on with
prob) confusing. So it'd help me if you explained those parts
of the JAGS code in words.

Finally, this model may not be easy to code in Stan because it
looks like the range of y[i,j] depends on K[j]. At least until
we get ragged arrays sorted out, which I'm working on now. Until
then, the Stan code would be ugly in different places.

- Bob

[Aaron Kaat]

I apologize if this posts out of order - it makes me reply to the first post not the recent additions.

Did getting rid of the loop in the parameters block make the code run? If not, note that the indexing variable in the loop was lower case - get rid of the loop and make the array of ordered logistic vectors size K and try again.

I am using the ordered_logistic statement for a GRM model as well and have no problems. Here is my code:

model.stereo <- "

  data{

    int<lower=2, upper=4> K; // four number of categoies

    int N; // 1803 number of individuals

    int D; // 7 number of items

    int<lower=1,upper=K> Y[N,D]; //array of responses

  }

  parameters {

    real<lower=0> slope[D]; // slope

    ordered[K-1] c[D]; //intercepts

    vector[N] theta; // latent trait

  }

  model{

    slope ~ cauchy(0,5);

    theta ~ normal(0,1);

    for (n in 1:N)

      for (d in 1:D)

        Y[n,d] ~ ordered_logistic(theta[n]*slope[d],c[d]);

  }

"

The differences from your code and mine is that I have no missing data (so I don't need to use the database-style coding), I use a cauchy prior on the slope instead of a normal prior with a hyperparameter variance (I'm confused why you use that type of a prior but to each his own), and I do not place a prior on the ordered intercepts (the data find them perfectly fine without them). 

One important note, though, is that the ordered_logistic distribution uses eta minus intercept, whereas other IRT programs (mirt-package, IRTPro, FlexMIRT) use eta plus intercept. To compare results with other programs you need to know this and flip the sign on all intercepts (order and values stay the same).

Hope that helps!

Aaron

[Jackson King]

My apologies -- you are correct. More information would be useful. Below I've included comments throughout the code. I ended up testing it on the Science data in the R ltm package. It consists of four five category items. I realize that Stan cannot yet accommodate ragged arrays, but the number of categories is constant by item here (and in my work), so the vector of response categories by item could just as easily be changed to a scalar where K=4. Other than that, I think the comments will help. If not, please let me know. And thanks in advance for the help. It's much appreciated.

require(ltm)

require(rjags)

data<-Science[c(1,3,4,7)]

data<-cbind(as.numeric(data[,1]),as.numeric(data[,2]), as.numeric(data[,3]), as.numeric(data[,4]) )

N<-length(data[,1]) #N rows

P<-length(data[1,]) #N items

K<-apply(data,2,max, na.rm=T) # This could be stored as a scalar, since K is the same for each item.

data.1 <- list(y=data,N=N,p=P,K=K)

modelstring <- "

model{

  for(i in 1:N){ #Loop over respondent, 

    for (j in 1:p){ #Loop over item

      y[i,j]~dcat(prob[i,j,1:K[j]]) 

    }

    theta[i]~dnorm(0, 1) #Latent ability

    for (j in 1:p){ #Loop over item

      for (k in 1:(K[j]-1)){ #Loop through K-1 response categories

        logit(P[i,j,k])<-kappa[j,k]-alpha[j]*theta[i] 

      }

      P[i,j,K[j]]<-1

    }

    for (j in 1:p){ #Loop over item

      prob[i,j,1]<-P[i,j,1]

      for (k in 2:K[j]){ #Loop through K-1 response categories

        prob[i,j,k]<-P[i,j,k]-P[i,j,k-1]

      }

    }

  }

  #Specify priors/hyperpriors for alpha

  for (j in 1:p){

    alpha[j]~dnorm(mu,prec) 

  }

  

  mu~dnorm(1,0.01)I(0,) #This censors the alpha parameters

  prec~dgamma(1,.1)

  

  for (j in 1:p){

    kappa0[j,1]~dunif(-10,10)

    for (k in 2:K[j]-1){

      kappa0[j,k]~dunif(kappa0[j,k-1],10)

    }

    kappa[j,1:K[j]-1]<-sort(kappa0[j,1:K[j]-1])

  }

}

"

model<-tempfile()

cat(modelstring,file=model)

rats.inits.2 <- function(){

  list("alpha"=rep(1, times=4), 

       "kappa0"=matrix(rep(c(1:3), each=4), nrow=4)

  )

}

m <- jags.model(file=model,  data=data.1, n.chains=3, n.adapt=1000, inits=rats.inits.2) 

draws <- coda.samples(m, 1000, variable.names=c("alpha", "kappa", "theta"), 10)

[Jackson King]

This seems to work. When I remove the loop from the parameters, the model produces estimates quite similar to grm in the ltm package. Thanks. I really appreciate all the help with this.

[Aaron Kaat]

so why can't you do (for the model block)

real<lower=0> mu;

real<lower=0> prec;

theta ~ normal(0,1);

alpha ~ normal(mu,var);

mu ~ normal(1,100);

var ~ //gamma(1,.01) was was you used for precision, new prior for the variance?

for(n in 1:N){ //loop over individuals

for(p in 1:P){ //loop over items

kappa[n,p] ~ ordered_logistic(alpha[p]*theta[n],cut[p]);

}

}

where "cut" is your ordered[P-1] array of vectors?

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/Kz6y8aIrKKY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.