updated rest of Ch17 ARM models with R files

Author: PeterLi2016

ARM/Ch.17/17.4_multilevel_logistic.stan

@@ -9,9 +9,10 @@ data {
   int<lower=0, upper=1> black[N];
   int<lower=0, upper=n_age> age[N];
   int<lower=0, upper=n_edu> edu[N];
-  int<lower=0, upper=n_state> region[N];
+  int<lower=0, upper=n_state> region[n_state];
   int<lower=0, upper=n_state> state[N];
   int<lower=0, upper=1> y[N];
+  vector[n_state] v_prev;
 }
 parameters {
   real<lower=0> sigma;
@@ -48,8 +49,10 @@ model {
   b_edu ~ normal(0, sigma_edu);
   b_region ~ normal(0, sigma_region);
 
-  for (j in 1:n_age)
-    b_age_edu[j,] ~ normal(0, sigma_age_edu);
+  for (j in 1:n_age) {
+    for (i in 1:n_edu)
+      b_age_edu[j,i] ~ normal(0, sigma_age_edu);
+  }
 
   b_v_prev ~ normal(0, 100);
 
@@ -59,10 +62,10 @@ model {
   b_hat ~ normal(b_state_hat, sigma_state);
 
   for (i in 1:N)
-    p[i] <- max(0, min(1, inv_logit(b_0 + b_female*female[i] 
+    p[i] <- fmax(0, fmin(1, inv_logit(b_0 + b_female*female[i] 
       + b_black*black[i] + b_female_black*female[i]*black[i] +
       b_age[age[i]] + b_edu[edu[i]] + b_age_edu[age[i],edu[i]] +
-      b_state[state[i]])));
+      b_hat[state[i]])));
 
-  y ~ binomial(1, p);
+  y ~ bernoulli(p);
 }

ARM/Ch.17/17.4_multilevel_logistic_regression.R

@@ -0,0 +1,43 @@
+library(rstan)
+
+polls.subset <- read.table ("polls.subset.dat")
+attach (polls.subset)
+
+ # Load in data for region indicators
+ # Use "state", an R data file (type ?state from the R command window for info)
+ #
+ # Regions:  1=northeast, 2=south, 3=north central, 4=west, 5=d.c.
+ # We have to insert d.c. (it is the 9th "state" in alphabetical order)
+
+data (state)                  # "state" is an R data file
+state.abbr <- c (state.abb[1:8], "DC", state.abb[9:50])
+dc <- 9
+not.dc <- c(1:8,10:51)
+region <- c(3,4,4,3,4,4,1,1,5,3,3,4,4,2,2,2,2,3,3,1,1,1,2,2,3,2,4,2,4,1,1,4,1,3,2,2,3,4,1,1,3,2,3,3,4,1,3,4,1,2,4)
+
+ # define other data summaries
+
+y <- bush                  # 1 if support bush, 0 if support dukakis
+n <- length(y)             # of survey respondents
+n.age <- max(age)          # of age categories
+n.edu <- max(edu)          # of education categories
+n.state <- max(state)      # of states
+n.region <- max(region)    # of regions
+
+ # also include a measure of previous vote as a state-level predictor
+
+library (foreign)
+presvote <- read.dta ("presvote.dta")
+attach (presvote)
+v.prev <- presvote$g76_84pr
+age.edu <- n.edu*(age-1) + edu
+
+ok <- !is.na(female+black+age+edu+state+y)
+ # election model
+dataList.1 <- list(N=length(y[ok]), n_age=n.age, n_edu=n.edu, n_region=n.region, n_state=n.state,
+                   female=female[ok], black=black[ok],age=age[ok], edu=edu[ok],
+                   region=region, state=state[ok],y=y[ok],v_prev=v.prev)
+multilevel_logistic.sf1 <- stan(file='17.4_multilevel_logistic.stan', data=dataList.1,
+                            iter=1000, chains=4)
+print(multilevel_logistic.sf1)
+

ARM/Ch.17/17.5_multilevel_poisson.stan

@@ -5,8 +5,8 @@ data {
 
   int<lower=0, upper=n_precint> precint[N];
   int<lower=0, upper=n_eth> eth[N];
-  vector offeset[N];
-  vector stops[N];
+  vector[N] offeset;
+  int<lower=0> stops[N];
 }
 parameters {
   real mu;
@@ -22,6 +22,7 @@ model {
   real mu_adj;
   vector[n_eth] b_eth_adj;
   vector[n_precint] b_precint_adj;
+  vector[N] lambda;
 
   mu ~ normal(0, 100);
   mu_adj <- mu + mean(b_eth) + mean(b_precint);

ARM/Ch.17/17.6_multilevel_ordered_categorical.stan

@@ -4,7 +4,8 @@ data {
   int<lower=0> n_player;
 
   int<lower=0, upper=n_player> player[N];
-  vector y[N];
+  int<lower=0, upper=n_cut> y[N];
+  vector[N] x;
 }
 parameters {
   vector[n_cut] mu_c;
@@ -41,6 +42,6 @@ model {
   }
 
   for (i in 1:N)
-    y[i] ~ categorical(P[i,]);
+    y[i] ~ categorical(transpose(row(P,i)));
 
 }

ARM/Ch.17/17.7_latent_glm.stan

@@ -9,10 +9,11 @@ data {
   int<lower=0, upper=1> black[N];
   int<lower=0, upper=n_age> age[N];
   int<lower=0, upper=n_edu> edu[N];
-  int<lower=0, upper=n_state> region[N];
+  int<lower=0, upper=n_state> region[n_state];
   int<lower=0, upper=n_state> state[N];
   int<lower=0, upper=1> y[N];
-  vector z[N]
+  vector[n_state] v_prev;
+  vector[N] z;
 }
 parameters {
   real<lower=0> sigma;
@@ -51,8 +52,10 @@ model {
   b_edu ~ normal(0, sigma_edu);
   b_region ~ normal(0, sigma_region);
 
-  for (j in 1:n_age)
-    b_age_edu[j,] ~ normal(0, sigma_age_edu);
+  for (j in 1:n_age) {
+    for (i in 1:n_edu)
+      b_age_edu[j,i] ~ normal(0, sigma_age_edu);
+  }
 
   b_v_prev ~ normal(0, 100);
 
@@ -61,18 +64,19 @@ model {
 
   b_state ~ normal(b_state_hat, sigma_state);
 
-  for (i in 1:N)
+  for (i in 1:N) {
     Xbeta[i] <- b_0 + b_female*female[i] 
       + b_black*black[i] + b_female_black*female[i]*black[i] +
       b_age[age[i]] + b_edu[edu[i]] + b_age_edu[age[i],edu[i]] +
-      b_state[state[i]]
-    p[i] <- max(0, min(1, inv_logit(Xbeta[i])));
-
-  y ~ binomial(1, p);
+      b_state[state[i]];
+    p[i] <- fmax(0, fmin(1, inv_logit(Xbeta[i])));
+  }
+  
+  y ~ bernoulli(p);
 
   for (i in 1:N) {
-    z_lo[i] <- 100 * equals(y[i],0);
-    z_hi[i] <- 100 * equals(y[i],1);
+    z_lo[i] <- 100 * (y[i] == 0);
+    z_hi[i] <- 100 * (y[i] == 1);
     z[i] ~ logistic(Xbeta[i],1) T[z_lo[i],z_hi[i]];  
   }
 

ARM/Ch.17/17.7_robit.stan

@@ -9,10 +9,11 @@ data {
   int<lower=0, upper=1> black[N];
   int<lower=0, upper=n_age> age[N];
   int<lower=0, upper=n_edu> edu[N];
-  int<lower=0, upper=n_state> region[N];
+  int<lower=0, upper=n_state> region[n_state];
   int<lower=0, upper=n_state> state[N];
   int<lower=0, upper=1> y[N];
-  vector z[N]
+  vector[N] z;
+  vector[N] v_prev;
 }
 parameters {
   real<lower=0> sigma_age;
@@ -27,9 +28,7 @@ parameters {
   real b_female_black;
 
   real b_v_prev;
-  real df_inv;
-  real df;
-  real<lower=0> sigma_z; 
+  real<lower=0> df_inv;
   vector[n_age] b_age;
   vector[n_edu] b_edu;
   vector[n_region] b_region;
@@ -42,6 +41,8 @@ model {
   vector[N] z_hi;
   vector[N] Xbeta;
   vector[n_state] b_state_hat;
+  real df;
+  real sigma_z; 
 
   b_0 ~ normal(0, 100);
   b_female ~ normal(0, 100);
@@ -52,8 +53,10 @@ model {
   b_edu ~ normal(0, sigma_edu);
   b_region ~ normal(0, sigma_region);
 
-  for (j in 1:n_age)
-    b_age_edu[j,] ~ normal(0, sigma_age_edu);
+  for (j in 1:n_age) {
+    for(i in 1:n_edu)
+      b_age_edu[j,i] ~ normal(0, sigma_age_edu);
+  }
 
   b_v_prev ~ normal(0, 100);
 
@@ -62,23 +65,24 @@ model {
 
   b_state ~ normal(b_state_hat, sigma_state);
 
-  for (i in 1:N)
+  for (i in 1:N) {
     Xbeta[i] <- b_0 + b_female*female[i] 
       + b_black*black[i] + b_female_black*female[i]*black[i] +
       b_age[age[i]] + b_edu[edu[i]] + b_age_edu[age[i],edu[i]] +
-      b_state[state[i]]
-    p[i] <- max(0, min(1, inv_logit(Xbeta[i])));
+      b_state[state[i]];
+    p[i] <- fmax(0, fmin(1, inv_logit(Xbeta[i])));
+  }
 
-  y ~ binomial(1, p);
+  y ~ bernoulli(p);
 
   df_inv ~ uniform(0, 0.5);
   df <- 1 / df_inv;
   sigma_z <- sqrt((df-2)/df);
 
   for (i in 1:N) {
-    z_lo[i] <- 100 * equals(y[i],0);
-    z_hi[i] <- 100 * equals(y[i],1);
-    z[i] ~ t(df,Xbeta[i],sigma_z) T[z_lo[i],z_hi[i]];  
+    z_lo[i] <- 100 * (y[i] == 0);
+    z_hi[i] <- 100 * (y[i] == 1);
+    z[i] ~ student_t(df,Xbeta[i],sigma_z) T[z_lo[i],z_hi[i]];  
   }
 
 }