Stan 'i' error reproduction

model.stan

data {

  int N; //number of obs
  int i; //number of legislators THIS CREATES THE ERROR
  int J; //number of bills
  int K; //legislator covariates

  matrix[i,K] Z; //predictors
  int<lower = 0, upper = 1> vote[N]; //observations
  int bills[N]; //bil indicator
  int legis[N]; //legislator indicator

}

parameters {

  vector[J] beta; //discrimination
  vector[J] alpha; //difficulty
  vector[i] theta; //ideology

  vector[K] gamma; //senator-level covariate coefficients

  real mu_beta; //prior for mean of discrimination
  vector<lower = 0.0001>[2] sigmas; //variances of priors for diff and disc


}

model {
  vector[i] mu_theta;
  mu_theta = Z * gamma;

  //likelihood
  vote ~ bernoulli_logit(beta[bills] .* theta[legis] - alpha[bills]);

  //priors
  theta ~ normal(mu_theta, 1);
  gamma ~ student_t(4,0,1);
  sigmas ~ cauchy(0,5);
  beta ~ normal(mu_beta, sigmas[1]);
  alpha ~ normal(0, sigmas[2]);
  mu_beta ~ normal(0,25);

}

r_file.R

## LAB 6
## SEAN NORTON

##LIBRARY-----------------------------------------------------------------------
library(tidyverse)
library(rstan)
options(mc.cores = 2)

#READ IN------------------------------------------------------------------------
s114 <- read.csv("https://voteview.polisci.ucla.edu/static/data/ord/senate_114.csv")
s114 <- s114[-1,] #get rid of obama

vote_mat <- as.matrix(s114[,10:ncol(s114)])
vote_mat[vote_mat>6] <- NA
vote_mat[vote_mat==6] <- 0

#flatten
vote_vec <- c(vote_mat)
miss <- is.na(vote_vec)
vote_vec <- vote_vec[!miss]

#senator indicator vector
senators <- rep(1:100, times = ncol(vote_mat))
senators <- senators[!miss]

#bills indicator vector
bills <- rep(1:ncol(vote_mat), each = 100)
bills <- bills[!miss]

#covariates
Z <- cbind(ifelse(s114$party==200,1,0))

#dataset
irt_data <- list(N = length(vote_vec),
                 K = 1,
                 i = 100,
                 J = ncol(vote_mat),
                 vote = vote_vec,
                 bills = bills,
                 legis= senators,
                 Z=Z
)

#model
init_vals <- list(list(theta=ifelse(Z[,1]==1, 1, -1)),
                                 list(theta=ifelse(Z[,1]==1, 1, -1)))

irt.fit <- stan('norton_lab6.stan',
                data = irt_data,
                chains = 2,
                init = init_vals,
                iter = 500, 
                seed = 1017)

#convergence checks
(id_res <- summary(irt.fit,"theta")$summary[,1])

summary(summary(irt.fit)$summary[,"Rhat"])

traceplot(irt.fit)

summary(summary(irt.fit)$summary[,"n_eff"])

#comparison to ideal
library(pscl)
ideal.irt <- ideal(rollcall(vote_mat)
                   ,normalize=TRUE
                   ,startvals=list(x=init_vals[[1]]$theta))
ideal_res <- summary(ideal.irt)$xm

#plot
comp_df <- cbind.data.frame(id_res, ideal_res)
colnames(comp_df)[2] <- 'ideal'

ggplot(comp_df, aes(x = id_res, y = ideal)) + 
  geom_point() +
  labs(x = 'Stan', y  = 'Ideal') +
  theme_bw()