CBPS: Covariate Balancing Propensity Score (CBPS) Estimation
In CBPS: Covariate Balancing Propensity Score

Description Usage Arguments Value Author(s) References See Also Examples

CBPS estimates propensity scores such that both covariate balance and prediction of treatment assignment are maximized. The method, therefore, avoids an iterative process between model fitting and balance checking and implements both simultaneously. For cross-sectional data, the method can take continuous treatments and treatments with a control (baseline) condition and either 1, 2, or 3 distinct treatment conditions.

Fits covariate balancing propensity scores.

### @aliases CBPS CBPS.fit print.CBPS

CBPS(
  formula,
  data,
  na.action,
  ATT = 1,
  iterations = 1000,
  standardize = TRUE,
  method = "over",
  twostep = TRUE,
  sample.weights = NULL,
  baseline.formula = NULL,
  diff.formula = NULL,
  ...
)

`formula`	An object of class `formula` (or one that can be coerced to that class): a symbolic description of the model to be fitted.
`data`	An optional data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. If not found in data, the variables are taken from `environment(formula)`, typically the environment from which `CBPS` is called.
`na.action`	A function which indicates what should happen when the data contain NAs. The default is set by the na.action setting of options, and is na.fail if that is unset.
`ATT`	Default is 1, which finds the average treatment effect on the treated interpreting the second level of the treatment factor as the treatment. Set to 2 to find the ATT interpreting the first level of the treatment factor as the treatment. Set to 0 to find the average treatment effect. For non-binary treatments, only the ATE is available.
`iterations`	An optional parameter for the maximum number of iterations for the optimization. Default is 1000.
`standardize`	Default is `TRUE`, which normalizes weights to sum to 1 within each treatment group. For continuous treatments, normalizes weights to sum up to 1 for the entire sample. Set to `FALSE` to return Horvitz-Thompson weights.
`method`	Choose "over" to fit an over-identified model that combines the propensity score and covariate balancing conditions; choose "exact" to fit a model that only contains the covariate balancing conditions.
`twostep`	Default is `TRUE` for a two-step estimator, which will run substantially faster than continuous-updating. Set to `FALSE` to use the continuous-updating estimator described by Imai and Ratkovic (2014).
`sample.weights`	Survey sampling weights for the observations, if applicable. When left NULL, defaults to a sampling weight of 1 for each observation.
`baseline.formula`	Used only to fit iCBPS (see Fan et al). Currently only works with binary treatments. A formula specifying the balancing covariates in the baseline outcome model, i.e., E(Y(0)\|X).
`diff.formula`	Used only to fit iCBPS (see Fan et al). Currently only works with binary treatments. A formula specifying the balancing covariates in the difference between the treatment and baseline outcome model, i.e., E(Y(1)-Y(0)\|X).
`...`	Other parameters to be passed through to `optim()`.

`fitted.values`	The fitted propensity score
`linear.predictor`	X * beta
`deviance`	Minus twice the log-likelihood of the CBPS fit
`weights`	The optimal weights. Let π_i = f(T_i \| X_i). For binary ATE, these are given by T_i/π_i + (1 - T_i)/(1 - π_i). For binary ATT, these are given by n/n_t (T_i - π_i)/(1 - π_i). For multi_valued treatments, these are given by ∑_{j=0}^{J-1} T_i,j / π_i,j. For continuous treatments, these are given by f(T_i) / f(T_i \| X_i)* . These expressions for weights are all before standardization (i.e. with standardize=`FALSE`). Standardization will make weights sum to 1 within each treatment group. For continuous treatment, standardization will make all weights sum to 1. If sampling weights are used, the weight for each observation is multiplied by the survey sampling weight.
`y`	The treatment vector used
`x`	The covariate matrix
`model`	The model frame
`converged`	Convergence value. Returned from the call to `optim()`.
`call`	The matched call
`formula`	The formula supplied
`data`	The data argument
`coefficients`	A named vector of coefficients
`sigmasq`	The sigma-squared value, for continuous treatments only
`J`	The J-statistic at convergence
`mle.J`	The J-statistic for the parameters from maximum likelihood estimation
`var`	The covariance matrix for the coefficients.
`Ttilde`	For internal use only.
`Xtilde`	For internal use only.
`beta.tilde`	For internal use only.
`simgasq.tilde`	For internal use only.

Christian Fong, Marc Ratkovic, Kosuke Imai, and Xiaolin Yang; The CBPS function is based on the code for version 2.15.0 of the glm function implemented in the stats package, originally written by Simon Davies. This documentation is likewise modeled on the documentation for glm and borrows its language where the arguments and values are the same.

Imai, Kosuke and Marc Ratkovic. 2014. “Covariate Balancing Propensity Score.” Journal of the Royal Statistical Society, Series B (Statistical Methodology). http://imai.princeton.edu/research/CBPS.html
Fong, Christian, Chad Hazlett, and Kosuke Imai. 2018. “Covariate Balancing Propensity Score for a Continuous Treatment.” The Annals of Applied Statistics. http://imai.princeton.edu/research/files/CBGPS.pdf
Fan, Jianqing and Imai, Kosuke and Liu, Han and Ning, Yang and Yang, Xiaolin. “Improving Covariate Balancing Propensity Score: A Doubly Robust and Efficient Approach.” Unpublished Manuscript. http://imai.princeton.edu/research/CBPStheory.html

summary.CBPS

###
### Example: propensity score matching
###

##Load the LaLonde data
data(LaLonde)
## Estimate CBPS
fit <- CBPS(treat ~ age + educ + re75 + re74 + 
			I(re75==0) + I(re74==0), 
			data = LaLonde, ATT = TRUE)
summary(fit)
## Not run: 
## matching via MatchIt: one to one nearest neighbor with replacement
library(MatchIt)
m.out <- matchit(treat ~ fitted(fit), method = "nearest", 
				 data = LaLonde, replace = TRUE)

### Example: propensity score weighting 
###
## Simulation from Kang and Shafer (2007).
set.seed(123456)
n <- 500
X <- mvrnorm(n, mu = rep(0, 4), Sigma = diag(4))
prop <- 1 / (1 + exp(X[,1] - 0.5 * X[,2] + 
			 0.25*X[,3] + 0.1 * X[,4]))
treat <- rbinom(n, 1, prop)
y <- 210 + 27.4*X[,1] + 13.7*X[,2] + 13.7*X[,3] + 13.7*X[,4] + rnorm(n)

##Estimate CBPS with a misspecified model
X.mis <- cbind(exp(X[,1]/2), X[,2]*(1+exp(X[,1]))^(-1)+10, 
			  (X[,1]*X[,3]/25+.6)^3, (X[,2]+X[,4]+20)^2)
fit1 <- CBPS(treat ~ X.mis, ATT = 0)
summary(fit1)
	
## Horwitz-Thompson estimate
mean(treat*y/fit1$fitted.values)
## Inverse propensity score weighting
sum(treat*y/fit1$fitted.values)/sum(treat/fit1$fitted.values)

rm(list=c("y","X","prop","treat","n","X.mis","fit1"))

### Example: Continuous Treatment as in Fong, Hazlett, 
### and Imai (2018).  See 
### https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/AIF4PI
### for a real data example.
set.seed(123456)
n <- 1000
X <- mvrnorm(n, mu = rep(0,2), Sigma = diag(2))
beta <- rnorm(ncol(X)+1, sd = 1)
treat <- cbind(1,X)%*%beta + rnorm(n, sd = 5)

treat.effect <- 1
effect.beta <- rnorm(ncol(X))
y <- rbinom(n, 1, (1 + exp(-treat.effect*treat - 
				   X%*%effect.beta))^-1)

fit2 <- CBPS(treat ~ X)
summary(fit2)
summary(glm(y ~ treat + X, weights = fit2$weights, 
			family = "quasibinomial"))

rm(list=c("n", "X", "beta", "treat", "treat.effect",
		  "effect.beta", "y", "fit2"))

### Simulation example: Improved CBPS (or iCBPS) from Fan et al
set.seed(123456)
n <- 500
X <- mvrnorm(n, mu = rep(0, 4), Sigma = diag(4))
prop <- 1 / (1 + exp(X[,1] - 0.5 * X[,2] + 0.25*X[,3] + 0.1 * X[,4]))
treat <- rbinom(n, 1, prop)
y1 <- 210 + 27.4*X[,1] + 13.7*X[,2] + 13.7*X[,3] + 13.7*X[,4] + rnorm(n)
y0 <- 210 + 13.7*X[,2] + 13.7*X[,3] + 13.7*X[,4] + rnorm(n)
##Estimate iCBPS with a misspecificied model
X.mis <- cbind(exp(X[,1]/2), X[,2]*(1+exp(X[,1]))^(-1)+10, 
			   (X[,1]*X[,3]/25+.6)^3, (X[,2]+X[,4]+20)^2)
fit1 <- CBPS(treat ~ X.mis, baseline.formula=~X.mis[,2:4], 
			 diff.formula=~X.mis[,1], ATT = FALSE)
summary(fit1)

## End(Not run)