Fixing horrible hacks in a hierarchical IRT model

[Christopher Gandrud]

Hi everyone

I have a hierarchical IRT model that runs ok. There are just two spots where I have made horrible hacks to accomplish my goals. I'm wondering if anyone could help me clean these up.

The full model code is below, but I'll just highlight the hacks first.

Hack 1: I am measuring ability (what I call transparency and denote alpha) for J countries for T years. For all years T > 1, the alpha Tn - 1 forms the prior for year T. I want to recenter the ability for T = 1 by subtracting the mean ability and dividing by the standard deviation.

I do this in the transformed parameters block:

transformed parameters {           
        // recenters transparency for t = 1
        vector[J] recentered_alpha1;
        real mean_alpha1;
        real<lower=0> sd_alpha1;
 
        mean_alpha1 <- mean(alpha1);
        sd_alpha1 <- sd(alpha1);
        for (j in 1:J)
            recentered_alpha1[j] <- ( alpha1[j] - mean_alpha1 ) / sd_alpha1;
}

The problem is that I want to put recentered_alpha[j] into an alpha matrix at alpha[j,1]. I can't seem to do this in the transformed parameters block so do the horrible hack:

alpha[j,1] ~ normal(recentered_alpha1[j], 0.001);

in the model block.

Hack 2: I would like to have a half-Cauchy prior for each country but I am having trouble vectorising this. I think STAN doesn't support truncation with Cauchy distributions and I can't seem to set the lower bound in all items of a J length list. My solution so far basically kills a lot of the benefits of vectorisation:

for (j in 1:J) {
            alpha[j,1] ~ normal(recentered_alpha1[j], 0.001); 
            sigma_alpha ~ cauchy(0,0.25);                    
            for (t in 2:T) 
                alpha[j,t] ~ normal(alpha[j,t-1], sigma_alpha); 
}

Any help would be much appreciated!

