Triad sou.

survivalパッケージのtt関数

はじめに

ちょっと色々とやっていて R の survival パッケージの tt 関数の挙動が気になるという話になりました。
survival パッケージでは時間変動共変量の解析が実施可能であり、時間変動共変量と時間変動効果(回帰係数)についての vignette には詳細な解説が書かれています。
vignette 中のコードには色々な事例が書かれているのですが、tt 関数を使っているものがあります。
tt 関数は、時間変換関数 を組み込むための機能です。
なんでこんなものが必要になるかというと、モデル中に共変量と時間の交互作用項(例えば、$\beta(t)x=\beta_x x+\beta_{tx}\log(t)x$ のようなもの)が含まれる場合には、部分尤度関数を計算する際にリスク集合毎に分母も含めてイチから全計算をする必要があり、単純な Cox 回帰モデルの場合と計算方法が変わるためです(おまけでも触れていますが、時間交互作用項がある場合でもデータをイベント時間毎に層を持つ形に再構成して、層化統合解析をするのと同じ計算になります)。
モデル中に共変量と時間の交互作用項がないときは、分母の計算にはリスク集合に含まれる部分までの和 $\sum_{j \in R(t_i)} \exp(\mathbf{X}_j^T \boldsymbol{\beta})$ が必要なのでソートして累積和を1セット計算するだけでよいのですが、交互作用項がある場合は $\sum_{j \in R(t_i)} \exp\{\mathbf{X}_j(t_i)^T \boldsymbol{\beta}\}$ になって毎回全計算が必要です(これを自動化するために tt 関数が用意されているということかと思います)。
この違いは、自力で Cox 回帰の部分尤度関数の計算を実装してみるとよく分かるので、やってみると良いでしょう(簡単です)。

この tt 関数は、時間関数形を

coxph(Surv(time, status) ~ ph.ecog + tt(age), data = lung,
     tt = function(x, t, ...) pspline(x + t/365.25))

こんな感じに定義して使うのですが、関数の定義を書かずに

coxph(Surv(time, status) ~ ph.ecog + tt(age), data = lung)

のようにしても結果が出てきます。
このように定義を書かない場合には一体何が出てきているのか、という話題になり、これを確認した結果をまとめてみようかなと思いました。

tt 関数のデフォルトの挙動

これについてはヘルプファイルには書いていないです(確認した限り、そのような記載は見つけられませんでした)。
正解は vignette をよく読むと書いてありました。

There are other interesting uses for the time-transform capability. One example is O'Brien's logit-rank test procedure [6]. He proposed replacing the covariate at each event time with a logit transform of its ranks. This removes the influence of any outliers in the predictor x. For this case we ignore the event time argument and concentrate on the groupings to create the following tt function.

> function(x, t, riskset, weights){
obrien <- function(x) {
r <- rank(x)
(r-.5)/(.5+length(r)-r)
}
unlist(tapply(x, riskset, obrien))
}

This relies on the fact that the input arguments to tt() are ordered by the event number or risk set. This function is used as a default if no tt argument is present in the coxph call, but there are tt terms in the model formula. (Doing so allowed me to depreciate the survobrien function).

ざっくりいうと「tt 関数は時間関数形の指定とは別の面白い使い方ができる、例えば、O'Brien's logit-rank testの実装にも使えて、これは coxphtt 引数に何も指定しない場合の tt 関数のデフォルトになっている(これにより、私は survobrien 関数を非推奨にできた)」という感じです。
どうやら、tt 引数の指定をしないで実行すると、リスク集合毎に共変量の順位を変換した変数を生成するため、通常よくやられているような時間関数形を指定した解析とは目的が違う結果が出てくるということが分かります(Therneau 先生...、これはちょっとやめてほしいかも...)。
これで計算される回帰係数は共変量の順位に基づいた結果になっていて時間は無視されるので、時間関数形については評価できていません!

tt (time-transform) 関数のデフォルトが、全然 time-transform と関係がない仕様というのは若干思うところがあります。

通常の tt 関数の使い方

これも念のため紹介しておきます。
一番よくある使い方は、比例ハザード性の検定だったり、比較的簡単な構造を持つ時間変動効果のモデリングではないかと思います。
例えば、対数時間のスケールにおける時間変動効果 $\beta(t)x=\beta_x x+\beta_{tx}\log(t)x$ を持つモデルは tt 関数を使って当てはめることが可能です。
また、これと同じモデルにおいて、時間変動する部分の回帰係数のスコア検定を実施できる cox.zph という関数もあります。

fit1 <- coxph(Surv(time, status) ~ karno + tt(karno), data = veteran, tt = function(x, t, ...) x*log(t))
fit2 <- coxph(Surv(time, status) ~ karno, data = veteran)
summary(fit1)
cox.zph(fit2, transform = "log")

両者の実行結果は、

> summary(fit1)
Call:
coxph(formula = Surv(time, status) ~ karno + tt(karno), data = veteran, 
    tt = function(x, t, ...) x * log(t))

  n= 137, number of events= 128 

               coef exp(coef)  se(coef)      z Pr(>|z|)    
karno     -0.083723  0.919686  0.016783 -4.989 6.08e-07 ***
tt(karno)  0.013408  1.013498  0.004196  3.195   0.0014 ** 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

          exp(coef) exp(-coef) lower .95 upper .95
karno        0.9197     1.0873    0.8899    0.9504
tt(karno)    1.0135     0.9867    1.0052    1.0219

Concordance= 0.709  (se = 0.024 )
Likelihood ratio test= 52.99  on 2 df,   p=3e-12
Wald test            = 49.81  on 2 df,   p=2e-11
Score (logrank) test = 56.97  on 2 df,   p=4e-13

> cox.zph(fit2, transform = "log")
       chisq df      p
karno   10.5  1 0.0012
GLOBAL  10.5  1 0.0012

となりました。
tt 関数を使った検定は Wald 検定($P = 0.0014$)で、zph 検定はスコア検定($P = 0.0012$)のため若干値は異なります。

参考文献

  • Therneau T, Crowson C, Atkinson E. Using time dependent covariates and time dependent coefficients in the Cox model. [Link]
  • O'Brien P. A nonparametric test for association with censored data. Biometrics 1978;34:243–250. DOI: 10.2307/2530014.

おまけ

一応 tt 関数の中身を見てみて、何を計算しているのかを確認しました。

library(survival)

