if (FALSE) {
# data generating process
p = 3 #two control variables and one moderator
n = 250
set.seed(1)
x = matrix(rnorm(n*p), nrow=n)
# create targeted selection
q = -1*(x[,1]>(x[,2])) + 1*(x[,1]<(x[,2]))
# generate treatment variable
pi = pnorm(q)
z = rbinom(n,1,pi)
# tau is the true (homogeneous) treatment effect
tau = (0.5*(x[,3] > -3/4) + 0.25*(x[,3] > 0) + 0.25*(x[,3]>3/4))
# generate the response using q, tau and z
mu = (q + tau*z)
# set the noise level relative to the expected mean function of Y
sigma = diff(range(q + tau*pi))/8
# draw the response variable with additive error
y = mu + sigma*rnorm(n)
# If you didn't know pi, you would estimate it here
pihat = pnorm(q)
bcf_fit = bcf(y, z, x, x, pihat, nburn=2000, nsim=2000)
# Get posterior of treatment effects
tau_post = bcf_fit$tau
tauhat = colMeans(tau_post)
plot(tau, tauhat); abline(0,1)
}
if (FALSE) {
# data generating process
p = 3 #two control variables and one moderator
n = 250
#
set.seed(1)
x = matrix(rnorm(n*p), nrow=n)
# create targeted selection
q = -1*(x[,1]>(x[,2])) + 1*(x[,1]<(x[,2]))
# generate treatment variable
pi = pnorm(q)
z = rbinom(n,1,pi)
# tau is the true (homogeneous) treatment effect
tau = (0.5*(x[,3] > -3/4) + 0.25*(x[,3] > 0) + 0.25*(x[,3]>3/4))
# generate the response using q, tau and z
mu = (q + tau*z)
# set the noise level relative to the expected mean function of Y
sigma = diff(range(q + tau*pi))/8
# draw the response variable with additive error
y = mu + sigma*rnorm(n)
pihat = pnorm(q)
# nburn and nsim should be much larger, at least a few thousand each
# The low values below are for CRAN.
bcf_fit = bcf(y, z, x, x, pihat, nburn=100, nsim=10)
# Get posterior of treatment effects
tau_post = bcf_fit$tau
tauhat = colMeans(tau_post)
plot(tau, tauhat); abline(0,1)
}
Run the code above in your browser using DataLab