2
$\begingroup$

I’m doing a registry based study with almost 200 000 observations and I want to perform a competing risk analysis. My problem is that the crr() in the cmprsk package is exponentially increasing with increasing number of observations. I therefore wrote a simulation for trying different approaches; check how factors, data frames and matrixes affect the performance so that I could choose the most efficient combination.

I have a 1 year old computer with 8 Gb of RAM and still it didn’t finish 70 000 observations when I left the computer overnight.

My main questions:

  • Is there a faster way of performing competing risk analysis?
  • Why does Win7 4 times perform better? Is it the 64-bit version that improves the performance?
  • Can I do something to speed things up?
  • Is the simulation similarly slow on your computer? (see simulation code at the end)

Thanks Max


My output with some of the errors that some of the alternatives cause

* Start crr performance check *
Using 10000 observations and 3 - 10000 covariates
..
Error: NA/NaN/Inf in foreign function call (arg 4)
In addition: Warning message:
NAs introduced by coercion
...
Error in drop(.Call("La_dgesv", a, as.matrix(b), tol, PACKAGE = "base")) :
  Lapack routine dgesv: system is exactly singular
..

The output in Linux

- Linux – Ubuntu – R32-bit 2.13.1
 *** Test results ***

Test 1 passed - using no factors
It took: 4.2 minutes

Test 2 passed - using factored status
Benchmark test 2: 91 %

Test 3 failed - factored covariates fail

Test 4 passed - using factored status & data frame covariates
Benchmark test 4: 78 %

Test 5 passed - using factored status & twice the covariates
                where the new covariates differ
Benchmark test 5: 76 %

Test 6 failed - exact copies of the covariates fail

Test 7 passed - using regular status & twice the covariates
                where the new covariates differ
Benchmark test 7: 80 %

Test 8 passed - using regular status & thrice the covariates
                where the new covariates differ
Benchmark test 8: 83 %

The output in Windows

- Windows 7 – R64-bit 2.13.1

 *** Test results ***


 *** Test results ***

Test 1 passed - using no factors
It took: 1.1 minutes

Test 2 passed - using factored status
Benchmark test 2: 99 %

Test 3 failed - factored covariates fail

Test 4 passed - using factored status & data frame covariates
Benchmark test 4: 104 %

Test 5 passed - using factored status & twice the covariates 
                where the new covariates differ
Benchmark test 5: 122 %

Test 6 failed - exact copies of the covariates fail

Test 7 passed - using regular status & twice the covariates 
                where the new covariates differ
Benchmark test 7: 122 %

Test 8 passed - using regular status & thrice the covariates 
                where the new covariates differ
Benchmark test 8: 144 %

Test with 20 000 (twice the data) on Windows 7 R64 2.13.1 (removed the failures)

*** Test results ***
Test 1 passed - using no factors
It took: 8.5 minutes

Test 2 passed - using factored status
Benchmark test 2: 100 %

Test 4 passed - using factored status & data frame covariates
Benchmark test 4: 100 %

Test 5 passed - using factored status & twice the covariates 
                where the new covariates differ
Benchmark test 5: 123 %

Test 7 passed - using regular status & twice the covariates 
                where the new covariates differ
Benchmark test 7: 122 %

Test 8 passed - using regular status & thrice the covariates 
                where the new covariates differ
Benchmark test 8: 147 %

The performance test code

library(cmprsk)
# Define a set size
my_set_size <- 10000

# Create the covariables
cov <- cbind(rbinom(my_set_size, 1, .5), 
  rbinom(my_set_size, 1, .05),
  rbinom(my_set_size, 1, .1))
dimnames(cov)[[2]] <- c('gender','risk factor 1','risk factor 2')
factored_cov <- data.frame(
  gender = factor(cov[,1], 
  levels=c(0,1), 
  labels=c("male", "female")),
  rf1 = factor(cov[,2], 
  levels=c(0,1), 
  labels=c("no", "yes")),
  rf2 = factor(cov[,3], 
  levels=c(0,1),
  labels=c("no", "yes")))

framed_cov <- data.frame(
  gender = cov[,1],
  rf1 = cov[,2], 
  rf2 = cov[,3])

cov_doubled <- cbind(cov[,1], 
  cov[,2],
  cov[,3],
  rbinom(my_set_size, 1, .5),
  rbinom(my_set_size, 1, .05),
  rbinom(my_set_size, 1, .1))
dimnames(cov_doubled)[[2]] <- 
  c('gender','risk factor 1','risk factor 2',
    'risk factor 3', 'risk factor 4', 'risk factor 5')

