Problem setting up Hierarchical IRT model with 'student' and 'test' levels

[Christopher Gandrud]

I have a question that I think has a simple answer, but I can't seem to see it.

Set up: I'm building an item-response model based directly on the multilevel 2PL example from the docs. My data is in country-year format with binary variables for whether or not country j reported a given data point to the World Bank in year t. I'm primarily interested in estimating each country's 'transparency' for each year. This is equivalent to 'student's ability' in the Stan IRT example.

In the model the transparency variables are denoted with alpha.

Problem: I can easily run the model for just one year. You can fully reproduce this with rstan using code at:

https://raw.githubusercontent.com/FGCH/FRTIndex/stanTest/source/FRT_Stan.R

What I'm stuck on is how to estimate alpha for each country year j,t (I guess this would be equivalent to having 'student' and 'number times the student took the test'). 

This is what I have so far:

    data {
        int<lower=1> J;                // number of countries
        int<lower=1> T;                // number of years
        int<lower=1> K;                // number of items
        int<lower=1> N;                // number of observations
        int<lower=1> jj[N];            // country for observation n
        int<lower=1> tt[N];            // time for observation n
        int<lower=1,upper=K> kk[N];    // question for observation n
        int<lower=0,upper=1> y[N];     // correctness for observation n
    }

    parameters {
        real delta;                    // mean transparency
        real alpha[J,T];               // transparency for j,t - mean
        real beta[K];                  // difficulty of item k
        real log_gamma[K];             // discrimination of k
        real<lower=0> sigma_alpha;     // scale of abilities
        real<lower=0> sigma_beta;      // scale of difficulties
        real<lower=0> sigma_gamma;     // scale of log discrimiation
    }

    model {
        alpha ~ normal(0,sigma_alpha);
        beta ~ normal(0,sigma_beta);
        log_gamma ~ normal(0,sigma_gamma);
        delta ~ cauchy(0,5);
        sigma_alpha ~ cauchy(0,5);
        sigma_beta ~ cauchy(0,5);
        sigma_gamma ~ cauchy(0,5);
        for (n in 1:N)
            y[n] ~ bernoulli_logit( exp(log_gamma[kk[n]])
                                * (alpha[jj[n],tt[n]] - beta[kk[n]] + delta) );
    }

However, using this model in stan returns the error no matches for function name="normal_log". The full replication code is available at:

https://raw.githubusercontent.com/FGCH/FRTIndex/stanTest/source/FRT_Stan_v0_2.R

Any help with this would be much appreciated.

[Daniel Lee]

Hi Christopher,

For future reference, it helps to provide the actual error message. My guess is that the error has to do with this line:

       alpha ~ normal(0,sigma_alpha);

The "no matches for function..." message probably lists a line and what arguments you provided with a list of available arguments.

The normal distribution (currently) does not accept matrix arguments, so you'll have to loop over the years or countries.

Daniel

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

[Christopher Gandrud]

Sorry about that. Here is the full error message:

no matches for function name="normal_log"
    arg 0 type=real[2]
    arg 1 type=int
    arg 2 type=real
available function signatures for normal_log:
0.  normal_log(real, real, real) : real
1.  normal_log(real, real, real[1]) : real
2.  normal_log(real, real, vector) : real
3.  normal_log(real, real, row vector) : real
4.  normal_log(real, real[1], real) : real
5.  normal_log(real, real[1], real[1]) : real
6.  normal_log(real, real[1], vector) : real
7.  normal_log(real, real[1], row vector) : real
8.  normal_log(real, vector, real) : real
9.  normal_log(real, vector, real[1]) : real
10.  normal_log(real, vector, vector) : real
11.  normal_log(real, vector, row vector) : real
12.  normal_log(real, row vector, real) : real
13.  normal_log(real, row vector, real[1]) : real
14.  normal_log(real, row vector, vector) : real
15.  normal_log(real, row vector, row vector) : rea

[Soga Shota]

probably,

alpha ~ normal(0,sigma_alpha);

to

alpha ~ normal(0.0,sigma_alpha);

may work? 0 is considered as an integer whereas 0.0 as considered as real.

