r - 加速手工制作的 Cox 模型拟合(v.s. `survival::coxph`)

Question

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

---

这是一些 survival data可以安装 Cox model .为此，需要最大化部分对数似然(维基百科链接中的第三个等式)。

在 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) 
+   #!!!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 更难编写。