Guidance for estimating multilevel survival models with "mi" (multiple imputation)

[Michael Harhay]

Dear Statalist members,

I am using Stata 15.1 (last updated 15 oct 2018). I am asking for general and specific guidance for estimating multilevel survival models with "mi" (multiple imputation).

I have used "mi impute chained" to impute clinical data from a transplant registry. I would like to then estimate a multilevel survival model. To the best I can understand, "mestreg" does not work with the mi suite based on a prior post from Feb 2016, so I am using a Cox shared frailty model. Below is the output for the following commands.

mi stset ptime_update, failure(death) scale(365.25)
mi estimate: stcox i.race_granular $covariates
mi estimate: stcox i.race_granular, strata(centre)
mi estimate: stcox i.race_granular $covariates, shared(centre)
stcox i.race_granular, shared(centre)

My specific questions to the community are the following:
1. Can Stata estimate stcox, shared(cluster)with mi as I did below without concern?
2. If yes, can someone explain me how to confirm and evaluate the estimation? The "strata" command works clearly from the output, and the shared estimates look like my non mi estimates, so it seems to be working, but I'd like to have some better confirmation.
3. If no to #1, is there an alternative for a time-to-event model in Stata other than using strata(cluster)?

Many thanks for your guidance.

-Michael

*** Output for Statalist
. mi stset ptime_update, failure(death) scale(365.25)
 
     failure event:  death != 0 & death < .
obs. time interval:  (0, ptime_update]
 exit on or before:  failure
    t for analysis:  time/365.25
 
------------------------------------------------------------------------------
     21,217 total observations
          0 exclusions
------------------------------------------------------------------------------
     21,217 observations remaining, representing
      8,743 failures in single-record/single-failure data
 71,196.252 total analysis time at risk and under observation
                                               at risk from t =         0
                                     earliest observed entry t =         0

mi estimate, cmdok: mestreg i.race_granular || centre:, distribution(exponential) time tratio
command prefix mestreg i.race_granular || centre: not allowed
 
** covariates only: mi
. mi estimate: stcox i.race_granular $nmiss $miss
 
Multiple-imputation estimates                  Imputations       =         10
Cox regression: Breslow method for ties        Number of obs     =     21,216
                                               Average RVI       =     0.0035
                                               Largest FMI       =     0.0394
DF adjustment:   Large sample                   DF:     min      =   5,894.53
                                                       avg       =   1.80e+09
                                                       max       =   3.70e+10
Model F test:       Equal FMI                   F(  27, 1.9e+07) =      21.34
Within VCE type:          OIM                   Prob > F          =    0.0000
 
---------------------------------------------------------------------------------------
                   _t |      Coef.  Std. Err.      t    P>|t|    [95% Conf. Interval]
----------------------+----------------------------------------------------------------
        race_granular |
  Non-Hispanic Black  |  .0041659   .0390737     0.11  0.915    -.0724172     .080749
            Hispanic  | -.1450088    .050287    -2.88  0.004    -.2435694   -.0464481
               Other  | -.1068894   .0760553    -1.41  0.160    -.2559551    .0421763
                      |
          tx_type_des |   .2523575  .0246634    10.23   0.000    .2040182    .3006969
                  age |   .0082351  .0010614     7.76   0.000    .0061548    .0103154
               female |  -.0854246  .0246798    -3.46   0.001   -.1337961   -.0370531
         grouping_des |   .0042175  .0098095     0.43  0.667    -.0150087    .0234437
        insurance_trr |   .0349039  .0104485     3.34   0.001    .0144253    .0553825
        recipient_edu |  -.0110632  .0082076    -1.35   0.178   -.0271497    .0050234
              abo_cat |   .0067155  .0077683     0.86   0.387   -.0085101    .0219411
        end_match_las |   .0010782  .0009063     1.19   0.234   -.0006981    .0028545
              age_don |   .0032827  .0008307     3.95   0.000    .0016545    .0049109
         female_donor |   .0530582  .0249816     2.12   0.034    .0040952    .1020212
             race_don |   .1408148  .0226417     6.22   0.000    .0964379    .1851917
         pulm_inf_don |  -.0314945  .0219263    -1.44   0.151   -.0744693    .0114803
          cod_cad_don |   -.000082  .0000711    -1.15   0.248   -.0002213    .0000572
            end_creat |   .0942102  .0161607     5.83   0.000    .0625358    .1258846
                 vent |   .2152456  .0522129     4.12   0.000    .1129103     .317581
                 ecmo |   .113633   .1279242     0.89  0.374    -.1370939    .3643598
          vent_p_ecmo |   .1212142  .1574985     0.77   0.442   -.1874772    .4299056
     hist_cig_don_num |   .1203965  .0341459     3.53   0.000    .0534698    .1873232
