R の survival package には、条件付きロジスティック回帰分析を行うための clogit() 関数なるものがあります。
最近この関数を使っていて、ちょっと気づいた点がありましたので、メモを残しておこうと思います。
method="approximate" は危ない
clogit() 関数はデフォルトでは method="exact" というオプションが設定されていて、この方法は局外パラメータの十分統計量を与えた条件付き尤度関数に基づいて、パラメータを推定してくれます。
しかし、この方法はうまくコーディングしないと計算量がとても大変なことになります。
マニュアルには、
The default is to use the exact conditional likelihood, a commonly used approximate conditional likelihood is provided for compatibility with older software.
とは書いてありますが、ぜんぜん答えが返ってこないし method="approximate" でいいかなぁ〜、という人もいないとは限りません (もちろん単に互換性のためにオプションが残っているということは認識しておくべきでしょう)。
そんなノリで近似してしまうと、なんだかよく分からない条件付き尤度関数が使われて、全然違う結果が返ってきてしまいます。
内部で何が行われているのか分からなかったので、Cox の比例ハザードモデルの結果 (条件付きロジスティック回帰分析と密接な関係があるので) と比較してみました。
path <- "C:/"; library(survival) data1 <- read.table(paste(path, "data1.csv", sep=""), header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE) data1$strata <- factor(data1$strata) data1$sy <- Surv(rep(1, length(data1$y)), data1$y)
clogit() 関数では、data1$sy <- Surv(rep(1, length(data1$y)), data1$y) と扱って、つまり全てを同じ時点とし、Case をイベント、 Control を打ち切りとして扱います。前回までは、ここが間違っていたらしい。
手法としては method="breslow" で計算しているみたいです (追記: 2010/05/05 以下も大幅に修正 & 追記)。
> # 正しい条件付き尤度関数 > clogit(y ~ x + strata(strata), data = data1) coef exp(coef) se(coef) z p x 3.09 22.0 1.15 2.68 0.0073 > # 近似 > clogit(y ~ x + strata(strata), data = data1, method = "approximate") coef exp(coef) se(coef) z p x 0.507 1.66 0.342 1.48 0.14 > # 無条件の尤度関数 > glm(formula = y ~ x + strata, family = binomial(), data = data1) (Intercept) x strata2 strata3 strata4 strata5 0.1192 3.5232 -1.3823 1.3451 -2.0330 0.4448 > # Cox の比例ハザードモデル (method = "efron") > coxph(sy ~ x + strata(strata), data = data1) coef exp(coef) se(coef) z p x 0.905 2.47 0.348 2.6 0.0094 > # Cox の比例ハザードモデル (method = "breslow") > coxph(sy ~ x + strata(strata), data = data1, method = "breslow") coef exp(coef) se(coef) z p x 0.507 1.66 0.342 1.48 0.14 > # Cox の比例ハザードモデル (method = "exact") > coxph(sy ~ x + strata(strata), data = data1, method = "exact") coef exp(coef) se(coef) z p x 3.09 22.0 1.15 2.68 0.0073
method="exact" が猛烈に遅い
これは何故なのかがよく分かりませんでした、Cox の比例ハザードモデル用に作られているためなのか (tphreg の結果を見ると、同じじゃまずい？中身をきちんと見ておかないといけないですね)、アルゴリズムがあまり良くないのか・・・、ソースまで見るのはちょっと止めておこうと思います。
ちなみに SAS で使われている方法 [1] を実装してみたところ、層が 6 つで全体の N = 975、層のサンプルサイズが最大で n = 242 のデータ [2] も 14 秒程度で計算できました (R で組んだので、コア部分を C で組むともっと早いはず)。
↓
C で dll 書いてみました。
$sampling.time [1] 14.34
↓
$sampling.time [1] 0.42
(・∀・)
SASの結果と比較してみよう
SAS の場合は logistic procedure の exact statement 以外はかなり高速です、n が数百ぐらいなら一瞬で計算が終わります。
exact statement では conditional exact な信頼区間の計算の部分で、局外パラメータの十分統計量の値を固定した標本空間上で確率変数の値を色々動かして条件付き分布をたくさん生成するので時間がかかってしまいます、条件付き尤度関数に基づいて点推定するだけならそれほど時間がかかるわけではないです。後は、Complete separation や Quasi-complete separation が起こったときに、median unbiased estimate を求めてくれます。
比較した結果を見てみます。
beta (x) の推定値 真値は 2 に設定したが、乱数を使ったので目安程度に SAS: Logistic-Conditional 3.0929 Logistic-Unconditional 3.5232 Logistic-Exact 3.0929 Cox-Conditional-Discrete 3.0929 Cox-Conditional-Breslow 0.5069 Cox-Conditional-Efron 0.9049 Cox-Conditional-Exact 1.5368 R: Logistic-Conditional 3.09 Logistic-Conditional-Approx 0.507 Logistic-Unconditional 3.5232 Cox-Conditional-Breslow 0.507 Cox-Conditional-Efron 0.905 Cox-Conditional-Exact 3.09
気づいた点
使ったデータ (SAS)
data Data1;
call streaminit(975861);
beta = 2;
do rep = 1 to 1;
do strata = 1 to 5;
do x = 0 to 1;
alpha = strata/10;
logit = alpha + beta * x;
theta = 1 / (1 + exp(-logit));
do n = 1 to floor(rand('Uniform')*10);
y = rand('Bernoulli', theta);
/* censored: censor = 1, not-censored: censor = 0 */
if y = 1 then do; sy = 1; censor = 0; end;
else do; sy = 1; censor = 1; end;
output;
end; end; end; end;
proc phreg data = Data1;
class x / param = ref ref = first;
strata strata;
model sy * censor(1) = x / ties = breslow;
run;
参考文献
[1] Gail MH, Lubin JH, Rubinstein LV. Likelihood calculations for matched case-control studies and survival studies with tied death times. Biometrika 1981; 68(3): 703–707.
[2] Breslow NE, Day NE. Statistical methods in cancer research. Volume I - The analysis of case-control studies. IARC Scientific Publications, 1980: p.151.
