Triad sou.

Docker DesktopでローカルにOverleaf環境を構築した時のメモ

Docker Desktopのインストールが簡単にできるようになったという記事を拝見したため、Overleaf Community Editionを設定してみたメモを作りました。

Docker Desktopをインストールします

以下から最新版をダウンロードしてインストーラーを実行しました。
インストール後にWindowsからのサインアウトが必要でした。

つまづいたときは以下の記事を参考にすると良いと思います。

docker-compose.ymlをダウンロード

OverleafのGitHubにおいてあるやつをダウンロードし、適当なディレクトリに置いておきます。

Windows PowerShellで設定をしていく

Windows PowerShellを起動して、docker-compose.ymlを配置したディレクトリに移動し、`docker-compose up -d`を実行します。

cd [docker-compose.ymlを配置したディレクトリ]
docker-compose up -d

必要なものが自動で`docker pull`され、コンテナを起動してくれます。
ファイヤーウォールの許可が必要と出たので許可しました。

コンテナの設定

以下の記事の通りにやります、大変参考になりました。

ほとんど同じなのですが、若干使ったコマンドが違うため、以下に自分がやったものをメモします。
`docker ps`で起動中のコンテナの確認をします。

docker ps

起動中のコンテナが表示されます。

CONTAINER ID   IMAGE                   COMMAND                  CREATED              STATUS                        PORTS                NAMES
ec9906335c40   sharelatex/sharelatex   "/sbin/my_init"          About a minute ago   Up About a minute             0.0.0.0:80->80/tcp   sharelatex
b9902eb5b9cf   mongo:4.0               "docker-entrypoint.s…"   About a minute ago   Up About a minute (healthy)   27017/tcp            mongo
1e6d413b422b   redis:5                 "docker-entrypoint.s…"   About a minute ago   Up About a minute             6379/tcp             redis

sharelatexのコンテナID(上ではec9906335c40)を確認し、コンテナ内の`bash`を起動します。

docker exec -it ec9906335c40 bash

起動すると以下のようになるはずです。

root@ec9906335c40:/# 

TeX環境が最小構成で組まれているそうなので、全部入れます(かなり時間がかかります、自分の環境では2時間かかりました)。
タイミングが良ければ以下で出来ます。

tlmgr update --self
tlmgr install scheme-full

が、TeX Liveのバージョン関係でエラーが出る場合があります。

tlmgr: Local TeX Live (2020) is older than remote repository (2021).
Cross release updates are only supported with
  update-tlmgr-latest(.sh/.exe) --update
See https://tug.org/texlive/upgrade.html for details.

となる場合は、`update-tlmgr-latest.sh`を持ってくる必要がありました。以下にある最新のやつを取ってきた方が良いでしょう。

wget http://mirror.ctan.org/systems/texlive/tlnet/update-tlmgr-latest.sh
chmod u+x update-tlmgr-latest.sh
./update-tlmgr-latest.sh -- --upgrade
tlmgr install scheme-full


TeXがインストールできたら、コンテナから出ます。

exit

そしてTeXが全部入ったコンテナを保存します。

docker commit -m "installing all latex packages" ec9906335c40 sharelatex/sharelatex:v1

コンテナの適用

これも以下の記事の通りにやります、大変参考になりました。

Docker Desktopからコンテナを一度停止します。
f:id:triadsou:20210408001916p:plain
そして、最初に配置した`docker-compose.yml`の中身を編集します。

version: '2.2'
services:
    sharelatex:
        restart: always
        # Server Pro users:
        # image: quay.io/sharelatex/sharelatex-pro
        image: sharelatex/sharelatex:v1         # ここを保存したコンテナに変更

日本語化したらちょっと見た目に違和感があったので、私はやりませんでした。

adminの作成

これも以下の記事の通りにやります、大変参考になりました。

Docker Desktopからコンテナを再起動して(これ以降はローカルのOverleaf不使用時はDocker Desktopから停止しておくことができます)、
f:id:triadsou:20210408022246p:plain
以下を実行し、表示されたアドレスに飛んでパスワードを設定します。

docker exec sharelatex /bin/bash -c "cd /var/www/sharelatex; grunt user:create-admin --email=[your e-mail address]"