hist_oth_drug_don_num |   .0597952   .0237674    2.52   0.012     .0132093     .106381
             bmi_calc |   .0021126  .0026663     0.79   0.428   -.0031132    .0073385
         bmi_don_calc |  -.0015305  .0021156    -0.72   0.469   -.0056771     .002616
         six_min_walk |  -.0001698  .0000279    -6.10   0.000   -.0002244   -.0001152
                   pa |   .0011132  .0007483     1.49   0.137   -.0003538    .0025802
                 po2  |  -.0000418   .0000725   -0.58   0.564    -.0001839   .0001002
 ---------------------------------------------------------------------------------------

 
*** strata works
. mi estimate: stcox i.race_granular, strata(centre)
 
Multiple-imputation estimates                  Imputations       =         10
Stratified Cox regr.: Breslow method for ties   Number of obs     =     21,217
                                               Average RVI       =     0.0000
                                               Largest FMI       =     0.0000
DF adjustment:   Large sample                   DF:     min      =          .
                                                       avg       =          .
                                                       max       =          .
Model F test:       Equal FMI                   F(   3,     .)   =       3.67
Within VCE type:          OIM                   Prob > F          =    0.0117
 
-------------------------------------------------------------------------------------
                 _t |     Coef.   Std. Err.      t   P>|t|     [95% Conf. Interval]
--------------------+----------------------------------------------------------------
      race_granular |
Non-Hispanic Black  | -.0517968   .0389612    -1.33  0.184    -.1281594    .0245658
          Hispanic  | -.1552613   .0511899    -3.03  0.002    -.2555916   -.0549309
             Other  | -.0756223   .0763356    -0.99  0.322    -.2252372    .0739927
-------------------------------------------------------------------------------------
 
*** trying a mi shared frailty model. It estimates, but no indication that it is multilevel. Estimates are different than covariate only model. 
mi estimate: stcox i.race_granular $nmiss $miss, shared(centre)
 
Multiple-imputation estimates                  Imputations       =         10
Cox regression: Breslow method for ties        Number of obs     =     21,216
                                               Average RVI       =     0.0037
                                                Largest FMI       =    0.0400
DF adjustment:   Large sample                   DF:     min      =   5,709.24
                                                       avg       =   1.49e+09
                                                       max       =  3.08e+10
Model F test:       Equal FMI                   F(  27, 1.8e+07) =      19.97
Within VCE type:          OIM                   Prob > F          =    0.0000
 
---------------------------------------------------------------------------------------
                   _t |      Coef.  Std. Err.      t    P>|t|    [95% Conf. Interval]
