added rest of ch17 prelim models

Author: PeterLi2016

ARM/Ch.17/17.4_multilevel_logistic.stan

@@ -11,6 +11,7 @@ data {
   int<lower=0, upper=n_edu> edu[N];
   int<lower=0, upper=n_state> region[N];
   int<lower=0, upper=n_state> state[N];
+  int<lower=0, upper=1> y[N];
 }
 parameters {
   real<lower=0> sigma;

ARM/Ch.17/17.7_latent_glm.stan

@@ -0,0 +1,80 @@
+data {
+  int<lower=0> N;
+  int<lower=0> n_age;
+  int<lower=0> n_edu;
+  int<lower=0> n_region;
+  int<lower=0> n_state;
+
+  int<lower=0, upper=1> female[N];
+  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> state[N];
+  int<lower=0, upper=1> y[N];
+  vector z[N]
+}
+parameters {
+  real<lower=0> sigma;
+  real<lower=0> sigma_age;
+  real<lower=0> sigma_edu;
+  real<lower=0> sigma_state;
+  real<lower=0> sigma_region;
+  real<lower=0> sigma_age_edu;
+
+  real b_0;
+  real b_female;
+  real b_black;
+  real b_female_black;
+
+  real b_v_prev;
+
+  vector[n_age] b_age;
+  vector[n_edu] b_edu;
+  vector[n_region] b_region;
+  matrix[n_age,n_edu] b_age_edu;
+  vector[n_state] b_state;
+}
+model {
+  vector[N] p;
+  vector[N] z_lo;
+  vector[N] z_hi;
+  vector[N] Xbeta;
+  vector[n_state] b_state_hat;
+
+  b_0 ~ normal(0, 100);
+  b_female ~ normal(0, 100);
+  b_black ~ normal(0, 100);
+  b_female_black ~ normal(0, 100);
+
+  b_age ~ normal(0, sigma_age);
+  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);
+
+  b_v_prev ~ normal(0, 100);
+
+  for (j in 1:n_state)
+    b_state_hat[j] <- b_region[region[j]] + b_v_prev * v_prev[j];
+
+  b_state ~ normal(b_state_hat, sigma_state);
+
+  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);
+
+  for (i in 1:N) {
+    z_lo[i] <- 100 * equals(y[i],0);
+    z_hi[i] <- 100 * equals(y[i],1);
+    z[i] ~ logistic(Xbeta[i],1) T[z_lo[i],z_hi[i]];  
+  }
+
+  
+}

ARM/Ch.17/17.7_robit.stan

@@ -0,0 +1,84 @@
+data {
+  int<lower=0> N;
+  int<lower=0> n_age;
+  int<lower=0> n_edu;
+  int<lower=0> n_region;
+  int<lower=0> n_state;
+
+  int<lower=0, upper=1> female[N];
+  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> state[N];
+  int<lower=0, upper=1> y[N];
+  vector z[N]
+}
+parameters {
+  real<lower=0> sigma_age;
+  real<lower=0> sigma_edu;
+  real<lower=0> sigma_state;
+  real<lower=0> sigma_region;
+  real<lower=0> sigma_age_edu;
+
+  real b_0;
+  real b_female;
+  real b_black;
+  real b_female_black;
+
+  real b_v_prev;
+  real df_inv;
+  real df;
+  real<lower=0> sigma_z; 
+  vector[n_age] b_age;
+  vector[n_edu] b_edu;
+  vector[n_region] b_region;
+  matrix[n_age,n_edu] b_age_edu;
+  vector[n_state] b_state;
+}
+model {
+  vector[N] p;
+  vector[N] z_lo;
+  vector[N] z_hi;
+  vector[N] Xbeta;
+  vector[n_state] b_state_hat;
+
+  b_0 ~ normal(0, 100);
+  b_female ~ normal(0, 100);
+  b_black ~ normal(0, 100);
+  b_female_black ~ normal(0, 100);
+
+  b_age ~ normal(0, sigma_age);
+  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);
+
+  b_v_prev ~ normal(0, 100);
+
+  for (j in 1:n_state)
+    b_state_hat[j] <- b_region[region[j]] + b_v_prev * v_prev[j];
+
+  b_state ~ normal(b_state_hat, sigma_state);
+
+  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);
+
+  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]];  
+  }
+  
+}