##Alternative pscore models estimated and matched using Imai's R program (with programming errors corrected by us).


psboot <- function(y, treat, pscores, sims=500, n.match=3, grp=NULL){ 
  uG<-unique(grp)
  nG<-length(uG)
  ate <- matrix(0,nrow=sims,ncol=nG+1)
  data<-cbind(y,treat,pscores,grp)
  labels <- seq(1:nrow(data))
  for(i in 1:sims){
    bootrows <- sample(labels,nrow(data),replace=TRUE)
    bootsample <- data[bootrows,]
    y.sample <- bootsample[,1]
    treat.sample <- bootsample[,2]
    grp.sample <- bootsample[,4]
    match <- psmatch(treat.sample, bootsample[,3],
                     n.match=n.match)
    diff <- y.sample[treat.sample==1]-y.sample[c(match$match.treat)]
    ate[i,1] <- mean(diff)
    for(j in 1:nG)
      ate[i,j+1] <- mean(diff[grp.sample[treat.sample==1]==uG[j]])
    if(i==sims){
      cat(i,"iter.(ate=",mean(ate[1:i,1]),"sd=",sd(ate[1:i,1]),")...\n")
      for(j in 1:nG)
        cat("**",uG[j],"grp.(ate=",mean(ate[1:i,j+1]),"sd=",sd(ate[1:i,j+1]),")...\n")
    }
  }
}

psmatch <- function(treat, pscores, n.match=3){
  n <- length(treat)
  n1 <- length(treat[treat==1])
  labels <- seq(1,n)
  clabels <- labels[treat==0]
  pscore1 <- pscores[treat==1]
  pscore0 <- pscores[treat==0]
  matchedT <- matrix(0, nrow=length(pscore1), ncol=n.match)
  p0<-pscore0
  for(i in 1:n1){
    needed<-n.match
    matched<-NULL
    for(j in 1:n.match)
      if(needed>0){
        deviation <- abs(p0-pscore1[i])
        mindev <- clabels[min(deviation)==deviation]
        if(length(mindev)>needed)
          temp <- sample(mindev,needed,replace=F)
        else 
          temp <- mindev
        for(k in 1:n.match)
          p0[clabels==temp[k]] <- 1000
        matched<-c(matched, temp)
        needed<-n.match-length(matched)
      }
    matchedT[i,] <- matched
  }
  list(match.treat = matchedT, match.pscore=pscores)
}

data2<-read.table("C:/Documents and Settings/Kevin Arceneaux/Desktop/apsr_corrected_latest.dat",header=T)
data2$WARD<-as.factor(data2$WARD)
data2$NEW<-(data2$VOTE96.0==0 & data2$VOTE96.1==0)*1
data2$phone<-((data2$PHONEGRP+data2$MAILCALL)>0)*1
data2$AGE2<-data2$AGE^2/100
D<-(data2$APPEAL>1 & data2$PHN.C1==1 & data2$MAILCALL==0 &
    data2$PERSNGRP==0 & data2$MAILINGS==0)*1 
C<-(data2$PHONEGRP==0 & data2$MAILCALL==0 & data2$PERSNGRP==0 &
    data2$MAILINGS==0)*1
dataP<-subset(data2, D+C>0)

#Linear model
pscore.glm<-glm(PHN.C1 ~ PERSONS + VOTE96.1 + NEW + MAJORPTY + AGE +
                WARD, 
                family=binomial(logit), data=dataP)
D<-dataP$PHN.C1
Y<-dataP$VOTED98
grp<-dataP$PERSONS
psboot(Y, D, fitted(pscore.glm), sims=500, n.match=5, grp=grp) 


#Everything Interacted model
pscore.glm<-glm(PHN.C1 ~ PERSONS + VOTE96.1 + NEW + MAJORPTY + AGE +
                WARD + PERSONS:VOTE96.1 + PERSONS:NEW + PERSONS:MAJORPTY +
		PERSONS:AGE + PERSONS:WARD + VOTE96.1:MAJORPTY + VOTE96.1:AGE +
		VOTE96.1:WARD + NEW:MAJORPTY + NEW:AGE + NEW:WARD + MAJORPTY:AGE +
		MAJORPTY:WARD + AGE:WARD, 
                family=binomial(logit), data=dataP)
D<-dataP$PHN.C1
Y<-dataP$VOTED98
grp<-dataP$PERSONS
psboot(Y, D, fitted(pscore.glm), sims=500, n.match=5, grp=grp)


#NXA added to Imai's model
pscore.glm<-glm(PHN.C1 ~ PERSONS + VOTE96.1 + NEW + MAJORPTY + AGE +
                WARD + PERSONS:VOTE96.1 + PERSONS:NEW + NEW:AGE + AGE2, 
                family=binomial(logit), data=dataP)
psboot(Y, D, fitted(pscore.glm), sims=500, n.match=5, grp=grp)