Weights from inverse probability treatment weighting in regression model in R? I obtained the weights for the data from inverse probability weighting for exposure variable.
Can I use these weights in rms::lrm model or glm binomial model to weight data?
Definig model as:
lrm(outcome ~ age + ethnicity + LTC_count + exposure,    weights = inv_prop_treatment_weights, data = data )  
The weights in documentation are described as "as a vector (same length as y) of possibly fractional case weights". It is unclear to me if these are the inverse probablitiy of treatment weights. If not, how can weights be included in a regression model?
The same model can be defined with weights with binomial family glm, instead of rms::lrm as glm in r also has weights parameter.
Alternatively - if I use the lrm or binomial glm model without weights, but use propensity scores as a variable, definig model as
lrm(outcome ~ age + ethnicity + LTC_count + exposure + prop_score_of_exposure, data = data)
Would this be a doubly robust regression?
I have read a lot on propensity scores and matching, but found litte on how to practically use weigths or propensity scores in regression models, so the help will be very appreciated.
 A: You need to be clear about the quantity you want to estimate. Most causal inference applications are concerned with the average marginal effect of the treatment on the outcome. This does not correspond to the coefficient on treatment in a logistic regression model. The way to use regression to estimate causal effects is to use g-computation, which involves fitting a model for the outcome given the treatment and covariates (and their interaction), then using this model to predict the potential outcomes under treatment and under control for all units, then taking the difference in means between those predict potential outcomes. G-computation is consistent if the outcome model is consistent.
Below is how you would do g-computation in R:
#Fit the outcome model with an interaction between the treatment and covariates
fit <- glm(y ~ t * (x1 + x2), data = data, family = quasibinomial)

#Estimate potential outcomes under treatment
data$t <- 1
pred1 <- predict(fit, newdata = data, type = "response")

#Estimate potential outcomes under control
data$t <- 0
pred0 <- predict(fit, newdata = data, type = "response")

#Compute risk difference
mean(pred1) - mean(pred0)

The coefficient on treatment in a logistic regression of the outcome on the treatment and covariates is an estimate of the conditional effect of the treatment assuming no effect heterogeneity on the odds ratio scale. This is an extremely strong and unnecessary assumption to make. The reason the effect is the conditional rather than the marginal effect is because of the noncollapsibility of the odds ratio. This is not the case for linear regression; it is possible (by correctly parameterizing the linear model) to have the coefficient on treatment be a valid estimate of the marginal treatment effect. There must be an interaction between treatment and the mean-centered covariates for this to work.
I will assume that you are now using g-computation and want to know how you can incorporate weights into it. You can include the propensity score weights into an outcome model and then perform g-computation using that model. This method is doubly robust. Kang and Schafer (2007) call this method "Regression Estimation with Inverse-Propensity Weighted Coefficients" (regression estimation is another name for g-computation). Simply add weights to the glm() model exactly as you specified, and insert the code below into its appropriate place in the g-computation code above.
#G-computation with a weighted outcome model
fit <- glm(y ~ t * (x1 + x2), data = data, family = quasibinomial,
           weights = weights)

There are many ways to produce doubly-robust estimators that are consistent if either the propensity score model or outcome model are correctly specified. The way you proposed is one, but you also have to get the conditional relationship between the outcome and propensity score correct, so simply adding it in as a covariate will not do. You need to allow it to be flexibly modeled. Using rcs() from rms is one way to do that.
library(rms)

#G-computation with PS as an additional covariate, modeled with a 
#restricted cubic spline (rcs)
fit <- glm(y ~ t * (x1 + x2 + rcs(ps)), data = data, family = quasibinomial)

This model has a chance of being doubly robust.
Estimating confidence intervals is a little challenging with g-computation; in some cases, you can use the delta method/estimating equations, but the most straightforward and recommended way is to bootstrap. You need to include the propensity score estimation in the bootstrap. I demonstrate this below using the boot package.
library(boot)
boot_fun <- function(d, i) {
  d_boot <- d[i,]
  
  #Estimate the PS
  ps_fit <- glm(t ~ x1 + x2, data = d_boot, family = binomial)
  
  d_boot$ps <- fitted(ps_fit)
  d_boot$ps_w <- with(d_boot, t/ps + (1-t)/(1-ps))
  
  #Fit the weighted outcome model with an interaction between the treatment and covariates
  fit <- glm(y ~ t * (x1 + x2), data = d_boot, family = quasibinomial,
             weights = ps_w)
  
  #Estimate potential outcomes under treatment
  d_boot$t <- 1
  pred1 <- predict(fit, newdata = d_boot, type = "response")
  
  #Estimate potential outcomes under control
  d_boot$t <- 0
  pred0 <- predict(fit, newdata = d_boot, type = "response")
  
  #Compute risk difference
  mean(pred1) - mean(pred0)
}
#Effect estimate ("original" column)
(boot_fit <- boot(data, boot_fun, R = 999))

#Confidence interval (percentile)
boot.ci(boot_fit, type = "perc")

This is for the average treatment effect in the population. Things are a little different for the average treatment effect in the treated. If you're using matching, see the guide on estimating effects after matching with MatchIt here. I also recommend you check out the WeightIt package, which makes estimating weights very straightforward. Both packages have extensive documentation explaining how to estimate effects.
