Counterfactual Bayesian survival analysis in pymc I am trying to determine mortality rates for untreated patients from an observational dataset where treatment has occurred (thus blocking the possibility of further untreated mortality). You can't just model this with a CoxPH model and consider treatment as a censoring event due to the coxph prohibition on informative censoring.  The situation is also complicated by the fact that the covariates that drive mortality also drive probability of treatment.  I've had success using inverse probability censor weighting, but now I'd like to try a Bayesian approach.
I can successfully model in pymc observed survival and observed rates of treatment using Poisson's regression.  I'm struggling to figure out how to tie those models together to figure out the unobserved counterfactual "what would mortality have been if the patients weren't treated".  Any suggestions?

What I ended up doing with (some success, but I'm still fiddling to improve) was use  a cure rate model where:
SurvivalObsevered(t) = TreatedPortion(t) + (1-TreatedPortion(t))*SurvivalUntreated(t)
In effect what I'm doing is considering treatment to be a "cure" for untreated mortality.  This allows me to model and ultimately solve for the quantity of interest which is untreated survival versus time which is not directly observed.