cov_doubled_exact <- cbind(cov[,1], 
  cov[,2],
  cov[,3],
  cov[,1], 
  cov[,2],
  cov[,3])
dimnames(cov_doubled_exact)[[2]] <- 
  c('1: gender','1: risk factor 1','1: risk factor 2',
    '2: gender','2: risk factor 1','2: risk factor 2')

cov_tripled <- cbind(cov[,1], 
  cov[,2],
  cov[,3],
  cov_doubled[,4],
  cov_doubled[,5],
  cov_doubled[,6],
  rbinom(my_set_size, 1, .5),
  rbinom(my_set_size, 1, .05),
  rbinom(my_set_size, 1, .1))
dimnames(cov_tripled)[[2]] <- 
  c('gender','risk factor 1','risk factor 2',
    'risk factor 3', 'risk factor 4', 'risk factor 5',
    'risk factor 6', 'risk factor 7', 'risk factor 8')

# Create random time to failure/cens periods
ftime <- rexp(my_set_size)

# Create events
my_event1 <- rbinom(my_set_size, 1, .04)
my_event2 <- rbinom(my_set_size, 1, .20)
# The competing event can't happen if 1 has already occurred
my_event2[my_event1 > 0] <- 0
fstatus <- my_event1 + my_event2*2

# Factor the censor variable
factored_status <- factor(fstatus, levels=c(0,1,2), labels=c("censored", "re-operation", "death"))

# Check that it seems Ok
#table(fstatus)

# The Sys.time() is a crude but I found it efficient enough 
# for my simple analysis
cat("\n * Start crr performance check *\nUsing", my_set_size, "observations and", 
  length(cov[1,]), "-", length(cov_doubled[,1]), "covariates\n")

# Initialize vars for later if() checks
my_test1 <- my_test2 <- my_test3 <- NULL
my_test4 <- my_test5 <- my_test6 <- NULL
my_test7 <-NULL

my_t1.a <- Sys.time()
my_test1 <- crr(ftime, fstatus, cov)
my_t1.b <- Sys.time()
cat(".")
my_t2.a <- Sys.time()
my_test2 <- crr(ftime, factored_status, cov, failcode="re-operation", cencode="censored")
my_t2.b <- Sys.time()
cat(".")
my_t3.a <- Sys.time()
my_test3 <- crr(ftime, factored_status, factored_cov, failcode="re-operation", cencode="censored")
my_t3.b <- Sys.time()
cat(".")
my_t4.a <- Sys.time()
my_test4 <- crr(ftime, factored_status, framed_cov, failcode="re-operation", cencode="censored")
my_t4.b <- Sys.time()
cat(".")
my_t5.a <- Sys.time()
my_test5 <- crr(ftime, factored_status, cov_doubled, failcode="re-operation", cencode="censored")
my_t5.b <- Sys.time()
cat(".")
my_t6.a <- Sys.time()
my_test6 <- crr(ftime, factored_status, cov_doubled_exact, failcode="re-operation", cencode="censored")
my_t6.b <- Sys.time()
cat(".")
my_t7.a <- Sys.time()
my_test7 <- crr(ftime, fstatus, cov_doubled)
my_t7.b <- Sys.time()
cat(".")
my_t8.a <- Sys.time()
my_test8 <- crr(ftime, factored_status, cov_tripled, failcode="re-operation", cencode="censored")
my_t8.b <- Sys.time()
#summary(my_test)

cat("\n *** Test results ***")

# Don't check if test 1 passed since this is the benchmark
cat("\n\nTest 1 passed - using no factors")
cat("\nIt took:", round(difftime(my_t1.b, my_t1.a, units="mins"), 1), "minutes")
benchmark_time <- as.numeric(my_t1.b - my_t1.a)

if (length(my_test2)){
  cat("\n\nTest 2 passed - using factored status")
  bench.t2 <- as.numeric(my_t2.b - my_t2.a)/benchmark_time
  cat("\nBenchmark test 2:", round(bench.t2*100), "%")
}else{
  cat("\n\nTest 2 failed")
}

if (length(my_test3)){
  cat("\n\nTest 3 passed - using factored status & factored covariates")
  bench.t3 <- as.numeric(my_t3.b - my_t3.a)/benchmark_time
  cat("\nBenchmark test 3:", round(bench.t3*100), "%")
}else{
  cat("\n\nTest 3 failed - factored covariates fail")
}

