加快手工制作的Cox模型拟合(v.s. `存活：：coxph`‘)

Question

下面，我将R-函数的结果与我自己的代码进行比较。该算法简单地包括最大化多个参数的函数(这里，19)。我的代码定义了函数并使用nlm进行优化。幸运的是，两者都返回相同的结果。然而，R-函数是惊人的快速.因此，我怀疑我可以比使用nlm (或者R中类似的优化例程)做得更好。有什么想法吗？

这是一些生存数据，可以安装一个Cox模型。要做到这一点，我们需要最大化部分日志的可能性(维基百科链接中的第三个方程)。

在R中，这可以通过coxph() ( survival包的一部分)来完成：

> library(survival)
> fmla <- as.formula(paste("Surv(time, event) ~ ", 
+                          paste(names(data)[-(1:3)], collapse=" +")))
> mod <- coxph(formula=fmla, data=data)
> round(mod$coef, 3)
    x1     x2     x3     x4     x5     x6     x7     x8     x9    x10    x11    x12    x13    x14    x15 
-0.246 -0.760  0.089 -0.033 -0.138 -0.051 -0.484 -0.537 -0.620 -0.446 -0.204 -0.112 -0.089 -0.451  0.043 
   x16    x17    x18    x19 
 0.106 -0.015 -0.245 -0.653

这可以通过显式写入部分日志可能性和使用一些数值优化例程来检查。这是一些粗略的代码来完成这项工作。

代码是根据我收到的评论进行编辑的

> #------ minus partial log-lik ------
> Mpll <- function(beta, data) 
+   #!must be ordered by increasing time! 
+   #--> data <- data[order(data$time), ]
+ {
+   #preparation
+   N <- nrow(data)  
+   linpred <- as.matrix(data[, -(1:3)]) %*% beta
+   
+   #pll
+   pll <- sum(sapply(X=which(data$event == 1), FUN=function(j) 
+     linpred[j] - log(sum(exp(linpred[j:N])))))
+   
+   #output
+   return(- pll)
+ }
> #-----------------------------------
> 
> data <- data[order(data$time), ]
> round(nlm(f=Mpll, p=rep(0, 19), data=data)$estimate, 3)
 [1] -0.246 -0.760  0.089 -0.033 -0.138 -0.051 -0.484 -0.537 -0.620 -0.446 -0.204 -0.112 -0.089 -0.451
[15]  0.043  0.106 -0.015 -0.245 -0.653

好吧，这很管用..。但是它要慢得多！

有没有人知道在coxph()中做了什么才能让它变得如此之快？

Answer 1

下面是代码的向量化版本。

Mpll2 <- function(beta, data) {
  X <- as.matrix(data[, -(1:3)])
  a <- X %*% beta
  b <- log(rev(cumsum(rev(exp(a)))))
  -sum((a - b)[data$event==1])
}

下面是运行时间的简单测试。

data <- data[order(data$time), ] # No reason to order every time

# Yours
system.time(round(nlm(f=Mpll, p=rep(0, 19), data=data)$estimate, 3))
#    user  system elapsed 
#    2.77    0.01    2.79 

# Vectorized
system.time(round(nlm(f=Mpll2, p=rep(0, 19), data=data)$estimate, 3))
#    user  system elapsed 
#    0.28    0.00    0.28 

# Optimized C code
fmla <- as.formula(paste("Surv(time, event) ~ ", 
                          paste(names(data)[-(1:3)], collapse=" +")))
system.time(round(coxph(formula=fmla, data=data)$coef,3)) 
#    user  system elapsed 
#    0.02    0.00    0.03 

所以，每种类型之间大约有一个数量级的差异。C是非常快的，你永远不会接近R中的速度，但是C很难写。