The full model:

    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];    // item for observation n
        int<lower=0,upper=1> y[N];     // response for observation n
    }
 
    parameters {
        real delta;                    // mean transparency
        vector[J] alpha1;              // initial alpha for t = 1 before recentering
        matrix[J,T] alpha;             // transparency for j,t - mean
        vector[K] beta;                // difficulty of item k
        vector[K] log_gamma;           // discrimination of k
 
        // all scale parameters have an implicit half Cauchy prior
        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 {           
        // recenters transparency for t = 1
        vector[J] recentered_alpha1;
        real mean_alpha1;
        real<lower=0> sd_alpha1;
 
        mean_alpha1 <- mean(alpha1);
        sd_alpha1 <- sd(alpha1);
        for (j in 1:J)
            recentered_alpha1[j] <- ( alpha1[j] - mean_alpha1 ) / sd_alpha1;
    }
 
    model {
        alpha1 ~ normal(0,100);                                 // diffuse prior
        for (j in 1:J) {
            alpha[j,1] ~ normal(recentered_alpha1[j], 0.001);   // horrible hack
            sigma_alpha ~ cauchy(0,0.25);                    
            for (t in 2:T) 
                alpha[j,t] ~ normal(alpha[j,t-1], sigma_alpha); 
        }
        beta ~ normal(0,sigma_beta);
        log_gamma ~ normal(0,sigma_gamma);
        delta ~ cauchy(0,0.25);
        sigma_beta ~ cauchy(0,0.25);
        sigma_gamma ~ cauchy(0,0.25);
        for (n in 1:N)
            y[n] ~ bernoulli_logit(
                    exp(log_gamma[kk[n]])
                    * (alpha[jj[n],tt[n]] - beta[kk[n]] + delta) );
}

[Marcus Brubaker]

On Tue, Sep 30, 2014 at 10:37 AM, Christopher Gandrud <christoph...@gmail.com> wrote:

Hi everyone

I have a hierarchical IRT model that runs ok. There are just two spots where I have made horrible hacks to accomplish my goals. I'm wondering if anyone could help me clean these up.

The full model code is below, but I'll just highlight the hacks first.

Hack 1: I am measuring ability (what I call transparency and denote alpha) for J countries for T years. For all years T > 1, the alpha Tn - 1 forms the prior for year T. I want to recenter the ability for T = 1 by subtracting the mean ability and dividing by the standard deviation.

I do this in the transformed parameters block:

transformed parameters {           
        // recenters transparency for t = 1
        vector[J] recentered_alpha1;
        real mean_alpha1;
        real<lower=0> sd_alpha1;
 
        mean_alpha1 <- mean(alpha1);
        sd_alpha1 <- sd(alpha1);
        for (j in 1:J)
            recentered_alpha1[j] <- ( alpha1[j] - mean_alpha1 ) / sd_alpha1;
}

The problem is that I want to put recentered_alpha[j] into an alpha matrix at alpha[j,1]. I can't seem to do this in the transformed parameters block so do the horrible hack:

alpha[j,1] ~ normal(recentered_alpha1[j], 0.001);

in the model block.

What you want to do is declare alpha as a transformed parameter, and create a variable alpha_raw which is a matrix[J,T-1].  E.G.,

parameters {

  ...

  matrix[J,T-1] alpha_raw;

  ...

}

transformed parameters {

  matrix[J,T] alpha;

  ...

  for (j in 1:J) {

    alpha[j,1] <- recentered_alpha1[j];

    for (t in 2:T) {

        alpha[j,t] <- alpha_raw[j,t-1];

    }

  }

}

  

 

Hack 2: I would like to have a half-Cauchy prior for each country but I am having trouble vectorising this. I think STAN doesn't support truncation with Cauchy distributions and I can't seem to set the lower bound in all items of a J length list. My solution so far basically kills a lot of the benefits of vectorisation:

for (j in 1:J) {
            alpha[j,1] ~ normal(recentered_alpha1[j], 0.001); 
            sigma_alpha ~ cauchy(0,0.25);                    
            for (t in 2:T) 
                alpha[j,t] ~ normal(alpha[j,t-1], sigma_alpha); 
}

Any help would be much appreciated!

Not sure I really follow what you're asking here, but you almost certainly don't want that sigma_alpha line inside the loop.  Doing so multiplies that prior into the posterior J times!

Also, you should be able to vectorize the normal distributions if you are willing to sacrifice a little readibility for (potentially) a lot of speed.  See the block() command which extracts a submatrix and that entire double loop can collapse down into a single line.

Cheers,
Marcus

 

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

Hi Marcus

Thanks for the insight. 

I actually do want to create a country (j) specific sigma_alpha. So sigma_alpha should be the same for all all alpha[t+1] for each j. This is for a couple of reasons, the simple on being that not doing so seems to lead the alpha parameters to grow dramatically. 

I would like to vectorize this though so that I only need to call sigma_alpha ~ cauchy(0,0.25) and it will have J values.

Thanks for the block suggestion. Any chance you could point me to an example?

Best

Christopher

--
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/77LB7iW_5fc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Marcus Brubaker]

On Tue, Sep 30, 2014 at 11:39 AM, Christopher Gandrud <christoph...@gmail.com> wrote:

Hi Marcus

Thanks for the insight. 

I actually do want to create a country (j) specific sigma_alpha. So sigma_alpha should be the same for all all alpha[t+1] for each j. This is for a couple of reasons, the simple on being that not doing so seems to lead the alpha parameters to grow dramatically. 

I would like to vectorize this though so that I only need to call sigma_alpha ~ cauchy(0,0.25) and it will have J values.

I just checked and the Cauchy distribution should be vectorized so if you had:

parameters {

  ...

  real<lower=0> sigma_alpha[J];

  ...

}

Then:

model {

  ...

  sigma_alpha ~ cauchy(0,0.25);

  ...

}

should work just fine.  If it doesn't, what error does it produce?

 

Thanks for the block suggestion. Any chance you could point me to an example?

For a single sigma_alpha:

block(alpha,1,2,J,T-1) ~ normal(block(alpha,1,1,J,T-1),sigma_alpha);

should work.

For different sigma_alphas you'll need to loop over J:

for (j in 1:J) {

  block(alpha,j,2,1,T-1) ~ normal(block(alpha,j,1,1,T-1),sigma_alpha);

[Christopher Gandrud]

The sigma_alpha suggestion works perfectly. Thank!

When I try the block suggestion I get:

EXPECTATION FAILURE - DIAGNOSTIC(S) FROM PARSER:

no matches for function name="normal_log"

    arg 0 type=matrix

    arg 1 type=matrix

    arg 2 type=real[1]

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) 

Any thoughts?

[Marcus Brubaker]

Maybe something like:

for (j in 1:J) {

  sub_row(alpha,j,2,T-1) ~ normal(sub_row(alpha,j,1,T-1),sigma_alpha[j]);

}

or

for (j in 1:J) {

  0.0 ~ normal(sub_row(alpha,j,1,T-1)-sub_row(alpha,j,2,T-1),sigma_alpha[j]);

}

Cheers,

Marcus

[Bob Carpenter]

What the message is telling you is that you are
passing

matrix ~ normal(matrix,real[]),

or equivalently,

normal_log(matrix, matrix, real[]).

This won't work. Only the signatures listed are available.
I'm not sure what your intent was with the statement with
matrix/matrix/real[] in terms of how the standard deviations
in the array line up with the means and outcomes in the two
matrices (and this is why we didn't choose to implement these
vectorizations!).

Marcus made a suggestion on how to deal with it, but I
thought I'd explain his reasoning.

- Bob

[Christopher Gandrud]

Thanks for these comments. They are really helpful.

I guess the last thing (and the thing I’m most concerned about substantively) is how to enter my recentered alphas for time 1 into the main alpha matrix.

Currently I create a vector[J] of recentered alpha in the transformed parameters block. I basically want these to be put into alpha[j,1].

At the moment I do this with a horrible hack in the model block like this:

for (j in 1:J)
            alpha[j,1] ~ normal(recentered_alpha1[j], 0.001);

But I clearly want to do something like:

for (j in 1:J)
            alpha[j,1] <- recentered_alpha1[j];

Any thoughts?

​

[Bob Carpenter]

There are two ways I can see to do this.

My preferred suggestion would be to alpha to start at t+1,
replacing

> for (j in 1:J) {
> > alpha[j,1] ~ normal(recentered_alpha1[j], 0.001);
> > sigma_alpha ~ cauchy(0,0.25);
> > for (t in 2:T)
> > alpha[j,t] ~ normal(alpha[j,t-1], sigma_alpha);
> > }

with

parameters {
matrix[J,T-1] alpha;
...
model {
sigma_alpha ~ cauchy(0,0.25); // this should go on outside, not be repeated

for (j in 1:J)

alpha[j,1] ~ normal(recentered_alpha[j], sigma_alpha);

for (t in 2:T)

for (j in 1:J)

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

...

Then, I'd suggest transposing alpha in the parameterization and
writing the above in vector form as:

alpha[j] ~ normal(recentered_alpha, sigma_alpha);

for (t in 2:T)

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

It'd be most efficient if alpha were declared as an array of vectors:

parameters {
vector[J] alpha[T-1];

If you really don't want to do that, you can use the col() extractor
to pull columns out.

If you really want to keep alpha all together (I'm not looking through the
rest of your model), then you can do roughly the same thing, change the
names a bit, and then paste alpha back together in the transformed paramters
block.

parameters {

matrix[J,T-1] alpha_raw;
...
transformed parameters {
matrix[J,T] alpha;

for (j in 1:J)

alpha[j,1] <- recentered_alpha[j];

for (t in 2:T)

alpha[j,t] <- alpha[j,t-1];
}

It'll create a lot of redundant output if you save alpha_raw (only
option in CmdStan, but you can delete variables on the fly in RStan
and presumably in PyStan).

- Bob