##Original code written by Imai. We modified his code to conduct matching on the actual sample rater than ##the bootstrapped sample. psboot <- function(y, treat, pscores, sims=1, n.match=3, grp=NULL){ uG<-unique(grp) nG<-length(uG) ate <- matrix(0,nrow=sims,ncol=nG+1) data<-cbind(y,treat,pscores) labels <- seq(1:nrow(data)) for(i in 1:sims){ #Do not bootstrap bootrows <- labels bootsample <- data[bootrows,] y.sample <- bootsample[,1] treat.sample <- bootsample[,2] 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[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) } #Phone data2<-read.table("C:/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) pscore.glm<-glm(PHN.C1 ~ PERSONS + VOTE96.1 + NEW + MAJORPTY + AGE + WARD + PERSONS:VOTE96.1 + PERSONS:NEW + AGE2, family=binomial(logit), data=dataP) D<-dataP$PHN.C1 Y<-dataP$VOTED98 grp<-dataP$PERSONS #sims refers to the number of times the actual sample is matched (it does not refer to the bootstrap sample anymore) psboot(Y, D, fitted(pscore.glm), sims=100, n.match=5, grp=grp)