if (length(my_test4)){
  cat("\n\nTest 4 passed - using factored status & data frame covariates")
  bench.t4 <- as.numeric(my_t4.b - my_t4.a)/benchmark_time
  cat("\nBenchmark test 4:", round(bench.t4*100), "%")
}else{
  cat("\n\nTest 4 failed")
}

if (length(my_test5)){
  cat("\n\nTest 5 passed - using factored status & twice the covariates",
  "\n                where the new covariates differ")
  bench.t5 <- as.numeric(my_t5.b - my_t5.a)/benchmark_time
  cat("\nBenchmark test 5:", round(bench.t5*100), "%")
}else{
  cat("\n\nTest 5 failed")
}

if (length(my_test6)){
  cat("\n\nTest 6 passed - using factored status & twice the covariates",
  "\n                where the new covariates are exactly a copy of first covariates")
  bench.t6 <- as.numeric(my_t6.b - my_t6.a)/benchmark_time
  cat("\nBenchmark test 6:", round(bench.t6*100), "%")
}else{
  cat("\n\nTest 6 failed - exact copies of the covariates fail")
}

if (length(my_test7)){
  cat("\n\nTest 7 passed - using regular status & twice the covariates",
  "\n                where the new covariates differ")
  bench.t7 <- as.numeric(my_t7.b - my_t7.a)/benchmark_time
  cat("\nBenchmark test 7:", round(bench.t7*100), "%")
}else{
  cat("\n\nTest 7 failed")
}

if (length(my_test8)){
  cat("\n\nTest 8 passed - using regular status & thrice the covariates",
  "\n                where the new covariates differ")
  bench.t8 <- as.numeric(my_t8.b - my_t8.a)/benchmark_time
  cat("\nBenchmark test 8:", round(bench.t8*100), "%")
}else{
  cat("\n\nTest 8 failed")
}

#summary(my_test1)
#summary(my_test2)
#summary(my_test4)
#summary(my_test5)
$\endgroup$
1
$\begingroup$

This is not an answer in terms of computational speed, but rather a practical solution based on my experience with other types of regression models: Have you tried with a much smaller sample? I'm guessing that even with n=10,000 you might get a decent model (in terms of precision). I'd also try estimating a few such small models and see how widely they vary to better capture precision.

$\endgroup$
  • 1
    $\begingroup$ Thank you for your answer. I guess that's my fall-back strategy, but right now I find it hard to believe that there is a regression that takes that long to compute. I also find it surprisingly hard to find examples on different competing risk analysis. I've read Melanie Pintilie's book on Competing Risks with a few helpful examples but not really the depth I was hoping for. I also stumbled upon Geskus article in Biometrics 2011 where he proposes weighting coxph() data but since I'm an MD I found the article a little hard to penetrate. There is a striking lack of examples... $\endgroup$ – Max Gordon Jul 17 '11 at 18:40
3
$\begingroup$

Variance estimation with the R crr() function from the cmprsk package accounts for slowing down the run time by perhaps an order of magnitude. If you call crr() and disable the variance calculation as follows:

crr(ftime, fstatus, cov, variance=FALSE)

then the run time will be dramatically reduced, but the returned model gives only the regression coefficients, and not the se(coef) or the p-values. It turns out that if you are doing stepwise competing risks regression, that you do not need these variance estimates for variable selection, only the log likelihoods, and can perform the variance calculation once at the end when you have finalized the variables in your model.

A single crr() run on ~270k observations without variance calculations can take <30 minutes, versus many hours with the default setting of variance=TRUE. We have been implementing a multicore version of the R crrstep package for stepwise competing risks regression, and discovered this variance bottleneck. We haven't been able to do anything about the slowness of the variance estimate, but have been able to run stepwise competing risks regression in a ~1 day run with ~270k observations and ~100 covariates on a 24-core machine, excluding the final variance estimates.

Some time ago I did some speed tests with the Stata competing risks regression and found it ran about an order of magnitude slower than the R cmprsk implementation. SAS also implements competing risks regression in its PROC PHREG, but I have no experience with its performance.

$\endgroup$
1
$\begingroup$

Stata's stcrreg uses the method of Fine and Gray (1999), I don't know if CRR uses that methodology or not. But if it does, there is an alternative methodology which does not estimate subhazard ratios but instead relies on stratifying by failure type and replicating records. For a description see in Rosthoj S. Andersen P.K., Abildstrom S.Z. (2004) which uses SAS macros; or Putter, Fiocco, and Geskus (2007). See also Stata's stcompadj.ado. Interpreting affects is largely done graphically, however. In Stata, stcrreg (Fine and Gray) is much slower than stcompadj (which uses additive hazards).

$\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.