# how to calculate manually propensity score weights for multinomial treatments where one of them is baseline

I want to get intuition into the calculation of propensity scores (PS) and inverse probability of treatment weights (IPTW) for a multinomial treatment using multinomial regression. One of the treatments is the baseline treatment.

I am aware that powerful packages automate this but for my analyses, I need to be able to create my own function that allows me to:

1. calculate manually the PS using a multinomial regression model and
2. calculate manually the IPTW from the estimated PS using the relevant formulas.

Here I use the data from here where each of the three categories of the variable group represents a treatment, and the group == 1 represents the baseline treatment (control group) to which each of the other treatments are compared. The approach discussed there gives several PS per row, which is not what I want - all I need is only one PS and one IPTW per row as I get when I use twang and any other similar package.

# data

library("survival")
require("survival")
library("nnet")
require("nnet")

set.seed(42)
days <- rpois(100, 3)
group <- sample(c(1,2,3), 100, replace=TRUE)
status <- rbinom(100,1,0.65)
demo1 <- rnorm(100,100,25)
demo2 <- rpois(100,10)
demo3 <- rbinom(100,1,0.67)

df <- data.frame(days, status, group, demo1, demo2,demo3)


Thanks in advance for any help.

For the ATE, the IPW is the inverse of the model-predicted probability of being in the treatment actually received. Using multinom() in nnet, we can generate a multinomial regression model and then extract the predicted probabilities for each treatment group for each individual using predict(). Then we find the predicted probability for each individual that corresponds to the treatment actually received, and take its inverse. Finally, I use bal.tab() in cobalt to assess balance, demonstrating that balance has been achieved.

library("survival")
library("nnet")
library("cobalt")

set.seed(42)
days <- rpois(100, 3)
group <- factor(sample(c(1,2,3), 100, #needs to be a factor
replace=TRUE))
status <- rbinom(100,1,0.65)
demo1 <- rnorm(100,100,25)
demo2 <- rpois(100,10)
demo3 <- rbinom(100,1,0.67)

df <- data.frame(days, status, group, demo1, demo2,demo3)

fit <- multinom(group ~ demo1 + demo2 + demo3, data = df)
#> # weights:  15 (8 variable)
#> initial  value 109.861229
#> iter  10 value 107.354113
#> final  value 107.353911
#> converged

ps.mat <- predict(fit, type = 'probs')

w <- rep(0, nrow(df)) #inititalize weights

for (i in levels(group)) {
w[group == i] <- 1/ps.mat[group == i, i]
}

bal.tab(group ~ demo1 + demo2 + demo3, data = df,
weights = w, un = TRUE)
#> Assuming "weighting". If not, specify with an argument to method.
#> Note: estimand and s.d.denom not specified; assuming ATE and pooled.
#> Balance summary across all treatment pairs
#> demo1 Contin.      0.2034       0.0153
#> demo2 Contin.      0.3131       0.2557
#> demo3  Binary      0.0753       0.0306
#>
#> Effective sample sizes
#>                 1      2      3


Created on 2019-01-23 by the reprex package (v0.2.1)

For ATT weights, you simply multiply everyone's ATE weight by their probability of being in the focal (i.e., treated) group. Those in the focal group will now have a weight of 1, just like in the binary case. For example, if the focal group was group 2, you would simply do the following after running the above code:

w_att <- w*ps.mat[,"2"]


This multiplies each unit's weight by the probability of being in the focal group. For those in group 2, their original weight was equal to 1/ps.mat[,"2"] so multiplying by ps.mat[,"2"] will yield weights of 1 for those in group 2.

• Just one clarification if you don't mind: there is also this formula weights = treatment/psvector + (1-treatment)/(1-psvector) in addition to the formula used in your answer w[group == i] <- 1/ps.mat[group == i, i]. How would this other formula be used in your answer? Thanks in advance. Jan 23, 2019 at 20:40
• That's a formula for weights when you only have 2 groups, and psvector is the probability of being in the treated group. It doesn't apply here. With only two groups, though, the formula I presented works when you have a matrix with two columns, one for the probability of being in each treatment. In this case, the second column would be 1 minus the first column.
– Noah
Jan 23, 2019 at 21:26
• What my line of code says in English is, "for each person in group i, take the probability of being in group i, which lives in column i of ps.mat, take the inverse of that probability, and then assign that value to their weight".
– Noah
Jan 23, 2019 at 21:29
• That explanation is very helpful. Thanks a lot for that, @Noah. Jan 23, 2019 at 21:49
• The estimand (i.e., the ATE vs. ATT) is simply a matter of the target population, the population to which you want to generalize the causal effect estimate. If your sample is a random sample from the target population, the ATE will generalize to that population. You need to reweight each group to have the same covariate distribution as the overall sample. The ATT generalizes to a population with a similar covariate distribution to the "treated" units, so you need to reweight the other groups to appear like the treated group.
– Noah
Jan 24, 2019 at 6:56