The error massage is sufficient for this case. I guess it is difficult for non C++ people to understand... :(

Shota

 
2014年7月11日金曜日 0時09分00秒 UTC+9 Christopher Gandrud:

[Daniel Lee]

Hi Chris,

The issue is:

no matches for function name="normal_log"
    arg 0 type=real[2]
    arg 1 type=int
    arg 2 type=real

For the left hand side of this expression:

alpha ~ normal(0,sigma_alpha);

you have alpha defined as:

real alpha[J,T];               // transparency for j,t - mean

If you want to use alpha as is, you'll need to replace "alpha ~ normal(0, sigma_alpha);" with:

    for (j in 1:J)

      for (t in 1:T)

        alpha[j, t] ~ normal(0, sigma_alpha);

or

    for (j in 1:J)
       alpha[j] ~ normal(0, sigma_alpha);

The second one uses the vectorized form for the left side.

So... the message is telling you that Stan doesn't have the normal distribution defined for a 2-d array of reals. We need to add better error messages, but that's how to interpret "arg 0 type=real[2]".

Daniel

--

[Bob Carpenter]

On Jul 10, 2014, at 7:52 PM, Soga Shota <soga....@gmail.com> wrote:

> probably,
>
> alpha ~ normal(0,sigma_alpha);
> to
>
> alpha ~ normal(0.0,sigma_alpha);
> may work? 0 is considered as an integer whereas 0.0 as considered as real.
>
> The error massage is sufficient for this case. I guess it is difficult for non C++ people to understand... :(

If there's ever a place where 0.0 works but 0 doesn't,
that's a bug and we'd appreciate hearing about it.
C++ doesn't have this problem with non-templated
functions, because int can be promoted to double. Ditto
for Java and C. But in C++ templates, it's possible to
specialize so that there's a definition for double inputs
but not int inputs; and if the templates can't resolve,
no go.

The problem here, as Daniel pointed out, is that the
normal_log function doesn't accept real[2] (i.e., 2D
array of real) arguments.

I want to replace the error output "real[2]" with
"real[,]" because that's what we're using in the function
syntax and I think it'll be clearer than just reporting the
number of arguments as we currently do:

> no matches for function name="normal_log"
> arg 0 type=real[2]
> arg 1 type=int
> arg 2 type=real

- Bob

[Soga Shota]

Hi Bob,

It was my bad. stan properly treats an integer argument as real.

Shota

2014年7月11日金曜日 12時45分09秒 UTC+9 Bob Carpenter:

[Christopher Gandrud]

Thanks! This discussion has been really useful. It's really helping me convert my IRT models to Stan from Jags (so far I've definitely had the advertised convergence speed improvements).

This worked perfectly:

for (j in 1:J)
    for (t in 1:T)
        alpha[j,t] ~ normal(0, sigma_alpha);

FYI for anyone else that might be interested: I've updated the priors slightly such that at time 1 alpha[j,t] ~ normal(0, sigma_alpha), but then for all t_1+n the previous year's alpha informs the prior for the following year, i.e.

for (j in 1:J)
    alpha[j,1] ~ normal(0, sigma_alpha);
for (j in 1:J)
    for (t in 2:T)
        alpha[j,t] ~ normal(alpha[j,t-1], sigma_alpha);

Any thoughts on further tweaks would be of course be welcome.

[Bob Carpenter]

You can probably continue to vectorize.
I'm not sure what your data types are, though.

If alpha[j,t] is a scalar, then it might help to
transpose alpha and use alpha_t[t,j], at which

point you could write:

for (t in 1:T)

alpha_t[t] ~ normal(alpha_t[t-1], sigma_alpha);

You'd also be able to write just

alpha_t[1] ~ normal(0, sigma_alpha);

Even better if you declare alpha_t as an array of vectors.
There's a chapter in the manual now on arrays versus
vectors and speed and where they are required.

This can be a really big speedup if J >> 1.

- Bob

[Christopher Gandrud]

Thanks for the suggestion!

This should be pretty easy, but I'm stuck on how to transpose alpha in the way that Bob recommends in the previous post. alpha[j,t] is scalar so I have:

    parameters {

        real delta;                    // mean transparency

        real alpha[J,T];               // transparency for j,t - mean

        vector[K] beta;                // difficulty of item k

        vector<lower=0>[K] gamma;      // discrimination of k

        real<lower=0> sigma_alpha;     // scale of abilities

        real<lower=0> sigma_beta;      // scale of difficulties

        real<lower=0> sigma_gamma;     // scale of log discrimination

    }

    transformed parameters {

        row_vector[J] alpha_t;

        alpha_t <- alpha';

    }

But get the basic error:

EXPECTATION FAILURE - DIAGNOSTIC(S) FROM PARSER:

no matches for function name="transpose"

    arg 0 type=real[2]

available function signatures for transpose:

0.  transpose(vector) : row vector

1.  transpose(row vector) : vector

2.  transpose(matrix) : matrix

no matches for function name="transpose"

    arg 0 type=real[2]

available function signatures for transpose:

0.  transpose(vector) : row vector

1.  transpose(row vector) : vector

2.  transpose(matrix) : matrix

expression is ill formed

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

PARSED:

