$ \newcommand\op{\operatorname} $
Mantel, Nathan and Haenszel, William. Statistical Aspects of the Analysis of Data From Retrospective Studies of Disease. Journal of the National Cancer Institute, Volume 22, Issue 4, April 1959, Pages 719–748. [link] [pdf]
Robins, J., S Greenland, S., and Breslow, N. E. A general estimator for the variance of the Mantel-Haenszel odds ratio. Am. J. Epidemiol
. 1986 Nov; 124(5):719-23. [link] [pdf]
共通オッズ比のMantel-Haenszelの推定量と最尤推定量の関係については次の文献に書いてある:
この文献ではMantel-Haenszel以外の共通オッズ比の推定量も扱っている. 共通オッズ比のMantel-Haenszelの推定量の対数の分散のRobins-Breslow-Greenlandの推定量の解説が次の文献に書いてある:
簡単のため $K$ 個の独立な2×2の分割表の各々は2つの二項分布の積分布に従っているという統計モデルを考える. $k$ 番目の分割表を行列
$$ A_k = \begin{bmatrix} a_k & b_k \\ c_k & d_k \\ \end{bmatrix} $$と書き,
$$ N_k = a_k + b_k + c_k + d_k $$とおく. $a_k + b_k = m_k$ と $c_k + d_k = n_k$ は定数であり, $a_k$, $c_k$ はそれぞれ独立な二項分布 $\op{Binomial}(m_k, p_k)$, $\op{Binomial}(n_k, q_k)$ に従っているとする. $k$ 番目のオッズ比 $\omega_k$ を
$$ \omega_k = \frac{p_k(1 - q_k)}{(1 - p_k)q_k} $$と定める. $0<p_k<1$ かつ $0<q_k<1$ のとき $\omega_k = 1$ と $p_k = q_k$ は同値である.
以下では $K$ が固定されており, $m_k$, $n_k$ 達が十分に大きい場合を扱う.
このとき, $K$ 番目のオッズ比の推定量
$$ \hat\omega_k = \frac{a_k d_k}{b_k c_k} $$は平均が $\omega_k$ で分散が
$$ \begin{aligned} v_k &= \omega^2\left( \frac{1}{m_k p_k} + \frac{1}{m_k (1 - p_k)} + \frac{1}{n_k q_k} + \frac{1}{n_k (1 - q_k)} \right) \\ &= \omega^2 \left( \frac{1}{m_k p_k (1 - p_k)} + \frac{1}{n_k q_k (1 - q_k)} \right) \end{aligned} $$の正規分布に近似的に従うことを示せる. (二項分布に関する中心極限定理と所謂デルタ法(逆数を取る操作の一次近似)の簡単な応用で示せる.)
以下では $\omega_1 = \cdots = \omega_K = \omega$ が成立していると仮定する. $\omega$ を 共通オッズ比 と呼ぶ.
この仮定のもとで以上で扱っている統計モデルの独立なパラメータは $q_1,\ldots,q_K, \omega$ の $K+1$ 個になる.
共通オッズ比の Mantel-Haenszelの推定量 $\hat\omega_{\op{MH}}$ を次のように定める:
$$ \hat\omega_{\op{MH}} = \frac {\sum_{k=1}^K a_k d_k/N_k} {\sum_{k=1}^K b_k c_k/N_k}. $$これは次のように書き直される:
$$ \hat\omega_{\op{MH}} = \frac {\sum_{k=1}^K \hat{w}_k \hat\omega_k} {\sum_{k=1}^K \hat{w}_k}. $$ここで
$$ \hat\omega_k = \frac{a_k d_k}{b_k c_k}, \quad \hat{w}_k = \frac{b_k c_k}{N_k} $$大数の法則より, $m_k, n_k$ 達が大きなとき, 近似
$$ \hat\omega_k = \frac{(a_k/m_k)(d_k/n_k)}{(b_k/m_k)(c_k/n_k)}\approx \frac{p_k(1-q_k)}{(1-p_k)q_k} = \omega_k = \omega $$が成立しているので, 近似
$$ \hat\omega_{\op{MH}} \approx \frac {\sum_{k=1}^K \hat{w}_k \omega} {\sum_{k=1}^K \hat{w}_k} = \omega. $$が成立する. すなわち大数の法則による近似が有効な場合にはMantel-Haenszelの推定量 $\hat\omega_{\op{MH}}$ は共通オッズ比 $\omega$ を近似している.
これはMantel-Haenszelの推定量 $\hat\omega_{\op{MH}}$ が共通オッズ比 $\omega$ の一致推定量になっていることを意味している.
共通オッズ比が $1$ の特殊な場合 ($\omega = 1$ の場合)にはさらに良いことを言える.
大数の法則より, 近似
$$ \frac{1}{\hat{w}_k} = \frac{m_k + n_k}{b_k c_k} = \frac{1}{(b_k/m_k)(c_k/n_k)}\left(\frac{1}{m_k} + \frac{1}{n_k}\right) \approx \frac{1}{p_k (1 - q_k)}\left(\frac{1}{m_k} + \frac{1}{n_k}\right) $$が成立しているので, $\omega_k = \omega = 1$ (すなわち $p_k = q_k$) ならば,
$$ \frac{1}{\hat{w}_k} \approx \frac{1}{p_k (1 - p_k)}\left(\frac{1}{m_k} + \frac{1}{n_k}\right) = \omega^2\left(\frac{1}{m_k p_k (1 - p_k)} + \frac{1}{n_k q_k (1 - q_k)}\right) = v_k. $$すなわち $\hat{w}_k$ の逆数は $\hat\omega_k = (a_k d_k)/(b_k c_k)$ の漸近的な分散 $v_k$ を近似する.
一般に, 正の実数達 $\sigma_1^2,\ldots,\sigma_K^2 > 0$ が与えられていて, 統計モデルが共通平均 $\mu$ と既知の分散 $\sigma_k^2$ を持つ $K$ 個の正規分布の積であるとき, 共通平均 $\mu$ の共通平均 $\mu$ の標本 $X_1,\ldots,X_n$ に対する最尤推定量が
$$ \hat\mu = \frac {\sum_{k=1}^K X_k/\sigma_k^2} {\sum_{k=1}^K 1/\sigma_k^2} $$になることを簡単な計算で示せる. これは分散の逆数を重みとする荷重平均になっている. 分散が小さな $X_k$ の方が確率的に共通平均の真の値に近くなる傾向があるはずなので, そのような $X_k$ を大きな重みで足し上げた方が共通平均の良い推定量になりそうである. この荷重平均は実際にそのようになっている. (この荷重平均は $\sigma_1^2 = \cdots = \sigma_K^2 = \sigma^2$ の場合の標本平均 $\bar{X}=\frac{1}{K}\sum_{k=1}^K X_k$ の一般化になっている.)
上で示した共通オッズ比が $1$ の場合 ($\omega_k = \omega = 1$) におけるMantel-Haenszelの推定量 $\hat\omega_{\op{MH}}$ に関する結果は, 共通オッズ比が $1$ の特殊な場合には, Mantel-Haenszelの推定量 $\hat\omega_{\op{MH}}$ が $\mu=\omega=1$, $\sigma_k^2=v_k$ の場合の 漸近的最尤推定量 $\hat\mu$ を近似していることを意味している.
ただし, これは共通のオッズ比が $1$ の特殊な場合におけるMantel-Haenszelの推定量 $\hat\omega_{\op{MH}}$ と最尤推定量の関係に過ぎない.
$K$ 個の独立な2×2の分割表の各々が2つの二項分布の積分布に従っているという統計モデルのパラメータ空間を共通オッズ比 $\omega$ を持つ場合に制限して得られる統計モデルを考えている. 独立なパラメータは $q_1,\ldots,q_K,\omega$ の $K+1$ 個になる. (ロジスティックモデルで書き直すこともできる.)
このモデルにおける標本
$$ A_k = \begin{bmatrix} a_k & b_k \\ c_k & d_k \\ \end{bmatrix} \quad (k = 1,\ldots,K) $$の統計モデルのパラメータ達の最尤推定量 $\hat{q}_1, \ldots, \hat{q}_K, \hat\omega$ を以下のようにして計算できることを示せる. (対数尤度函数を偏微分したものが $0$ になるという方程式を整理すればよい.)
(1) 次の条件を満たす $\delta_k = \delta_k(\omega)$ を求める:
$$ \frac{(a_k - \delta_k)(d_k - \delta_k)}{(b_k + \delta_k)(c_k + \delta_k)} = \omega, \quad -\min(b_k, c_k) < \delta_k < \min(a_k, d_k) $$具体的に $\delta_k = \delta_k(\omega)$ は2次方程式を解いて次のように表される:
$$ \delta_k = \delta_k(\omega) = \frac{2C_k}{B_k + \sqrt{B_k^2 - 4A_k C_k}}. $$ここで
$$ A_k = 1 - \omega, \quad B_k = a_k + d_k + \omega(b_k + c_k), \quad C_k = a_k d_k - \omega b_k c_k. $$このとき, $\delta_k = \delta_k(\omega)$ は $\omega$ の単調減少函数になり, 次を満たしている:
$$ \delta_k(0) = \min(a_k, d_k), \quad \delta_k(1) = \frac{a_k d_k - b_k c_k}{N_k}, \quad \delta_k(\hat\omega_k) = 0, \quad \delta_k(\infty) = -\min(b_k, c_k). $$$\hat\omega_k = (a_k d_k)/(b_k c_k)$ と定めたのであった.
(2) $x > 0$ に関する次の方程式の解を $\hat\omega$ とする:
$$ \sum_{k=1}^K \delta_k(x) = 0. \tag{$*$} $$これを共通オッズ比の 最尤方程式 と呼ぶことにする. さらに
$$ \begin{bmatrix} \hat{a}_k & \hat{b}_k \\ \hat{c}_k & \hat{d}_k \\ \end{bmatrix} = \begin{bmatrix} a_k - \delta_k(\hat\omega) & b_k + \delta_k(\hat\omega) \\ c_k + \delta_k(\hat\omega) & d_k - \delta_k(\hat\omega) \\ \end{bmatrix} $$とおき, $\hat{q}_k = \hat{c}_k/n_k$ と定める. このとき, $p_k$ の最尤推定量は $\hat{p}_k = \hat{a}_k/m_k$ になる.
この節の内容は上で引用した Yamagimoto-Yamamoto (1985) の Proposition 1 の紹介である.
$y = \delta_k(x)$ は $x = \hat\omega_k = (a_k d_k)/(b_k c_k)$ で $y = 0$ となり, $x=1$ で $y = (a_k d_k - b_k c_k)/N_k$ となる. その2点を通る直線
$$ y = \frac{(a_k d_k - b_k c_k)/N_k}{1 - (a_k d_k)/(b_k c_k)}\left(x - \frac{a_k d_k}{b_k c_k}\right) = -\frac{b_k c_k}{N_k}x + \frac{a_k d_k}{N_k} $$で $y = \delta_k(x)$ を近似することにしよう(これはかなり大胆な近似である).
この近似を採用すると, 最尤方程式 ($*$) を次の一次方程式 ($**$) で近似することになる:
$$ -\left(\sum_{k=1}^K \frac{b_k c_k}{N_k}\right)x + \sum_{k=1}^K \frac{a_k d_k}{N_k} = 0. \tag{$**$} $$これを共通オッズ比の 線形近似最尤方程式 と呼ぶことにする. この線形近似最尤方程式の解はMantel-Haenszel推定量
$$ \hat\omega_{\op{MH}} = \frac {\sum_{k=1}^K a_k d_k/N_k} {\sum_{k=1}^K b_k c_k/N_k}. $$に一致する. このように共通オッズ比のMantel-Haenszel推定量は最尤方程式($*$)を線形近似最尤方程式($**$)で近似することによって得られる.
この節の内容は上で引用した Silcocks (2005) の紹介になっている.
Mantel-Haenszelの推定量は次のように表される:
$$ \hat\omega_{\op{MH}} = \frac{\sum_{k=1}^K a_k d_k/N_k}{\sum_{k=1}^K b_k c_k/N_k} = \frac{\sum_{k=1}^K \hat{w}_k \hat\omega_k}{\sum_{k=1}^K \hat{w}_k}, \quad \hat\omega_k = \frac{a_k d_k}{b_k c_k}, \quad \hat{w}_k = \frac{b_k c_k}{N_k} $$デルタ法(対数の一次近似)によって,
$$ \op{var}(\log \hat\omega_{\op{MH}}) \approx \frac{\op{var}(\hat\omega_{\op{MH}})}{\hat\omega_{\op{MH}}^2}. $$$k$ 番目の分割表のオッズ比の最尤推定量 $\hat\omega_k$ は近似的に平均が共通オッズ比 $\omega$ で分散が
$$ \op{var}(\hat\omega_k) \approx \omega^2\left( \frac{1}{m_k p_k} + \frac{1}{m_k (1 - p_k)} + \frac{1}{n_k q_k} + \frac{1}{n_k (1 - q_k)} \right) $$の正規分布に近似的に従うことを示せるのであった. さらに,
$$ \op{var}(\hat\omega_k) \approx \hat\omega_{\op{MH}}^2\left(\frac{1}{a_k}+\frac{1}{b_k}+\frac{1}{c_k}+\frac{1}{d_k}\right) $$という近似も使えると仮定し, 重み達 $\hat{w}_k$ が定数であるかのような近似が成立していると仮定する. このとき,
$$ \op{var}(\hat\omega_{\op{MH}}) \approx \frac{\sum_{k=1}^K \hat{w}_k^2 \op{var}(\hat\omega_k)}{\left(\sum_{k=1}^K \hat{w}_k\right)^2}. $$以上の近似をまとめると次の近似が得られる:
$$ \op{var}(\log\hat\omega_{\op{MH}}) \approx \frac{\op{var}(\hat\omega_{\op{MH}})}{\hat\omega_{\op{MH}}^2} \approx \frac {\sum_{k=1}^K \hat{w}_k^2\left(\frac{1}{a_k}+\frac{1}{b_k}+\frac{1}{c_k}+\frac{1}{d_k}\right)} {\left(\sum_{k=1}^K \hat{w}_k\right)^2}. $$さらに, $\hat\omega_{\op{MH}} \approx (a_k d_k)/(b_k c_k)$ という近似も使えると仮定する. このとき
$$ \frac{1}{a_k}+\frac{1}{b_k}+\frac{1}{c_k}+\frac{1}{d_k} = \frac{a_k + d_k}{a_k d_k}+\frac{b_k + c_k}{b_k c_k} \approx \frac{1}{b_k c_k}\left(\frac{a_k + d_k}{\hat\omega_{\op{MH}}}+ b_k + c_k\right). $$これと $\hat{w}_k = b_k c_k/N_k$ を使うと, 上で示した近似式より,
$$ \op{var}(\log\hat\omega_{\op{MH}}) \approx \frac {\sum_{k=1}^K b_k c_k (\hat\omega_{\op{MH}}^{-1}(a_k + d_k)+ b_k + c_k) / N_k^2} {\left(\sum_{k=1}^K b_k c_k/N_k\right)^2} = \frac {\sum_{k=1}^K \hat\omega_{\op{MH}}^{-1} b_k c_k (a_k + d_k + \hat\omega_{\op{MH}}(b_k + c_k)) / N_k^2} {\left(\sum_{k=1}^K b_k c_k/N_k\right)^2}. $$$(a_k, b_k)$ と $(c_k, d_k)$ を交換すると, $\hat\omega_{\op{MH}}$ は逆数になり, $\log\hat\omega_{\op{MH}}$ は $-1$ 倍になって分散は変わらないので,
$$ \op{var}(\log\hat\omega_{\op{MH}}) \approx \frac {\sum_{k=1}^K a_k d_k (\hat\omega_{\op{MH}}(b_k + c_k) + a_k + d_k) / N_k^2} {\left(\sum_{k=1}^K a_k d_k/N_k \right)^2}. $$記号の簡単のため, ここで次のようにおく:
$$ P_k = \frac{a_k + d_k}{N_k}, \quad Q_k = \frac{b_k + c_k}{N_k}, \quad R_k = \frac{a_k d_k}{N_k}, \quad S_k = \frac{b_k c_k}{N_k}, \quad R = \sum_{k=1}^K R_k, \quad S = \sum_{k=1}^K S_k. $$このとき, 上で示した2つの近似式の平均を取ると,
$$ \begin{aligned} \op{var}(\log\hat\omega_{\op{MH}}) & \approx \frac {\sum_{k=1}^K (R^2\hat\omega_{\op{MH}}^{-1} b_k c_k + S^2 a_k d_k) (a_k + d_k + \hat\omega_{\op{MH}}(b_k + c_k)) / N_k^2} {2 R^2 S^2} \\ & = \frac {\sum_{k=1}^K (R^2\hat\omega_{\op{MH}}^{-1} S_k + S^2 R_k) (P_k + \hat\omega_{\op{MH}}Q_k)} {2 R^2 S^2} \end{aligned} $$さらに分子分母を $S^2$ で割って, $\hat\omega_{\op{MH}} = R/S$ を使うと,
$$ \begin{aligned} \op{var}(\log\hat\omega_{\op{MH}}) & \approx \frac {\sum_{k=1}^K (\hat\omega_{\op{MH}} S_k + R_k) (P_k + \hat\omega_{\op{MH}}Q_k)} {2 R^2} \\ & = \frac {\sum_{k=1}^K (P_k R_k + \hat\omega_{\op{MH}}(P_k S_k + Q_k R_k) + \hat\omega_{\op{MH}}^2 Q_k S_k)} {2 R^2} \\ & = \frac{\sum_{k=1}^K P_k R_k}{2R^2} + \frac{\sum_{k=1}^K (P_k S_k + Q_k R_k)}{2RS} + \frac{\sum_{k=1}^K Q_k S_k}{2S^2}. \end{aligned} $$この最後の量を Mantel-Haenszelの推定量の対数の分散のRobins-Breslow-Greenlandの推定量 と呼び, $V_{\op{RBG}}$ と書くことにする.
Mantel-Haenszelの推定量の対数 $\log\hat\omega_{\op{MH}}$ は平均 $\log \omega$ (共通オッズ比の対数), 分散 $V_{\op{RBG}}$ の正規分布に近似的に従う.
この事実を使って共通オッズ比に関するP値函数や信頼区間を構成することができる.
以下の数値例を見ればわかるように, このRBG公式は非常に良い公式であり, 最尤法(スコア検定)によって構成したP値函数をRBG公式を使って構成したP値函数は近似的によく一致する.
using RCall
using Distributions
using StatsPlots
default(fmt = :png)
using Roots
safediv(x, y) = x == 0 ? x/one(y) : x/y
safesqrt(x) = √max(0, x)
oddsratio(a, b, c, d) = safediv(a*d, b*c)
oddsratio(A::AbstractVecOrMat) = oddsratio(A...)
function delta(a, b, c, d, ω)
A, B, C = 1 - ω, a + d + ω*(b + c), a*d - ω*b*c
safediv(2C, B + safesqrt(B^2 - 4A*C))
end
delta(A, ω) = delta(A..., ω)
function linear_approx_delta(a, b, c, d, ω)
N = a + b + c + d
-(b*c/N)*ω + a*d/N
end
linear_approx_delta(A::AbstractVecOrMat, ω) = linear_approx_delta(A..., ω)
eachmatrix(A::AbstractArray{<:Any, 3}) = eachslice(A; dims=3)
function maximum_likelihood_equation(A::AbstractArray{<:Any, 3}, x)
sum(A -> delta(A, x), eachmatrix(A))
end
function linear_approx_maximum_likelihood_equation(A::AbstractArray{<:Any, 3}, x)
sum(A -> linear_approx_delta(A, x), eachmatrix(A))
end
function mantel_haenszel_estimator(A::AbstractArray{<:Any, 3})
num = sum(A -> A[1,1]*A[2,2]/sum(A), eachmatrix(A))
den = sum(A -> A[1,2]*A[2,1]/sum(A), eachmatrix(A))
num/den
end
function maximum_linkelihood_estimator(A::AbstractArray{<:Any, 3})
f(t) = maximum_likelihood_equation(A, exp(t))
logomegahat = find_zero(f, 0.0, Order2())
omegahat = exp(logomegahat)
end
maximum_linkelihood_estimator (generic function with 1 method)
function vardelta(a, b, c, d, ω)
N = a + b + c + d
δ = delta(a, b, c, d, ω)
(N - 1)/N/(1/(a - δ) + 1/(b + δ) + 1/(c + δ) + 1/(d - δ))
end
vardelta(A::AbstractVecOrMat, ω) = vardelta(A..., ω)
function chisq_mantel_haenszel(A::AbstractArray{<:Any, 3}, ω = 1.0)
num = sum(A -> delta(A, ω), eachmatrix(A))^2
den = sum(A -> vardelta(A, ω), eachmatrix(A))
num/den
end
function pvalue_chisq(A::AbstractArray{<:Any, 3}, ω = 1.0)
chisq = chisq_mantel_haenszel(A, ω)
ccdf(Chisq(1), chisq)
end
function ci_chisq(A::AbstractArray{<:Any, 3}, α = 0.05)
f(t) = pvalue_chisq(A, exp(t)) - α
logci = find_zeros(f, -1e1, 1e1)
exp(first(logci)), exp(last(logci))
end
function mh_chisq(A::AbstractArray{<:Any, 3}; ω₀ = 1.0, α = 0.05)
common_odds_ratio = maximum_linkelihood_estimator(A)
chisq = chisq_mantel_haenszel(A, ω₀)
df = 1
p_value = ccdf(Chisq(df), chisq)
conf_int = ci_chisq(A, α)
(; common_odds_ratio, ω₀, p_value, α, conf_int, chisq, df)
end
mh_chisq (generic function with 1 method)
function mantel_haenszel_robins_breslow_greenland(A::AbstractArray{<:Any, 3})
@views a, b, c, d = A[1,1,:], A[1,2,:], A[2,1,:], A[2,2,:]
abcd = zip(a, b, c, d)
R = sum(((a, b, c, d),) -> a*d/(a+b+c+d), abcd)
S = sum(((a, b, c, d),) -> b*c/(a+b+c+d), abcd)
logOR = log(R) - log(S)
sumPR = sum(((a, b, c, d),) -> (a+d)*a*d/(a+b+c+d)^2, abcd)
sumPSplusQR = sum(((a, b, c, d),) -> ((a+d)*b*c + (b+c)*a*d)/(a+b+c+d)^2, abcd)
sumQS = sum(((a, b, c, d),) -> (b+c)*b*c/(a+b+c+d)^2, abcd)
SE² = sumPR/(2R^2) + sumPSplusQR/(2(R*S)) + sumQS/(2S^2)
SE = √SE²
(; logOR, SE)
end
function pvalue_mhrbg(A::AbstractArray{<:Any, 3}, ω = 1.0)
(; logOR, SE) = mantel_haenszel_robins_breslow_greenland(A)
normal = Normal(logOR, SE)
min(1, 2cdf(normal, log(ω)), 2ccdf(normal, log(ω)))
end
function ci_mhrbg(A::AbstractArray{<:Any, 3}, α = 0.05)
(; logOR, SE) = mantel_haenszel_robins_breslow_greenland(A)
normal = Normal(logOR, SE)
exp.(quantile.(normal, (α/2, 1 - α/2)))
end
function mh_rbg(A::AbstractArray{<:Any, 3}; ω₀ = 1.0, α = 0.05)
(; logOR, SE) = mantel_haenszel_robins_breslow_greenland(A)
common_odds_ratio = exp(logOR)
normal = Normal(logOR, SE)
p_value = min(1, 2cdf(normal, log(ω₀)), 2ccdf(normal, log(ω₀)))
conf_int = exp.(quantile.(normal, (α/2, 1 - α/2)))
(; common_odds_ratio, ω₀, p_value, α, conf_int)
end
mh_rbg (generic function with 1 method)
R"""
library(tidyverse)
df.forMH <- tribble(
~X, ~L, ~Y, ~n,
0, 0, 0, 325,
0, 0, 1, 273,
0, 1, 0, 324,
0, 1, 1, 363,
1, 0, 0, 292,
1, 0, 1, 321,
1, 1, 0, 278,
1, 1, 1, 323
)
df.MH1 <- uncount(df.forMH, weights = n)
"""
R"""
library(samplesizeCMH)
partial_tables <- table(df.MH1[,c("X","Y","L")])
partial_tables
"""
┌ Warning: RCall.jl: -- Attaching packages --------------------------------------- tidyverse 1.3.1 -- │ v ggplot2 3.3.5 v purrr 0.3.4 │ v tibble 3.1.6 v dplyr 1.0.8 │ v tidyr 1.2.0 v stringr 1.4.0 │ v readr 2.1.2 v forcats 0.5.1 │ -- Conflicts ------------------------------------------ tidyverse_conflicts() -- │ x dplyr::filter() masks stats::filter() │ x dplyr::lag() masks stats::lag() └ @ RCall D:\.julia\packages\RCall\6kphM\src\io.jl:172
RObject{IntSxp} , , L = 0 Y X 0 1 0 325 273 1 292 321 , , L = 1 Y X 0 1 0 324 363 1 278 323
R"""
result_mh <- mantelhaen.test(partial_tables)
"""
RObject{VecSxp} Mantel-Haenszel chi-squared test with continuity correction data: partial_tables Mantel-Haenszel X-squared = 3.319, df = 1, p-value = 0.06848 alternative hypothesis: true common odds ratio is not equal to 1 95 percent confidence interval: 0.9921162 1.3588118 sample estimates: common odds ratio 1.161077
R"""
mantelhaen.test(partial_tables)
library(jtools)
log.model <- glm(Y ~ X + L, data=df.MH1, family = binomial)
result_logistic <- summ(log.model, exp=T, digits=5)
"""
RObject{VecSxp} MODEL INFO: Observations: 2499 Dependent Variable: Y Type: Generalized linear model Family: binomial Link function: logit MODEL FIT: χ²(2) = 7.90056, p = 0.01925 Pseudo-R² (Cragg-Uhler) = 0.00421 Pseudo-R² (McFadden) = 0.00228 AIC = 3460.95990, BIC = 3478.43084 Standard errors: MLE ---------------------------------------------------------------------- exp(Est.) 2.5% 97.5% z val. p ----------------- ----------- --------- --------- ---------- --------- (Intercept) 0.89251 0.77739 1.02469 -1.61387 0.10656 X 1.16124 0.99218 1.35911 1.86222 0.06257 L 1.19097 1.01759 1.39389 2.17714 0.02947 ----------------------------------------------------------------------
R"df_MH1 <- df.MH1"
@rget df_MH1
@rget partial_tables
2×2×2 Array{Int64, 3}: [:, :, 1] = 325 273 292 321 [:, :, 2] = 324 363 278 323
A = partial_tables
@show A
@show mh_chisq(A)
@show mh_rbg(A)
f(x) = maximum_likelihood_equation(A, x)
g(x) = linear_approx_maximum_likelihood_equation(A, x)
plot(f, 0.5, 2; label="maximum likelihood equation")
plot!(g; label="Mantel-Haenszel linear approx.", ls=:dash)
hline!([0]; label="", c=:black, lw=0.5, ls=:dot)
plot!(; xlabel="common odds ratio", ylabel="score")
plot!(; xtick=0.1:0.1:3)
A = [325 273; 292 321;;; 324 363; 278 323] mh_chisq(A) = (common_odds_ratio = 1.1612448334670686, ω₀ = 1.0, p_value = 0.062406142852502865, α = 0.05, conf_int = (0.9922868538703246, 1.3589713679005824), chisq = 3.472256039291442, df = 1) mh_rbg(A) = (common_odds_ratio = 1.1610767473163213, ω₀ = 1.0, p_value = 0.06269910907099763, α = 0.05, conf_int = (0.992116222048988, 1.3588117835372746))
このグラフのMantel-Haenszel linear approx.が線形近似最尤方程式であり, その零点が共通オッズ比のMantel-Haenszel推定量になる. 最尤方程式(maximum likelihood equation)の零点が(unconditionalな場合の)最尤法で求めた共通オッズ比の推定量になる.
plot(title="P-value functions")
plot!(ω -> pvalue_chisq(A, ω), 0.7, 1.7; label="Mantel-Haenszel χ²")
plot!(ω -> pvalue_mhrbg(A, ω); label="Robins-Breslow-Greenland", ls=:dash)
plot!(; xtick=0.1:0.1:3, ytick=0:0.05:1)
vline!([1]; label="", c=:black, lw=0.5)
この場合に, 最尤法の場合のスコア検定のP値函数(上のグラフのMantel-Haenszel χ²)と共通オッズ比のMantel-Haenszel推定量の対数の分散のRobins-Breslow-Greenlandの推定量を使った正規分布近似で作ったP値函数(上のグラフのRobins-Breslow-Greenland)はほぼぴったり一致している.
TenStudies = [
215 229 311-215 306-229
38 33 59-38 51-33
161 174 293-161 293-174
76 88 164-76 163-88
103 105 129-103 133-105
65 67 120-65 125-67
81 75 113-81 110-75
48 63 160-48 159-63
22 21 60-22 62-21
56 51 137-56 140-51
]
ElevenStudies = [
TenStudies
468 480 697-468 685-480
]
11×4 Matrix{Int64}: 215 229 96 77 38 33 21 18 161 174 132 119 76 88 88 75 103 105 26 28 65 67 55 58 81 75 32 35 48 63 112 96 22 21 38 41 56 51 81 89 468 480 229 205
oddsratio.(eachrow(ElevenStudies))
11-element Vector{Float64}: 0.7530476710334789 0.987012987012987 0.8341605712295368 0.7360537190082644 1.0564102564102564 1.023066485753053 1.18125 0.6530612244897959 1.1303258145363408 1.206487533284919 0.8728165938864629
A = A10 = reshape(TenStudies', 2, 2, :) |> collect
@show A
@show mh_chisq(A)
@show mh_rbg(A)
f(x) = maximum_likelihood_equation(A, x)
g(x) = linear_approx_maximum_likelihood_equation(A, x)
plot(f, 0.5, 1.5; label="maximum likelihood equation")
plot!(g; label="Mantel-Haenszel linear approx.", ls=:dash)
hline!([0]; label="", c=:black, lw=0.5, ls=:dot)
plot!(; xlabel="common odds ratio", ylabel="score")
plot!(; xtick=0.1:0.1:3)
A = [215 96; 229 77;;; 38 21; 33 18;;; 161 132; 174 119;;; 76 88; 88 75;;; 103 26; 105 28;;; 65 55; 67 58;;; 81 32; 75 35;;; 48 112; 63 96;;; 22 38; 21 41;;; 56 81; 51 89] mh_chisq(A) = (common_odds_ratio = 0.8779119042043689, ω₀ = 1.0, p_value = 0.08661968813950241, α = 0.05, conf_int = (0.7564370173217948, 1.018894296145525), chisq = 2.936099126901802, df = 1) mh_rbg(A) = (common_odds_ratio = 0.8781804164991955, ω₀ = 1.0, p_value = 0.08757101848883686, α = 0.05, conf_int = (0.7565860364363416, 1.0193167819422062))
plot(title="P-value functions of Ten Studies")
plot!(ω -> pvalue_chisq(A10, ω), 0.7, 1.1; label="Mantel-Haenszel χ²")
plot!(ω -> pvalue_mhrbg(A10, ω); label="Robins-Breslow-Greenland", ls=:dash)
plot!(; xtick=0.1:0.05:3, ytick=0:0.05:1)
vline!([1]; label="", c=:black, lw=0.5)
plot!(; xlabel="common odds ratio", ylabel="P-value")
A = A11 = reshape(ElevenStudies', 2, 2, :) |> collect
@show A
@show mh_chisq(A)
@show mh_rbg(A)
f(x) = maximum_likelihood_equation(A, x)
g(x) = linear_approx_maximum_likelihood_equation(A, x)
plot(f, 0.5, 1.5; label="maximum likelihood equation")
plot!(g; label="Mantel-Haenszel linear approx.", ls=:dash)
hline!([0]; label="", c=:black, lw=0.5, ls=:dot)
plot!(; xlabel="common odds ratio", ylabel="score")
plot!(; xtick=0.1:0.1:3)
A = [215 96; 229 77;;; 38 21; 33 18;;; 161 132; 174 119;;; 76 88; 88 75;;; 103 26; 105 28;;; 65 55; 67 58;;; 81 32; 75 35;;; 48 112; 63 96;;; 22 38; 21 41;;; 56 81; 51 89;;; 468 229; 480 205] mh_chisq(A) = (common_odds_ratio = 0.8763754526900092, ω₀ = 1.0, p_value = 0.037876956707385036, α = 0.05, conf_int = (0.7737258422531788, 0.9926435960137233), chisq = 4.310539094380579, df = 1) mh_rbg(A) = (common_odds_ratio = 0.876565283202828, ω₀ = 1.0, p_value = 0.03831919505862805, α = 0.05, conf_int = (0.7738329501187073, 0.9929361312394169))
plot(title="P-value functions")
plot!(ω -> pvalue_chisq(A11, ω), 0.7, 1.1; label="Mantel-Haenszel χ²")
plot!(ω -> pvalue_mhrbg(A11, ω); label="Robins-Breslow-Greenland", ls=:dash)
plot!(; xtick=0.1:0.05:3, ytick=0:0.05:1)
vline!([1]; label="", c=:black, lw=0.5)
plot!(; xlabel="common odds ratio", ylabel="P-value")
plot(title="P-value functions")
plot!(ω -> pvalue_chisq(A10, ω), 0.7, 1.1; label="Ten Studies")
plot!(ω -> pvalue_chisq(A11, ω); label="Eleven Studies", ls=:dash)
plot!(; xtick=0.1:0.05:3, ytick=0:0.05:1)
vline!([1]; label="", c=:black, lw=0.5)
plot!(; xlabel="common odds ratio", ylabel="P-value")
TenStudies = [
274 276 311-274 306-276
45 39 59-45 51-39
215 222 293-215 293-222
134 139 164-134 163-139
119 127 129-119 133-127
86 92 120-86 125-92
81 75 113-81 110-75
95 92 168-95 162-92
51 54 60-54 62-54
56 51 137-56 140-51
]
ElevenStudies = [
TenStudies
468 480 698-468 687-480
]
11×4 Matrix{Int64}: 274 276 37 30 45 39 14 12 215 222 78 71 134 139 30 24 119 127 10 6 86 92 34 33 81 75 32 35 95 92 73 70 51 54 6 8 56 51 81 89 468 480 230 207
oddsratio.(eachrow(ElevenStudies))
11-element Vector{Float64}: 0.8049353701527615 0.989010989010989 0.8815546315546315 0.7712230215827338 0.5622047244094488 0.9072890025575447 1.18125 0.990172721858249 1.2592592592592593 1.206487533284919 0.8775
A = reshape(ElevenStudies', 2, 2, :) |> collect
@show A
@show mh_chisq(A)
@show mh_rbg(A)
f(x) = maximum_likelihood_equation(A, x)
g(x) = linear_approx_maximum_likelihood_equation(A, x)
plot(f, 0.5, 1.5; label="maximum likelihood equation")
plot!(g; label="Mantel-Haenszel linear approx.", ls=:dash)
hline!([0]; label="", c=:black, lw=0.5, ls=:dot)
plot!(; xlabel="common odds ratio", ylabel="score")
plot!(; xtick=0.1:0.1:3)
A = [274 37; 276 30;;; 45 14; 39 12;;; 215 78; 222 71;;; 134 30; 139 24;;; 119 10; 127 6;;; 86 34; 92 33;;; 81 32; 75 35;;; 95 73; 92 70;;; 51 6; 54 8;;; 56 81; 51 89;;; 468 230; 480 207] mh_chisq(A) = (common_odds_ratio = 0.917829906301257, ω₀ = 1.0, p_value = 0.2236611785030908, α = 0.05, conf_int = (0.7994583071772404, 1.053728254953918), chisq = 1.480724650002029, df = 1) mh_rbg(A) = (common_odds_ratio = 0.917886640779304, ω₀ = 1.0, p_value = 0.22446117848630437, α = 0.05, conf_int = (0.7993745034991397, 1.0539689239938608))
plot(title="P-value functions")
plot!(ω -> pvalue_chisq(A, ω), 0.7, 1.2; label="Mantel-Haenszel χ²")
plot!(ω -> pvalue_mhrbg(A, ω); label="Robins-Breslow-Greenland", ls=:dash)
plot!(; xtick=0.1:0.05:3, ytick=0:0.05:1)
vline!([1]; label="", c=:black, lw=0.5)
plot!(; xlabel="common odds ratio", ylabel="P-value")
data = [
10 20 50 100
150 50 50 20
50 50 50 50
200 160 150 180
]
4×4 Matrix{Int64}: 10 20 50 100 150 50 50 20 50 50 50 50 200 160 150 180
oddsratio.(eachrow(data))
4-element Vector{Float64}: 1.0 1.2 1.0 1.5
A = reshape(data', 2, 2, :) |> collect
@show A
@show mh_chisq(A)
@show mh_rbg(A)
f(x) = maximum_likelihood_equation(A, x)
g(x) = linear_approx_maximum_likelihood_equation(A, x)
plot(f, 0.5, 2.2; label="maximum likelihood equation")
plot!(g; label="Mantel-Haenszel linear approx.", ls=:dash)
hline!([0]; label="", c=:black, lw=0.5, ls=:dot)
plot!(; xlabel="common odds ratio", ylabel="score")
plot!(; xtick=0.1:0.1:3)
A = [10 50; 20 100;;; 150 50; 50 20;;; 50 50; 50 50;;; 200 150; 160 180] mh_chisq(A) = (common_odds_ratio = 1.3093105722597338, ω₀ = 1.0, p_value = 0.022249891715306905, α = 0.05, conf_int = (1.0391363562380977, 1.6497246250786697), chisq = 5.226137306969316, df = 1) mh_rbg(A) = (common_odds_ratio = 1.309886547811994, ω₀ = 1.0, p_value = 0.02242202026521213, α = 0.05, conf_int = (1.0389493572639923, 1.651478732955071))
plot(title="P-value functions")
plot!(ω -> pvalue_chisq(A, ω), 0.7, 2; label="Mantel-Haenszel χ²")
plot!(ω -> pvalue_mhrbg(A, ω); label="Robins-Breslow-Greenland", ls=:dash)
plot!(; xtick=0.1:0.1:3, ytick=0:0.05:1)
vline!([1]; label="", c=:black, lw=0.5)
plot!(; xlabel="common odds ratio", ylabel="P-value")
data = [
10 5 5 10
10 5 10 20
10 10 5 10
10 20 5 10
]
4×4 Matrix{Int64}: 10 5 5 10 10 5 10 20 10 10 5 10 10 20 5 10
oddsratio.(eachrow(data))
4-element Vector{Float64}: 4.0 4.0 2.0 1.0
A = reshape(data', 2, 2, :) |> collect
@show A
@show mh_chisq(A)
@show mh_rbg(A)
f(x) = maximum_likelihood_equation(A, x)
g(x) = linear_approx_maximum_likelihood_equation(A, x)
plot(f, 0.8, 3.6; label="maximum likelihood equation")
plot!(g; label="Mantel-Haenszel linear approx.", ls=:dash)
hline!([0]; label="", c=:black, lw=0.5, ls=:dot)
plot!(; xlabel="common odds ratio", ylabel="score")
plot!(; xtick=0.2:0.2:5)
A = [10 5; 5 10;;; 10 10; 5 20;;; 10 5; 10 10;;; 10 5; 20 10] mh_chisq(A) = (common_odds_ratio = 2.3699386690760393, ω₀ = 1.0, p_value = 0.012425118136135197, α = 0.05, conf_int = (1.2017576744162632, 4.671608629234173), chisq = 6.249174748913902, df = 1) mh_rbg(A) = (common_odds_ratio = 2.297872340425532, ω₀ = 1.0, p_value = 0.016369155956235767, α = 0.05, conf_int = (1.1649832375762472, 4.532440573031957))
このグラフのMantel-Haenszel linear approx.が線形近似最尤方程式であり, その零点が共通オッズ比のMantel-Haenszel推定量になる. 最尤方程式(maximum likelihood equation)の零点が(unconditionalな場合の)最尤法で求めた共通オッズ比の推定量になる. この場合にそれらは互いに少しずれている.
plot(title="P-value functions")
plot!(ω -> pvalue_chisq(A, ω), 0.5, 7; label="Mantel-Haenszel χ²")
plot!(ω -> pvalue_mhrbg(A, ω); label="Robins-Breslow-Greenland", ls=:dash)
plot!(; xtick=0.5:0.5:10, ytick=0:0.05:1)
vline!([1]; label="", c=:black, lw=0.5)
plot!(; xlabel="common odds ratio", ylabel="P-value")
この場合には, 最尤法の場合のスコア検定のP値函数(上のグラフのMantel-Haenszel χ²)と共通オッズ比のMantel-Haenszel推定量の対数の分散のRobins-Breslow-Greenlandの推定量を使った正規分布近似で作ったP値函数(上のグラフのRobins-Breslow-Greenland)は互いに少しずれている.
data = [
3 9 104 1059
1 3 5 86
]
2×4 Matrix{Int64}: 3 9 104 1059 1 3 5 86
oddsratio.(eachrow(data))
2-element Vector{Float64}: 3.394230769230769 5.733333333333333
A = reshape(data', 2, 2, :) |> collect
@show A
@show mh_chisq(A)
@show mh_rbg(A)
f(x) = maximum_likelihood_equation(A, x)
g(x) = linear_approx_maximum_likelihood_equation(A, x)
plot(f, 0.8, 5; label="maximum likelihood equation")
plot!(g; label="Mantel-Haenszel linear approx.", ls=:dash)
hline!([0]; label="", c=:black, lw=0.5, ls=:dot)
plot!(; xlabel="common odds ratio", ylabel="score")
plot!(; xtick=0.2:0.2:10)
A = [3 104; 9 1059;;; 1 5; 3 86] mh_chisq(A) = (common_odds_ratio = 3.7877789915115923, ω₀ = 1.0, p_value = 0.015663013464546467, α = 0.05, conf_int = (1.25569855189032, 11.471349137718786), chisq = 5.840315991419105, df = 1) mh_rbg(A) = (common_odds_ratio = 3.781172274625746, ω₀ = 1.0, p_value = 0.024415938786306244, α = 0.05, conf_int = (1.1873402722478363, 12.04142073218091))
plot(title="P-value functions")
plot!(ω -> pvalue_chisq(A, ω), 0.5, 16; label="Mantel-Haenszel χ²")
plot!(ω -> pvalue_mhrbg(A, ω); label="Robins-Breslow-Greenland", ls=:dash)
plot!(; xtick=0:30, ytick=0:0.05:1)
vline!([1]; label="", c=:black, lw=0.5)
plot!(; xlabel="common odds ratio", ylabel="P-value")