2
$\begingroup$

I need your help with my problem. So, after step of imputation missing data through MICE method, I got multiple imputed dataset. Then, I pooled the estimates and coefficients with mixed effect cox model. But in the final result, I got NA for the random effect terms. Is it normal ?

For more infos, I have 27 860 patients in the database which has a multilevel structure. In the database, I have a character variable called "Area" that represents an area of living of people. For the variable Area, there are 17 000 values distincts.

My goal is to examine the role spatial accessibility to a medical event. I intended to use the mixed effects cox model with Area as random effect.

Here is the steps of my analysis:

watt_final_imp <- mice(watt_final, m = 5, maxit = 10, seed = 123, method = "rf") #I used rf method because all of the missing data variables are factors.
cox.fit <- with(data_imputed, coxph(Surv(time, event_status) ~ age + sex + ...  + smoke + spatial_accessibility +  frailty(Area)))
cox.fit
pool.fit <- pool(cox.fit)
summary(pool.fit)

Below are the results of one imputed dataset:

                             coef      se(coef)    se2      Chisq      DF       p
age                         -4.29e-02  1.10e-03  1.01e-03  1.52e+03    1     < 2e-16
....
SmokeYes                     1.46e-01  3.32e-02  3.08e-02  1.94e+01    1     1.0e-05
Spatial_accessibilityHigh   -2.43e-02  9.06e-02  8.17e-02  7.16e-02    1     0.78908
frailty(Area)                                              2.30e+03  1890    2.0e-10

Iterations: 6 outer, 34 Newton-Raphson
     Variance of random effect= 0.34   I-likelihood = -63326.5 
Degrees of freedom for terms=    0.8    0.9    0.8    0.8    1.7    0.9    3.4    0.9    2.8    4.3    1.9    0.9    0.9    0.9    0.9    0.9    2.5    0.8    0.8    0.9    0.9    1.7    0.8    0.8 1889.6 
Likelihood ratio test=13559  on 1923 df, p=<2e-16
n= 27860, number of events= 6791

Below are the results after pooling estimates and coeff:

                         term      estimate    std.error    statistic         df      p.value
1                        age    -0.042260183 0.001401088 -30.16240416   25.05626 3.248679e-21
.....
37                  SmokeYes     0.152341302 0.036133085   4.21611667  134.33406 4.536397e-05
38    Spatial_accessibilityHigh -0.021317365 0.089576125  -0.23798043 5174.21229 8.119057e-01
39        frailty(Area)             NA          NA           NA         NA           NA

Thanks in advance

Improve this question
$\endgroup$
1
  • 2
    $\begingroup$ I am not an expert in this area but since you do not get a coefficient for frailty in the individual mode why did you expect pooling to give you one? $\endgroup$
    – mdewey
    Commented Feb 26 at 13:59

1 Answer 1

Reset to default
3
$\begingroup$

First, I'm a bit worried about a frailty term having 17000 distinct categorical values when you only have 27860 individuals and 6791 events. At first glance at least, that would seem to make it impossible to get unique estimates for frailties for most of the levels of Area, those that lack an event. It's not clear how much that type of frailty term can help here. A frailty term can help to take unmodeled outcome-associated predictors into account, but with fewer than 2 individuals per Area on average (and less than 1 event per Area) it's not clear how well that will work with your data.

Second, to answer your question, so far as I know the pool() function in the mice package does not pool random effects across the multiple models. Some other packages do that for the Gaussian random effects used in typical mixed models. Those include the merTools and mitml packages; I don't know whether they can handle Cox mixed models directly. It's also not clear that those methods would apply to your default choice of a (skewed) gamma-distributed frailty term. If you used a Gaussian frailty term instead, you might consider using the code in those packages to help you get corresponding random-effect estimates for your coxph() models.

Finally, if you're mainly interested in the fixed effects, it might be adequate just to display the list of Variance of random effect values for the models.

Improve this answer
$\endgroup$
2
  • $\begingroup$ thanks for your answer. Well, I know that putting a random effect here is dumb because Area with only 1 patient will get a variance as 0. But my supervisor is forced me to do it bcz the database is such a multilevel. Anw, do you know any model can fit these type of data ? Or just a normal cox is ok. I'm interested in the survival time.between Area but it seems that I don't have enough data to do that. $\endgroup$
    – Hoang-Giang Pham
    Commented Mar 1 at 9:14
  • $\begingroup$ @Hoang-GiangPham if Areas are themselves grouped at a higher level, say into Provinces, then you could consider using frailties at that higher level, or move to the coxme package that (unlike coxph()) allows for multiple levels of random effects. Alternatively, use coxph() with a cluster(Area) term. That won't model the variance among Area values directly, but it will give a "robust" error estimate that takes the Area groupings into account. Or stick with your current model and just report the list of Variance of random effect values. $\endgroup$
    – EdM
    Commented Mar 1 at 9:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.