Survival analysis with inverse probability weighting after multiple imputation

[Ester Olthof]

Dear Statalist-users,

Hopefully, there is someone who can help me with the following:
I'm trying to perform a survival analysis with inverse probability weighting after multiple imputation.
I've tried to apply the approach as described previously: https://www.statalist.org/forums/for...opensity-score
But after (successful) multiple data imputation (m=10), I don't know how to proceed from here:

 * Running my propensity score model for M1-M10 and save complete-data estimation results to miest.ster using mi estimate's saving() option.
 mi estimate, saving(miest, replace): logit thergr leeft i.FIGO_2009 i.ni_pb i.ni_loc newnode_dm tumsize i.morf_cat i.invasiedieptegr i.diffgrad newbmi i.cci i.lvsi  

* Obtain multiple-imputation linear predictions and store them as variable xb mi in the original data (m=0).
mi predict xb_mi using miest  

* Apply the inverse-logit transformation to obtain the probabilities.  
quietly mi xeq: generate phat = invlogit(xb_mi)
mi xeq: generate phat = invlogit(xb_mi)  

* Generate IPTW's  
mi xeq: gen psweight=.
mi xeq: replace psweight = (1/phat) if thergr==1
mi xeq: replace psweight = (1/(1-phat)) if thergr==0  

* Set survival time
mi stset vitfup_years [pweight=psweight], failure(vit_stat) id(rn)

If I understand correctly, the propensity scores (phat), and thereby the IPTW (psweight), are already pooled estimates of M1-M10.
They are stored in the original data (m=0) and empty/missing for m=1-10 in the browser.
Should I, therefore, proceed with:

 stcox i.thergr phat

(And is it necessary to add the propensity score (phat) in the cox-regression while the weights (psweight) are already added during the declaration of survival data?)

Because if I proceed with mi estimate, I end up with only 59 observations which were the complete cases in M0.

 mi estimate , eform: stcox i.thergr phat

Lastly, this method does not seem to align with the MIte approach to me, or is it?
(Leyrat, C, et al. (available online; in press), "Propensity score analysis with partially observed covariates: How should multiple imputation be used?" Statistical Methods in Medical Research)

Thanks in advance for helping me figure this out!

Kind regards, Ester

[Ester Olthof]

For those who are interested, I might have figured out a way to calculate survival estimates using inverse probability treatment weighting (IPTW) after multiple imputation.

As mentioned in the previous post, I’m interested in overall survival, comparing two treatment groups, indicated by variable thergr.
Because of missing data, I performed multiple imputation with 20 datasets.

I. First, I've calculated the propensity scores and IPW per imputed dataset (M1-20)

 * Calculate propensity scores by pscore
mi xeq: pscore thergr ln_leeft figo ni_pb ni_loc ln_newnode_dm ln_newklin_grootte morf_cat diffgrad ln_newbmi cci lvsi, pscore(ps) blockid(myblock) logit detail

* Compute inverse probability weights
gen ipw=1
replace ipw=1/ps if thergr==1
replace ipw=1/(1-ps) if thergr==0
egen sumofweights = total(ipw)
gen norm_weights  = ipw/sumofweights

* Balance check of included variables
xi: mpbalchk thergr leeft i.figo i.ni_pb i.ni_loc newnode_dm newklin_grootte i.morf_cat i.diffgrad newbmi i.cci i.lvsi, wt(norm_weights) graph

I believe this is necessary to include the propensity score in my cox-regression analysis later on, and to check for covariate balance.

II. Second, calculate pooled propensity scores and IPW which are stored in M0.

 * Calculate pooled propensity scores, to be saved in M0:
mi estimate, saving(miest, replace): logit thergr ln_leeft i.figo i.ni_pb i.ni_loc ln_newnode_dm ln_newklin_grootte i.morf_cat i.diffgrad ln_newbmi i.cci i.lvsi
mi predict xb_mi using miest
quietly mi xeq: generate psM0 = invlogit(xb_mi)

* Generate IPTW's 
mi xeq: gen psweight=.
mi xeq: replace psweight = (1/psM0) if thergr==1
mi xeq: replace psweight = (1/(1-psM0)) if thergr==0

* Combine PS en IPW for M0 with M1-20
replace ps=psM0 if _mi_m==0
replace psweight=ipw if _mi_m!=0

I believe this is necessary to include the IPW in the set for survival time.

III. Third, calculate the cox-proportional hazard ratio

mi stset vitfup_years [pweight=psweight], failure(vit_stat) id(_mi_id)
mi estimate, hr baselevels: stcox i.thergr ps

IV. Fourth, calculate the 5-years adjusted survival function graph and rate of pooled estimates

 save results, emptyok replace // file to hold results

local isets = r(M) /* number of imputed data sets */
mi stset vitfup_years [pweight=psweight], failure(vit_stat) exit(time 5)
return list
mi estimate, hr baselevels: stcox i.thergr ps
return list

local isets = 20
tempfile d0
save `d0', replace
quietly {
forvalues i = 1/`isets'{
  tempfile d`i'
  use `d0', clear
  mi extract `i' 
  qui stcox i.thergr ps
  stcurve , surv at1(thergr==0) at2(thergr==1) ///
  outfile(`d`i'', replace) ylabel(0.00(0.25)1.00)
  use results, clear
  append using `d`i''
  save, replace
}
}

use results, clear
collapse (mean) surv1 (mean) surv2, by(_t)
label var surv1 "Chemoradiation" 
label var surv2 "Radical hysterectomy"
sort _t
twoway scatter surv1 _t, c(stairstep) ms(i) ///
 || scatter surv2 _t, c(stairstep) ms(i) ///
  ti("Overall-survival") ytitle("Survival probability") xtitle("Years from diagnosis") ylabel(0.00(0.25)1.00)

rename surv1 CRT
rename surv2 RH
list CRT RH in -1

However, I’m not an expert on this field, so hopefully, there is someone who can verify or comment on this method.


Kind regards, Ester Olthof

[Ester Olthof]

Anyone who can confirm?
I would like to use this approach but I'm not sure about it...

[Dana MS]

Hi Esther,
Just checking whether you had confirmation as this code is exactly what I am trying to do

[Darina Peycheva]

Hi, has anyone been successful with this? I have the same problem.