3

以下では、R 関数の結果を自分のコードと比較します。アルゴリズムは単純に、多くのパラメーター (ここでは 19) の関数を最大化することで構成されています。私のコードは関数を定義し、nlm最適化に使用します。幸いなことに、どちらも同じ結果を返します。しかし、R 関数は驚くほど高速です。nlmしたがって、 （またはRの同様の最適化ルーチン）を使用するよりもうまくいくと思います。何か案が？

Cox モデルに当てはめることができるいくつかの生存データを次に示します。そのためには、部分対数尤度を最大化する必要があります (ウィキペディア リンクの 3 番目の方程式)。

では、これは (パッケージの一部) で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

OK、動作します...しかし、はるかに遅いです!

内部で何が行われているの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

したがって、各タイプ間で約1桁の違いがあります。Cは非常に高速であり、Rでこれらの速度に近づくことは決してありません。しかし、Cは書くのが難しいです。

于 2013-01-04T18:55:33.933 に答える