Merge pull request #709 from stan-dev/feature/ARM

Fixes #707. Feature/arm models

Author: syclik

ARM/Ch.10/10.3_Matching.R

@@ -0,0 +1,13 @@
+## Computation of propensity score matches
+#FIXME: MISSING DATA
+
+ # starting point
+
+ # second model
+
+ # predicted values
+
+ # regression on matched data
+
+## Propensity score as a one-number summary (Figure 10.7)
+

ARM/Ch.10/10.4_LackOfOverlapWhenTreat.AssignmentIsUnknown.R

@@ -0,0 +1,122 @@
+library(rstan)
+library(ggplot2)
+
+ideo2 <- read.table ("ideo2.dat")
+candidate.effects <- read.table ("candidate_effects.dat", row.names=1)
+
+# Simple correction for incumbency advantage
+incadv <- function (years){
+  ifelse (years<46, .02,
+    ifelse (years<66, .02 + .08*(years-46)/(66-46), .10))
+}
+
+# Regions of the country
+region <- c(3,4,4,3,4,4,1,1,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)
+south <- region==3
+
+
+# add stuff to ideo2
+
+normvote <- rep (NA, nrow(ideo2))
+dum1 <- ideo2[,"dpvote"]
+dum2 <- ideo2[,"stgj"]
+dum3 <- (ideo2[,"year"]-44)/4
+for (i in 1:length(normvote)){
+  normvote[i] <- dum1[i] - candidate.effects[dum2[i],dum3[i]]
+}
+dvfix <- ideo2[,"dv"] - incadv(ideo2[,"year"])*ideo2[,"i2"]
+dvpfix <- ideo2[,"dvp"] - incadv(ideo2[,"year"]-2)*ideo2[,"i1"]
+ideo2 <- cbind(ideo2[,1:19], normvote, dvfix, dvpfix)
+names(ideo2) <- c(names(ideo2)[1:19], "normvote", "dvfix", "dvpfix")
+
+# impute normal vote from 2 years earlier:  ideo2
+
+for (year in seq(62,94,4)){
+  yr.cond <- ideo2[, "year"] == year
+  normvote <- rep (NA, sum(yr.cond))
+  dpvote <- rep (NA, sum(yr.cond))
+  stgj <- ideo2[yr.cond,"stgj"]
+  cd <- ideo2[yr.cond,"cd"]
+  indexes <- (1:nrow(ideo2))[yr.cond]
+  yr2.cond <- ideo2[, "year"] == year-2
+  data2 <- ideo2[yr2.cond,c("stgj","cd","normvote","dpvote")]
+  for (i in 1:sum(yr.cond)){
+    cond <- data2[,"stgj"]==stgj[i] & data2[,"cd"]==cd[i]
+    if (sum(cond)==1){
+      normvote[i] <- data2[cond,c("normvote")]
+      dpvote[i] <- data2[cond,c("dpvote")]
+    }
+  }
+  ideo2[yr.cond,c("normvote")] <- normvote
+  ideo2[yr.cond,c("dpvote")] <- dpvote
+}
+
+year <- 94
+yr.cond <- ideo2[, "year"] == year
+data <- ideo2[yr.cond,  ]
+deminc.cond <- (data[, "dvp"] > 0.5) & (abs(data[,"i2"]) == 1)
+repinc.cond <- (data[, "dvp"] < 0.5) & (abs(data[,"i2"]) == 1)
+#
+# fudge for 1992, 1994
+if (year>=92) data[,"occup"] <- rep(0,nrow(data))
+#
+dum <- apply(is.na(data),1,sum)
+ok <- dum==0 & !south[data[,"stgj"]]
+attach(data)
+
+ # Plot figure 10.8
+frame1 = data.frame(x1=1-dvp[deminc.cond],y1=score1[deminc.cond])
+frame2 = data.frame(x2=1-dvp[repinc.cond],y2=score1[repinc.cond])
+p1 <- ggplot() +
+      geom_point(data=frame1,aes(x=x1,y=y1),shape="x") +
+      geom_point(data=frame2,aes(x=x2,y=y2),shape="o") +
+      scale_y_continuous("(liberal)           ideology score       (conservative)") +
+     scale_x_continuous("Republican's vote share") +
+     theme_bw()
+print(p1)
+
+ # regression discontinuity analysis
+
+x <- 1 - dvp
+party <- ifelse (dvp<.5, 1, 0)
+
+## Regression in the area near the discontinuity (ideo_two_pred.stan)
+## lm (score1 ~ party + x, subset=overlap)
+
+overlap <- (deminc.cond | repinc.cond) & dvp>.45 & dvp<.55 & !is.na(score1+party+x)
+sc1 <- score1[overlap]
+p1 <- party[overlap]
+x1 <- x[overlap]
+dataList.1 <- list(N=length(sc1), score1=sc1,party=p1,x=x1)
+ideo_two_pred.sf1 <- stan(file='ideo_two_pred.stan', data=dataList.1,
+                          iter=1000, chains=4)
+print(ideo_two_pred.sf1)
+
+## Regression fit to all data (ideo_two_pred.stan)
+## lm (score1 ~ party + x, subset=incs)
+incs <- (deminc.cond | repinc.cond) & !is.na(score1+party+x)
+sc2 <- score1[incs]
+p2 <- party[incs]
+x2 <- x[incs]
+dataList.2 <- list(N=length(sc2), score1=sc2,party=p2,x=x2)
+ideo_two_pred.sf2 <- stan(file='ideo_two_pred.stan', data=dataList.2,
+                          iter=1000, chains=4)
+print(ideo_two_pred.sf2)
+
+## Regression with interactions (ideo_interactions.stan)
+## lm (score1 ~ party + x + party:x, subset=incs)
+
+ideo_interactions.sf1 <- stan(file='ideo_interactions.stan', data=dataList.2,
+                              iter=1000, chains=4)
+print(ideo_interactions.sf1)
+
+## Reparametrized regression (ideo_reparam.stan)
+## lm (score1 ~ party + I(z*(party==0)) + I(z*(party==1)), subset=incs)
+
+z <- x2 - 0.5
+z.1 <- I(z*(p2==0))
+z.2 <- I(z*(p2==1))
+dataList.3 <- list(N=length(sc2), score1=sc2,party=p2,z1=z.1,z2=z.2)
+ideo_reparam.sf1 <- stan(file='ideo_reparam.stan', data=dataList.3,
+                         iter=1000, chains=4)
+print(ideo_reparam.sf1)

