R/EstPS.R

In CrystalXuR/ITRinCEA: Estimating the Optimal ITR from a Cost-effectiveness Perspective

#------------------------------------------------------------------------------------#######
#----------- 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  
}