Multinomial IRT Model Error : attempt to assign variable

[John]

Hi Stan users,

I have a question regarding a multinomial item response theory model. I am trying to transform the following BUGS code to a STAN code.

BUGS code here.

for(i in 1:lengthy){
      yy[i] ~ dcat(prob[nidx[i],midx[i],1:K[midx[i]]])
    }
    for (i in 1:lengthy){
        for (k in 1:K[midx[i]]){
           p[nidx[i],midx[i],k] <- exp(beta[midx[i],k]*x[nidx[i]] + beta2[midx[i],k]*x2[nidx[i]] - alpha[midx[i],k])
        }
        prob[nidx[i],midx[i],1:K[midx[i]] ] <- p[nidx[i],midx[i],1:K[midx[i]] ]/sum(p[nidx[i],midx[i],1:K[midx[i]]])
    } 

Through the second for loop, I will obtain an array (number of individuals by number of questions by number of categories in questions) of probabilities. These probabilities will go through dcat to the response yy.  

I tried to code this in Stan because it is faster to implement this IRT model. But I got the following error message:

attempt to assign variable in wrong block. left-hand-side variable origin=parameter

Here is the Stan code that I wrote.

library(rstan)
stan.code <- "data {
  int<lower=1> J; // number of Individuals
  int<lower=1> M; // number of Questions
  int<lower=1> N; // number of observations
  int<lower=1> K; // maximum number of categories
  int<lower=1, upper=J> j[N]; // Individual for observation n
  int<lower=1, upper=M> m[N]; // Question for observation n
  int<lower=1, upper=K> y[N]; // Answer of observation n
  int<lower=1, upper=K> k[M]; // Maximum categories of question m
  int<lower=1, upper=K> q[N]; // Maximum categories of observation n
  int<lower=1> D; //no. of dimensions
}
parameters {
  matrix[M,K] alpha;
  matrix[M,K] beta;
  matrix[M,K] beta2;
  real theta[J];
  real theta2[J];
  real<lower=0> p[J, M, K];
  real<lower=0> prob[J, M, K];
}
model {
  to_vector(alpha) ~ normal(0,1); 
  for (i in 1:M){
    alpha[i,1] ~ normal(0,0.000001);
    for (g in (k[i]+1):K)
      alpha[i,g] ~ normal(999,0.000001);} 
  
  to_vector(beta) ~ normal(0,1); 
  for (i in 1:M){
    beta[i,1] ~ normal(0,0.000001);
    for (g in (k[i]+1):K)
      beta[i,g] ~ normal(999,0.000001);}
  
  to_vector(beta2) ~ normal(0,1); 
  for (i in 1:M){
    beta2[i,1] ~ normal(0,0.000001);
    for (g in (k[i]+1):K)
      beta2[i,g] ~ normal(999,0.000001);}
  
  theta ~ normal(0,0.1);
  theta2 ~ normal(0,0.1);   // priors
  
  for (n in 1:N){
    for (h in 1:q[n]){
      p[j[n], m[n], h] = exp(theta[j[n]] * beta[m[n],h] + theta2[j[n]] * beta2[m[n],h] - alpha[m[n],h])}
      prob[j[n], m[n], 1:q[n]] = p[j[n], m[n], 1:q[n]]/sum(p[j[n], m[n], 1:q[n]])}
   for (n in 1:N){
      y[n] ~ categorical(prob[j[n], m[n], 1:q[n]])
   }
}
"  
  
  J <- dim(jdata$y)[1]      # number of Individuals
  M <- dim(jdata$y)[2]      # number of Questions
  N <- length(jdata$y)      # number of observations
  K <- max(jdata$K)         # maximum number of categories
  
  j <- rep(1:J,M)           # individual for observation n
  m <- rep(1:M, each=J)     # question for observation n
  y <- as.vector(jdata$y)   # answer of observation n
  k <- jdata$K              # maximum categories of question m
  q <- rep(jdata$K, each=J) # maximum categories of observation n
  D <- 2                    # no of dimensions
    
  # deleting missing values
  miss <- which(is.na(y))
  N <- N - length(miss)
  j <- j[-miss]
  m <- m[-miss]
  y <- y[-miss]
  q <- q[-miss]
  
  ## data and initial values
  
  stan.data <- list(J=J, M=M, N=N, K=K, j=j, m=m, y=y, k=k, q=q, D=D)
  starting=list(theta=startvalues+rnorm(length(startvalues)),
                theta2=rnorm(length(startvalues)))
  
  stan.fit <- stan(model_code=stan.code, data=stan.data, iter=500, warmup=200,
                   chains=1, thin=2, init=starting)

Why is my code wrong? How should I fix the code to implement the same thing as BUGS?

Thank you in advance!

John

[Stephen Martin]

That error just means you are trying to overwrite a parameter (as defined in parameters block), which you can't do.

From just a quick glance, that seems to be p[j[n], m[n], h], which you're defining as a parameter, but then assigning to in the model block.

I'm guessing this is due to some confusion about what the parameters block is for. Parameters{} isn't for listing "things I want estimated", it's for things you want "sampled" (really, for defining ficticious particles that are simulated in multidimensional space). If p[] is a COMPUTED parameter, then you should compute it in the transformed parameters block. If it's just a utility variable that you don't need to save or do inference on, but is merely used in computing the posterior, then declare and define it directly in the model block.