Skip to content

Instantly share code, notes, and snippets.

@stnorton
Last active Apr 11, 2019
Embed
What would you like to do?
Stan 'i' error reproduction
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);
}
## 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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment