« 読了：「高い窓」「プレイバック」「マルタの鷹」「隠し剣孤影抄」「隠し剣秋風抄」「沿岸 頼むから静かに死んでくれ」 | メイン | 読了：Richarson, et al. (2018) 共変量のある多項回帰モデルを、共変量で予測した曝露確率の逆数でウェイティングして共変量なしで済ませる »
2018年5月 5日 (土)
対応のある二値データについて考える。たとえば、n人の人に製品PとQをみせ、それぞれについて買ってみたいかどうかを判断してもらった、というような場合ですね。両方の製品について「はい」と答えた人数を$a$, 製品Qについてのみ「はい」と答えた人数を$b$, 製品Pについてのみ「はい」と答えた人数を$c$, 両方の製品について「いいえ」と答えた人数を$d$とする(当然ながら$a+b+c+d=n$)。
製品PとQのあいだで「はい」率は異なるか。これを調べる古典的な方法が、ご存じMcNemar検定である。母集団における「はい」率が製品P,Q間で同一であるという帰無仮説の下で、$(b-c)^2/(b+c)$が自由度1のカイ二乗分布に従うことを利用する。
分野の違いというのは面白いもので、市場調査の参考書ではこれを「対応のあるZ検定」と呼び、検定統計量$(b-c)/\sqrt{b+c}$を$N(0,1)$と比べよ、と書いてあることが多い。まあどっちでも同じことである。
ところがこのMcNemar検定、自分で手計算してみるとわかるけど、ちょっと不思議な面がある。
たとえば、サンプルサイズは$n=100$、うちQについてのみ「はい」と答えた人数が$b=15$, Pについてのみ「はい」と答えた人数が$c=10$だったとしよう。このときの検定統計量の値は$(b-c)^2/(b+c)=25/25=1$。
今度は、サンプルサイズは$n=100000$, うち$b=15, c=10$だったとしよう。このときも$(b-c)^2/(b+c)=25/25=1$。
つまり、サンプルサイズがどんなに変わろうと、2x2クロス表の対角セル$a, d$がどんなに変わろうと、非対角セル$b$と$c$さえ同じなら、検定統計量の値は同じなのだ。
でも、そもそもここで比較したいのは、Pに「はい」と答えた人の割合$(a+c)/n$と、Qに「はい」と答えた人の割合$(a+b)/n$でしょう？ サンプルサイズ$n$が大きいとき、これらの割合の信頼性は高くなるはずでしょう？
...という素朴な疑問は、私自身がかつて感じたことでもあり、前職の社内研修のときに訊ねられたことでもある。このたび仕事のなかで、ふとこの話を思い出す機会があり、試しにAgrestiの分厚い教科書をめくってみたら、その話はAgresti & Min (2003)に書いたから読んでね、とのことであった。
Agresti, A., & Min, Y. (2003) Effects and non-effects of paired identical observations in comparing proportions with binary matched-pairs data. Statistics in Medicine, 22. 65-75.
というわけで、読んでみました。難しくてよくわかんない箇所が多かったんだけど...
あるペア(たとえばあるヒト)$i=1, \ldots, n$が提供する$t=1,2$個目の値を$y_{it}$とする。値$1$を成功、$0$を失敗と呼ぶ。
このデータを表す形式を2種類考えよう。
- subject-specific形式。表側を$t$, 表頭を(成功, 失敗)とした2x2の二値行列を$n$枚縦積みする。2列のどっちかに1がはいっているわけだ。
- population-averaged形式。表側を$t=1$の(成功, 失敗), 表頭を$t=2$の(成功, 失敗)とした2x2のクロス表。セルのなかには頻度$a,b,c,d$がはいっていて、これらは多項分布$(n, \{\pi_{jk}\})$の実現値と捉えることができる。
subject-specific形式のデータに当てはめられる標準的モデルは、
$link[P(y_{i1})=1] = \alpha_i$
$link[P(y_{i2})=1] = \alpha_i+\beta_c$
であろう。これをconditional modelという(効果$\beta_c$の定義がペアに条件づけられているという意味)。linkというのはたとえばlogit。モデル・フィッティングにおいては、$(y_{i1}, y_{i2})$が$\{\alpha_i\}$の下で独立だと仮定するのが普通である。
いっぽう、ランダムに取り出したあるヒトの最初の観察を$y_1$, ランダムに取り出した別のヒトの二番目の観察を$y_2$として
$link[P(y_1)=1] = \alpha$
$link[P(y_2)=1] = \alpha+\beta_m$
とモデル化するという手もある。これをmarginal modelという。ペアの母集団を通じた平均効果のモデル化である。population-averaged形式の2x2クロス表の周辺分布のモデル化だといえる。
ここでちょっとこの論文のメモから離れて、Agrestiの"An Introduction to Categorical Data Analysis"からメモしておくと；
一般化線形モデルのパラメータ$\beta$について$H_0: \beta = 0$を検定する方法には３つある。
- ワルド検定。$\beta$の最尤推定量を$\hat{\beta}$とする。$\hat{\beta}$は大標本の下で正規性を持つと仮定できるので、漸近標準誤差をASEと略記して
$z^2 = (\hat{\beta}/ASE)^2$
は$H_0$の下で自由度1のカイ二乗分布に従う。 - 尤度比検定。$H_0$のもとで最大化した尤度を$\ell_0$(その対数を$L_0$), 対立仮説のもとで最大化した尤度を$\ell_1$(その対数を$L_1$)とすると、
$-2 \log(\ell_0/\ell_1) = -2(L_0 - L_1)$
は$H_0$の下で自由度1のカイ二乗分布に従う。 - スコア検定。対数尤度関数は上に凸で、その頂点が$\hat{\beta}$である。$\beta_0$における傾き(つまり導関数の値)とそのASEの比を二乗した値をスコア統計量という。これもまた$H_0$の下で自由度1のカイ二乗分布に従う。モデルをあてはめなくても得られる、というのがポイント。
論文に戻って。。。
$\pi_1=P(y_1=1)$と$\pi_2=P(y_2=1)$が等しいという帰無仮説の検定について考えよう。この帰無仮説は$\beta_m=0$と等価であり、$\pi_{12}=\pi_{21}$と等価である。
marginal modelに基づいて考えると、$\beta_m=0$という条件の下で最大化された尤度と、条件なしで最大化された尤度の比は、$b,c$だけに依存している。尤度比統計量(=-2*対数尤度比)は$2b \log(2b/(b+c)) + 2c \log(2c/(b+c))$となる。
$b+c$が小さい場合、$c$を試行数$b+c$, パラメータ$1/2$の二項分布と比べる正確検定が可能である。大標本の場合、これは正規近似できて、
$\displaystyle z=\frac{c-(1/2)(b+c)}{\sqrt{(b+c)(1/2)(1/2)}} = \frac{c-b}{\sqrt{b+c}}$
これを二乗したのがMcNemar統計量である。これは周辺分布が等しいという検定のためのスコア統計量になっている。
ワルド統計量はどうか。ワルド統計量は、推定された効果をSEで割ったものである。だからリンク関数の定義によって変わってくる。でも後述するように、$a,c$への依存性は無視できる程度である。
このように、この検定において情報を提供するのは$b$と$c$だけである。
ここからは、割合の差、オッズ比、相対リスクの推定、という３つの観点から考える。
まず割合の差の推定について。ここ、短いんだけど途中からいきなりわかりにくくなるので、思い余って全訳した。
identity linkの場合、subject-specificな効果とpopulation-averagedな効果は同一である。たとえばconditional modelでは、すべての$i$について$\beta_c = P(Y_{i2}=1)-P(Y_{i1}=1)$である。これを母集団内のペアを通じて平均すると、marginal modelにおける周辺確率の差$\beta_m=\pi_2-\pi_1$と等しくなる。
別の見方として次の結果に着目することもできる。indentity linkのconditional modelを$2 \times 2 \times n$表に適用したとき、層別変数が説明変数と周辺独立であるときにcollapsibilityが成立する。このとき、これら2変数の標本クロス表は全セルが1となる。[←先生すいません、この段落、なにいってんだかさっぱりわかんないです]
標本における割合の差は
$p_2 = p_1 = p_{21} - p_{12} = (c-b)/n$
となる。これはindentity linkのmarginal modelにおける$\beta_m$の最尤推定値である。conditional modelの観点からいえば、これはヒト別の表のMantel-Haenszel型加重平均からも得られる。またこれは、Chen(1996)が示しているように、$\{\alpha_i\}$がなんらかのparametric familyからランダムに抽出されたものだと仮定した、ランダム効果版conditional modelでの最尤推定値になっている。
周辺表の多項抽出において
$Var(p_2-p_1) = [(\pi_{12}+\pi{21}) - (\pi_{21}-\pi_{12})^2]/n$
その標本推定値は
$\hat{Var}(p_2-p_1) = \frac{(b+c)-(c-b)^2/n}{n^2}$
周辺等質性($\pi_1=\pi_2$)という帰無仮説の下で、分散は$(\pi_{12}+\pi_{21})/n$に、その推定値は$(b+c)/n$に、それぞれ帰着する。ここからMcNemar検定が導かれる。$n=a+b+c+d$であるから、帰無仮説が真でない時、効果のサイズの推定値とその標準誤差は、$(a+d)$が増大するにつれて減少する。
ここでも、また以下のモデルでも、$(a+d)$の効果について考える際には、$b$と$c$を固定して$(a+d)$だけを大きくすると考える。$n$自体も大きくなるので、帰無仮説のもとでの漸近性において$a$, $d$が果たす役割をチェックすることができる。
$(p_2-p_1)$の分散の推定値は、$n^{-2}$のオーダーで、$(b+c)/n^2$すなわち帰無仮説の下での分散の推定値になる。同様に、$n$が大きければ、$(p_2-p_1)$とその標準誤差の推定値との比は近似的に$(c-b)/\sqrt{b+c}$となる(これはMcNemar統計量の標準正規形式である)。このように、とりわけ$n$が大きいとき、割合の差のワルド検定における$(a+d)$の寄与はごくわずかである。
では、リンクlogitである場合(パラメータは対数オッズ比になる)は...[パス]
リンクがlogである場合(対数相対リスク)は...[パス]
ベイズ推論の観点から見るとどうか。
marginal modelの場合。$\{\pi_{jk}\}$の事前分布としてパラメータ$\{\mu_{jk}\}$のディリクレ分布を与えるとしよう。$\pi_2 > \pi_1$の事後確率は、$\pi_{21} / (\pi_{12}+\pi_{21}) > 1/2$の事後確率に等しい。左辺の事後分布はパラメータ$(c+\mu_{21}, b+\mu_{12})$のベータ分布になる。だから$a, d$は効かない。ところが、conditional modelになると話が変わってきて...[めんどくさいので省略]
云々、云々。。。というわけで、$a, d$は効果のサイズやそのSEに全然効かないわけじゃないけど、たいして効かない、ということなのだそうである。すいませんAgresti先生、途中で心が折れました。
論文：データ解析(2018-) - 読了：Agresti & Min (2003) 対応のある二条件間で割合を比べるとき、「両方1」「両方0」の事例を無視してよいのか