#IDSS analysis

library(foreign)
library(car)
library(memisc)
library(effects)
library(stargazer)

load("IDSS.jeps.rep.Rdata")

#TABLE 1, FIGURES 1 AND 2

#All samples combined
model1 <- lm(welfexp ~ welfrace * welfimm, data=IDSS[IDSS$citborn==1,])
linearHypothesis(model1, "welfrace1 = welfimm1") #for footnote 3
eff1 <- effect("welfrace * welfimm", model1, typical="mean")
eff1r <- cbind(eff1$x,fit=eff1$fit,lower=eff1$lower,higher=eff1$upper)
eff1r <- eff1r[order(eff1r$welfrace,eff1r$welfimm),]
rownames(eff1r) <- NULL
eff1r$treat <- as.numeric(rownames(eff1r)) 

#Canada, ROC
model2 <- lm(welfexp ~ welfrace * welfimm, data=IDSS[IDSS$sample=="Canada, ROC" & IDSS$citborn==1,])
linearHypothesis(model2, "welfrace1 = welfimm1") #for footnote 3
eff2 <- effect("welfrace * welfimm", model2, typical="mean")
eff2r <- cbind(eff2$x,fit=eff2$fit,lower=eff2$lower,higher=eff2$upper)
eff2r <- eff2r[order(eff2r$welfrace,eff2r$welfimm),]
rownames(eff2r) <- NULL
eff2r$treat <- as.numeric(rownames(eff2r)) 

#Canada, QC
model3 <- lm(welfexp ~ welfrace * welfimm, data=IDSS[IDSS$sample=="Canada, QC"  & IDSS$citborn==1,])
linearHypothesis(model3, "welfrace1 = welfimm1") #for footnote 3
eff3 <- effect("welfrace * welfimm", model3, typical="mean")
eff3r <- cbind(eff3$x,fit=eff3$fit,lower=eff3$lower,higher=eff3$upper)
eff3r <- eff3r[order(eff3r$welfrace,eff3r$welfimm),]
rownames(eff3r) <- NULL
eff3r$treat <- as.numeric(rownames(eff3r)) 

#US
model4 <- lm(welfexp ~ welfrace * welfimm, data=IDSS[IDSS$sample=="US"  & IDSS$citborn==1,])
linearHypothesis(model4, "welfrace1 = welfimm1") #for footnote 3
eff4 <- effect("welfrace * welfimm", model4, typical="mean")
eff4r <- cbind(eff4$x,fit=eff4$fit,lower=eff4$lower,higher=eff4$upper)
eff4r <- eff4r[order(eff4r$welfrace,eff4r$welfimm),]
rownames(eff4r) <- NULL
eff4r$treat <- as.numeric(rownames(eff4r)) 

#table 1.
stargazer(model1,model2,model3,model4,type="text",no.space=T)
stargazer(model1,model2,model3,model4,type="html",no.space=T,out="table1.html")

#figure 2.
pdf("Figure2.jeps.pdf", width = 6, height = 7)
par(mar=c(5, 5, 1, 2) + 0.1)
plot(eff1r$treat,eff1r$fit, type="p", pch=15, cex=1.5, ylim=c(.55,.85),xlim=c(.6,4.4),axes=F,ann=F)
arrows(eff1r$treat,eff1r$low,eff1r$treat,eff1r$high, length=.05, code=3, angle=90, col="black",lwd=2,lty=3)
axis(side=2, at=c(.55,.6,.65,.7,.75,.8,.85), las=1, cex.axis=.8, col="gray")
axis(side=4, at=c(.55,.6,.65,.7,.75,.8,.85), labels=F, col="gray")
mtext(side=2,"Support For Assistance",cex=1,line=3)
axis(side=1, at=c(1,2,3,4), labels=c("John\nNative","John\nImmigrant","Sulaiman\nNative","Sulaiman\nImmigrant"), las=1, cex.axis=1, col="gray", padj=.8)
dev.off()

#appendix table 4
model7 <- lm(welfexp ~ welfrace * welfimm * therm_muslim, data=IDSS[IDSS$citborn==1,])
stargazer(model7, type="text")

#figure 3
eff7 <- effect("welfrace * welfimm * therm_muslim", model7, typical="mean")
eff7r <- cbind(eff7$x,fit=eff7$fit,lower=eff7$lower,higher=eff7$upper)
eff7r <- eff7r[eff7r$therm_muslim==0 | eff7r$therm_muslim==100, ]
eff7r <- eff7r[order(eff7r$therm_muslim,eff7r$welfrace,eff7r$welfimm),]
rownames(eff7r) <- NULL
eff7r$treat <- c(6,7,8,9,1,2,3,4) 
pdf("Figure3.jeps.pdf", width = 9, height = 7)
par(mar=c(5, 5, 1, 2) + 0.1)
plot(eff7r$treat[eff7r$therm_muslim==0],eff7r$fit[eff7r$therm_muslim==0], type="p", pch=15, cex=1.5, ylim=c(.4,.9),xlim=c(0,10),axes=F,ann=F,col="red")
arrows(eff7r$treat[eff7r$therm_muslim==0],eff7r$low[eff7r$therm_muslim==0],eff7r$treat[eff7r$therm_muslim==0],eff7r$high[eff7r$therm_muslim==0], length=.05, code=3, angle=90, col="red",lwd=2,lty=3)
points(eff7r$treat[eff7r$therm_muslim==100],eff7r$fit[eff7r$therm_muslim==100], type="p", pch=15, cex=1.5,col="blue")
arrows(eff7r$treat[eff7r$therm_muslim==100],eff7r$low[eff7r$therm_muslim==100],eff7r$treat[eff7r$therm_muslim==100],eff7r$high[eff7r$therm_muslim==100], length=.05, code=3, angle=90, col="blue",lwd=2,lty=3)
axis(side=2, at=c(.4,.5,.6,.7,.8,.9), las=1, cex.axis=.8, col="gray")
axis(side=4, at=c(.4,.5,.6,.7,.8,.9), labels=F, col="gray")
mtext(side=2,"Support For Assistance",cex=1,line=3)
axis(side=1, at=c(1,2,3,4), labels=c("John\nNat","John\nImm","Sulmn\nNat","Sulmn\nImm"), las=1, cex.axis=1, col="gray", padj=.8)
axis(side=1, at=c(6,7,8,9), labels=c("John\nNat","John\nImm","Sulmn\nNat","Sulmn\nImm"), las=1, cex.axis=1, col="gray", padj=.8)
text(1,.4,pos=4,"Therm=100",col="blue")
text(6,.4,pos=4,"Therm=0",col="red")
dev.off()