CBPS (version 0.1)

CBPSlogit: Covariate Balancing Propensity Score (CBPS) Estimation via Logistic Regression

Description

CBPSlogit and CBPSlogit.fit estimate propensity score 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.

Usage

CBPSlogit(formula, data, na.action, ATT, k = 0, bayes = FALSE, ...)
	  CBPSlogit.fit(treat, X, ATT, X.bal = X, k, bayes, ...)

Arguments

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

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

Set to TRUE to find the average treatment effect on the treated, and FALSE to find the average treatment effect. Default is TRUE.

treat

A vector of treatment assignments.

X

A covariate matrix.

X.bal

A matrix to be balanced. Currently not accepted.

k

The model dimensionality. 0 is the default, and is estimated from X. Useful with ill-conditioned matrices.

bayes

Set to TRUE if using bayesglm() in arm package for initial values, and FALSE otherwise. Currently not accepted.

...

Other parameters to be passed through to optim()

Value

• coefficients
  A named vector of coefficients
• residuals
  The working residuals from the final iteration of CBPS
• fitted.values
  The fitted mean values, obtained by transforming the linear predictors by the inverse of the binomial function
• rank
  The numeric rank of the fitted model
• family
  The family object used.
• deviance
  Minus twice the log-likelihood of the CBPS fit. Will be lower than the maximum likelihood deviance.
• null.deviance
  The deviance for the null model, comparable with deviance. Not implemented yet.
• weights
  The CBPS balancing weights calculated from the propensity scores.
• df.residual
  The residual degrees of freedom.
• df.null
  The residual degrees of freedom for the null model.
• 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.
• J
  The J-statistic at convergence.
• df
  The degrees of freedom.
• var
  The covariance matrix, evaluated numerically from optim().
• sd.coef
  Coefficient standard deviations.

Details

Fits covariate balancing propensity scores.

References

Imai, Kosuke and Marc Ratkovic. ``Covariate Balancing Propensity Score.'' Unpublished Manuscript, Princeton University. http://imai.princeton.edu/research/CBPS.html

See Also

summary.CBPSlogit

Examples

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

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

### Example: propensity score weighting 
###
## Simulation from Kang and Shafer (2007).
set.seed(123456)
n <- 1000
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 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 <- CBPSlogit(treat ~ X.mis, ATT = TRUE)
	
## Horwitz-Thompson estimate
mean(treat*y/fit1$fitted.values)
## Inverse probability weighting
sum(treat*y/fit1$fitted.values)/sum(treat/fit1$fitted.values)