#------------------------------------------------------------------------------------#######
#----------- Estimate propensity score using logistic regression --------------------#######
#----------- Outcome: treatment A --------------------#######
#----------- Predictor: covariate Xs --------------------#######
#------------------------------------------------------------------------------------#######
EstPS = function(A,Xs){
if(length(A)!= nrow(as.matrix(Xs))) stop("Treatment and covariates have different length!")
if(ncol(as.matrix(A))!=1 | length(unique(A))!= 2) stop("Multiple treatments present, or treatment levels are not 2!")
if(!(unique(A)[1] %in% c(0,1)) | !(unique(A)[2] %in% c(0,1))) stop("Treatment levels are not 0 and 1")
dat = data.frame(A,Xs)
fit = glm(A ~., family="binomial", data=dat)
e = predict(fit,dat,type = "response")
PS = data.frame((1-e),e);names(PS) = c("Untreated","Treated");PS
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.