simulating_crap.R

# Simulating shit
library(tidyverse)
library(magrittr)
# Model 

# 10 parties
np <- 10
# uniform distribution for this parties seq(10) - number of legislators
nl <- sample(10:40, np, replace = T)

# parties ideal points
repo <- Vectorize(rep)
rollcalls <- data_frame(id_party=1:np, phi=rnorm(np), nl=nl) # priors for parties N(0,1)
for (i in 1:nrow(rollcalls)) {
  rollcalls$id_party[i] <- paste0(rep(as.character(rollcalls$id_party[i]), rollcalls$nl[i]), collapse="#")
}
rollcalls %<>%
  unnest(id_party = strsplit(id_party, "#")) %>%
  select(-nl)
# \theta_i <- \gamma_i + \phi_k


rollcalls$id_party <- as.numeric(rollcalls$id_party)
rollcalls$id_legis <- 1:nrow(rollcalls)
rollcalls$gamma <- rnorm(nrow(rollcalls))
rollcalls$theta <- rollcalls$phi + rollcalls$gamma

nr <- 200 # number rollcalls
matrix_rollcalls <- cbind(rnorm(nr), rnorm(nr))
vector_rollcalls <- paste0(apply(matrix_rollcalls, 1, function(x) paste0(x, collapse=",")), collapse='#')

rollcalls$alpha_beta <- vector_rollcalls

rollcalls %<>%
  unnest(alpha_beta = strsplit(alpha_beta, "#")) %>%
  separate(alpha_beta, c("alpha","beta"), sep=",")

rollcalls$alpha <- as.numeric(rollcalls$alpha)
rollcalls$beta <- as.numeric(rollcalls$beta)

rollcalls$id_rollcalls <- rep(1:nr, sum(nl))

logit_f <- function(z) {
  return(1/(1+exp(-z)))
}
rollcalls$y <- rbinom(nrow(rollcalls), 1, logit_f(-rollcalls$alpha + rollcalls$beta * rollcalls$theta))
# Afterwards, we could test a second model: # \theta_i <- \gamma_i + \rho_i * \phi_k + \epsilon_i 

# Legislator 1 is left or right?
if ( rollcalls$gamma[1] < 0 )  prior_orient <- -1
if ( rollcalls$gamma[1] > 0 )    prior_orient <- 1

########

# Running first model without hierarchy

stan_data <- list(nrow=nrow(rollcalls), 
                  L=max(rollcalls$id_legis), 
                  R=max(rollcalls$id_rollcalls),
                  l=rollcalls$id_legis, r=rollcalls$id_rollcalls,
                  prior_orient=prior_orient,y=rollcalls$y)

stan_str <- "
data {
  int nrow;
  int L;
  int R;
  int<lower=0> l[nrow];
  int<lower=0> r[nrow];
  real prior_orient;
  int y[nrow];
}
parameters {
  vector[R] alpha;
  vector[R] beta;
  vector[L] theta;
}
model {
  for (i in 1:nrow) {
    y[i] ~ bernoulli_logit(-alpha[r[i]] - beta[r[i]] * theta[l[i]]);
  }
  alpha ~ normal(0,1);
  beta ~ normal(0,1);
  theta[1] ~ normal(prior_orient, 1);
  for (i in 2:L) {
    theta[i] ~ normal(0, 1);
  }
}
"
library(rstan)
fit1 <- stan(model_code = stan_str, data=stan_data, iter=3000, warmup=1500, thin=10, chains=1)

thetas <- extract(fit1, pars="theta")
thetas <- apply(thetas$theta, 2, mean)

thetas_originais <- rollcalls %>%
  select(theta)
thetas_originais <- thetas_originais[!duplicated(thetas_originais),]
gamas_originais <- rollcalls %>%
  select(gamma)
gamas_originais <- gamas_originais[!duplicated(gamas_originais),]

library(ggplot2)

theta_plot <- qplot(x=thetas, y=thetas_originais$theta, geom="point") + geom_smooth(method="lm")
gama_plot <- qplot(x=thetas, y=gamas_originais$gamma, geom="point") + geom_smooth(method="lm")

library(gridExtra)
grid.arrange(theta_plot, gama_plot)