3

下面，我将 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) 
+   #!!!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()才能让它如此之快？

4

1 回答 1

6

这是您的代码的矢量化版本。

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 更难编写。

于 2013-01-04T18:55:33.933 回答