transformed parameters {

        real alpha_t[T,J];

        alpha_t <- transpose(alpha)

EXPECTED: <eps> BUT FOUND: 

;

    }

-------

Any help would be much appreciated.

[Daniel Lee]

Hi Chris,

If all you're trying to do is transpose alpha, wouldn't it just be easier to define alpha in the parameters block as "row_vector[J] alpha[T]" or something similar?

Daniel

--

[Christopher Gandrud]

Thanks for the idea Daniel. 

Probably again missing something simple. So I do:

    parameters {

        real delta;                    // mean transparency

        row_vector[J] alpha[T];        // transparency for j,t - mean

        vector[K] beta;                // difficulty of item k

        vector<lower=0>[K] gamma;      // discrimination of k

        real<lower=0> sigma_alpha;     // scale of abilities

        real<lower=0> sigma_beta;      // scale of difficulties

        real<lower=0> sigma_gamma;     // scale of log discrimination

    }

    model {

        alpha[1] ~ normal(0, sigma_alpha); 

        for (t in 1:T) 

            alpha[t] ~ normal(alpha[t-1], sigma_alpha); 

The model compiles but I get the error:

SAMPLING FOR MODEL 'frt_code' NOW (CHAIN 1).

Error : INDEX OPERATOR [] OUT OF BOUNDS; index=0; lower bound=1; upper bound=14; index position=1; alpha

error occurred during calling the sampler; sampling not done

[Daniel Lee]

In this loop:

        for (t in 1:T) 

            alpha[t] ~ normal(alpha[t-1], sigma_alpha);

when t == 1, you're asking for alpha[0] (alpha[t-1]), which is not allowed.

--

[Christopher Gandrud]

Ha, yeah that was simple.  

Though with:

for (t in 2:T) 

            alpha[t] ~ normal(alpha[t-1], sigma_alpha);

I get:

SAMPLING FOR MODEL 'frt_code' NOW (CHAIN 1).

Error : INDEX OPERATOR [] OUT OF BOUNDS; index=15; lower bound=1; upper bound=14; index position=1; alpha

error occurred during calling the sampler; sampling not done

Not sure where the 15 comes from as T = 14

[Daniel Lee]

I'd check your other uses of alpha. Can't remember your model off the top of my head.

--

[Bob Carpenter]

Each data structure has a dimensionality, for instance, both of

matrix[M,N] a;

vector[N] a[M];

define an (M x N) structure. To use Daniel's suggestion, you need
to read the data in transposed.

To just transpose in the model, all you need to replace the declaration

real alpha[J,T];

with

matrix[J,T] alpha;

and then define:

transformed data {
matrix[T,J] alpha_t;
alpha_t <- alpha';
}

even more efficient would be:

transformed data {
vector[J] alpha_t[T];

for (t in 1:T)

alpha_t[t] <- col(alpha,t);
}

Then just use alpha_t in the code as suggested earlier.

But I wouldn't worry about any of this if the model's fast
enough without it.

- Bob

[Christopher Gandrud]

Just using:

matrix[J,T] alpha;

(with the corresponding change in the priors setup) rather than:

real alpha[J,T]; 

seems to have halved the computation time!

[Bob Carpenter]

I wouldn't have expected that. Typically, almost all
of Stan's compute time is taken up computing derivatives
and matrix vs. 2D array doesn't change that.

There's a chapter in the manual that explains what's going on
with efficiency for arrays and vectors.

- Bob

[Avraham Adler]

Something to do with row vs. column order memory perhaps?

Avi

On Monday, July 14, 2014 2:10:36 PM UTC-4, Bob Carpenter wrote:

I wouldn't have expected that.  Typically, almost all
of Stan's compute time is taken up computing derivatives
and matrix vs. 2D array doesn't change that.

There's a chapter in the manual that explains what's going on
with efficiency for arrays and vectors.  

- Bob

On Jul 14, 2014, at 1:31 PM, Christopher Gandrud <christoph...@gmail.com> wrote:

> Just using:
>
> matrix[J,T] alpha;
>
> (with the corresponding change in the priors setup) rather than:
>
> real alpha[J,T];
>
> seems to have halved the computation time!
>
>
> --
> 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+unsubscribe@googlegroups.com.

[Bob Carpenter]

The manual chapter explains all of that.

real[,] : std::vector<std:vector<double> >.

matrix : Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>.

No copying should be needed because indexing returns a reference.

We haven't seen these 2D array vs. matrix differences elsewhere, but
maybe we just haven't profiled carefully enough in the context of
very simple models where it might make a difference.

- Bob

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

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

[Christopher Gandrud]

I suspect it doesn't have to do with going from real to matrix. Instead, doing this allowed me to get rid of a couple of loops.