st: Cox regression using a shared frailty model in multiply imputed data

[Justin Schaffer <justin.schaffer@gmail.com>]

Hello, first time poster here (please excuse my ignorance).

Stata (using Stata 13) will allow me to create a Cox regression model
with shared frailty on a multiply imputed dataset. However, it does
not give me an estimate of theta after running the command as it does
when I run the regression model with shared frailty on non-imputed
data. Can someone explain why this is, and whether I am violating some
obscure law of statistics when I create a Cox regression model with
shared frailty on my imputed dataset? I assume that estimating the
standard error of theta is not statistically valid on an imputed
dataset which is why the theta value is not shown, but to be honest, I
am out of my league here. One option I'm considering is running the
regression on each of my imputed data sets after creating a separate
file for each imputed dataset (using mi set flongsep) and then
averaging the theta values and estimating the standard error using
Rubin's rules (although I have no clue if this is statistically
valid). Notably, the variable that I regress on as well as the
"random" variable of my shared frailty model need not be imputed (i.e.
in the example below I run a Cox regression with shared frailty on the
variable "age" with the variable "hosp_id" as my random effect
variable, although neither "age" nor "hosp_id" were actually imputed
in my dataset).

Example:

stset daysAlive, fail(dead)
sts generate nelsonAalen = na

mi set wide model01
mi register regular age hosp_id nelsonAalen dead
mi register impute imputed_variable
mi impute chained (pmm) imputed_variable = age hosp_id dead
nelsonAalen, add(20) augment chaindots burnin(10) rseed(1234)

mi stset daysAlive, fail(dead==1)

//Cox model with shared frailty on imputed dataset
mi estimate, hr dots: stcox age, shared(hosp_id)
> Imputations (20):
>   .........10.........20 done
>
> Multiple-imputation estimates                     Imputations     =   20
> Cox regression: Breslow method for ties           Number of obs   = 3174
>                                                   Average RVI     = 0.0000
>                                                   Largest FMI     = 0.0000
> DF adjustment:   Large sample                     DF:     min     =    .
>                                                           avg     =    .
>                                                           max     =    .
> Model F test:       Equal FMI                     F(   1,      .) =  31.39
> Within VCE type:          OIM                     Prob > F        = 0.0000
>
> ------------------------------
------------------------------------------------
>           _t | Haz. Ratio   Std. Err.      t    P>|t|     [95% Conf. Interval]
> -------------+----------------------------------------------------------------
>          age |   1.030194   .0054693     5.60   0.000      1.01953  1.040969
> ------------------------------------------------------------------------------

//Cox model with shared frailty on non-imputed dataset
stcox age, shared(hosp_id)
>          failure _d:  dead == 1
>    analysis time _t:  daysAlive
>
> Fitting comparison Cox model:
>
> Estimating frailty variance:
>
> Iteration 0:   log profile likelihood = -8334.2065
> Iteration 1:   log profile likelihood = -8334.2065  (backed up)
> Iteration 2:   log profile likelihood = -8333.3181
> Iteration 3:   log profile likelihood = -8333.2928
> Iteration 4:   log profile likelihood = -8333.2927
>
> Fitting final Cox model:
>
> Iteration 0:   log likelihood = -8369.2282
> Iteration 1:   log likelihood = -8333.5078
> Iteration 2:   log likelihood = -8333.2927
> Iteration 3:   log likelihood = -8333.2927
> Refining estimates:
> Iteration 0:   log likelihood = -8333.2927
>
> Cox regression --
>          Breslow method for ties                Number of obs      = 3174
>          Gamma shared frailty                   Number of groups=        60
> Group variable: hosp_id
>
> No. of subjects =         3174                  Obs per group: min=         1
> No. of failures =         1138                                 avg = 52.9
> Time at risk    =      2951413                                 max=       283
>
>                                                 Wald chi2(1)       =31.39
> Log likelihood  =   -8333.2927                  Prob > chi2        =0.0000
>
> ------------------------------------------------------------------------------
>           _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf.Interval]
> -------------+----------------------------------------------------------------
>          age |   1.030194   .0054693     5.60   0.000      1.01953  1.040969
> -------------+----------------------------------------------------------------
>        theta |   .0391666   .0190069
> ------------------------------------------------------------------------------
> Likelihood-ratio test of theta=0: chibar2(01) =    11.50 Prob>=chibar2 =0.000
>
> Note: standard errors of hazard ratios are conditional on theta.

Many thanks in advance for any statistical advice you gurus have to offer.

Sincerely,

JMS
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/