あとは設定したメールアドレスとパスワードでログインできます。
f:id:triadsou:20210408022701p:plain


以上で完成です。
f:id:triadsou:20210410080121p:plain

\RequirePackage{plautopatch}
\documentclass[uplatex,dvipdfmx]{jsarticle}

\title{サンプル文書}
\date{\today}
\author{著者}

\begin{document}
\maketitle

\begin{abstract}
日本語もテストしてみます。
\end{abstract}

\section{はじめに}
ローカルからテスト。

\[
\frac{1}{\sqrt{2\pi\sigma^2}}
\exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}
\]

\end{document}

ある仮定のもとでのKaplan–Meier推定量の正確な分布について

はじめに

この記事はあまり実用的ではないです。
調べた内容を忘れてしまうともったいないような気がしたので、一応メモしておくことにしました。

生存関数のKaplan–Meier推定量$\hat{S}(t)$は非常によく使われていて、$\sqrt{n}\{\hat{S}(t)-S(t)\} \to^d -S(t)W(t)$が成り立つことはよく知られています。
ただし、$W(t)$は$\mathrm{E}[W(t)]=0, \mathrm{Var}[W(t)]=\int_0^t \{\Pr(U>s)S(s)\}^{-1} \mathrm{d}\Lambda(s)$なるブラウン運動、$U$は打切り時間の確率変数です。
この性質から、Kaplan–Meier推定量のpoint-wiseの信頼区間は正規近似を使って構成することができます。
ただし、そのまま正規近似すると$n$が小さかったり生存関数が0や1に近いところで近似精度が悪く、変換を行った統計量の方が正規近似の精度が良いことが知られています(Borgan & Liestøl, 1990など)。
シミュレーションの結果などから、より保守的な二重対数変換($\log\{-\log S(t)\}$)と分散安定化変換っぽい逆正弦平方根変換($\mathrm{arc sin}\sqrt{S(t)}$)を利用した近似信頼区間が推奨されています。
普通はこの変換推定量を利用した近似信頼区間が用いられますが、正確な分布に関する研究もいくつかあり、その一つを調べた結果をメモしています。

Chang (1996)

この論文では、Coxの比例ハザードモデルとは別の比例ハザードモデルを考え、その仮定のもとでKaplan–Meier推定量の正確な分布を導いたものでした。

生存時間モデルの基本の仮定

$T_i, U_i$(for $i=1, \ldots, n$)をそれぞれ生存時間と打切り時間の確率変数とし、$T_i$は分布$F$にしたがい$U_i$は分布$G$にしたがいこれらがi.i.d.であるという仮定と、$T_i$と$U_i$の間の独立性を仮定します。
そして$X_i=\min(T_i, U_i)$, $\delta_i=I(T_i \leq U_i)$が観測可能である(生存時間$T_i$は常に観測できるわけではない)というモデルを考えます。
また、確率$K(t)=\Pr(X \leq t)$を考え、この経験推定量を$K_n(t)=n^{-1}\sum_{i=1}^n I(X_i \leq t)$とします。
Kaplan–Meier推定量は生存関数$\bar{F}(t)=1-F(t)$の推定量であり、
\[
\bar{F}_n(t)=\prod_{i=1}^{nK_n(t)} c_{in}^{\delta_{(i)}}I(X_{(n)} \geq t),
\] と表わすことができます。
$X_{(i)}$は$X_{(1)} < X_{(2)} < \ldots < X_{(n)}$となるように$X_i$を並べ替えたもので、$\delta_{(i)}$は$\delta_i$を同じ順番に並べ替えたもので、$c_{in}=(n-i)/(n-i+1)$とします。

比例ハザードモデル

上の設定のもとで、比例ハザードモデル
\[
\bar{G}=1-G=(\bar{F})^{\beta}, \beta>0
\] を仮定します。
この仮定のもとでは、$\delta_1,\ldots,\delta_n$と$X_1,\ldots,X_n$の間の独立性が成り立ちます(Allen, 1963)。

正確な分布