temp <- subset(pbc, id <= 312, select = c(id:sex, stage))
pbc2 <- tmerge(temp, temp, id = id, death = event(time, status))
pbc2 <- tmerge(pbc2, pbcseq, id = id, lbili = tdc(day, log(bili)))
cox_model0 <- coxph(Surv(tstart, tstop, death == 2) ~ lbili, pbc2)
cox_model1 <- coxph(Surv(tstart, tstop, death == 2) ~ tt(lbili), pbc2)

tt <- function(x, riskset) {
  obrien <- function(x) {
    r <- rank(x)
    (r - 0.5)/(0.5 + length(r) - r)
  }
  unlist(tapply(x, riskset, obrien))
}
Y <- cbind(pbc2$tstart, pbc2$tstop, as.numeric(pbc2$death == 2))
sort.end <- order(-Y[, 2], Y[, 3])
sort.start <- order(-Y[, 1])
newstrat <- c(1L, rep(0, nrow(Y) - 1))
counts <- .Call(survival:::Ccoxcount2, Y, as.integer(sort.start - 1L), as.integer(sort.end - 1L), as.integer(newstrat))
tindex <- counts$index
istrat <- rep(1:length(counts$nrisk), counts$nrisk)
Y <- Surv(rep(counts$time, counts$nrisk), counts$status)
tv <- pbc2$lbili[tindex]
ttv <- tt(tv, istrat)
pbc3 <- data.frame(time = Y[, 1], status = counts$status, ttv, istrat)
cox_modelttv <- coxph(Surv(time, status) ~ ttv + strata(istrat), data = pbc3)

cox_model1
cox_modelttv

これが tt 関数をそのまま実行した結果で、

> cox_model1
Call:
coxph(formula = Surv(tstart, tstop, death == 2) ~ tt(lbili), 
    data = pbc2)

               coef exp(coef)  se(coef)    z      p
tt(lbili) 0.0073115 1.0073383 0.0005712 12.8 <2e-16

Likelihood ratio test=74.29  on 1 df, p=< 2.2e-16
n= 1807, number of events= 125 

こっちは tt 関数が作っているデータを確認するために検算した方です。

> cox_modelttv
Call:
coxph(formula = Surv(time, status) ~ ttv + strata(istrat), data = pbc3)

         coef exp(coef)  se(coef)    z      p
ttv 0.0073115 1.0073383 0.0005712 12.8 <2e-16

Likelihood ratio test=74.29  on 1 df, p=< 2.2e-16
n= 24422, number of events= 125

データ数以外は一致していることが分かります。
検算用のデータ(pbc3)を見ると、リスク集合毎に層が作られていることが確認できます。

> pbc3
    time status         ttv istrat
1   4191      0  1.36363636      1
2   4191      0  0.36842105      1
3   4191      0  0.23809524      1
4   4191      0  2.71428571      1
5   4191      0  0.08333333      1
6   4191      0  0.62500000      1
7   4191      0  0.62500000      1
8   4191      0  1.88888889      1
9   4191      0  4.20000000      1
10  4191      0  1.00000000      1
11  4191      0  0.08333333      1
12  4191      0 25.00000000      1
13  4191      1  7.66666667      1
14  4079      0  0.78947368      2
15  4079      0  0.13333333      2
16  4079      0  0.54545455      2
17  4079      0  3.85714286      2
18  4079      0  0.30769231      2
19  4079      0  1.00000000      2
20  4079      0  1.61538462      2
21  4079      0  2.09090909      2
22  4079      0  5.80000000      2
23  4079      0  1.26666667      2
24  4079      0  0.30769231      2
25  4079      0 10.33333333      2
26  4079      0  0.54545455      2
27  4079      0  0.03030303      2
28  4079      0  0.13333333      2
29  4079      0  2.77777778      2
30  4079      1 33.00000000      2
31  3853      0  0.61290323      3

『宇宙怪人しまりす統計よりも重要なことを学ぶ』を読みました

楽しく読ませていただきました。

いくつか印象に残った点についての感想です。

第1話(としまりすコラム1)、報告時に探索的研究と検証的研究を区別しないといけない、という話題が出てきていました。やっぱり、これは大事だなと思いました。完全に切り分けるのも難しいのかもしれませんが、論文に読者の判断の助けになるような記号をつけられるような仕組みがほしいと思いました(P値を使うことを禁止しないんだったらこれぐらいはやってもいいのでは)。事前に主要評価項目が定義されている、サンプルサイズ設計がちゃんとされている、解析計画が事前に公表されている等。某系列のようにsupplementaryにチェックリストへの報告義務を課すところもありますが、これはどのぐらい読まれているんだろうといつも思います。

第2話、裁判記録のお話は聞いたことがなかったので個人的に衝撃を受けました、みんな読んだ方が良いなと思いました。誠実であり続けることを誓います。

第3話、合流の話も出てきて非常に参考になりました。

第4話、層別解析とサブグループ解析の話はあるあるですね。自分も授業では触れるようにしています。確かに層化統合解析とか別の呼び方にしたほうが間違われにくそうです。ターゲット集団の話が出てきて、授業で話すんだったらどこまで話すのがいいんだろう、とふと思いました。

第5話、審査会でしまりすくんと先生はさらっと正しいコメントしていますね(カッコいい)。現状では同様の指摘に該当する研究がどのぐらいあるんだろうかと思います。このようなコメントが来ないことが常識になるような世界を整備できるようにしないといけないな、と思いました。

自分が関わる医学研究の範囲に限っての感想ですが、因果推論を導入するメリットは因果関係を示すことができるからではなく、研究計画の段階で効果の定義をめちゃくちゃ意識できるところかな、とより強く感じました。

追記

COI開示です。

割合の差のスコア検定

割合の差のスコア検定

以下のツイートをしましたところ、情報提供をいただきましたので、まとめた文章を作成しました。


非劣性仮説の割合の差のWald検定でマージンを0にした時の結果と、分割表のカイ二乗検定(スコア検定に一致)の結果がずれるという状況だったので、それならば非劣性仮説の割合の差のスコア検定でマージンを0にした時の結果を出せば、分割表のカイ二乗検定と一致すると思い、Farrington–Manning検定の結果が同じになることを確認したという内容です。


