# NOT RUN {
set.seed(1234)
## number of units in sample
n <- 2000
## measured potential confounders
z1 <- rnorm(n)
z2 <- rnorm(n)
z3 <- rnorm(n)
z4 <- rnorm(n)
## treatment assignment
prob.treated <- pnorm(-0.5 + 0.75*z2)
x <- rbinom(n, 1, prob.treated)
## potential outcomes
y0 <- z4 + rnorm(n)
y1 <- z1 + z2 + z3 + cos(z3*2) + rnorm(n)
## observed outcomes
y <- y0
y[x==1] <- y1[x==1]
## put everything in a data frame
examp.data <- data.frame(z1, z2, z3, z4, x, y)
## augment data with interactions and powers of covariates
examp.data$z1z1 <- examp.data$z1^2
examp.data$z2z2 <- examp.data$z2^2
examp.data$z3z3 <- examp.data$z3^2
examp.data$z4z4 <- examp.data$z4^2
examp.data$z1z2 <- examp.data$z1 * examp.data$z2
examp.data$z1z3 <- examp.data$z1 * examp.data$z3
examp.data$z1z4 <- examp.data$z1 * examp.data$z4
examp.data$z2z3 <- examp.data$z2 * examp.data$z3
examp.data$z2z4 <- examp.data$z2 * examp.data$z4
examp.data$z3z4 <- examp.data$z3 * examp.data$z4
## check balance of a propensity score model that is not sufficient to
## control confounding bias
bal.1 <- balance.IPW(pscore.formula=x~s(z3)+s(z4),
pscore.family=binomial(probit),
treatment.var="x",
outcome.var="y",
data=examp.data,
nboot=250)
print(bal.1) ## some big z-statistics here indicating balance not so great
## try again
bal.2 <- balance.IPW(pscore.formula=x~z1+z2+z3+z4,
pscore.family=binomial(probit),
treatment.var="x",
outcome.var="y",
data=examp.data,
nboot=250)
print(bal.2) ## balance looks much better--
## only 1 out of 14 zs > 2.0 in absval
# }
Run the code above in your browser using DataLab