ARM/Ch.10/10.5_CasualEffectsUsingIV.R

@@ -0,0 +1,37 @@
+library(rstan)
+library(ggplot2)
+library(foreign)
+
+sesame <- read.dta("sesame.dta")
+attach(sesame)
+
+## Rename variables of interest
+watched <- regular
+encouraged <- encour
+y <- postlet
+
+## Instrumental variables estimate (sesame_one_pred_a.stan)
+## lm (watched ~ encouraged)
+
+dataList.1 <- list(N=length(watched), watched=watched,encouraged=encouraged)
+sesame_one_pred_a.sf1 <- stan(file='sesame_one_pred_a.stan', data=dataList.1,
+                              iter=1000, chains=4)
+print(sesame_one_pred_a.sf1)
+
+beta.post <- extract(sesame_one_pred_a.sf1, "beta")$beta
+beta.mean1 <- colMeans(beta.post)
+
+## (sesame_one_pred_b.stan)
+## lm (y ~ encouraged)
+
+dataList.2 <- list(N=length(y), watched=y,encouraged=encouraged)
+sesame_one_pred_b.sf1 <- stan(file='sesame_one_pred_a.stan', data=dataList.2,
+                              iter=1000, chains=4)
+print(sesame_one_pred_b.sf1)
+
+beta.post <- extract(sesame_one_pred_b.sf1, "beta")$beta
+beta.mean2 <- colMeans(beta.post)
+
+
+iv.est.1 <- beta.mean2[2] / beta.mean1[2]
+print(iv.est.1)

ARM/Ch.10/10.6_IVinaRegressionFramework.R

