Unfolding models in Stan; problems when used transformed parameters block

unfolding.R

library(rstan)
library(mudfold)
library(reshape2)

## Simulate data
K <- 12
J <- 100
simulation1 <- mudfoldsim(N=K, n=J, seed = 1)
dat <- simulation1$dat

## Transform data from wide to long format
dat$respid <- 1:nrow(dat)
dat <- melt(dat, id.var = "respid")

### Number things
dat$variable <- factor(dat$variable,
                       levels = LETTERS[1:K],
                       ordered = TRUE)
dat$variable.num <- as.numeric(dat$variable)

## Turn into a list for Stan
stan_data <- list(J = J,
                  K = K,
                  N = nrow(dat),
                  y = dat$value,
                  jj = dat$respid,
                  kk = dat$variable.num)

### #######################################################
### Model 1: unfolding model
### #######################################################

stan_code <- '
data {
  int<lower=1> J;              // number of respondents
  int<lower=1> K;              // number of questions
  int<lower=1> N;              // number of observations
  int<lower=1,upper=J> jj[N];  // legislator for observation n
  int<lower=1,upper=K> kk[N];  // question for observation n
  int<lower=0,upper=1> y[N];   // vote for observation n
}
parameters {    
  real theta[J];                  // respondent position
  real alpha[K];               // random question intercept
  real omega[K];                // question position
  real <lower=0> beta; // weight on resp.-question distance
}
transformed parameters {
	    real mu[N];
	    for (n in 1:N) {
	    	mu[n] = alpha[kk[n]] - beta * (omega[kk[n]] - (theta[jj[n]]))^2;
	    }
}

model {
  theta ~ normal(0, 1);
  omega ~ normal(0, 1);
  alpha ~ normal(0, 10);
  beta ~ normal(0, 1);

  for (n in 1:N) {
    y[n] ~ bernoulli_logit(mu[n]);
  }


}
'

writeLines(stan_code, con = "stan_model_temp.stan")

m <- stan_model("stan_model_temp.stan")

f <- sampling(m,
              data = stan_data,
              cores = 1, chains = 1,
              iter = 1000,
              control = list(adapt_delta = 0.95,
                             max_treedepth = 20))

## Did the model recover the correct ordering?

### #######################################################
### Model 2: unfolding model, scale theta by 2/3rds
### Do this in the transformed parameters block
### DOES NOT WORK
### #######################################################

stan_code_2 <- '
data {
  int<lower=1> J;              // number of respondents
  int<lower=1> K;              // number of questions
  int<lower=1> N;              // number of observations
  int<lower=1,upper=J> jj[N];  // legislator for observation n
  int<lower=1,upper=K> kk[N];  // question for observation n
  int<lower=0,upper=1> y[N];   // vote for observation n
}
parameters {    
  real theta[J];                  // respondent position
  real alpha[K];               // random question intercept
  real omega[K];                // question position
  real <lower=0> beta; // weight on resp.-question distance
}
transformed parameters {
	    real mu[N];
	    real thetastar[J];
	    for (n in 1:N) {
	    	mu[n] = alpha[kk[n]] - beta * (omega[kk[n]] - (thetastar[jj[n]]))^2;
	    }
	    for (j in 1:J) {
	    	thetastar[j] = 2.0/3.0 * theta[j];
	    }
}

model {
  theta ~ normal(0, 1);
  omega ~ normal(0, 1);
  alpha ~ normal(0, 10);
  beta ~ normal(0, 1);

  for (n in 1:N) {
    y[n] ~ bernoulli_logit(mu[n]);
  }

}
'

writeLines(stan_code_2, con = "stan_model_temp_2.stan")

m2 <- stan_model("stan_model_temp_2.stan")

f2 <- sampling(m2,
              data = stan_data,
              cores = 1, chains = 1,
              iter = 1000,
              control = list(adapt_delta = 0.95,
                             max_treedepth = 20))


### #######################################################
### Model 2b: unfolding model, scale theta by 2/3rds
### Do this inline whilst creating variable mu
### WORKS
### #######################################################

stan_code_2b <- '
data {
  int<lower=1> J;              // number of respondents
  int<lower=1> K;              // number of questions
  int<lower=1> N;              // number of observations
  int<lower=1,upper=J> jj[N];  // legislator for observation n
  int<lower=1,upper=K> kk[N];  // question for observation n
  int<lower=0,upper=1> y[N];   // vote for observation n
}
parameters {    
  real theta[J];                  // respondent position
  real alpha[K];               // random question intercept
  real omega[K];                // question position
  real <lower=0> beta; // weight on resp.-question distance
}
transformed parameters {
	    real mu[N];
	    for (n in 1:N) {
	    	mu[n] = alpha[kk[n]] - beta * (omega[kk[n]] - (2.0/3.0 * theta[jj[n]]))^2;
	    }
}

model {
  theta ~ normal(0, 1);
  omega ~ normal(0, 1);
  alpha ~ normal(0, 10);
  beta ~ normal(0, 1);

  for (n in 1:N) {
    y[n] ~ bernoulli_logit(mu[n]);
  }
}
'

writeLines(stan_code_2b, con = "stan_model_temp_2b.stan")

m2b <- stan_model("stan_model_temp_2b.stan")

f2b <- sampling(m2b,
              data = stan_data,
              cores = 1, chains = 1,
              iter = 1000,
              control = list(adapt_delta = 0.95,
                             max_treedepth = 20))