You can use any causal inference studies data to test the package.
The database needs to have the following attributes:
The package can generate synthetic data that can be used to test
different features of the package. At the current implementation, the
code only generates synthetic numerical data; however, with small
innovation, one can add categorical data. In the following, we present
some reproducible examples that you can copy and build upon them.
Estimating GPS Values for the Dataset
In the following example, first, we generate 10000 synthetic data
samples, then we feed them into estimate_gps function to
estimate their GPS values. Please note that the Y variable is not used
inside the code, and it is only provided to cbind with the generated GPS
values. You can read more about different arguments in the documentation
(?estimate_gps). In summary, - We want to run a code using
SuperLearner (sl) prediction model. Inside the SuperLearner
package, we want to use the XGBoost package. As you know, we internally
generate a customized wrapper for the XGBoost package; thus, the correct
way to activate that is by passing m_xgboost, which stands
for modified xgboost. - We requested for one thread
(nthread = 1); you can use as much as you want; the package
will use the available one and will ignore the rest of them. XGBoost
uses OpenMP backend to use all cores. Sometimes it becomes a really big
challenge on Mac systems. So if you are using Mac and do not see any
performance improvement, be aware of that. - In params, we
passed the list of parameters; the function will choose one from each
list at random and will generate a customize wrapper based on them. If
you want to use a specific value, just give a list one number (e.g.,
xgb_max_depth = c(3)). All parameters that start with xgb_
will only change XGBoost hyperparameters.
library("CausalGPS")
mydata <- generate_syn_data(sample_size = 10000)
data_with_gps <- estimate_gps(mydata$Y,
mydata$treat,
mydata[c("cf1","cf2","cf3","cf4","cf5","cf6")],
pred_model = "sl",
internal_use = FALSE,
params = list(xgb_max_depth = c(3,4,5),
xgb_nrounds=c(10,20,30,40,50,60)),
nthread = 1,
sl_lib = c("m_xgboost")
)
Now, let’s add some categorical data. Let’s say our data belongs to 5
different years (2000 observations per year), and each year we have data
from 4 different regions, including North, South,
East, and West.
library("CausalGPS")
mydata <- generate_syn_data(sample_size = 10000)
year <- c(rep(c("2001"), each=2000),
rep(c("2002"), each=2000),
rep(c("2003"), each=2000),
rep(c("2004"), each=2000),
rep(c("2005"), each=2000))
region <- rep(c(rep("North",each=500),
rep("South",each=500),
rep("East",each=500),
rep("West",each=500)), each=5)
mydata$year <- as.factor(year)
mydata$region <- as.factor(region)
data_with_gps <- estimate_gps(mydata$Y,
mydata$treat,
mydata[c("cf1","cf2","cf3","cf4","cf5","cf6", "year", "region")],
pred_model = "sl",
internal_use = FALSE,
params = list(xgb_max_depth = c(3,4,5),
xgb_nrounds=c(10,20,30,40,50,60)),
nthread = 1,
sl_lib = c("m_xgboost")
)
Generating Pseudo Population
Generating pseudo population (gen_pseudo_pop())is one of
the important parts of the package. It internally uses
estimate_gps() package. An acceptable pseudo population
should pass the covariate balance test. Users choose which test and what
threshold, and how many attempts. If the covariate balance test is met,
the function stops trying and returns the generated population. However,
if it cannot pass the test, it still returns what is generated with a
message to the user indicating the population did not pass the test. At
each iteration, we change the hyperparameters to modify GPS values (to
better or worse) to pass the covariate balance test. After estimating
the GPS values, we need to compile the population. There are three major
approaches to compile pseudo population, including:
- Matching
- Adjusting (not implemented)
- Weighting (not implemented)
Among different methods for testing covariate balance test, only
absolute approach is implemented.
library("CausalGPS")
mydata <- generate_syn_data(sample_size = 10000)
pseudo_pop <- generate_pseudo_pop(mydata$Y,
mydata$treat,
mydata[c("cf1","cf2","cf3","cf4","cf5","cf6")],
ci_appr = "matching",
pred_model = "sl",
sl_lib = c("m_xgboost"),
params = list(xgb_nrounds=c(10,20,30),
xgb_eta=c(0.1,0.2,0.3)),
nthread = 1,
covar_bl_method = "absolute",
covar_bl_trs = 0.1,
covar_bl_trs_type= "mean",
max_attempt = 1,
matching_fun = "matching_l1",
delta_n = 1,
scale = 0.5)
The second causal inference approach is weighting. Here
is an example to generate pseudo population using weighting
approach.
mydata <- generate_syn_data(sample_size = 10000)
year <- c(rep(c("2001"), each=2000),
rep(c("2002"), each=2000),
rep(c("2003"), each=2000),
rep(c("2004"), each=2000),
rep(c("2005"), each=2000))
region <- rep(c(rep("North",each=500),
rep("South",each=500),
rep("East",each=500),
rep("West",each=500)), each=5)
mydata$year <- as.factor(year)
mydata$region <- as.factor(region)
pseudo_pop <- generate_pseudo_pop(mydata$Y,
mydata$treat,
mydata[c("cf1","cf2","cf3","cf4","cf5","cf6","year","region")],
ci_appr = "weighting",
pred_model = "sl",
sl_lib = c("m_xgboost"),
params = list(xgb_nrounds=c(10,20,30),
xgb_eta=c(0.1,0.2,0.3)),
nthread = 1,
covar_bl_method = "absolute",
covar_bl_trs = 0.1,
covar_bl_trs_type = "mean",
max_attempt = 1
)
Further Processing
After generating a pseudo population, we can process the data for
different purposes. So far, estimating exposure rate function
(estimate_erf) is implemented.
# library("CausalGPS")
mydata <- generate_syn_data(sample_size = 10000)
pseudo_pop <- generate_pseudo_pop(mydata$Y,
mydata$treat,
mydata[c("cf1","cf2","cf3","cf4","cf5","cf6")],
ci_appr = "matching",
pred_model = "sl",
sl_lib = c("m_xgboost"),
params = list(xgb_nrounds=c(10,20,30),
xgb_eta=c(0.1,0.2,0.3)),
nthread = 1,
covar_bl_method = "absolute",
covar_bl_trs = 0.1,
covar_bl_trs_type= "mean",
max_attempt = 1,
matching_fun = "matching_l1",
delta_n = 1,
scale = 0.5)
erf_val <- estimate_npmetric_erf(pseudo_pop$pseudo_pop$Y,
pseudo_pop$pseudo_pop$w,
bw_seq=seq(0.2,2,0.2),
w_vals = seq(2,20,0.5))