比例ハザードモデルのもとでは$\gamma=\Pr(T \leq U)=1/(1+\beta)$が成り立ちます。
$T$と$U$の独立性から、$\bar{K}(t)=1-K(t)=\bar{F}(t)\bar{G}(t)$が成り立ちます。
また、$X\leq t$を満たす人数を$A_n(t)=n K_n(t)$とし、$t$以下の死亡の人数を$D_n(t)=\sum_{i=1}^{A_n(t)}\delta_{(i)}$とします。
このとき、
\[
\Pr(A_n(t)= q)=\binom{n}{q} \{K(t)\}^q\{\bar{K}(t)\}^{n-q},
\] であり(独立に確率$K(t)$で起こる事象が$n$人中$q$人で観測される確率より)、
\[
\Pr(D_n(t)= i \mid A_n(t)= q)=\binom{q}{i} \gamma^i (1-\gamma)^{q-i},
\] が成り立ちます($\delta_1,\ldots,\delta_n$と$X_1,\ldots,X_n$の間の独立性より、独立に確率$\gamma$で起こる事象が$q$人中$i$人で観測される確率より)。
すると、Kaplan–Meier推定量が時点$t$で$x$以下となる確率$H_n^t(x)=\Pr(\bar{F}_n(t) \leq x)$($0 < x \leq 1$)は
\[
\begin{aligned}
H_n^t(x)&=
\sum_{q=0}^n \Pr(A_n(t)= q) \Pr(\bar{F}_n(t) \leq x \mid A_n(t)= q)\\&=
\sum_{q=0}^n \Pr(A_n(t)= q) \left[\sum_{i=0}^q \Pr(D_n(t)= i \mid A_n(t)= q)\Pr(\bar{F}_n(t) \leq x \mid A_n(t)= q, D_n(t)=i) \right]
\end{aligned}
\] となります(周辺分布の定義より)。
ここで、$\Pr(\bar{F}_n(t) \leq x \mid A_n(t)= q, D_n(t)=i)$は
\[
\Pr\left(\prod_{j=1}^q c_{jn}^{\delta_{(j)}} \leq x \mid A_n(t)= q, D_n(t)=i\right)=
\Pr\left(\sum_{j=1}^q \delta_{(j)} \log c_{jn} \leq \log x \mid A_n(t)= q, D_n(t)=i\right),
\] から計算ができます。
具体的には、$q$と$i$が固定されているため、$\sum_{j=1}^q \delta_{(j)}=i$($j=1,\ldots,q$)を満たすような$\delta_{(j)}$の並べ替え分布を作り、$\sum_{j=1}^q \delta_{(j)} \log c_{jn} \leq \log x$を満たす割合を計算すれば良いです。
自分は以下のコードで計算しました。

hnt <- function(t, x, n, gamma, kt) {
  ghat <- gamma
  prob <- kt^n
  for (q in 0:(n - 1)) {
    if (q == 0) {
      gx <- 1
    } else {
      gx <- 0
      bs <- e1071::bincombinations(q)
      grp <- rowSums(bs)
      bs[,1] <- bs[,1]*log((n - 1)/n)
      if (q > 1) {
        for (j in 2:q) {
          bs[,1] <- bs[,1] + bs[,j]*log((n - j)/(n - j + 1))
        }
      }
      hx <- aggregate(x = data.frame(prob = as.numeric(bs[,1] <= log(x))),
                      by = data.frame(grp), FUN = mean)
      for (i in 0:q) {
        gx <- gx + choose(q, i)*ghat^i*(1 - ghat)^(q - i)*hx[i + 1, 2]
      }
    }
    prob <- prob + gx*choose(n, q)*kt^q*(1 - kt)^(n - q)
  }
  prob
}

信頼区間

正確な分布に基づく信頼区間についてはこの論文には書かれていません。
$H_n^t(x)$を求めるためには$K(t)$と$\gamma$が必要ですが、実際にはデータから求めた推定値$K_n(t)=n^{-1}\sum_{i=1}^n I(X_i \leq t)$と$\hat{\gamma}=n^{-1}\sum_{i=1}^n \delta_i$しか得られないため、これらをプラグインしたものを使ってシミュレーションをしてみました。
また、時点$t$における両側$100(1-\alpha/2)$%信頼区間は$\alpha/2=H_n^t(x_l) < H_n^t(x) < H_n^t(x_u)=1-\alpha/2$から求められるので、これを満たす$x_l$と$x_u$を数値的に探さなくてはなりません。
以下は、下側信頼区間を計算する目的で$x_l$を二分探索する適当なコードです(計算が重くなりすぎるので、精度は適当です)。

