How to improve structural bias correction using propensity scores

Given a control group and a treatment group, which are not equivalent (different distributions on a number of variables), I want to estimate the probability of an event given that the observation was treated or not. However, since my groups are nonequivalent, I cannot simply compare the computed probability of my control group with that of my treatment group.

I want to find a way to "correct" for the structural bias between my two groups. I know that one method for this is propensity score matching. I have tried to implement the methodology described in Austin, 2011, A Tutorial and Case Study in Propensity Score Analysis.

I run a simulation in order to know my expectations and so that all of my event's variance can be explained by the included covariates.

There are two datasets. The "Equivalent" dataset's treatment and control group are equivalent, so the effect of treatment can be measured directly. The "Non-Equivalent" dataset's T and C groups are non-equivalent (a large portion of group C is skewed, with a much lower probability of having score=1)

The treatment group (20000 observations) has the same properties for each dataset. For the treatment group:

• The "pre-score" binary variable follows a binomial distribution with p = 0.01
• The "activity" continuous variable is normally distributed with mean 200 and std 100
• The "score" binary variable has the following property: P(score|pre-score = 1) = 0.12 P(score|pre-score = 0) = 0.035

In the "Equivalent" dataset, the control group (20000 observations) has the same properties for "pre-score" and "activity", but P(score=1) is lower:

P(score|pre-score = 1) = 0.08
P(score|pre-score = 0) = 0.03

In the "Non-Equivalent" dataset, the control group is built from two sub-groups: - Subgroup 1 (n=5000) has the same properties as the Control group in the "Equivalent" dataset - Subgroup 2's (n=15000) variables follow different distributions: "pre-score" follows a binomial distribution with p = 0.0001 "activity" follows a normal distribution with means 100 and stdev 100

P(score=1|pre-score=1) = 0.06
P(score=1|pre-score=0) = 0.02

For the Equivalent dataset, I measure and compare P(score=1|treatment=1) with P(score=1|treatment=0) directly.

For the Non-Equivalent dataset, I attempt to balance my groups based on the propensity score matching method:

First, calculate propensity scores for each observation by fitting a logistic regression:

treatment ~ pre-score + activity

The propensity score is equal to the predicted probability, for each observation, of being in the treatment group.

Second, match each observation of the smallest group (for instance, treatment=1) with one observation of the other group, using shortest distance as the matching criterion, with a caliper of 0.05.

Third, drop all observations that contain NA values.

Finally, compute P(score=1|treatment=1) and P(score=1|treatment=0) and the risk ratio for each dataset and see how the propensity score method affects my results.

In my Equivalent dataset, I find p1/p0 = 1.15.

In my Non-Equivalent dataset, p1/p0 = 1.55 for the whole dataset and p1/p0 = 1.51 for the matched observations. Of course, I would have expected to find something around 1.15...

This seems to show that my method is not working very well. My matched data does not seem very balanced. This is particularly clear upon looking at the means of the covariates (pre-score and activity) in the treatment and control group, remain very different albeit more similar than they are in the Non-Equivalent Group.

I would like to know how to improve the efficiency of this method (or how to use another method to get what I want). I haven't included code but I can include code if you think it can help or that there may be a mistake in my code or something I've omitted.