Farrington–Manning検定 (Farrington & Manning, 1990) は、非劣性仮説の割合の差のスコア検定と(割合の比の方もですが)として知られていて、SASやRのパッケージではこのスコア検定はFarrington–Manning検定と称されています。
しかし、実際にはこれと同値な割合の差のスコア検定統計量がMee (1984) で最初に与えられており、佐藤(1988)ではMeeの方法がスコア検定の式に一致することが示されています。
Farrington & Manning(1990)は非劣性仮説の検定とサンプルサイズ設計を主題にしたもので、なんらかの理由でこっちの方が有名になり(Mee(1984)はレターだったことも関連しているでしょうか?)、Farrington–Manning検定という呼び方が定着してしまったのかと思います。


Meeの方法と佐藤俊哉先生の方法の詳細については、佐藤(1988)で導出も含めて非常に丁寧に書かれていますので、そちらをご参照ください。

割合の差のスコア検定統計量

佐藤(1988)の(3)式
\[
\frac{x - n(\delta + \hat{p}_2)}{(\delta + \hat{p}_2)(1 - \delta - \hat{p}_2)} + \frac{y - m \hat{p}_2}{\hat{p}_2(1 - \hat{p}_2)} = 0
\] の下の
\[
T(\delta)=\left\{ \frac{x - n \hat{p}_1}{\hat{p}_1(1 - \hat{p}_1)} \right\}^2
\left\{ \frac{\hat{p}_1(1 - \hat{p}_1)}{n} + \frac{\hat{p}_2(1 - \hat{p}_2)}{m} \right\}
\] がスコア検定統計量で、$\hat{p}_1 = \delta + \hat{p}_2$ となっています。
(3)式は3次方程式となるので $\hat{p}_2$ を解けば $P$ 値が簡単に計算できます(という話も教えていただきました)。

ということなので、(3)式を整理してみると
\[
a \hat{p}_2^3 + b \hat{p}_2^2 + c \hat{p}_2 + d = 0
\] となり、
\[
\begin{aligned}
a &= 1 + \theta, \\
b &= -\{1 + \theta + \theta\tilde{p}_1 + \tilde{p}_2 - \delta(2 + \theta) \}, \\
c &= \delta^2 - \delta(\theta + 2\tilde{p}_2 + 1) + \theta\tilde{p}_1 + \tilde{p}_2, \\
d &= \delta \tilde{p}_2 (1 - \delta),
\end{aligned}
\] となります、ただし、$\theta=n / {m}$, $\tilde{p}_1=x/n$, $\tilde{p}_2 = y / {m}$ です。
これを解くと
\[
\begin{aligned}
p &= \frac{b^2}{3a^2} - \frac{c}{a}, \\
q &= \frac{2b^3}{27a^3} - \frac{bc}{3a^2} + \frac{d}{a}, \\
\hat{p}_2 &= 2\sqrt{\frac{p}{3}} \cos\left\{ \frac{1}{3}\cos^{-1}\left( -\frac{3q}{2p}\sqrt{\frac{3}{p}} \right)+\frac{4\pi}{3} \right\}-\frac{b}{3a},
\end{aligned}
\] となります(導出は飛ばしていますが、3次方程式の解の公式を使って普通に出しています)。
あとは、$\hat{p}_2$ を $T(\delta)$ に代入し、$T(\delta)$ が帰無仮説のもとで漸近的に自由度1の $\chi^2$ 分布に従うことを利用して $P$ 値を得ることができます。


なお、しまりすくんの本で使用しているのは計算が簡単なWald検定です(あそこの説明で上の式を使うのはちょっと厳しいですね)。
Wald統計量は
\[
T_W(\delta)=\frac{(\tilde{p}_1-\tilde{p}_2-\delta)^2}{\frac{\tilde{p}_1(1-\tilde{p}_1)}{n}+\frac{\tilde{p}_2(1-\tilde{p}_2)}{m}}
\] となります。



数値例

ツイートの数値例を使って検算してみます。
仮説$H: p_1-p_2= 0$を検定すると、

x <- 58
n <- 58 + 22
y <- 62
m <- 62 + 38
d <- 0

a3 <- 1 + n/m
a2 <- -1*(1 + n/m + x/m + y/m - d*(2 + n/m))
a1 <- d^2 - d*(n/m + 2*y/m + 1) + x/m + y/m
a0 <- d*y/m*(1 - d)
p <- a2^2/(3*a3^2) - a1/a3
q <- 2*a2^3/(27*a3^3) - a2*a1/(3*a3^2) + a0/a3
p2h <- 2*sqrt(p/3)*cos(1/3*acos(-3*q/(2*p)*sqrt(3/p)) + 4*pi/3) - a2/(3*a3)
p1h <- d + p2h
Td <- ((x - n*p1h)/(p1h*(1 - p1h)))^2*(p1h*(1 - p1h)/n + p2h*(1 - p2h)/m)
1 - pchisq(Td, 1)

実行結果は

> 1 - pchisq(Td, 1)
[1] 0.1375639

でした。

Farrington–Manning検定の結果も出して確認してみます。
Farrington–Manning検定はMeeタイプの統計量を二乗しないものを使っていて非劣性の検定として実装されているので、片側検定の結果が出てきます。
そのため、$P$ 値を2倍した値で確認する必要があります。

x <- c(rep(TRUE, 58), rep(FALSE, 22))
y <- c(rep(TRUE, 62), rep(FALSE, 38))
fm <- DescrTab2::farrington.manning(x, y, 0)
fm$p.value*2

実行すると

> fm$p.value*2
[1] 0.1375639

となってめでたく一致しました。

$\chi^2$ 検定の方の $P$ 値も確認してみます(載せ忘れていたので追記しました [2024/3/1 10:07])。

X <- matrix(c(58, 62, 22, 38), ncol = 2)
s <- c(colSums(X)[1]/sum(X), colSums(X)[2]/sum(X))
E <- rbind(rowSums(X)[1]*s, rowSums(X)[2]*s)
1 - pchisq(sum((X-E)^2/E), 1)

実行すると

> 1 - pchisq(sum((X-E)^2/E), 1)
[1] 0.1375639

こちらも一致します。


各方法と仮説で計算した $P$ 値がどうなるかをまとめると以下になります。

検定 $\delta=0$ $\delta=0.2$
Wald型 0.132 0.172
スコア型 0.138 0.164

Meeの方法とスコア統計量