sch <- function(prob, t, n, gamma, kt) {
  if (hnt(t, 0.5, n, gamma, kt) < prob) {
    if (hnt(t, 0.75, n, gamma, kt) < prob) {
      a <- 1
      b <- 0.75
    } else {
      a <- 0.75
      b <- 0.5
    }    
  } else {
    if (hnt(t, 0.25, n, gamma, kt) < prob) {
      a <- 0.5
      b <- 0.25
    } else {
      a <- 0.25
      b <- 0
    }
  }
  while (abs(a - b) >= 0.002) {
    if (hnt(t, round((a + b)/2, 3), n, gamma, kt) < prob) {
      b <- round((a + b)/2, 3)
    } else {
      a <- round((a + b)/2, 3)
    }
  }
  b
}

シミュレーションの結果は載せませんが、比例ハザードモデルが成り立つ条件下では、近似信頼区間よりもかなり保守的ではあるものの正確にType I errorが$\alpha$以下になっているような結果にはなりました。
計算方法も適当であり、大きい$q$についてe1071::bincombinations(q)を呼ぶとメモリが足りなくなるので、$n=25$ぐらいまでが計算の限界でした。
特に何の参考にもなりませんが、こういう方法もあるのだなということが分かって勉強になりました。

参考文献

  • Borgan Ø, Liestøl K. A note on confidence intervals and bands for the survival function based on transformations. Scandinavian Journal of Statistics 1990; 17: 35–41.
  • Chang MN. Exact distribution of the Kaplan–Meier estimator under the proportional hazards model. Statistics & Probability Letters 1996; 28: 153–157.
  • Allen WR. A note on conditional probability of failure when hazards are proportiona. Operations Research 1963; 11: 658–659.

WSLをインストールしてRStudio Serverを動かしてみた

ということで、インストールしてみた。


WSLのインストールとUbuntuの入手。
管理者権限で起動したPowerShell上でWSLをインストール。

Enable-WindowsOptionalFeature -Online -FeatureName Microsoft-Windows-Subsystem-Linux

再起動してMicrosoft StoreからUbuntu 18.04 LTSをインストール。


Ubuntuのターミナルを起動&初期設定し、Rの最新版をCRANからインストールするように設定(2020/04/25現在ではまだR-3.6.3だった、2020/04/29にR-4.0.0のバイナリが配置されたので更新)。

echo -e "\n## For R package"  | sudo tee -a /etc/apt/sources.list
echo "deb https://cloud.r-project.org/bin/linux/ubuntu $(lsb_release -cs)-cran40/" | sudo tee -a /etc/apt/sources.list
sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys E298A3A825C0D65DFD57CBB651716619E084DAB9

3行目が失敗する場合は以下もご参照のこと。


make, g++, r-baseをインストール。

sudo apt update
sudo apt install make g++ r-base


Rstudio Serverをインストール(適宜最新のファイルに置き換えると良いかも; 2020/5/22に1.3.959に更新されていた)。

sudo apt-get install gdebi-core
wget https://download2.rstudio.org/server/bionic/amd64/rstudio-server-1.3.959-amd64.deb
sudo gdebi rstudio-server-1.3.959-amd64.deb


RStudio Severを起動。

sudo rstudio-server start

ブラウザで開く。

http://localhost:8787/


WindowsのCドライブは/mnt/c/にマウントされている。

パッケージ作成環境構築

パッケージ作成環境構築のため、devtoolsなどを入れておきたかったが、WSL上ではxml2がうまくインストールできない事があるようなので、対処法をメモしておく。
ターミナル上で、libxml2-dev, libssl-dev, libcurl4-openssl-devをインストールしておく。

sudo apt install libxml2-dev libssl-dev libcurl4-openssl-dev

RStudioでdevtoolsと必要パッケージを入れる。

install.package("devtools")

devtoolsの関連パッケージのインストールに失敗する場合は、エラーメッセージ