@@ -0,0 +1,58 @@
+library(rstan)
+library(ggplot2)
+
+source("10.5_CasualEffectsUsingIV.R") # where data was cleaned
+
+## Rename variables of interest
+
+pretest <- prelet
+
+## 2 stage least squares (sesame_one_pred_a.stan)
+## lm (watched ~ encouraged)
+
+dataList.1 <- list(N=length(watched), watched=watched,encouraged=encouraged)
+sesame_one_pred_2a.sf1 <- stan(file='sesame_one_pred_a.stan', data=dataList.1,
+                               iter=1000, chains=4)
+print(sesame_one_pred_2a.sf1)
+beta.post <- extract(sesame_one_pred_2a.sf1, "beta")$beta
+beta.mean2a <- colMeans(beta.post)
+
+watched.hat <- beta.mean2a[1] + beta.mean2a[2] * encouraged
+
+## (sesame_one_pred_2b.stan)
+## lm (y ~ watched.hat)
+
+dataList.2 <- list(N=length(y), watched=y,encouraged=watched.hat)
+sesame_one_pred_2b.sf1 <- stan(file='sesame_one_pred_a.stan', data=dataList.2,
+                               iter=1000, chains=4)
+print(sesame_one_pred_2b.sf1)
+
+## Adjusting for covariates in a IV framework (sesame_multi_preds_3a.stan)
+## lm (watched ~ encouraged + pretest + as.factor(site) + setting)
+
+dataList.3 <- list(N=length(watched), watched=watched,encouraged=encouraged,pretest=pretest, site=site,setting=setting)
+sesame_multi_pred_3a.sf1 <- stan(file='sesame_multi_preds_3a.stan',
+                                 data=dataList.3,
+                                 iter=1000, chains=4)
+print(sesame_multi_pred_3a.sf1)
+
+beta.post <- extract(sesame_multi_pred_3a.sf1, "beta")$beta
+beta.mean3a <- colMeans(beta.post)
+
+watched.hat <- beta.mean3a[1] + beta.mean3a[2] * encouraged + beta.mean3a[3] * pretest + beta.mean3a[4] * (site==2) + beta.mean3a[5] * (site==3) + beta.mean3a[6] * (site==4) + beta.mean3a[7] * (site==5) + beta.mean3a[8] * setting
+
+## (sesame_multi_preds_3b.stan)
+## lm (y ~ watched.hat + pretest + as.factor(site) + setting)
+dataList.4 <- list(N=length(watched.hat), watched=y,encouraged=watched.hat,pretest=pretest, site=site,setting=setting)
+sesame_multi_pred_3b.sf1 <- stan(file='sesame_multi_preds_3b.stan',
+                                 data=dataList.4,
+                                 iter=1000, chains=4)
+print(sesame_multi_pred_3b.sf1)
+
+## Se for IV estimates (FIXME)
+
+## Performing 2sls automatically
+
+ # regression without pre-treatment variables
+
+ # regression controlling for pre-treatment variables

ARM/Ch.10/README

@@ -0,0 +1,62 @@
+Data
+====
+
+1. ideo_incs.data.R
+  - N     : number of observations
+  - party : party association? 1: Republican, 0: Democrat
+  - score1: effect of election results
+  - x     : pretreatment variable
+
+2. ideo_overlap.data.R
+  - N     : number of observations
+  - party : party association? 1: Republican, 0: Democrat
+  - score1: effect of election results
+  - x     : pretreatment variable
+
+3. ideo_reparam.data.R
+  - N     : number of observations
+  - party : party association? 1: Republican, 0: Democrat
+  - score1: effect of election results
+  - z1    : reparameterized term 1
+  - z2	: reparameterized term 2
+
+4. sesame.data.R
+  - N          : number of observations
+  - encouraged : encouraged to watch? 1: Yes, 0: No
+  - watched    : watched Sesame Street? 1: Yes, 0: No
+  - watched_hat: estimated casual effect of watching sesame street
+  - y          : post test scores
+
+5. sesame_one_pred_2b.data.R
+  - N          : number of observations
+  - encouraged : encouraged to watch? 1: Yes, 0: No
+  - watched    : watched Sesame Street? 1: Yes, 0: No
+  - watched_hat: estimated casual effect of watching sesame street
+  - y          : post test scores
+
+6. sesame_multi_preds_3b.data.R
+  - N          : number of observations
+  - encouraged : encouraged to watch? 1: Yes, 0: No
+  - pretest    : pre test scores
+  - setting    : setting category
+  - site       : site category
+  - watched    : watched Sesame Street? 1: Yes, 0: No
+  - watched_hat: estimated casual effect of watching sesame street
+  - y          : post test scores
+
+Models
+======
+
+1. One predictor
+  sesame_one_pred_2b.stan: lm(y ~ watched_hat)
+  sesame_one_pred_a.stan : lm(watched ~ encouraged)
+  sesame_one_pred_b.stan : lm(y ~ encouraged)
+
+3. Multiple predictors without interaction
+  ideo_reparam.stan        : lm(score1 ~ party + z1 + z2)
+  ideo_two_pred.stan       : lm(score1 ~ party + x)
+  sesame_mult_preds_3a.stan: lm(y ~ encouraged + pretest + factor(site) + setting)
+  sesame_mult_preds_3b.stan: lm(y ~ watched_hat + pretest + factor(site) + setting)
+
+3. Multiple predictors with interaction
+  ideo_interactions.stan: lm(score1 ~ party + x + party:x)

