条件付きロジスティック回帰でクラスターデータを解析する（clogistic関数, clogit関数）[R]

条件付きロジスティック回帰（CLR）
通常のロジスティック回帰を行う
Epiパッケージのclogistic( )を使ってCLR
survivalパッケージのclogit( )を使ってCLR
おわりに
参考資料

対象集団をいくつかの異なったグループに分けられるデータをクラスターデータ（clustered data）と呼ぶ。

例えば、

学校 - 生徒
個体 - 複数の患部（眼など）
個体 - 毎月の測定（反復測定）

などが挙げられる。

クラスター内の類似性を考慮しないとデータが持っている情報量を過大評価することになるほか、クラスターの異質性を考慮することで検出力を向上させうる。

同一対象者の前後比較やマッチングは、同一対象者やマッチさせたペア同士の類似性を差し引くことで検出力を向上させることが目的であり、この点でクラスターデータと同じと考えられる。

条件付きロジスティック回帰（CLR）

2値アウトカムの解析方法としてロジスティック回帰が行われることが多いが、クラスターを考慮するために取られる方法としては以下のいずれかの方法が用いられることが多い。

一般化線形混合効果モデル（generalized linear mixed model, GLMM）
一般化推定方程式（generalized estimating equations, GEE）
条件付きロジスティック回帰（conditional logistic regression, CLR）

（他にalternating logistic regressionというものがあるみたい。）

条件付きロジスティック回帰は別名fixed effect logistic regressionと呼ばれる。通常のロジスティック回帰（standard logistic regression, SLR）の推定ではunconditional likelihoodを用いた最尤法が使われるが、CLRではconditional likelihoodが使われる。
重要な違いは、conditional likelihoodではモデルの切片を推定しないこと。各クラスター（あるいはマッチングの各ペア）の固定効果を推定しなくていいので、その分パラメータ数が減ってお得。
しかし、切片が推定されないので予測確率が算出できない。

Conditional likelihoodを使うべきときにunconditional likelihoodを使ってしまうと回帰係数βの推定にバイアスを生じるが、unconditional likelihoodを使うべき時にconditional likelihoodを使ってもバイアスは生じない。よって、迷ったらconditional likelihoodを使う解析をしていいとのこと。

ここではサンプルデータとして、Epiパッケージのbdendoデータセットを使う。
bdendoデータは子宮内膜癌患者に対して4人のmatched controlを設定したケースコントロール研究のもの。

> library(Epi)
> data(bdendo)
> head(bdendo,3)
  set d gall hyp   ob est dur non duration age cest agegrp  age3
1   1 1   No  No  Yes Yes   4 Yes       96  74    3  70-74 65-74
2   1 0   No  No <NA>  No   0  No        0  75    0  70-74 65-74
3   1 0   No  No <NA>  No   0  No        0  74    0  70-74 65-74

ここで使う変数は以下の3つ。

d：1=case, 0=control
cest：結合型エストロゲン用量（0,1,2,3のカテゴリー変数）
dur：結合型エストロゲン治療の期間（0,1,2,3,4のカテゴリー変数）
set：マッチさせたグループID（case1人とcontrol4人から成る）

通常のロジスティック回帰を行う

マッチングを無視してSLRを行ってみる。

> fit <- glm(d~cest+dur, data=bdendo, family=binomial(link="logit"))
> library(broom)
> tidy(fit, exponentiate = TRUE, conf.int = TRUE)
# A tibble: 8 x 7
  term        estimate std.error statistic   p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
1 (Intercept)   0.0839     0.301    -8.24   1.65e-16   0.0441     0.145
2 cest1         5.21       0.526     3.14   1.70e- 3   1.83      14.6  
3 cest2         5.24       0.550     3.01   2.61e- 3   1.76      15.4  
4 cest3        15.7        0.532     5.17   2.30e- 7   5.60      45.7  
5 dur1          0.358      0.579    -1.77   7.60e- 2   0.108      1.08 
6 dur2          0.739      0.496    -0.611  5.41e- 1   0.277      1.96 
7 dur3          0.914      0.556    -0.162  8.71e- 1   0.302      2.71 
8 dur4         NA         NA        NA     NA         NA         NA

Interceptが推定されている。
（durのカテゴリー4が推定不可なのは何故だろう...）

Epiパッケージのclogistic( )を使ってCLR

クラスター/マッチングIDをstrata=で指定する。

fit <- clogistic(d ~ cest+dur, strata=set, data=bdendo)

> fit
Call: 
clogistic(formula = d ~ cest + dur, strata = set, data = bdendo)

        coef exp(coef) se(coef)      z       p
cest1  1.890     6.621    0.580  3.258 1.1e-03
cest2  1.846     6.333    0.596  3.095 2.0e-03
cest3  2.944    18.991    0.609  4.835 1.3e-06
dur1  -1.247     0.287    0.620 -2.011 4.4e-02
dur2  -0.521     0.594    0.542 -0.961 3.4e-01
dur3  -0.177     0.838    0.565 -0.313 7.5e-01
dur4      NA        NA    0.000     NA      NA

Likelihood ratio test=35.3  on 6 df, p=3.8e-06, n=254

今回は切片が推定されていない。

survivalパッケージのclogit( )を使ってCLR

こちらはsummary( )やbroom( )で整形して表示可。

> library(survival)
> fit <- clogit(d~cest+dur+strata(set), data=bdendo, method="exact")
> summary(fit)
Call:
coxph(formula = Surv(rep(1, 315L), d) ~ cest + dur + strata(set), 
    data = bdendo, method = "exact")

  n= 295, number of events= 54 
   (20 observations deleted due to missingness)

         coef exp(coef) se(coef)      z Pr(>|z|)    
cest1  1.8902    6.6207   0.5801  3.258  0.00112 ** 
cest2  1.8458    6.3334   0.5964  3.095  0.00197 ** 
cest3  2.9439   18.9906   0.6089  4.835 1.33e-06 ***
dur1  -1.2471    0.2873   0.6201 -2.011  0.04431 *  
dur2  -0.5206    0.5942   0.5419 -0.961  0.33674    
dur3  -0.1767    0.8381   0.5651 -0.313  0.75457    
dur4       NA        NA   0.0000     NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

      exp(coef) exp(-coef) lower .95 upper .95
cest1    6.6207    0.15104   2.12387   20.6384
cest2    6.3334    0.15789   1.96796   20.3828
cest3   18.9906    0.05266   5.75731   62.6411
dur1     0.2873    3.48022   0.08523    0.9688
dur2     0.5942    1.68298   0.20542    1.7186
dur3     0.8381    1.19323   0.27685    2.5369
dur4         NA         NA        NA        NA

Concordance= 0.72  (se = 0.051 )
Likelihood ratio test= 35.28  on 6 df,   p=4e-06
Wald test            = 27.73  on 6 df,   p=1e-04
Score (logrank) test = 36.83  on 6 df,   p=2e-06