佐藤(1988)ではMeeの方法がスコア検定の式に一致することが示されていますが、導出が省略されているので検算した結果をメモしておきます。
Meeの統計量は
\[
T_M(\delta)=\frac{(\tilde{p}_1-\tilde{p}_2-\delta)^2}{\frac{\hat{p}_1(1-\hat{p}_1)}{n}+\frac{\hat{p}_2(1-\hat{p}_2)}{m}}
\] で、これが $T(\delta)$ に等しいことを示します。
$T_M(\delta)$ の分子は、
\[
(\tilde{p}_1-\tilde{p}_2-\delta)^2=
\left( \frac{x}{n} - \hat{p}_1 + \hat{p}_1 - \frac{y}{m} - \delta \right)^2=
\left\{ \left(\frac{x}{n} - \hat{p}_1\right) - \left( \frac{y}{m} - \hat{p}_2 \right) \right\}^2=
\left(\frac{x}{n} - \hat{p}_1\right)^2 -2\left(\frac{x}{n} - \hat{p}_1\right) \left(\frac{y}{m} - \hat{p}_2\right) + \left( \frac{y}{m} - \hat{p}_2 \right)^2
\] と変形できます。
また、(3)式から得られる
\[
\frac{y}{m}-\hat{p}_2=-\frac{\hat{p}_2(1-\hat{p}_2)}{m}\frac{x-n\hat{p}_1}{\hat{p}_1(1-\hat{p}_1)}
\] という関係式を用いれば、$T_M(\delta)$ の分子は、
\[
\left\{ \frac{\hat{p}_1(1-\hat{p}_1)}{n} \right\}^2 \left\{ \frac{x-n\hat{p}_1}{\hat{p}_1(1-\hat{p}_1)} \right\}^2+
2\frac{\hat{p}_1(1-\hat{p}_1)}{n} \frac{\hat{p}_2(1-\hat{p}_2)}{m} \left\{ \frac{x-n\hat{p}_1}{\hat{p}_1(1-\hat{p}_1)} \right\}^2+
\left\{ \frac{\hat{p}_2(1-\hat{p}_2)}{m} \right\}^2 \left\{ \frac{x-n\hat{p}_1}{\hat{p}_1(1-\hat{p}_1)} \right\}^2 =
\left\{ \frac{x-n\hat{p}_1}{\hat{p}_1(1-\hat{p}_1)} \right\}^2 \left\{ \frac{\hat{p}_1(1-\hat{p}_1)}{n} + \frac{\hat{p}_2(1-\hat{p}_2)}{m}\right\}^2
\] となることが分かるので、Meeの統計量とスコア統計量が同じになることが分かります。

参考文献

  • Farrington CP, Manning G. Test statistics and sample size formulae for comparative binomial trials with null hypothesis of non-zero risk difference or non-unity relative risk. Statistics in Medicine 1990; 9(12): 1447–1454. DOI: 10.1002/sim.4780091208.
  • Mee RW. Confidence bounds for the difference between two probabilities. Biometrics 1984; 40(4): 1175–1176. [Link].
  • 佐藤俊哉.コホート研究で用いる疫学指標のefficient scoreにもとづく信頼区間.応用統計学 1988; 17(1): 43–54. DOI: 10.5023/jappstat.17.43.
  • 佐藤俊哉.宇宙怪人しまりす 統計よりも重要なことを学ぶ.朝倉書店 2024.[Link].

おわりに

手抜かりがありました、正確な情報をお伝えできるよう精進致します。

forestplotパッケージのコードメモ

役に立ちそうだったので、メモを作成してみました。

library(forestplot)
library(dplyr)