----------------------+----------------------------------------------------------------
        race_granular |
  Non-Hispanic Black  | -.0223719   .0397676    -0.56  0.574    -.1003149    .0555711
            Hispanic  | -.1831117   .0516951    -3.54  0.000    -.2844321   -.0817912
               Other  | -.1199062   .0764192    -1.57  0.117    -.2696851    .0298728
          tx_type_des |   .2595752  .0263676     9.84   0.000    .2078957    .3112546
                  age |    .007673  .0010875     7.06   0.000    .0055416    .0098044
               female |  -.0880246  .0248277    -3.55   0.000    -.136686   -.0393632
         grouping_des |   .0065191  .0099099     0.66   0.511   -.0129039    .0259422
        insurance_trr |   .0368581  .0106699     3.45   0.001    .0159454    .0577708
        recipient_edu |  -.0020136  .0085693    -0.23   0.814   -.0188091    .0147818
              abo_cat |   .0050944  .0077766     0.66   0.512   -.0101475    .0203363
        end_match_las |   .0011873  .0009394     1.26   0.206    -.000654    .0030285
              age_don |   .0035487  .0008396     4.23   0.000    .0019031    .0051943
         female_donor |   .0560035   .025171     2.22   0.026    .0066693    .1053377
             race_don |   .1262233  .0231297     5.46   0.000    .0808898    .1715567
         pulm_inf_don |  -.0365005    .02236    -1.63   0.103   -.0803252    .0073243
          cod_cad_don |  -.0000799  .0000712    -1.12   0.262   -.0002194    .0000597
            end_creat |   .0894781  .0167006     5.36   0.000    .0567456    .1222107
                 vent |   .2415701  .0539368     4.48   0.000    .1358559    .3472844
                 ecmo |   .1160986   .128481     0.90   0.366   -.1357195    .3679167
          vent_p_ecmo |    .094155  .1582924     0.59   0.552   -.2160923    .4044024
     hist_cig_don_num |   .1158391   .034369     3.37   0.001    .0484747    .1832036
hist_oth_drug_don_num |   .0535818   .0238907    2.24   0.025     .0067546   .1004091
             bmi_calc |   .0015292  .0026879     0.57   0.569   -.0037389    .0067974
         bmi_don_calc |  -.0020372  .0021261    -0.96   0.338   -.0062043    .0021299
         six_min_walk |  -.0001853  .0000314    -5.90   0.000   -.0002468   -.0001238
                   pa |   .0013551  .0007552     1.79  0.073    -.0001255    .0028357
                  po2 |  -7.80e-06  .0000749    -0.10   0.917   -.0001545    .0001389
---------------------------------------------------------------------------------------
 
 
** complete case model, which shows center estimation 
. stcox i.race_granular, shared(centre)
 
         failure _d:  death
   analysis time _t:  ptime_update/365.25
 
Fitting comparison Cox model:
 
Estimating frailty variance:
 
Iteration 0:   log profile likelihood =  -79831.26 
Iteration 1:   log profile likelihood = -79827.878  
Iteration 2:   log profile likelihood = -79827.866  
Iteration 3:   log profile likelihood = -79827.866  
 
Fitting final Cox model:
 
Iteration 0:   log likelihood = -79962.988
Iteration 1:   log likelihood = -79829.421
Iteration 2:   log likelihood = -79827.867
Iteration 3:   log likelihood = -79827.866
Refining estimates:
Iteration 0:   log likelihood = -79827.866
 
Cox regression -- Breslow method for ties
 
Gamma shared frailty                           Number of obs     =     21,217
Group variable: centre                         Number of groups  =         78
                                               Obs per group:
No. of subjects =       21,217                                min =          1
No. of failures =        8,743                                avg =  272.01282
Time at risk    =  71196.25188                                max =      1,190
 
                                               Wald chi2(3)      =       9.69
Log likelihood  =   -79827.866                  Prob > chi2       =    0.0214
 
-------------------------------------------------------------------------------------
                 _t | Haz. Ratio   Std. Err.      z   P>|z|     [95% Conf. Interval]
--------------------+----------------------------------------------------------------
      race_granular |
Non-Hispanic Black  |  .9632818   .0371996    -0.97  0.333     .8930628    1.039022
          Hispanic  |  .8613009   .0436307    -2.95  0.003     .7798944    .9512047
             Other  |  .9370008   .0712056    -0.86  0.392     .8073364     1.08749
--------------------+----------------------------------------------------------------
              theta |   .0350027  .0087813
-------------------------------------------------------------------------------------
LR test of theta=0: chibar2(01) = 168.80               Prob >= chibar2 = 0.000

[Steve Samuels]

Although they do not appear in the output, theta and its standard error are in the e-returned results of the mi model with shared frailty. To see them, type:

di e(theta)
di e(se_theta)

[Michael Harhay]

Thank you Steve for finding the way to confirm this for me.
-Michael

[Michael Harhay]

For completeness, I wanted to add the response from Stata Technical support.

***
The command -mestreg- doesn't currently work with -mi estimate, cmdok-
just because of a parsing issue, but you can safely use it by
writing a wrapper, for example:


**** begin code
capture program drop mymestreg
program mymestreg
mestreg age female || patient:, distribution(weibull)
end

webuse catheter, clear
replace female = . in 35

mi set flong
mi register imputed female
mi impute logit female age time infect, add(5) rseed(123)
mi stset time, failure(infect)

mi estimate, saving(miest,replace) cmdok: mymestreg
**** end code

Concerning -stcox-, you can see the output for the individual
imputations by using the option -noisily- with the -mi estimate-
command. You will see that the individual calls are actually
accounting for the -shared()- option, so you will get the
coefficients that you want with -mi estimate-. It will not compute
the mi-combined estimate for the frailty parameter theta, though.
If you need it, you can obtain it also with a wrapper.

[Steve Samuels]

Thanks for sharing this information. To write the wrapper for mi-combined estimate of theta and its estimated standard error, see Isabel Cañette and Yulia Marchenko's FAQ
How can I combine results other than coefficients in e(b) with multiply imputed data?

Perhaps your source is mistaken and e(theta) and e(se_theta) are the mi-combined estimates. Else, what are they?

[Laura Beres]

Dear Statalist Members,
I am have a similar question to the thread above. I am not able to see the estimate of theta or the se of theta using e-returned results. It appears as if they are not being estimated. If the estimate of theta is not being returned, does that mean that the coefficients are not accounting for the shared option? Is there somewhere else I need to look? Noisily option with the MI estimate returns output similar to the MI estimate without noisily. Is there another command I need to appropriately use the shared option with MI?
When I run the complete case model without MI I see the estimated theta and SE of theta in the output, as well as the likelihood ratio test of theta =0. I am using STATA 15.1, updated in February 2019. Very grateful for your assistance.

Code:
mi set wide
mi register imputed lastcd4cat ill_at_enrollment
mi impute chained (mlogit) lastcd4cat ill_at_enrollment = Male ageatinterview_cat_2 facility_cat ///
b110_disclosed_hivstatus_r visittointv_weeks1 visittointv_weeks2 hourstoclinic_cat i_clinic_contact_for_return_c i_b19a_month_away ///
i_i_cat_anticipated i_b150_confronted_r i_b15h_any_structural i_b15h_any_psychosocial i_b15h_any_clinic i_b15h_any_medical ///
i_b15i_any_structural i_b15i_any_psychosocial i_b15i_any_clinic i_b15i_returning b04_clinic, add(10) rseed(11111)
mi stset days_return_or_censor, failure(returner_SC==1)
mi estimate, hr: stcox Male i.ageatinterview_cat_2 i.lastcd4cat ill_at_enrollment i.facility_cat ///
b110_disclosed_hivstatus_r visittointv_weeks1 visittointv_weeks2 i.hourstoclinic_cat i.i_clinic_contact_for_return_c i_b19a_month_away ///
i.i_i_cat_anticipated i.i_b150_confronted_r i_b15h_any_structural i_b15h_any_psychosocial i_b15h_any_clinic i_b15h_any_medical ///
i_b15i_any_structural i_b15i_any_psychosocial i_b15i_any_clinic i_b15i_returning, shared(b04_clinic)
ereturn list

*with noisily, following the guidance in the thread above to see how the individual imputations are called
noisily: mi estimate, hr: stcox Male i.ageatinterview_cat_2 i.lastcd4cat ill_at_enrollment i.facility_cat ///
b110_disclosed_hivstatus_r visittointv_weeks1 visittointv_weeks2 i.hourstoclinic_cat i.i_clinic_contact_for_return_c i_b19a_month_away ///
i.i_i_cat_anticipated i.i_b150_confronted_r i_b15h_any_structural i_b15h_any_psychosocial i_b15h_any_clinic i_b15h_any_medical ///
i_b15i_any_structural i_b15i_any_psychosocial i_b15i_any_clinic i_b15i_returning, shared(b04_clinic)
ereturn list