mv: cannot move '/usr/lib/R/library/00LOCK-ps/00new/xxxx' to
'/usr/lib/R/library/xxxx': Permission denied

が出ている場合は、

Sys.setenv(R_INSTALL_STAGED = FALSE)

を設定した上で、再度インストールする。
終わったら一応設定を戻しておく。

Sys.setenv(R_INSTALL_STAGED = TRUE)

または、Ubuntuに何かのパッケージが足りない場合は、エラーメッセージを探してsudo apt install package-nameする。

追記

加賀谷先生の解説も大変参考になります。ありがとうございます。



Rをソースからビルドしてみよう

多分そのままだとエラーが出たためsource.listを編集しておく。

sudo vi /etc/apt/sources.list

ファイルを開いた後に

:%s/# deb-src/deb-src/
:wq


ビルドに必要なものを入れておく。

sudo apt update
sudo apt build-dep r-base

バージョン番号を変数に代入しておき、ソースファイルをダウンロードして展開。

export R_VERSION=4.0.0
curl -O https://cran.rstudio.com/src/base/R-4/R-${R_VERSION}.tar.gz
tar -xzvf R-${R_VERSION}.tar.gz
cd R-${R_VERSION}

ビルドとインストール。

./configure \
    --prefix=/opt/R/${R_VERSION} \
    --enable-memory-profiling \
    --enable-R-shlib \
    --with-blas \
    --with-lapack
make
sudo make install

インストールの結果確認。

/opt/R/${R_VERSION}/bin/R --version

シンボリックリンクを貼っておく。

sudo ln -s /opt/R/${R_VERSION}/bin/R /usr/local/bin/R
sudo ln -s /opt/R/${R_VERSION}/bin/Rscript /usr/local/bin/Rscript


RStudio Serverを一旦止める。

sudo rstudio-server stop

設定ファイルを編集する。

sudo vi /etc/rstudio/rserver.conf

ファイルを開いた後にiで挿入モードにし

rsession-which-r=/usr/local/lib/R

を追記、escで挿入モードを抜けて、:wqで保存して閉じる。
RStudio Serverを起動してブラウザから開く。

sudo rstudio-server start


うまくいかない時はリセットしたりして対応した。

sudo rstudio-server force-suspend-all
sudo rstudio-server stop
rm -rf /home/user/.rstudio
sudo rstudio-server start

部分集団と集団全体での割合の差に対する推測

以下のような状況を想定する。

  • あるマーカーが陽性の時の真の有病率を$p_1$、あるマーカーが陰性の時の真の有病率を$p_2$とする
  • ある集団では、陽性の$n_1$人のうち病気は$X_1$人であり、陰性の$n_2$人のうち病気は$X_2$人であった
  • ある集団全体では$n=n_1+n_2$人中$X_1+X_2$人が病気であった

一般的には、陽性の有病率$p_1$と陰性の有病率$p_2$を比較したリスク差$p_2-p_1$を考える事が多いと思われるが、ここでは陽性の有病率$p_1$と全体の有病率$p$を比較したリスク差$p-p_1$の推測について考えてみたい。

リスク差$p-p_1$に対する推測

いわゆる二項分布の仮定$X_1 \sim \mathrm{Bin}(n_1, p_1)$および$X_2 \sim \mathrm{Bin}(n_2, p_2)$のもとでは、$j=1,2$について、$\mathrm{E}[X_j]=n_jp_j$、$\mathrm{Var}[X_j]=n_jp_j(1-p_j)$である。

各有病率$p_j$の推定量は$\hat{p}_j=\frac{X_j}{n_j}$である。$\mathrm{E}[\hat{p}_j]=\frac{n_jp_j}{n_j}=p_j$であり、これは不偏推定量である。

$p$の推定量は$\hat{p}=\frac{X_1+X_2}{n}$である。陽性割合を$w=n_1/n$とおけば、$\mathrm{E}[\hat{p}]=\frac{n_1p_1+n_2p_2}{n}=wp_1+(1-w)p_2$と書ける。

そして、リスク差$p-p_1$の点推定量は、
\[
\begin{split}
\hat{p}-\hat{p}_1
&=
\frac{X_1+X_2}{n}-\frac{X_1}{n_1}
\\&=
\frac{n_2}{n}\frac{X_2}{n_2}-\frac{n_2}{n}\frac{X_1}{n_1}
\\&=
(1-w)(\hat{p}_2-\hat{p}_1),
\end{split}
\]であり、期待値は$\mathrm{E}[\hat{p}-\hat{p}_1]=(1-w)(p_2-p_1)$となる。

また、$\hat{p}-\hat{p}_1$の分散は、
\[
\begin{split}
\mathrm{Var}[\hat{p}-\hat{p}_1]
&=
\mathrm{Var}[(1-w)(\hat{p}_2-\hat{p}_1)]
\\&=
(1-w)^2(\mathrm{Var}[\hat{p}_2]+\mathrm{Var}[\hat{p}_1]),
\end{split}
\]となり、$\hat{p}-\hat{p}_1$の分散推定量は、$\widehat{\mathrm{Var}}[\hat{p}-\hat{p}_1]=(1-w)^2(\widehat{\mathrm{Var}}[\hat{p}_2]+\widehat{\mathrm{Var}}[\hat{p}_1])$である。ただし、$\mathrm{Var}[\hat{p}_j]=\frac{p_j(1-p_j)}{n_j}$、$\widehat{\mathrm{Var}}[\hat{p}_j]=\frac{\hat{p}_j(1-\hat{p}_j)}{n_j}$である。

正規近似できる事を想定すれば、$p-p_1=0$を帰無仮説とした検定統計量
\[
\begin{split}
Z
&=
\frac{(1-w)(\hat{p}_2-\hat{p}_1)}{\sqrt{(1-w)^2(\widehat{\mathrm{Var}}[\hat{p}_2]+\widehat{\mathrm{Var}}[\hat{p}_1])}}
\\&=
\frac{\hat{p}_2-\hat{p}_1}{\sqrt{\widehat{\mathrm{Var}}[\hat{p}_2]+\widehat{\mathrm{Var}}[\hat{p}_1]}},
\end{split}
\]が得られるが、これは$p_2-p_1=0$を帰無仮説としたナイーブな検定統計量に完全に一致する。

リスク差$p-p_1$の信頼区間は、
\[
(1-w)\left[(\hat{p}_2-\hat{p}_1) \pm \Phi^{-1}(1-\alpha/2) \sqrt{\widehat{\mathrm{Var}}[\hat{p}_2]+\widehat{\mathrm{Var}}[\hat{p}_1]}\right]
\]と書け、リスク差$p_2-p_1$の信頼区間を$(1-w)$倍すれば求められる事が分かる。

まとめ

リスク差$p-p_1$の点推定量区間定量はリスク差$p_2-p_1$の点推定量区間定量から計算でき、$p-p_1=0$を帰無仮説としたナイーブな検定は$p_2-p_1=0$を帰無仮説としたナイーブな検定と同一の$P$値を与える事が分かる。

補足


篠崎先生、ありがとうございました。

層内に相関のあるK×2×2分割表の共通オッズ比の推測について

Liang KY. Odds ratio inference with dependent data. Biometrika 1985; 72(3): 678–682.


と言ったしたものの確認したことがなかったのでシミュレーションしてみました。
マルコフ連鎖よりもベータ二項分布の方が簡単そうだったので、積ベータ二項分布モデルについて検討しています(原典と微妙な違いがあるかもしれません)。

シミュレーション用コード