dat1 <- read.table(textConnection("
l0	l1	l2	l3	l4	l5	hr	lcl	ucl
    Variable 1	100	50	50	0.83 (0.70,  0.99)	0.041	0.833	0.698	0.992
    Variable 2	100	50	50	2.34 (0.58,  9.51)	0.233	2.344	0.577	9.510
    Variable 3	100	50	50	0.44 (0.12,  1.58)	0.207	0.440	0.123	1.577
    Variable 4	100	50	50	1.50 (0.77,  2.90)	0.230	1.497	0.774	2.896
    Variable 5	100	50	50	2.12 (0.54,  8.26)	0.280	2.118	0.543	8.263
    Variable 6	100	50	50	2.03 (1.04,  3.99)	0.039	2.035	1.037	3.992
    Variable 7	100	50	50	0.85 (0.70,  1.02)	0.077	0.846	0.702	1.018
    Variable 8	100	50	50	0.97 (0.94,  1.01)	0.092	0.970	0.935	1.005
    Variable 9	100	50	50	0.43 (0.15,  1.25)	0.123	0.434	0.150	1.254
    Variable 10	100	50	50	1.81 (0.61,  5.41)	0.286	1.813	0.608	5.406
    Variable 11	100	50	50	2.63 (0.66, 10.49)	0.171	2.628	0.659	10.488
"), stringsAsFactors = FALSE, header = TRUE, sep = "\t")
dat1

dat1 %>%
  forestplot(
    mean = hr,
    lower = lcl,
    upper = ucl,
    labeltext = c(l0, l1, l2, l3, l4, l5),
    xlog = TRUE,
    lwd.ci = 3,
    boxsize = 0.25,
    ci.vertices.height = 0.2,
    vertices = TRUE,
    xticks.digits = 1
  ) %>%
  fp_set_style(
    box = "black",
    line = "black",
    summary = "black",
    align = "lccccc",
    hrz_lines = "black",
    txt_gp = fpTxtGp(ticks = gpar(cex = 1))
  ) %>%
  fp_add_header(
    l0 = " Variable",
    l1 = "N",
    l2 = "Label1",
    l3 = "(Label2)",
    l4 = "HR (95%CI)",
    l5 = "P-value"
  ) %>%
  fp_decorate_graph(graph.pos = 5) %>%
  fp_set_zebra_style("#EFEFEF") -> p1

svg("figure.svg", width = 10, height = 3.85)
plot(p1)
dev.off()

感想

forestplotパッケージはかなり高機能で、使いそうだなと思われるオプションが一通り揃っており、非常にカスタムしやすいと思います。
ヘルプに事例があまりついていないため、若干わかりにくい所はありますが、大体雰囲気で使えそうだなと思いました。

少し工夫した点

1列目の変数名とラベルには適宜スペースを入れて調整すると見栄えが良いです。
X軸の文字の大きさを変更するオプション txt_gp = fpTxtGp(ticks = gpar(cex = 1)) が調べてもなかなか分からずにかなり困りました。
個人的に、forest plotは幅をいい感じに調整し、高さをピッタリに詰めると見栄えが良くなると思いますので、そのようにwidthとheightをいじってファイルに吐き出しています。

GEEの注意点の論文(Pepe & Anderson, 1994)を読んだ話

Pepe & Anderson (1994)

先日雑談をしていて、GEEで時間依存性共変量を使うときに注意が必要という話題になり、この論文を思い出して改めて読みなおしてみようかなと思った次第です。

A Key Condition

細かい話は論文を確認すれば良いとして、この論文で書かれているKey Conditionを見てみます(導出はセクション2に丁寧に書いてあります)。


When using an estimating equation of the form (2) to estimate $\beta$ from a marginal model $E(Y_{it} \mid X_{it}) = g(X_{it}\beta)$, either a diagonal working covariance matrix should be employed or one must validate that the marginal expectation $E(Y_{it} \mid X_{it}) = E(Y_{it} \mid X_{ij}, j=1,...,n_{i})$.
推定方程式$E(Y_{it} \mid X_{it}) = g(X_{it}\beta)$を使って$\beta$を推定するなら、対角の作業相関行列を使うか、$E(Y_{it} \mid X_{it}) = E(Y_{it} \mid X_{ij}, j=1,...,n_{i})$が成立することを確認せよということでした。
$X_{it}$が時間について固定されている場合には成立つことが自明ですが、$\mid X_{ij}$が時間変動する場合は微妙なケースが多そうです。

事例

論文ではシミュレーションを実施しており、そのうちの一つが、
\[
Y_{it}=\alpha Y_{i(t-1)}+\beta X_{it} + \epsilon_{it},~
t=1,\ldots,n_i,
\] というモデルでした。
ただし、$Y_{i0}=0$で、$(X_{it}, \epsilon_{it})$が平均$0$で互いに独立かつ$Y_{i(t-1)}$とも独立としていました。
これらの条件を設定しておくと、周辺平均に対しては
\[
E(Y_{it} \mid X_{it})=X_{it}\beta,
\] が成り立ちます。
一方で、条件付き周辺平均に対しては
\[
E(Y_{it} \mid X_{ij}, j=1,...,n_{i})=\beta \sum_{j \leq t} \alpha^{t-j}X_{ij}
\] が成り立ちます。
この例は、$E(Y_{it} \mid X_{it}) \ne E(Y_{it} \mid X_{ij}, j=1,...,n_{i})$となります。

検算

検算しやすいように、$X$は二値にして平均$0$になるように設定しました。

library(dplyr)
library(tidyr)

beta <- 1
alpha <- 1.5
df <- NULL
for(i in 1:1000) {
  x <- rep(NA, 3)
  for(t in 1:3) {
    x[t] <- sample(c(-1, 1), 1)
    if (t == 1) {
      ypre <- 0
    }
    y <- alpha*ypre + beta*x[t] + rnorm(1, 0, 0.1)
    pcm <- 0
    for(k in 1:t) {
      pcm <- pcm + alpha^(t - k)*x[k]
    }
    pcm <- beta*pcm
    ypre <- y
    df <- rbind(df, data.frame(i, t, y, x = x[t], pcm))
  } 
}

df %>%
  group_by(t, x) %>%
  summarize(mm = mean(y), mm_altc = mean(beta*x))

df2 <- df %>%
  pivot_wider(id_cols = i, names_from = t, values_from = x,
              names_prefix = "x") %>%
  mutate(patrn = paste0(x1, x2, x3))

df <- left_join(df, df2, by = "i")

df %>%
  group_by(patrn, t) %>%
  summarize(pcmm = mean(y), pcmm_altc = mean(pcm)) %>%
  print(n = 24)

周辺平均は以下の通り(mmはデータから計算、mm_altcは真値です)。

> df %>%
+   group_by(t, x) %>%
+   summarize(mm = mean(y), mm_altc = mean(beta*x))
# A tibble: 6 × 4
# Groups:   t [3]
      t     x    mm mm_altc
  <int> <dbl> <dbl>   <dbl>
1     1    -1 -1.00      -1
2     1     1  1.00       1
3     2    -1 -1.05      -1
4     2     1  1.06       1
5     3    -1 -1.23      -1
6     3     1  1.01       1

条件付き周辺平均は以下の通りです(pcmmはデータから計算、pcmm_altcは真値です)。

> df %>%
+   group_by(patrn, t) %>%
+   summarize(pcmm = mean(y), pcmm_altc = mean(pcm)) %>%
+   print(n = 24)
# A tibble: 24 × 4
# Groups:   patrn [8]
   patrn      t   pcmm pcmm_altc
   <chr>  <int>  <dbl>     <dbl>
 1 -1-1-1     1 -1.01      -1   
 2 -1-1-1     2 -2.52      -2.5 
 3 -1-1-1     3 -4.77      -4.75
 4 -1-11      1 -0.995     -1   
 5 -1-11      2 -2.50      -2.5 
 6 -1-11      3 -2.75      -2.75
 7 -11-1      1 -1.00      -1   
 8 -11-1      2 -0.515     -0.5 
 9 -11-1      3 -1.79      -1.75
10 -111       1 -0.993     -1   
11 -111       2 -0.489     -0.5 
12 -111       3  0.295      0.25
13 1-1-1      1  0.995      1   
14 1-1-1      2  0.483      0.5 
15 1-1-1      3 -0.274     -0.25
16 1-11       1  1.00       1   
17 1-11       2  0.486      0.5 
18 1-11       3  1.72       1.75
19 11-1       1  1.00       1   
20 11-1       2  2.51       2.5 
21 11-1       3  2.76       2.75
22 111        1  1.02       1   
23 111        2  2.54       2.5 
24 111        3  4.81       4.75

明らかに、$E(Y_{it} \mid X_{it}) \ne E(Y_{it} \mid X_{ij}, j=1,...,n_{i})$になっていることが確認できます。

上記のデータにGEEを適用してみます。

library(geepack)
summary(geeglm(y ~ x, id = i, data = df, corstr = "independence"))
summary(geeglm(y ~ x, id = i, data = df, corstr = "exchangeable"))
summary(geeglm(y ~ x, id = i, data = df, corstr = "unstructured"))

作業相関構造が"independence"以外では、$\hat{\beta}$が$1$に近い値になりませんでした。
あと、geepackではエラーが出ないですが、"unstructured"は推定に失敗していそうです。

> summary(geeglm(y ~ x, id = i, data = df, corstr = "independence"))

Call:
geeglm(formula = y ~ x, data = df, id = i, corstr = "independence")

 Coefficients:
            Estimate Std.err   Wald Pr(>|W|)    
(Intercept)  -0.0345  0.0428   0.65     0.42    
x             1.0584  0.0338 978.72   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Correlation structure = independence 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)     3.31  0.0736
Number of clusters:   1000  Maximum cluster size: 3 

> summary(geeglm(y ~ x, id = i, data = df, corstr = "exchangeable"))

Call:
geeglm(formula = y ~ x, data = df, id = i, corstr = "exchangeable")

 Coefficients:
            Estimate Std.err   Wald Pr(>|W|)    
(Intercept)  -0.0563  0.0535   1.11     0.29    
x             0.3477  0.0176 392.38   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Correlation structure = exchangeable 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)     3.82   0.105
  Link = identity 

Estimated Correlation Parameters:
      Estimate Std.err
alpha    0.624 0.00495
Number of clusters:   1000  Maximum cluster size: 3 

> summary(geeglm(y ~ x, id = i, data = df, corstr = "unstructured"))

Call:
geeglm(formula = y ~ x, data = df, id = i, corstr = "unstructured")

 Coefficients:
            Estimate Std.err Wald Pr(>|W|)   
(Intercept)  -0.0691  0.0651 1.13   0.2878   
x            -0.0860  0.0285 9.10   0.0026 **
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Correlation structure = unstructured 
Estimated Scale Parameters:

            Estimate Std.err
(Intercept)     4.62    0.18
  Link = identity 

Estimated Correlation Parameters:
          Estimate Std.err
alpha.1:2    0.366 0.00869
alpha.1:3    0.556 0.00329
alpha.2:3    1.149 0.00537
Number of clusters:   1000  Maximum cluster size: 3 

おまけ

(制限付き)最尤法ベースでもあまり変わらないですね。

library(nlme)
summary(glm(y ~ x, data = df))
summary(gls(y ~ x, data = df, correlation = corCompSymm(form = ~ 1 | i)))
summary(gls(y ~ x, data = df, correlation = corSymm(form = ~ 1 | i)))
> summary(glm(y ~ x, data = df))

Call:
glm(formula = y ~ x, data = df)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-4.437  -1.191   0.033   1.116   4.719  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.0345     0.0333   -1.04      0.3    
x             1.0584     0.0333   31.83   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 3.31)

    Null deviance: 13292.8  on 2999  degrees of freedom
Residual deviance:  9935.1  on 2998  degrees of freedom
AIC: 12112

Number of Fisher Scoring iterations: 2

> summary(gls(y ~ x, data = df, correlation = corCompSymm(form = ~ 1 | i)))
Generalized least squares fit by REML
  Model: y ~ x 
  Data: df 
    AIC   BIC logLik
  11401 11425  -5696

Correlation Structure: Compound symmetry
 Formula: ~1 | i 
 Parameter estimate(s):
  Rho 
0.625 

Coefficients:
             Value Std.Error t-value p-value
(Intercept) -0.056    0.0535   -1.05   0.293
x            0.347    0.0260   13.37   0.000

 Correlation: 
  (Intr)
x 0.015 

Standardized residuals:
    Min      Q1     Med      Q3     Max 
-2.6229 -0.4340  0.0394  0.4639  2.6845 

Residual standard error: 1.95 
Degrees of freedom: 3000 total; 2998 residual

> summary(gls(y ~ x, data = df, correlation = corSymm(form = ~ 1 | i)))
Generalized least squares fit by REML
  Model: y ~ x 
  Data: df 
    AIC   BIC logLik
  10602 10638  -5295

Correlation Structure: General
 Formula: ~1 | i 
 Parameter estimate(s):
 Correlation: 
  1      2     
2 -0.677       
3  0.062  0.596

Coefficients:
            Value Std.Error t-value p-value
(Intercept) 0.024    0.0192     1.3   0.206
x           0.577    0.0155    37.2   0.000

 Correlation: 
  (Intr)
x 0.04  

Standardized residuals:
      Min        Q1       Med        Q3       Max 
-2.52e+00 -5.56e-01 -3.64e-06  4.87e-01  2.50e+00 

Residual standard error: 1.98 
Degrees of freedom: 3000 total; 2998 residual

感想

完全に理解できたとは言い難いので、$E(Y_{it} \mid X_{it}) = E(Y_{it} \mid X_{ij}, j=1,...,n_{i})$の事例と$E(Y_{it} \mid X_{it}) \ne E(Y_{it} \mid X_{ij}, j=1,...,n_{i})$の事例をいっぱい作ってもうちょっと勉強したいなと思いました。

文献

  • Pepe MS, Anderson GL. A cautionary note on inference for marginal regression models with longitudinal data and general correlated response data. Communications in Statistics - Simulation and Computation 1994; 23(4): 939-951. doi: 10.1080/03610919408813210.

クロスオーバーデザインにおける順序カテゴリデータの解析 (Senn, 1993)

Senn (1993)

クロスオーバーデザインで順序カテゴリアウトカムを用いた場合の解析方法を知らなかったので、このレターを読んでみて、メモを作成しました。
このレターは、Ezzet & Whitehead (1991, Stat Med) に対してコメントをしたもので、Ezzet & Whitehead が提案した変量効果を含む比例オッズモデルが複雑で実装も大変なので、別の方法を提案するという内容でした。

事例のデータとchange score

事例のデータは、クロスオーバー試験において0, 1, 2, 3のような4水準の順序カテゴリデータが得られるケースでした。
治療AとBについて1期と2期がある試験で、ABの順で治療を受ける群と、BAの順で治療を受ける群がランダム化されたものでした。

まず、データをchange scaleに変換したchange scoreを用いることを提案しています。
例えば、1期から2期でスコアが何段階変化したかを求め、-3 を 'much worse', -2 or -1 を 'worse', 0 を 'same', +1 or +2 を 'better', +3 を 'much better' のように変換します。
この例では、4水準の順序カテゴリデータの変化量を5水準に変換していて、partial orderingを用いていると説明されていました。
Partial orderingを用いているので、+1 or +2は変化量が違うのに両方 'better' と変換していて、0から2への変化と2から3への変化を同じとみなすことが妥当かどうかは不明、ということです。
不完全であるということは認めた上での提案のようです。

