3
$\begingroup$

My magelastbirth_pctl is a categorical variable with 5 levels. I want to calculate the hazard ratio of the 4 last categories compared to the first categories using Fine and Gray competing risk model.

This is easy to do with normal cox model using the following code

elcs.cox_CoD = coxph(Surv(fuptime, status_CoD==17) ~ factor(magelastbirth_pctl) +
                                                     factor(magefirstbirth_pctl) +
                                                     totchildren +
                                                     prochildrendeceased_bef18 + byr_cat +
                                                     factor(npses_husb_high_qtl),
                     data = elcs6)

However, I don't know how to do it using crr command. It works if I treat the magelastbirth_pctl as a continuous variable per the following code

elcs.cmprsk = crr(ftime=elcs6$fuptime, fstatus=elcs6$status_CoD, failcode=1, cencode=0,
              cov1=elcs.cov)

I tried to use the FGR command of the riskRegressions package

test = FGR(Hist(fuptime,status_CoD)~factor(magelastbirth_pctl)+factor(magefirstbirth_pctl) + 
                                    totchildren + prochildrendeceased_bef18 + byr_cat,
           data=elcs6)

but it kept giving me the error "Argument cause missing. Analyse cause: 1"

Anybody knows what I'm doing wrong? There are surprisingly limited examples in cmprsk and none of them mentioned categorical data.

| cite | improve this question | |
$\endgroup$
3
$\begingroup$

Create a model.matrix:

cov1 <- model.matrix(~ factor(magelastbirth_pctl) + factor(magefirstbirth_pctl) +
                       totchildren + prochildrendeceased_bef18 + byr_cat +
                       factor(npses_husb_high_qtl),
                     data = elcs6)[, -1]
crr(ftime=elcs6$fuptime, fstatus=elcs6$status_CoD, failcode=1, cencode=0, cov1=cov1)

from ?cmprsk::crr

the model.matrix function can be used to generate suitable matrices of covariates from factors, eg model.matrix(~factor1+factor2)[,-1] will generate the variables for the factor coding of the factors factor1 and factor2. The final [,-1] removes the constant term from the output of model.matrix

For a formula interface, see survival::finegray

library('cmprsk')
fg <- finegray(Surv(futime, event) ~ age + sex + abo, na.omit(transplant),
               etype = 'death')
coef(
  coxph(Surv(fgstart, fgstop, fgstatus) ~ age + sex + abo, fg, weights = fgwt)
)
#        age        sexf        aboB       aboAB        aboO 
# 0.01890731 -0.42597601  0.39092581  0.13277852  0.38908033

Or cmprsk2 which is still under development but includes some useful functions for cmprsk including a formula interface, crr2

# devtools::install_github('raredd/cmprsk2')
library('cmprsk2')
c2 <- crr2(Surv(futime, event(censored) == death) ~ age + sex + abo,
           na.omit(transplant))
coef(c2$`CRR: death`)
#        age        sexf        aboB       aboAB        aboO 
# 0.01890467 -0.42585553  0.39078253  0.13283606  0.38898005

If you look at the call for the model, it will tell you what needed to be ran in order to get this fit (including how to construct the model matrix)

c2$`CRR: death`$call

# crr(na.omit(transplant)[, "futime"], na.omit(transplant)[, "event"], 
#     cov1 = model.matrix(~age + sex + abo, na.omit(transplant))[, 
#         -1L, drop = FALSE], cencode = "censored", failcode = "death", 
#     variance = TRUE)

eval(c2$`CRR: death`$call)

# convergence:  TRUE 
# coefficients:
#     age    sexf    aboB   aboAB    aboO 
#  0.0189 -0.4259  0.3908  0.1328  0.3890 
# standard errors:
# [1] 0.01116 0.25350 0.38170 0.62300 0.28140
# two-sided p-values:
#   age  sexf  aboB aboAB  aboO 
# 0.090 0.093 0.310 0.830 0.170
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ cmprsk2 seems great! However, why are the coefficients slightly different between the 2 function? $\endgroup$ – Dan Chaltiel Dec 22 '19 at 11:44
  • 1
    $\begingroup$ @DanChaltiel hmm, good point. I'm not sure, but it seems to do with how finegray handles discrete time (or times with ties). For example, if you add a tiny amount of noise to the time variable, you get the same coefficients: run transplant$futime <- transplant$futime + runif(nrow(transplant), 0, 0.1) first then fit the models. There is a line in finegray which is newtime <- matrix(findInterval(Y[, 1:2], utime), ncol = 2) where utime are unique times, so that may be the culprit $\endgroup$ – rawr Dec 25 '19 at 2:36
1
$\begingroup$

The abovementioned answer is great, just to mention one tiny detail in the matrix;

cov1 <- model.matrix(~factor(magelastbirth_pctl) + factor(magefirstbirth_pctl) +
                       totchildren + prochildrendeceased_bef18 + byr_cat +
                       factor(npses_husb_high_qtl),
                     data = elcs6)[, -1]

The ~ sign is essential to prevent errors like '+' is not meaningful for factors

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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