ARM/Ch.10/candidate_effects.dat

@@ -0,0 +1,50 @@
+AL -0.0074 0.0128 -0.0042 -0.0309 -0.1137 -0.0116 -0.0107 0.0601 0.0708 0.0125 0.0033 0.0595
+AK NA NA NA -0.011 0.0073 -0.0007 -0.0069 -0.0331 -0.0463 -0.0144 -0.0144 -0.0085
+AZ -0.006 -0.0035 -0.0029 0.0032 -0.0363 -0.0116 -0.0149 -0.0404 -0.0441 -0.0144 -0.0127 -0.0085
+AR 0.0074 0.0076 0.0068 -0.0319 -0.0233 -0.0046 -0.0092 0.0604 0.0794 0.0144 0.0035 0.094
+CA -0.0242 -0.0217 -0.0211 -0.0315 0.0053 -0.0368 -0.041 0.0057 -0.0277 -0.033 0.0022 0.0055
+CO -0.0074 -0.007 -0.0064 -0.0127 -0.0033 -0.0066 -0.012 -0.0364 -0.0396 -0.0093 -0.0073 -0.0006
+CT -0.006 -0.0035 -0.0028 0.0737 0.0399 0.0038 -0.0018 0.0058 0.0006 0.0015 0.001 0.003
+DE -0.0016 0.0009 0.0015 -0.014 0.0072 -0.0029 -0.0091 0.0025 0.0069 0.0094 0.0045 -0.0007
+FL -0.0001 0.0024 0.0031 -0.0223 -0.0363 -0.0116 -0.0149 0.0582 0.0665 0.0012 -0.0058 0.0565
+GA 0.0074 0.0076 0.0068 -0.0338 -0.1137 -0.0062 -0.0149 0.0898 0.1109 0.0144 0.0089 0.0596
+HI NA NA NA 0.0274 0.0159 0.0092 0.0105 0.0068 0.0041 0.011 0.0106 0.0085
+ID -0.0074 -0.0076 -0.0068 -0.0387 -0.0026 -0.0116 -0.0149 -0.0404 -0.0481 -0.0144 -0.0152 -0.0085
+IL -0.0043 0.0297 0.0303 0.012 0.0036 -0.0048 0.0222 0.0036 0.0026 0.0038 0.0069 0.0074
+IN -0.0074 -0.0062 -0.0056 -0.0222 -0.0039 -0.0108 -0.0149 -0.0056 -0.0042 -0.0037 -0.0202 -0.0241
+IA -0.0074 -0.0076 -0.0068 -0.0188 -0.0054 -0.0076 0.0181 0.0016 0.0035 0.0013 0.0086 0.0085
+KS -0.0074 -0.0076 -0.0068 -0.0256 -0.0144 -0.0116 0.0167 -0.0245 -0.0056 -0.0076 -0.0101 -0.0053
+KY 0.0224 0.0072 0.0068 -0.0191 0.0014 -0.0038 -0.0098 -0.0022 0.0062 0.0105 0.0043 -0.0005
+LA 0.0074 0.0076 0.0068 0.0287 -0.1137 -0.0074 -0.0089 0.0624 0.068 0.0079 0.002 0.0647
+ME -0.0074 -0.0076 -0.0068 0.0357 0.0212 0.0148 0.0021 0.0057 0.0024 0.0067 0.0028 -0.0008
+MD -0.0058 -0.0033 -0.0027 -0.006 0.0082 -0.0157 -0.0222 0.0024 0.0065 0.0126 0.0143 0.0085
+MA -0.001 0.0015 0.0021 0.1039 0.0495 0.0116 0.0324 0.0068 0.0144 0.0141 0.0464 0.0085
+MI -0.0058 -0.0033 -0.0027 0.0074 0.007 0.0012 0.0304 -0.0247 0.0013 0.0036 0.0037 0.0022
+MN -0.0009 0.0016 0.0022 0.0127 0.0254 0.0303 0.0323 0.0245 0.0321 0.0459 0.0152 0.0085
+MS 0.0074 0.0076 0.0068 -0.031 -0.1137 -0.0116 -0.0113 0.0582 0.0609 0.0089 0.002 0.058
+MO 0.0356 0.0066 0.0068 -0.0037 0.0092 -0.0014 -0.0067 0.0014 0.0049 0.0064 0.0032 0.0042
+MT -0.0072 -0.0047 -0.0041 -0.0024 0.0031 -0.0059 -0.0129 -0.0329 -0.0351 -0.0077 -0.0036 0.0023
+NE -0.0074 -0.0076 -0.0068 -0.0159 -0.0159 -0.0116 0.0167 -0.0068 -0.0144 -0.0144 -0.0152 -0.0085
+NV -0.0074 -0.0076 -0.0068 -0.0134 0.0063 -0.0046 -0.0129 -0.0356 -0.0346 -0.0133 -0.0138 -0.0085
+NH -0.0074 -0.0076 -0.0068 0.0555 0.0288 -0.0042 -0.0106 -0.0026 -0.0058 -0.0125 -0.0142 -0.0085
+NJ -0.0074 -0.0059 -0.0053 0.0281 0.0008 0.0001 -0.0065 0.0009 0.001 0.0006 0.0006 -0.0018
+NM -0.002 0.0005 0.0012 0.0333 0.0048 -0.0047 -0.0146 -0.0356 -0.0328 -0.0022 -0.0005 0.0034
+NY -0.0374 -0.0034 -0.0028 0.0205 -0.0115 0.0036 -0.0013 0.0065 0.007 0.0266 0.0112 0.0085
+NC 0.0074 0.0076 0.0068 -0.0351 -0.027 -0.0059 -0.0149 0.0582 0.0709 0.0114 0.0022 0.0609
+ND -0.0074 -0.0076 -0.0068 -0.0009 -0.0051 -0.0096 0.0167 -0.0041 -0.0354 -0.0144 -0.0125 -0.0044
+OH -0.0064 -0.0039 -0.0033 -0.0072 -0.0018 -0.0046 -0.0077 0.0022 0.0033 0.0031 0.0027 0.0005
+OK 0.0021 0.0046 0.0052 -0.0331 -0.0063 -0.0116 -0.0149 -0.0068 -0.0041 -0.0054 -0.0124 -0.0085
+OR -0.0074 -0.0076 -0.0068 -0.0216 0.0023 -0.0034 -0.0094 0.0063 0.0056 0.0031 0.0067 0.0085
+PA -0.0033 -0.0008 -0.0002 0.02 0.0067 0.0001 -0.0036 0.0047 0.0052 0.0058 0.0103 0.0085
+RI 0.0038 0.0063 0.0068 0.1113 0.0495 0.0116 0.0149 0.0068 0.0139 0.0144 0.0152 0.0085
+SC 0.0074 0.0076 0.0068 -0.0341 -0.1137 -0.0116 -0.0149 0.0582 0.0716 0.0125 -0 0.0565
+SD -0.0074 -0.0076 -0.0068 -0.0138 -0.0069 -0.0116 0.0482 0.0058 0.0052 -0.0078 -0.0069 0.0009
+TN 0.0042 0.0067 0.0245 -0.0326 -0.0354 -0.0103 -0.0149 0.0582 0.0727 0.013 0.0071 0.0807
+TX 0.0074 -0.0239 -0.0247 0.0148 -0.0005 -0.0013 -0.0054 0.0633 0.05 -0.0152 -0.0169 0.0305
+UT -0.0067 -0.0042 -0.0035 -0.042 -0.0058 -0.0116 -0.0149 -0.0404 -0.0481 -0.0144 -0.0152 -0.0085
+VT -0.0074 -0.0076 -0.0068 0.0459 0.0184 -0.0058 -0.0102 -0.0016 -0.0051 0.0016 0.0047 0.0053
+VA -0.0043 -0.0018 -0.0011 -0.0356 -0.039 -0.011 -0.0149 0.0582 0.0636 0.0014 -0.0016 0.058
+WA -0.0019 0.0006 0.0013 -0.0161 0.0048 -0.0036 -0.0055 0.0048 0.0012 0.0003 0.0053 0.0085
+WV 0.0067 0.0076 0.0068 -0.0268 0.0103 0.0032 0.0004 0.0031 0.0122 0.0144 0.0118 0.0085
+WI -0.0074 -0.0076 -0.0068 0.0144 -0.0003 -0.0041 0.0228 0.0068 0.0069 0.0068 0.0099 0.0085
+WY -0.0074 -0.0076 -0.0068 -0.022 -0.0036 -0.0103 -0.0149 -0.0404 -0.0467 -0.0144 -0.0152 -0.0085