ちなみに、アウトカムが二値のときの変化量は-1 (worse), 0 (same), 1 (better) に正確に変換できて、クロスオーバー試験において二値アウトカムの結果を変換してから仮説検定する Prescott's test と近い考え方の方法と言える、ということでした。

比例オッズモデルによるchange scoreの解析

データをchange scoreに変換し、比例オッズモデルを用いて解析すると、非常に簡単に順序AB vs BA間の比較ができることが紹介されていました。
McCullagh (1980) の経験ロジット変換の方法を使うと、分割表の数値から手計算で各change score間のロジットの推定値が得られ、それをもとに治療効果パラメータ $\phi$ の近似推定値が得られることが紹介されていました。
以下に示したレターで使われていた事例データに対して、経験ロジット変換の方法を用いた場合の推定結果は $\hat{\phi}=1.571, \text{SE}(\hat{\phi})=0.34$ と書かれていました。

Change score AB BA
Much worse 2 0
Worse 41 13
Same 86 79
Better 12 50
Much better 1 2

Rで比例オッズモデルを用いて解析すると、以下のようになり、レターに書かれていた比例オッズモデルによる推定結果と一致しました ($\hat{\phi}=1.604, \text{SE}(\hat{\phi})=0.26$)。
経験ロジット変換の結果とは、まあまあ近い結果になるようでした。
また、Ezzet & Whiteheadの方法をここで使ったモデルに近い形に修正した結果は、$\hat{\phi}=1.283, \text{SE}(\hat{\phi})=0.21$ であったと書かれていました。

library(MASS)

df <- data.frame(
  sequence = factor(
    c(rep(0, 142), rep(1, 144)),
    labels = c("AB", "BA")),
  change = factor(c(
    rep(0, 2), rep(1, 41), rep(2, 86), rep(3, 12), rep(4, 1),
    rep(1, 13), rep(2, 79), rep(3, 50), rep(4, 2)),
    labels = c("Much worse", "Worse", "Same", "Better", "Much better"))
)
table(df$change, df$sequence)
fit <- polr(change ~ sequence, data = df, Hess = TRUE)
summary(fit)
confint(fit)
> table(df$change, df$sequence)
             
              AB BA
  Much worse   2  0
  Worse       41 13
  Same        86 79
  Better      12 50
  Much better  1  2
> fit <- polr(change ~ sequence, data = df, Hess = TRUE)
> summary(fit)
Call:
polr(formula = change ~ sequence, data = df, Hess = TRUE)

Coefficients:
           Value Std. Error t value
sequenceBA 1.604     0.2637   6.085

Intercepts:
                   Value   Std. Error t value
Much worse|Worse   -4.4364  0.7132    -6.2208
Worse|Same         -0.8032  0.1740    -4.6153
Same|Better         2.2037  0.2331     9.4540
Better|Much better  5.6437  0.6156     9.1680

Residual Deviance: 556.4482 
AIC: 566.4482 
> confint(fit)
Waiting for profiling to be done...
   2.5 %   97.5 % 
1.099457 2.135862 

感想

1993年の論文ということで、計算機事情もかなり違ったようですが、シンプルな解析方法があるのだなと思いました。
このレターでは他にも代替案(Koch et al., 1977; Jones & Kenward, 1989; Koch, 1972)が挙げられており、そちらも少し調べてみようかなと思いました。
あと、治療効果パラメータは推定できるのですが、どういう解釈ができるのかについては悩ましいなと思っています。

文献

  • Senn S. A random effects model for ordinal responses from a crossover trial. by F. Ezzet and J. Whitehead, Statistics in Medicine, 10, 901–907 (1991). Statistics in Medicine 1993; 12(12): 2147–2151. doi: 10.1002/sim.4780122208.
  • Ezzet F, Whitehead J. A random effects model for ordinal responses from a crossover trial. Statistics in Medicine 1991; 10(6): 901–907. doi: 10.1002/sim.4780100611.
  • Prescott R. The comparison of success rates in cross-over trials in the presence of an order effect. Applied Statistics 1981; 30(1): 9–15. doi: 10.2307/2346652.
  • McCullagh P. Regression model for ordinal data (with discussion). Journal of the Royal Statistical Society, Series B 1980; 42(2): 109–142. doi: 10.1111/j.2517-6161.1980.tb01109.x.
  • Koch GG, Landis JR, Freeman JL, Freeman Jr DH, Lehnen RC. A general methodology for the analysis of experiments with repeated measurement of categorical data. Biometrics 1977; 33(1): 133–158. doi: 10.2307/2529309.
  • Jones B, Kenward M. Design and Analysis of Cross-over Trials, Chapman and Hall, London and New York, 1989.
  • Koch GG. The use of non-parametric methods in the statistical analysis of the two-period change-over design. Biometrics 1972; 28(2): 577–584. doi: 10.2307/2556170.

1:1マッチングの場合の2×2×K表の共通オッズ比の無条件の最尤推定量

$K$層の積二項分布モデルの無条件の尤度関数から導かれるスコア方程式

$k$番目の層のCase群の曝露ありの人数が二項分布$X_{11k} \sim Bin(n_{1+k}, p_{1k})$に従い、$k$番目の層のControl群の曝露ありの人数が二項分布$X_{01k} \sim Bin(n_{0+k}, p_{0k})$に従い、$k$番目の層のオッズ比を$\psi_k=\frac{p_{1k}/(1-p_{1k})}{p_{0k}/(1-p_{0k})}$で表したもとで、共通オッズ比 ($\psi_1=\ldots=\psi_K=\psi$) の仮定をおいたときの積二項分布モデルの対数尤度関数は、$\mathbf{p}_{0}=(p_{01},\ldots,p_{0K})^T$とおくと、
\[
\begin{aligned}
l(\psi, \mathbf{p}_{0})=
\sum_{k=1}^K &
\log \binom{n_{1+k}}{x_{11k}}+
x_{11k} \log \psi p_{0k} -
x_{11k} \log \{1-p_{0k}(1-\psi)\} +\\&
(n_{1+k}-x_{11k}) \log (1-p_{0k})-
(n_{1+k}-x_{11k}) \log \{1-p_{0k}(1-\psi)\}+\\&
\log \binom{n_{0+k}}{x_{01k}}+
x_{01k}\log p_{0k}+
(n_{0+k}-x_{01k})\log (1-p_{0k})
\end{aligned}
\] となります。