*complete case model
stcox Male i.ageatinterview_cat_2 i.lastcd4cat ill_at_enrollment i.facility_cat ///
b110_disclosed_hivstatus_r visittointv_weeks1 visittointv_weeks2 i.hourstoclinic_cat i.i_clinic_contact_for_return_c i_b19a_month_away ///
i.i_i_cat_anticipated i.i_b150_confronted_r i_b15h_any_structural i_b15h_any_psychosocial i_b15h_any_clinic i_b15h_any_medical ///
i_b15i_any_structural i_b15i_any_psychosocial i_b15i_any_clinic i_b15i_returning, shared(b04_clinic)

Partial Output (log file attached):

Multiple-imputation estimates	Imputations	= 10
Cox regression: Breslow method for ties	Number of obs	= 571
Average RVI	= 0.0249
Largest FMI	= 0.1569
DF adjustment: Large sample	DF: min	= 386.79
avg	= 3649206.19
max	= 5.19e+07
Model F test: Equal FMI	F( 30,442811.3)	= 1.54
Within VCE type: OIM	Prob > F	= 0.0300
_t Haz. Ratio	Std. Err. t	P>t [95% Conf.	Interval]
Male .807657	.1315787 -1.31	0.190 .5868842	1.11148
ageatinterview_cat_2
1=25-34 .9356844	.2515952 -0.25	0.805 .5523956	1.584925
2=35-44 1.105942	.3022218 0.37	0.713 .647324	1.889482
3=45+ 1.401485	.4357413 1.09	0.278 .761969	2.577743
lastcd4cat
101-350 1.301411	.5620558 0.61	0.542 .556818	3.041697
351-500 1.129187	.5014454 0.27	0.784 .4723363	2.69948
>500 1.046443	.4735666 0.10	0.920 .4298296	2.547622
ill_at_enrollment 1.040607	.1982324 0.21	0.835 .7160659	1.512238
facility_cat
urban .7028782	.1601835 -1.55	0.122 .4496698	1.098668
hospital .5145202	.149297 -2.29	0.022 .2913473	.9086445
b110_disclosed_hivstatus_r 1.100414	.2565641 0.41	0.682 .6967808	1.737864
visittointv_weeks1 .9891687	.0070779 -1.52	0.128 .9753931	1.003139
visittointv_weeks2 .9964407	.0038782 -0.92	0.360 .9888686	1.004071
hourstoclinic_cat
1 -	.1944398 0.30	0.761 .7375328	1.516327
2 hours + .6915834	.1325405 -1.92	0.054 .4750217	1.006875
i_clinic_contact_for_return_c
Contacted 1-3 times .8932873	.2166357 -0.47	0.642 .5553438	1.43688
Contacted 4-12 times 3.185982	1.395717 2.65	0.008 1.350045	7.518625
i_b19a_month_away .9338023	.1461428 -0.44	0.662 .6871305	1.269026
i_i_cat_anticipated
1=Medium 1.043006	.207763 0.21	0.833 .7058793	1.541144
2=High 1.299374	.2395313 1.42	0.155 .9053582	1.864866
i_b150_confronted_r
1=Once 1.607536	.440144 1.73	0.083 .9399425	2.749288
2=More than once .5432483	.175894 -1.88	0.059 .2880014	1.024713
i_b15h_any_structural 1.053098	.2002711 0.27	0.786 .7254262	1.528779
i_b15h_any_psychosocial .935808	.1600014 -0.39	0.698 .6693464	1.308346
i_b15h_any_clinic 1.082094	.191139 0.45	0.655 .7654376	1.52975
i_b15h_any_medical 1.143014	.2052839 0.74	0.457 .8038531	1.625273
i_b15i_any_structural 1.166824	.2484351 0.72	0.469 .7687249	1.771087
i_b15i_any_psychosocial .7468199	.1512369 -1.44	0.149 .5021585	1.110685
i_b15i_any_clinic .9965548	.1795507 -0.02	0.985 .7000659	1.418611
i_b15i_returning .9250287	.1654049 -0.44	0.663 .6515553	1.313285
. ereturn list
scalars:
e(df_m) = .
e(se_theta) = .
e(chi2_c) = .
e(ll) = .
e(g_min) = 5
e(risk) = 271374
e(converged) = 1
e(theta) = .
e(N) = 571