names(big.est.store)[1] <- paste("nsim=", nsim, "; CtoT ratio=", r, sep = "")
# summarize mse
store.mse <- matrix(NA,3,length(est))
colnames(store.mse) <- est
rownames(store.mse) <- paste(rep(c("Y1","Y2","Y3")))
p <- 1
for(k in 1:3){
store.mse[p,] <- apply((big.est.store[[length(big.para.store)]][[k]]-0)^2,2,mean,na.rm=T)
p <- p + 1
}
mse <- round(store.mse*100,1)
# summarize bias
store.bias <- matrix(NA,3,length(est))
colnames(store.bias) <- est
rownames(store.bias) <- paste(rep(c("Y1","Y2","Y3")))
p <- 1
for(k in 1:3){
store.bias[p,] <- apply(big.est.store[[length(big.para.store)]][[k]]-0,2,mean,na.rm=T)
p <- p + 1
}
bias <- round(store.bias*100,1)
cat("Simulation with no. of obs: ", nsim, "CtoT ratio of",r,"\n")
print(bias)
print(mse)
names(est.store) <- c("Y1","Y2","Y3")
big.est.store[[1]]  <- est.store
names(big.est.store)[1] <- paste("nsim=", nsim, "; CtoT ratio=", r, sep = "")
# summarize mse
store.mse <- matrix(NA,3,length(est))
colnames(store.mse) <- est
rownames(store.mse) <- paste(rep(c("Y1","Y2","Y3")))
p <- 1
for(k in 1:3){
store.mse[p,] <- apply((big.est.store[[length(big.est.store)]][[k]]-0)^2,2,mean,na.rm=T)
p <- p + 1
}
mse <- round(store.mse*100,1)
# summarize bias
store.bias <- matrix(NA,3,length(est))
colnames(store.bias) <- est
rownames(store.bias) <- paste(rep(c("Y1","Y2","Y3")))
p <- 1
for(k in 1:3){
store.bias[p,] <- apply(big.est.store[[length(big.est.store)]][[k]]-0,2,mean,na.rm=T)
p <- p + 1
}
bias <- round(store.bias*100,1)
cat("Simulation with no. of obs: ", nsim, "CtoT ratio of",r,"\n")
print(bias)
print(mse)
rm(list=ls())
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
n
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=2, ds=FALSE, exclude=list(c("educ", "nodegr")), print.level = -1)
summary(att(out.eb, method="lm_robust"))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=3, exclude=list(c("educ", "nodegr")))
summary(att(out.eb, method="lm_robust"))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=3, ds=FALSE, exclude=list(c("educ", "nodegr")), print.level = 1)
summary(att(out.eb, method="lm_robust"))
X.sq <- hbal::covarExpand(X, treatment=rep(0, nrow(X)), expand.degree=2)$mat
hbal::covarExpand
?hbal::covarExpand
X.sq <- hbal::covarExpand(X, treatment=rep(0, nrow(X)), exp.degree=2)$mat
hbal::covarExpand(X, treatment=rep(0, nrow(X)), exp.degree=2)
hbal::covarExpand(X, treatment=rep(0, nrow(X)))
hbal::covarExpand(as.matrix(X), treatment=rep(0, nrow(X)), exp.degree=2)
X.sq <- hbal::covarExpand(as.matrix(X), treatment=rep(0, nrow(X)), exp.degree=2)$mat
library(ebal)
?ebal
?ebalanace
?ebalance
X.sq <- hbal::covarExpand(as.matrix(X), treatment=rep(0, nrow(X)), exp.degree=2, exclude=list(c("educ", "nodegr")))$mat
out.eb <- ebalance(Treatment=treatment, X=X.sq)
rm(list=ls())
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
n
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=2, exclude=list(c("educ", "nodegr")))
summary(att(out.eb, method="lm_robust"))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=2, Y=Y, exclude=list(c("educ", "nodegr")), shuffle.treat=FALSE, print.level = -1)
summary(att(out.eb))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=3, Y=Y, exclude=list(c("educ", "nodegr")), shuffle.treat=FALSE, print.level = -1)
summary(att(out.eb))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=2, exclude=list(c("educ", "nodegr")))
summary(att(out.eb, method="lm_robust"))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=3, exclude=list(c("educ", "nodegr")), print.level = 1)
summary(att(out.eb, method="lm_robust"))
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=3, exclude=list(c("educ", "nodegr")), print.level = 2)
summary(att(out.eb, method="lm_robust"))
out.eb$weights
out.eb$weights[2000:2020]
treatment
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=3, exclude=list(c("educ", "nodegr")), print.level = -1, max.iterations = 500)
summary(att(out.eb, method="lm_robust"))
?hbal
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
n
#dat <- cbind(Y, treatment, X)
X.sq <- hbal::covarExpand(as.matrix(X), treatment=rep(0, nrow(X)), exp.degree=2, exclude=list(c("educ", "nodegr")))$mat
rank(X.sq)
qr(X.sq)
out.eb <- ebalance(Treatment=treatment, X=X.sq, print.level=2)
?ebalance
X.sq <- hbal::covarExpand(as.matrix(X), treatment=rep(0, nrow(X)), exp.degree=3, exclude=list(c("educ", "nodegr")))$mat
out.eb <- ebalance(Treatment=treatment, X=X.sq, print.level=2)
ncol(X.sq)
qr(X.sq)
hbal::covarExpand()
hbal::covarExpand
X.sq <- hbal::covarExpand(as.matrix(X), treatment=rep(0, nrow(X)), exp.degree=2, exclude=list(c("educ", "nodegr")))$mat
out.eb <- ebalance(Treatment=treatment, X=scale(X.sq))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=2, Y=Y, ds=TRUE, exclude=list(c("educ", "nodegr")), print.level = -1)
summary(att(out.eb))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=3, Y=Y, ds=TRUE, exclude=list(c("educ", "nodegr")), print.level = -1)
summary(att(out.eb))
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=2, Y=Y,  print.level = -1)
summary(att(out.eb))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=2, Y=Y,  print.level = -1)
summary(att(out.eb))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=2)
summary(att(out.eb, method="lm_robust"))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=3)
summary(att(out.eb, method="lm_robust"))
rm(list=ls())
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
n
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=2, exclude=list(c("educ", "nodegr")))
summary(att(out.eb, method="lm_robust"))
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=3, exclude=list(c("educ", "nodegr")), print.level = -1, max.iterations = 500)
summary(att(out.eb, method="lm_robust"))
out.eb$coefs
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=3, exclude=list(c("educ", "nodegr")), print.level = 2)
summary(att(out.eb, method="lm_robust"))
mena(lalonde$nsw)
mean(lalonde$nsw)
mean(Y[treatment==1])
sum(out.eb$weights)
sum(Y[treatment==0]*out.eb$weights)
6349.145 - 21553.92
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, Y=Y, cv=F, expand.degree=2, exclude=list(c("educ", "nodegr")))
mean(Y[treatment==1]) - sum(Y[treatment==0]*out.eb$weights)
rm(list=ls())
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
n
rm(list=ls())
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
n
dat <- data.frame(Y=Y, X)
dat <- data.frame(Y=Y, T=tratment, X)
?lm_robust
dat <- data.frame(Y=Y, T=treatment, X)
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=2, Y=Y,  print.level = -1)
mod1 <- lm_robust(Y~tratment, weights=c(n, out.eb$weights * c)
)
library(estimatr)
mod1 <- lm_robust(Y~tratment, weights=c(n, out.eb$weights * c))
c
c <- sum(treatment==0)
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
c <- sum(treatment==0)
dat <- data.frame(Y=Y, T=treatment, X)
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=2, Y=Y,  print.level = -1)
dat$w <- c(rep(1, n), out.eb$weights * c)
rm(list=ls())
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
c <- sum(treatment==0)
dat <- data.frame(Y=Y, T=treatment, X)
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=2, Y=Y,  print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod1 <- lm_robust(Y~tratment, weights=w)
mod1 <- lm_robust(Y~tratment, weights=w, data=dat)
mod1 <- lm_robust(Y~treatment, weights=w, data=dat)
summary(mod1)
rm(list=ls())
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
c <- sum(treatment==0)
dat <- data.frame(Y=Y, T=treatment, X)
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=2, Y=Y,  print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod1 <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=3, Y=Y,  print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod2 <- lm_robust(Y~treatment, weights=w, data=dat)
summary(mmod2)
summary(mod2)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=2, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod3 <- lm_robust(Y~treatment, weights=w, data=dat)
summary(mod2)
rm(list=ls())
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
c <- sum(treatment==0)
dat <- data.frame(Y=Y, T=treatment, X)
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=2, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod3 <- lm_robust(Y~treatment, weights=w, data=dat)
summary(mod3)
rm(list=ls())
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
c <- sum(treatment==0)
dat <- data.frame(Y=Y, T=treatment, X)
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=2, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod1 <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=3, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod2 <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=2, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod3 <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=3, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod4 <- lm_robust(Y~treatment, weights=w, data=dat)
summary(mod4)
summary(mod3)
summary(mod2)
summary(mod1)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=0, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
ori <- lm_robust(Y~treatment, weights=w, data=dat)
summary(ori)
library(tidyverse)
m1 <- tidy(mod1)
m2 <- tidy(mod2)
m3 <- tidy(mod3)
m4 <- tidy(mod4)
ori <- tidy(ori)
eb <- rbind(m2, m4)
reg <- rbind(m1, m3)
ori <- ori[c(2),]
eb <- eb[c(2, 4),]
reg <- reg[reg$term=="treat",]
ori[,1] <- c("Linear terms only")
eb[,1]<-c("Series expansion (degree=2)", "Series expansion (degree=3)")
reg[,1]<-c("Series expansion (degree=2)", "Series expansion (degree=3)")
ori <- mutate(ori, model = "Model 2")
eb <- mutate(eb, model = "Model 2")
reg <- mutate(reg, model = "Model 1")
tmod <- rbind(ori,reg,eb)
reg
reg <- rbind(m1, m3)
reg
ori <- tidy(ori)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=0, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
ori <- lm_robust(Y~treatment, weights=w, data=dat)
ori <- tidy(ori)
ori[c(2),]
rm(list=ls())
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
c <- sum(treatment==0)
dat <- data.frame(Y=Y, T=treatment, X)
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=0, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
ori <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=2, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod1 <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=3, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod2 <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=2, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod3 <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=3, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod4 <- lm_robust(Y~treatment, weights=w, data=dat)
m1 <- tidy(mod1)
m2 <- tidy(mod2)
m3 <- tidy(mod3)
m4 <- tidy(mod4)
ori <- tidy(ori)
eb <- rbind(m1, m2)
reg <- rbind(m3, m4)
ori <- ori[c(2),]
eb <- eb[c(2, 4),]
reg <- reg[c(2, 4),]
ori[,1] <- c("Linear terms only")
eb[,1]<-c("Series expansion (degree=2)", "Series expansion (degree=3)")
reg[,1]<-c("Series expansion (degree=2)", "Series expansion (degree=3)")
ori <- mutate(ori, model = "Model 2")
eb <- mutate(eb, model = "Model 2")
reg <- mutate(reg, model = "Model 1")
tmod <- rbind(ori,reg,eb)
dwplot(tmod,dodge_size = 0.15, dot_args = list(aes(shape = model), size = 6, colour="black"),
whisker_args = list(colour="black", size = 1.4),
vline = geom_vline(xintercept = 1794, linetype = 2, size=2)) + theme_bw() +
geom_text(aes(x=1900,
y=0.6),label="Benchmark=1794.00",hjust=1, size=8) +
xlab("Coefficient Estimate")+ ylab("") + theme(legend.title=element_blank()) +
theme(legend.position = c(0.1, 0.9),
legend.justification=c(0.6, 0.6), legend.text = element_text(size = 30),
legend.background = element_rect(colour="grey80"),
legend.title = element_blank(), legend.key.size = unit(2.0, "cm"), legend.key.width = unit(2.5,"cm")) +
scale_shape_discrete(breaks = c("Model 1", "Model 2"),
labels = c("hbal", "ebal")) +
theme(text = element_text(size=30), axis.text.y = element_text(size=30, color="black")) +
scale_x_continuous(breaks=seq(-2000,6000,2000), limits=c(-2500, 6000))
library(dotwhisker)
dwplot(tmod,dodge_size = 0.15, dot_args = list(aes(shape = model), size = 6, colour="black"),
whisker_args = list(colour="black", size = 1.4),
vline = geom_vline(xintercept = 1794, linetype = 2, size=2)) + theme_bw() +
geom_text(aes(x=1900,
y=0.6),label="Benchmark=1794.00",hjust=1, size=8) +
xlab("Coefficient Estimate")+ ylab("") + theme(legend.title=element_blank()) +
theme(legend.position = c(0.1, 0.9),
legend.justification=c(0.6, 0.6), legend.text = element_text(size = 30),
legend.background = element_rect(colour="grey80"),
legend.title = element_blank(), legend.key.size = unit(2.0, "cm"), legend.key.width = unit(2.5,"cm")) +
scale_shape_discrete(breaks = c("Model 1", "Model 2"),
labels = c("hbal", "ebal")) +
theme(text = element_text(size=30), axis.text.y = element_text(size=30, color="black")) +
scale_x_continuous(breaks=seq(-2000,6000,2000), limits=c(-2500, 6000))
dwplot(tmod,dodge_size = 0.15, dot_args = list(aes(shape = model), size = 6, colour="black"),
whisker_args = list(colour="black", size = 1.4),
vline = geom_vline(xintercept = 1794, linetype = 2, size=2)) + theme_bw() +
geom_text(aes(x=1900,
y=0.6),label="Benchmark=1794.00",hjust=1, size=8) +
xlab("Coefficient Estimate")+ ylab("") + theme(legend.title=element_blank()) +
theme(legend.position = c(0.1, 0.9),
legend.justification=c(0.6, 0.6), legend.text = element_text(size = 30),
legend.background = element_rect(colour="grey80"),
legend.title = element_blank(), legend.key.size = unit(2.0, "cm"), legend.key.width = unit(2.5,"cm")) +
scale_shape_discrete(breaks = c("Model 1", "Model 2"),
labels = c("hbal", "ebal")) +
theme(text = element_text(size=30), axis.text.y = element_text(size=30, color="black")) +
scale_x_continuous(breaks=seq(-15000,15000,5000), limits=c(-16500, 16500))
library(hbal)
library(estimatr)
library(tidyverse)
library(dotwhisker)
library(ggplot2)
setwd('/Users/eddie/Desktop/PA Replication/Appendix')
rm(list=ls())
data(lalonde, package = "kbal")
lalonde$nodegr=as.numeric(lalonde$educ<=11)
xvars=c("age","black","educ","hisp","married","re74","re75","nodegr","u74","u75")
W <- lalonde$nsw
X <- lalonde[order(lalonde$nsw, decreasing=T),]
treatment <- X$nsw
Y <- X$re78
X <- X[,xvars]
n <- sum(treatment)
c <- sum(treatment==0)
dat <- data.frame(Y=Y, T=treatment, X)
set.seed(92092)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=0, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
ori <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=2, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod1 <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, cv=F, expand.degree=3, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod2 <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=2, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod3 <- lm_robust(Y~treatment, weights=w, data=dat)
out.eb <- hbal(Treatment=treatment, X=X, expand.degree=3, Y=Y, print.level = -1, exclude=list(c("educ", "nodegr")))
dat$w <- c(rep(1, n), out.eb$weights * c)
mod4 <- lm_robust(Y~treatment, weights=w, data=dat)
m1 <- tidy(mod1)
m2 <- tidy(mod2)
m3 <- tidy(mod3)
m4 <- tidy(mod4)
ori <- tidy(ori)
eb <- rbind(m1, m2)
reg <- rbind(m3, m4)
ori <- ori[c(2),]
eb <- eb[c(2, 4),]
reg <- reg[c(2, 4),]
ori[,1] <- c("Linear terms only")
eb[,1]<-c("Series expansion (degree=2)", "Series expansion (degree=3)")
reg[,1]<-c("Series expansion (degree=2)", "Series expansion (degree=3)")
ori <- mutate(ori, model = "Model 2")
eb <- mutate(eb, model = "Model 2")
reg <- mutate(reg, model = "Model 1")
tmod <- rbind(ori,reg,eb)
png("lalonde_rep.png",  width = 1440, height = 1000)
dwplot(tmod,dodge_size = 0.15, dot_args = list(aes(shape = model), size = 6, colour="black"),
whisker_args = list(colour="black", size = 1.4),
vline = geom_vline(xintercept = 1794, linetype = 2, size=2)) + theme_bw() +
geom_text(aes(x=1900,
y=0.6),label="Benchmark=1794.00",hjust=1, size=8) +
xlab("Coefficient Estimate")+ ylab("") + theme(legend.title=element_blank()) +
theme(legend.position = c(0.1, 0.9),
legend.justification=c(0.6, 0.6), legend.text = element_text(size = 30),
legend.background = element_rect(colour="grey80"),
legend.title = element_blank(), legend.key.size = unit(2.0, "cm"), legend.key.width = unit(2.5,"cm")) +
scale_shape_discrete(breaks = c("Model 1", "Model 2"),
labels = c("hbal", "ebal")) +
theme(text = element_text(size=30), axis.text.y = element_text(size=30, color="black")) +
scale_x_continuous(breaks=seq(-15000,15000,5000), limits=c(-16500, 16500))
dev.off()