対数尤度関数を$\psi$で偏微分すると、
\[
\frac{\partial l(\psi, \mathbf{p}_{0})}{\partial \psi}=
\sum_{k=1}^K
\frac{x_{11k}}{\psi}-
n_{1+k}\frac{p_{0k}}{1-p_{0k}(1-\psi)}
\] が得られ、これを$=0$と置いて解くことで
\[
\sum_{k=1}^K x_{11k}=
\sum_{k=1}^K n_{1+k}\frac{\hat{\psi} \hat{p}_{0k}}{1-\hat{p}_{0k}(1-\hat{\psi})}
\] という関係式が得られます。
また、定義から$\frac{\hat{\psi} \hat{p}_{0k}}{1-\hat{p}_{0k}(1-\hat{\psi})}=\hat{p}_{1k}$なので、$\sum_{k=1}^K x_{11k}=\sum_{k=1}^K n_{1+k} \hat{p}_{1k}=\sum_{k=1}^K E^A[X_{11k}; \psi]$と書くこともできます。

次に、対数尤度関数を$p_{0k}$で偏微分すると、
\[
\begin{aligned}
\frac{\partial l(\psi, \mathbf{p}_{0})}{\partial p_{0k}}&=
\frac{x_{11k}}{p_{0k}}+
\frac{x_{11k}(1-\psi)}{1-p_{0k}(1-\psi)}-
\frac{n_{1+k}-x_{11k}}{1-p_{0k}}+
\frac{(n_{1+k}-x_{11k})(1-\psi)}{1-p_{0k}(1-\psi)}+
\frac{x_{01k}}{p_{0k}}-
\frac{n_{0+k}-x_{01k}}{1-p_{0k}}\\&=
\frac{x_{11k}+x_{01k}}{p_{0k}}-
\frac{n_{1+k}+n_{0+k}-x_{11k}-x_{01k}}{1-p_{0k}}+
\frac{n_{1+k}(1-\psi)}{1-p_{0k}(1-\psi)}\\&=
\frac{1}{1-p_{0k}}\left(
\frac{x_{11k}+x_{01k}}{p_{0k}}-
\frac{n_{1+k}\psi}{1-p_{0k}(1-\psi)}-
n_{0+k}
\right)
\end{aligned}
\] が得られます。

1:1マッチングの場合の$\hat{p}_{0k}$

1:1マッチングの場合、$n_{1+k}=n_{0+k}=1$となり、$x_{11k}+x_{01k}=0, 1, 2$の3パターンだけになります。
下で使うため、3パターンの場合について対数尤度関数を$p_{0k}$で偏微分して$=0$とおいて解き、それぞれの場合の$\hat{p}_{0k}$を導出しておきます。

$x_{11k}+x_{01k}=0$のとき

\[
\frac{0}{\hat{p}_{0k}}-
\frac{\hat{\psi}}{1-\hat{p}_{0k}(1-\hat{\psi})}-
1=0
\iff
\hat{p}_{0k}=-\frac{\hat{\psi}\hat{p}_{0k}}{1-\hat{p}_{0k}(1-\hat{\psi})}=-\hat{p}_{1k}
\] $\hat{p}_{0k}=-\hat{p}_{1k}$が成り立つのは$\hat{p}_{0k}=0$の時しかないため、$\hat{p}_{0k}=0$が得られます。

$x_{11k}+x_{01k}=1$のとき

\[
\frac{1}{\hat{p}_{0k}}-
\frac{\hat{\psi}}{1-\hat{p}_{0k}(1-\hat{\psi})}-
1=0
\iff
\hat{p}_{0k}=\frac{1}{1 \pm \sqrt{\hat{\psi}}}
\] $\hat{p}_{0k}=\frac{1}{1-\sqrt{\hat{\psi}}}$の方は$\hat{p}_{1k}$が負になるので、$\hat{p}_{0k}=\frac{1}{1+\sqrt{\hat{\psi}}}$が得られます。

$x_{11k}+x_{01k}=2$のとき

\[
\frac{2}{\hat{p}_{0k}}-
\frac{\hat{\psi}}{1-\hat{p}_{0k}(1-\hat{\psi})}-
1=0
\iff
\hat{p}_{0k}=1,\frac{2}{1-\hat{\psi}}
\] $\hat{p}_{0k}=\frac{2}{1-\hat{\psi}}$の方は$\hat{p}_{0k}>1$か$\hat{p}_{0k}<0$になってしまうので、$\hat{p}_{0k}=1$が得られます。

1:1マッチングの場合の$\hat{\psi}$

上で求めた対数尤度関数を$\psi$で偏微分して$=0$とおいて解いた式について、1:1マッチングの場合を考えることで$\hat{\psi}$を得ることができます。
全体で$K$個の層がありますが、

  • $x_{11k}=1, x_{01k}=1$となっている層の数を$n_{11}$
  • $x_{11k}=1, x_{01k}=0$となっている層の数を$n_{10}$
  • $x_{11k}=0, x_{01k}=1$となっている層の数を$n_{01}$
  • $x_{11k}=0, x_{01k}=0$となっている層の数を$n_{00}$

と書くことにします ($K=n_{11}+n_{10}+n_{01}+n_{00}$です)。
$\sum_{k=1}^K x_{11k}=\sum_{k=1}^K n_{1+k}\frac{\hat{\psi} \hat{p}_{0k}}{1-\hat{p}_{0k}(1-\hat{\psi})}$に対して、今度は4パターンについて、上で導出した$\hat{p}_{0k}$も代入してそれぞれ当てはめると、
\[
\sum_{k=1}^K x_{11k}=\sum_{k=1}^K n_{1+k}\frac{\hat{\psi} \hat{p}_{0k}}{1-\hat{p}_{0k}(1-\hat{\psi})}
\iff
n_{11}+n_{10}+0+0=
n_{11}+n_{10}\frac{\sqrt{\hat{\psi}}}{1+\sqrt{\hat{\psi}}}+n_{01}\frac{\sqrt{\hat{\psi}}}{1+\sqrt{\hat{\psi}}}+0
\] が得られます。
これを解くと、
\[
\hat{\psi}=\left(\frac{n_{10}}{n_{01}}\right)^2
\] が得られ、1:1マッチングの時の無条件の最尤推定量が条件付き最尤推定量$\hat{\psi}_C=\frac{n_{10}}{n_{01}}$を二乗したものになっていることがわかります。

参考文献

  • Breslow N. Odds ratio estimators when the data are sparse. Biometrika 1981; 68(1): 73–84.

ロジスティック回帰版

ironman さんが解いています。