以下では、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()
か、非常に高速にするためのアイデアを持っている人はいますか?