require("VGAM")
require("reshape2")
require("foreach")
require("doParallel")
mhsim <- function(k, phi, rho, n = 5, m = 5, r = 5000) {
  p2k0 <- 0.2
  cluster <- makeCluster(detectCores() - 1)
  registerDoParallel(cluster)
  res <- foreach(i = 1:r, .combine = "rbind") %dopar% {
    tab <- array(dim = c(2, 2, k))
    p2 <- p2k0
    for (l in 1:k) {
      p2 <- p2 + 0.6/k
      p1 <- (p2*phi) / (1 - p2 + p2*phi)
      n11 <- VGAM::rbetabinom(1, n, p2, rho)
      m11 <- VGAM::rbetabinom(1, m, p1, rho)
      tab[2, 1, l] <- n11
      tab[1, 1, l] <- m11
      tab[2, 2, l] <- n - n11
      tab[1, 2, l] <- m - m11
    }
    dat1 <- reshape2::melt(tab)
    dat1 <- data.frame(x=rep(dat1$Var1 - 1, dat1$value),
                       y=rep(dat1$Var2 - 1, dat1$value),
                       k=rep(dat1$Var3, dat1$value))
    mh <- mantelhaen.test(tab)
    cmle <- mantelhaen.test(tab, exact = TRUE)
    mle <- glm(y ~ factor(x) + factor(k),
                   family = binomial(link = "logit"),
                   data = dat1)
    data.frame(mh = mh$estimate, cmle = cmle$estimate,
               mle = exp(mle$coefficients[2]))
  }
  stopCluster(cluster)
  res
}
summary(mhsim(k =  25, phi = 5, rho = 0.2))
summary(mhsim(k =  50, phi = 5, rho = 0.2))
summary(mhsim(k = 100, phi = 5, rho = 0.2))
summary(mhsim(k = 500, phi = 5, rho = 0.2, r = 100))
summary(mhsim(k = 100, phi = 5, rho = 0.0))
summary(mhsim(k = 500, phi = 5, rho = 0.0, n = 1, m = 1, r = 100))

相関の影響についての結果

オッズ比$\phi=5$、相関のパラメータ$\rho=0.2$、各層内のサンプルサイズが$N_k=M_k=5$、$K=25, 50, 100$のときを検討します。
mh:Mantel–Haenszel推定量、cmle:条件付き最尤推定量、mle:無条件の最尤推定量です。
Mantel–Haenszel推定量は$K$が大きければ$\phi=5$に近づいていきます。
条件付き最尤推定量は$K$が大きくなっても$\phi=5$から離れたままです。
無条件の最尤推定量はもともと一致性がありません。

> summary(mhsim(k =  25, phi = 5, rho = 0.2))
       mh              cmle              mle         
 Min.   : 1.086   Min.   :  1.098   Min.   :  1.109  
 1st Qu.: 3.833   1st Qu.:  4.402   1st Qu.:  5.236  
 Median : 5.167   Median :  6.059   Median :  7.497  
 Mean   : 6.018   Mean   :  7.074   Mean   :  9.313  
 3rd Qu.: 7.125   3rd Qu.:  8.459   3rd Qu.: 11.038  
 Max.   :89.333   Max.   :115.649   Max.   :231.215  
> summary(mhsim(k =  50, phi = 5, rho = 0.2))
       mh              cmle             mle        
 Min.   : 2.042   Min.   : 2.148   Min.   : 2.340  
 1st Qu.: 4.130   1st Qu.: 4.786   1st Qu.: 5.764  
 Median : 5.064   Median : 5.968   Median : 7.395  
 Mean   : 5.410   Mean   : 6.381   Mean   : 8.096  
 3rd Qu.: 6.261   3rd Qu.: 7.511   3rd Qu.: 9.606  
 Max.   :35.167   Max.   :32.358   Max.   :54.375
> summary(mhsim(k = 100, phi = 5, rho = 0.2))
       mh              cmle             mle        
 Min.   : 2.458   Min.   : 2.770   Min.   : 3.109  
 1st Qu.: 4.365   1st Qu.: 5.101   1st Qu.: 6.186  
 Median : 5.041   Median : 5.961   Median : 7.392  
 Mean   : 5.218   Mean   : 6.167   Mean   : 7.730  
 3rd Qu.: 5.880   3rd Qu.: 7.012   3rd Qu.: 8.899  
 Max.   :14.608   Max.   :16.252   Max.   :24.016  
> summary(mhsim(k = 500, phi = 5, rho = 0.2, r = 100))
       mh             cmle            mle        
 Min.   :3.928   Min.   :4.636   Min.   : 5.555  
 1st Qu.:4.757   1st Qu.:5.628   1st Qu.: 6.930  
 Median :5.061   Median :6.040   Median : 7.504  
 Mean   :5.098   Mean   :6.042   Mean   : 7.514  
 3rd Qu.:5.433   3rd Qu.:6.390   3rd Qu.: 7.995  
 Max.   :6.756   Max.   :7.992   Max.   :10.339  

おまけ:無相関の積ベータ二項分布モデル

相関がない積ベータ二項分布モデルは積二項分布モデルに等しく、Mantel–Haenszel推定量も条件付き最尤推定量も一致性があります。

> summary(mhsim(k = 100, phi = 5, rho = 0.0))
       mh             cmle            mle        
 Min.   :2.945   Min.   :3.050   Min.   : 3.465  
 1st Qu.:4.498   1st Qu.:4.495   1st Qu.: 5.363  
 Median :5.005   Median :5.007   Median : 6.059  
 Mean   :5.090   Mean   :5.075   Mean   : 6.169  
 3rd Qu.:5.597   3rd Qu.:5.573   3rd Qu.: 6.842  
 Max.   :9.277   Max.   :9.877   Max.   :13.188  

おまけ:無条件のMLE

よく知られた結果ですが、無条件のMLEは1:1マッチングのときに$\phi^2$に収束します(Breslow, 1981)。

> summary(mhsim(k = 500, phi = 5, rho = 0.0, n = 1, m = 1, r = 100))
       mh             cmle            mle       
 Min.   :3.423   Min.   :3.423   Min.   :11.72  
 1st Qu.:4.311   1st Qu.:4.310   1st Qu.:18.58  
 Median :4.872   Median :4.872   Median :23.73  
 Mean   :5.013   Mean   :5.012   Mean   :25.99  
 3rd Qu.:5.457   3rd Qu.:5.457   3rd Qu.:29.78  
 Max.   :7.966   Max.   :7.965   Max.   :63.45  
  • Breslow N. Odds ratio estimators when the data are sparse. Biometrika 1981; 68(1): 73–84.

SAS/STAT 14.3

New Procedure

14.2に続いて因果推論系のプロシジャが追加されていますね。

  • CAUSALMED Procedure: Causal mediation analysis用のプロシジャですね (The CAUSALMED Procedure)。いくつかの方法で、直接効果・間接効果を分解した推定ができるようです。

Enhancements

コンパートメントモデルやcause-specific PHモデルなど実務で使いそうなものがサポートされています。

  • FREQ Procedure: COMMONRISKDIFFオプションでoverall risk differenceの推定値、信頼区間が出力されるようになった。METHOD=SCORE2オプションが追加され、リスク差やリスク比の正確な信頼区間の計算方法が追加された (Agresti and Min 2001)。Bowker’s symmetry testのexact版が追加された。(What’s New in SAS/STAT 14.3 FREQ Procedure)。
  • NLMIXED and MCMC Procedures: CMPTMODEL statementが追加され、コンパートメントモデルの解析が簡単に?できるようになった (CMPTMODEL Statement)。
  • PHREG Procedure: Cause-specific proportional hazards modelsのサポート (PHREG Procedure MODEL Statement)。
  • TTEST Procedure: ブートストラップ法に対応 (What’s New in SAS/STAT 14.3 TTEST Procedure)。

メモ

このページ に新情報が結構まとまっています。

SAS/STAT 14.2

New Procedures

最近は因果推論関係に力をいれているのか、procedure がふたつ追加されていました What's New in SAS/STAT 14.2 - New Procedures)。

  • CAUSALTRT Procedure: 二値 or 連続アウトカム、二値の治療変数の場合の average causal effect (ATT と ATE) を推定できる (The CAUSALTRT Procedure)。
  • PSMATCH Procedure: モデルを指定して簡単にマッチングが実施できる (The PSMATCH Procedure)。Exact マッチング、Mahalanobis distance マッチング、PS マッチングができる。

Enhancements

機能追加関連も、時間依存性ROC解析やNLMIXEDのマルチスレッド化など色々とあるようでした (What's New in SAS/STAT 14.2 - Highlights of Enhancements)。

  • PHREG Procedure: 時間依存性ROC解析に対応 (SUGI 2017 の paper に解説あり; Evaluating Predictive Accuracy of Survival Models with PROC PHREG)。Concordance Statistics (Harrell's C-statistic, Uno's C-statistic)、経時AUC・各時点のROC曲線の図示。
  • NLMIXED Procedure: 変量効果を指定している場合にマルチスレッド化して処理できるように (STAT14.1までは変量効果がある場合にNTHREADSを指定しても無視されていたらしいです)。