Kaplan–Meier推定量に基づく単群の生存関数の検定のサンプルサイズ設計
この文章について
本ページは、Nagashima et al. (Pharm Stat 2021) を概説したものです、詳細な内容は論文 (10.1002/pst.2090. arXiv:2012.03355) をご参照ください。
サンプルサイズ設計のWebアプリ (年次生存割合から計算, MSTから計算) も公開しています。
はじめに
生存時間をエンドポイントとした単群の臨床試験では、Kaplan–Meier推定量の信頼区間の下限が閾値を超えるか否かを判断基準とした評価が行われる場合がある。 例えば、がんの第II相試験では奏功割合をエンドポイントとすることが標準的ではあるものの、奏功割合がoverall survivalのサロゲートエンドポイントにならないなどの様々な理由で、生存時間を主要評価項目として上述の評価が行われている。 通常、単群の臨床試験は早期の相で実施されたり、希少疾患や倫理的な側面から対照群を置けない場合など、小標本で行われるためこの点の考慮が必須となる。
Kaplan–Meier推定量が漸近正規性を持つことはよく知られており (Breslow & Crowly, 1974)、小標本下では正規近似がうまくいかず特に生存関数が$0$や$1$に近い場合に第一種の過誤の確率を制御できなくなる事も知られている (Bie et al., 1987; Borgan & Liestøl, 1990)。 また、Kaplan–Meier推定量の変換を考えることで小標本下での正規近似を改善することができることもわかっている (Bie et al., 1987; Borgan & Liestøl, 1990)。 Bie et al. (1987) やBorgan & Liestøl (1990) では、無変換 (identity) と対数変換 (log) は正規近似がうまくいかず、逆正弦平方根変換 (arcsin) か二重対数変換 (log-log) を行うことで第一種の過誤の確率の制御や信頼区間の被覆確率が改善することが示されている。 なお、小標本に対応するための別のアプローチとして正確な分布 (Chang, 1996) を用いる手法も提案されているが、並べ替え分布に基づくことから計算量が多く ($n \lt 20$程度が限度)、検定の保守性が非常に強いため、本研究では検討の対象から外した。 また、Kaplan–Meier推定量以外の統計量に基づくサンプルサイズ設計の手法 (Wu, 2015) も提案されているが、非常に適用範囲が狭いため、こちらも本研究では検討の対象から外した。
本研究の開始時点では、Kaplan–Meier推定量をベースにした単群の生存関数の検定のサンプルサイズ設計に関する文献は存在せず、SWOGが公開しているWebアプリケーション (One Sample Nonparametric Survival) を引用してサンプルサイズ設計がなされている試験が複数確認された。 また、このアプリケーションは正規近似がうまくゆかない対数変換に基づいている上に、誤ったサンプルサイズ計算式が用いられていたことがわかった。 そのため、本研究では妥当なサンプルサイズ計算式を提案したうえで、どの変換が適切なのかを検討することを目的とした。
Kaplan–Meier推定量の漸近分布とその変換推定量
$T_i$ ($i=1, \ldots, n$)は分布$F(t)$に従う生存時間を表わす確率変数、$U_i$は特定の分布に従う打ち切り時間を表わす確率変数、$X_i=\min\{T_i, U_i\}$は観測時間、$\delta_i=I(T_i \lt U_i)$は打ち切りか否かの指示関数とする。 このとき、時点$t$におけるKaplan–Meier推定量 (Kaplan & Meier, 1958) は \[ \hat{S}(t)= \prod_{s \lt t}\left\{ 1-\frac{\Delta\bar{N}(s)}{\bar{Y}(s)} \right\}, \] で定義され、その分散は \[ \hat{\sigma}^2(t)= \hat{S}(t)^2 \int_{0}^{t} \frac{1}{\bar{Y}(s)\{\bar{Y}(s)-\Delta \bar{N}(s)\}} \mathrm{d}\bar{N}(s), \] となる。 ただし、$\bar{N}(t)=\sum_{i=1}^{n}N_i(t)$、$\bar{Y}(t)=\sum_{i=1}^{n}Y_i(t)$、$\Delta X(t)=X(t)-X(t-)$、$X(t-)=\lim_{h\downarrow 0}X(t-h)$であり、$N_i(t)=I(X_i \leq t, \delta_i=1)$は$i$番目の個体のイベントを数え上げる ($0$または$1$) 計数過程、$Y_i(t)=I(X_i \geq t)$は$i$番目の個体のリスク集合過程とする (時点$t$でリスク集合に含まれていれば$1$、そうでなければ$0$をとる)。 $\Lambda(t)=-\log S(t)=\int_0^t \lambda(s) \mathrm{d}s$を累積ハザード関数、$\lambda(t)=-\frac{1}{S(t)}\frac{\mathrm{d}S(t)}{\mathrm{d}t}$をハザード関数、$\hat{\Lambda}(t)=\int_0^t \{\bar{Y}(t)\}^{-1} \mathrm{d}\bar{N}(s)$をNelson–Aalen推定量 (Nelson, 1969; Aalen, 1978) とする。 $L(t)$を$\mathrm{E}[L(t)]=0$で$\mathrm{Var}[L(t)]=\lim_{n\to\infty} \langle L_n, L_n \rangle(t)=\int_0^t \{\Pr(U>s)S(s)\}^{-1} \mathrm{d}\Lambda(s)$なるブラウン運動とする。 マルチンゲール中心極限定理から、$L_n(t)=\sqrt{n}\{\hat{\Lambda}(t)-\Lambda(t)\}=\int_{0}^{t} \sqrt{n}\{\bar{Y}(s)\}^{-1} \mathrm{d}\bar{M}(t) \to_L L(t)$となることが知られている (Fleming & Harrington, 1991)。 ただし、$\bar{M}(t)=\bar{N}(t)-\int_0^t \bar{Y}(s)\lambda(s) \mathrm{d}s$である。 Kaplan–Meier推定量とNelson–Aalen推定量は次の関係を持つことが知られている:$\sqrt{n}\{\hat{S}(t)-S(t)\}=\sqrt{n}[\exp\{-\hat{\Lambda}(t)\}-\exp\{-\Lambda(t)\}]+o_P(1)$ (Breslow & Crowley, 1974)。 汎関数デルタ法を適用すると、Kaplan–Meier推定量の漸近分布$\sqrt{n}\{\hat{S}(t)-S(t)\} \to_L -S(t)L(t)$ (Gill & Johansen, 1990)と、漸近分散 \[ \sigma^2(t)=\mathrm{Var}[\sqrt{n}\{\hat{S}(t)-S(t)\}]=S^2(t)\int_0^t \{\Pr(U>s)S(s)\}^{-1} \mathrm{d}\Lambda(s), \] を得ることができる (Andersen et al., 1993). 小標本下 (特に$S(t)$が$0$や$1$に近い場合) において、Kaplan–Meier推定量の正規近似の精度が不十分であることが知られている (Bie et al., 1987; Borgan & Liestøl, 1990). 小標本下での正規近似の精度を改善するため、複数の変換推定量が提案されている。これらの定義を以下に示す:
- 恒等変換: $g\{S(t)\}=S(t)$.
- 対数変換 (Borgan & Liestøl, 1990; $\Lambda(t)$の推定量に対応): $g\{S(t)\}=\log S(t)$.
- 二重対数変換 (Kalbfleisch & Prentice, 1980): $g\{S(t)\}=\log\{-\log S(t)\}$.
- ロジット変換 (Meeker & Escobar, 1998): $g\{S(t)\}=\log[S(t)\{1-S(t)\}^{-1}]$.
- 逆正弦平方根変換 (Thomas & Grunkemeier, 1975): $g\{S(t)\}=\arcsin \sqrt{S(t)}$.
Kalbfleisch & Prentice (1980)は二重対数変換、Thomas & Grunkemeier (1975)は逆正弦平方根変換、Borgan & Liestøl (1990)やAalen et al. (2008)はそのどちらかの利用を推奨している。 これらの研究等では、小標本下でもKaplan–Meier推定量の二重対数変換や逆正弦平方根変換に基づく信頼区間が恒等変換や対数変換に基づくものよりも、名義の被覆確率に近いことが述べられている。 また、Borgan & Liestøl (1990)では、小標本下においても二重対数変換や逆正弦平方根変換に基づく変換推定量の性質がよくなる理由についても検討されており、変換推定量の分布の歪みが小さいことが明らかにされている。 さらに、Borgan & Liestøl (1990)は時点$t$に固定した (point-wiseの) ときの変換推定量の漸近正規性についても議論しており、$0 \leq t \leq T$において$S(t)$の変換$g$が微分$g'$を持ち、$g'$が非ゼロであれば、 \[ \sqrt{n}[g\{\hat{S}(t)\} - g\{S(t)\}] \to_L g'\{S(t)\}\{-S(t)L(t)\}, \] が成り立ち、$\tau^2(t)=\mathrm{Var}(\sqrt{n}[g\{\hat{S}(t)\} - g\{S(t)\}])=g'\{S(t)\}^2 \sigma^2(t)$であることが示されている。 ここで、各変換における微分$g'\{S(t)\}=\frac{\mathrm{d} g\{S(t)\}}{\mathrm{d}S(t)}$は以下となる:
- 恒等変換: $g'\{S(t)\}=1$.
- 対数変換: $g'\{S(t)\}=\{S(t)\}^{-1}$.
- 二重対数変換: $g'\{S(t)\}=\{S(t)\log S(t)\}^{-1}$.
- ロジット変換: $g'\{S(t)\}=[S(t) \{1-S(t)\}]^{-1}$.
- 逆正弦平方根変換: $g'\{S(t)\}=[4S(t) \{1-S(t)\}]^{-1/2}$.
変換推定量の分散$\tau^2(t)$は$g'\{S(t)\}^2$をそれぞれ上の通りに置換することで得ることができる。 ただし、対数変換では$0 \lt S(t) \leq 1$の仮定が必要であり、二重対数変換、ロジット変換、逆正弦平方根変換では$0 \lt S(t) \lt 1$の仮定が必要なことに注意が必要である。
提案するサンプルサイズ計算式
提案するサンプルサイズ計算式は一般的な仮説検定のサンプルサイズ計算式を直接適用しただけのものである。 まずは、片側仮説$H_0: \epsilon \leq 0$と$H_1: \epsilon \gt 0$を考える。 ただし、$\epsilon = g\{S_1(t)\} - g\{S_0(t)\}$は変換した生存関数の差とし、$S_0(t)$と$S_1(t)$は時点$t$における ($H_0$のときの) 閾値生存確率と ($H_1$のときの) 期待値生存確率とする。 第一種の過誤の確率を$\alpha$と設定し、 \[ Z=\frac{g\{\hat{S}(t)\} - g\{S_0(t)\}}{\hat{\tau}(t)/\sqrt{n}}>z_{1-\alpha}, \] である時に$H_0$を棄却する。 これは$H_0$のもとで$n\to \infty$のときに$Z \to_L N(0, 1)$が成り立つ事に基づいている。 ただし、$\hat{\tau}(t)=g'\{\hat{S}(t)\}\hat{\sigma}(t)$であり、$z_p=\Phi^{-1}(p)$は標準正規分布の$p$における分位点である。 $H_1$のもとでは、$Z$は$n$が大きい時に$N(\epsilon\sqrt{n}/\tau_1, 1)$に近づくことから、検出力は \[ \Pr\left(Z>z_{1-\alpha} \,\middle|\, H_1 \right) \approx \Phi\left(-z_{1-\alpha}+\frac{\epsilon\sqrt{n}}{\tau_1} \right), \] で近似することができる。 ただし、$\tau_j^2= g'\{S_j(t)\}^2 \sigma_j^2(t)$ ($j=0,1$)はそれぞれ$H_0$と$H_1$における漸近分散である。 このとき、検出力$1-\beta$を達成するには、 \[ 1-\beta = \mathrm{\Phi} \left(-z_{1-\alpha}+\frac{\epsilon\sqrt{n_1}}{\tau_1} \right), \] となる必要があり、以上に基づいたサンプルサイズ計算式 \[ n_1=\left\{ \frac{\tau_1(z_{1-\alpha}+z_{1-\beta})}{\epsilon} \right\}^2, \] を提案する。 この方法は前節で示した全ての変換$g\{S(t)\}$に適用でき、小標本下で変換推定量の正規近似が良ければ性能が期待できる。
SWOG のサンプルサイズ計算式の問題点
SWOG の Web アプリケーション (Cancer Research and Biostatistics, 2002)
は無料で公開されており、これを引用してサンプルサイズ設計がなされている試験が複数存在する。
SWOG アプリの解説が複数公開されていることや、個人的な経験からも、実務上よく使われていると言っても過言ではないと思われる。
公式の解説に書いてあるとおり、SWOG の One Sample Nonparametric Survival は対数変換 $g\{S(t)\}=\log S(t)$
に基づくサンプルサイズ計算式を用いている。
しかし、上述の通り、小標本下では対数変換の性質は良くないことがわかっている (Borgan & Liestøl, 1990)。
対数変換は累積ハザード関数の推定量に対応しており、累積ハザード関数は$0$から$+\infty$の値を取り、強い非対称性を持つためである。
ソースコードを読むと、このアプリは
\[
n_2 = \left( \frac{\tau_1 z_{1-\alpha}+\tau_0 z_{1-\beta}}{\epsilon} \right)^2,
\]
というサンプルサイズ計算式に基づいていることが分かる (SWOG アプリは JavaScript で実装されているため、計算式やアルゴリズムを確認できる)。
この式を逆算すると、検出力を
\[
1-\beta=\Phi\left( -\frac{\tau_1}{\tau_0}z_{1-\alpha} + \frac{\epsilon\sqrt{n_2}}{\tau_0} \right),
\]
として計算していることになる。
これを $H_1$ のもとでの統計量の形に逆算すると、
\[
\begin{split}
\Pr\left(Z>z_{1-\alpha} \,\middle|\, H_1 \right)
&\approx
\Phi\left(-\frac{\tau_1}{\tau_0}z_{1-\alpha} + \frac{\epsilon\sqrt{n}}{\tau_0} \right)
\\&=
\Pr\left(\frac{\tau_0 W + \epsilon \sqrt{n}}{\tau_1} >z_{1-\alpha} \right),
\end{split}
\]
としていることになる。
ただし、$W \sim N(0,1)$ なる確率変数とする。
しかし、$\tau_0 \ne \tau_1$ であるため、$H_1$ のもとで $(\tau_0 W + \epsilon
\sqrt{n})/\tau_1$は$N(\epsilon\sqrt{n}/\tau_1, 1)$ に従わず $N(\epsilon\sqrt{n}/\tau_1, \tau_0^2/\tau_1^2)$
に従う。
この近似は明らかに不適切であり、サンプルサイズの計算に大きな影響を与える。
※これは個人的な想像だが、製作者がシミュレーションをしてみたら対数変換のサンプルサイズ計算式の性能が悪く、とりあえず分散を入れ替えてみたら上手くいったため、この計算式が実装されたのではないかと思っている。標準的な考え方からいえば、この式は絶対に出てこないのである。単純なタイプミスかもしれないが、SWOGの研究者はきっとシミュレーションで性能の確認はしているはずだと思う。
また、対数変換の場合には分散を逆にした方がうまくいくという点は面白く、もしかすると理論的に何か示すことができるのかもしれない。
連続で微分可能な生存関数を仮定した時のサンプルサイズの計算
この節では、単純な例をもとに実際の数値計算について説明する。
以下では、解析時点を $t$ とし、生存関数 $S$ は連続微分可能なパラメトリックな関数であり、時間に対して一様な組入れが行われ、追跡不能例がない状況を考える。
組入については、組入期間 $a$ の間に一様な組入が行われ、追跡期間$b$の間は追跡されるとし、総期間を $c=a+b$ としておく。
このとき、打ち切りの確率変数 $U$ の生存関数は、
\[
\Pr(U \gt t)=
\begin{cases}
1 & 0 \lt t \leq b\\
a^{-1}(c-t) & b \lt t \leq c\\
\end{cases},
\]
で与えられ、Kaplan–Meier 推定量の漸近分散は
\[
\sigma^2(t)=
S^2(t)\int_0^t \frac{\mathrm{d}\Lambda(s)}{\Pr(U \gt s)S(s)}=
\begin{cases}
S^2(t) \int_0^t \frac{\mathrm{d}\Lambda(s)}{S(s)} & 0 \lt t \leq b\\
S^2(t) \int_0^b \frac{\mathrm{d}\Lambda(s)}{S(s)}+S^2(t) \int_b^t \frac{a}{(c-s)S(s)}
\mathrm{d}\Lambda(s) & b \lt t \leq c\\
\end{cases},
\]
となる。
$S(t)$ が連続で微分可能であることを仮定しているため、$0 \lt t \leq b$ (すなわち $\Pr(U \gt t)=1$) であれば
\[
\begin{split}
\int_0^t & \frac{\mathrm{d}\Lambda(s)}{S(s)}
=
\int_0^t -\frac{\mathrm{d}S(s)}{\mathrm{d}s} \frac{1}{\{S(s)\}^2} \mathrm{d}s
=
\left[-\{S(s)\}^{-1} \right]_0^t+
\int_0^t -2\frac{\mathrm{d}S(s)}{\mathrm{d}s} \frac{1}{\{S(s)\}^2} \mathrm{d}s \\
\Leftrightarrow&
\int_0^t -\frac{\mathrm{d}S(s)}{\mathrm{d}s} \frac{1}{\{S(s)\}^2} \mathrm{d}s
=
-\left[-\{S(s)\}^{-1} \right]_0^t
=
-\left\{ -\frac{1}{S(t)}+\frac{1}{S(0)} \right\}=
S^{-1}(t)-1,
\end{split}
\]
が成り立つ (生存時間の定義から $S(0)=1$)。
このとき、
\[
\sigma^2(t)=
S^2(t)\{S^{-1}(t)-1\}=
S(t)\{1-S(t)\},
\]
が成り立つ。
つまり、打ち切り確率変数の生存関数 $\Pr(U \gt t)$ が$t$に依存しない場合 (例えば、$0 \lt t \leq
b$や定数になる場合なども考えられる)、漸近分散は $p=S(t)$ とした二項分布の分散に一致する。
提案法は二項分布のサンプルサイズ設計の拡張になっていると言える。
さらに、逆正弦平方根変換を用いた場合、漸近分散は $\tau^2(t)=g'\{S(t)\}^2
\sigma^2(t)=[4S(t)\{1-S(t)\}]^{-1}S(t)\{1-S(t)\}=1/4$ となり、分散安定化変換になっている事が分かる。
これは逆正弦平方根変換を用いた場合の利点と言えるだろう。
次に、$b \lt t \leq c$ の場合を考えていく。
このとき、$\Pr(U \gt t)$ は $t$ の関数 ($\Pr(U \gt t)=a^{-1}(c-t)$) になっているため、漸近分散を求めるためには数値積分を行う必要がある。
例えば、生存時間の確率変数 $T$ が指数分布に従うと仮定した場合、$S(t)=\exp(-\lambda t)$ となり、漸近分散は
\[
\sigma^2(t)=\exp(-2 \lambda t)
\left[\{\exp(\lambda b)-1\}+
\int_b^t \frac{a \lambda}{(c-s) \exp(-\lambda s)} \mathrm{d}s
\right],
\]
を計算することで得ることができる。
生存時間の確率変数に指数分布を仮定すると $S(t)=\exp(-\lambda t)$ であり、$\lambda=-\log
S(t)/t$ の関係から時間に対して一定なハザード $\lambda$ を見積もることができる。
Brookmeyer & Crowley (1982)の方法を用いれば、$\lambda=-\log(0.5)/M$ の関係から生存期間中央値 $M$ を見積もることができる。
提案法を実装したサンプルサイズ設計の Web アプリ (年次生存割合から計算, MST から計算)
では、生存時間の確率変数に指数分布を仮定し、組入れが一様で、追跡不能例がない状況を想定した実装になっている。
若干の脱落などが想定される場合は、脱落率 $r$ とした一般的な方法 ($n'= \lceil n/(1-r) \rceil$ とする) を用いれば、保守的なサンプルサイズを得ることができる
(なお、このケースは $\Pr(U \gt t)$ が $t$ に依存しない場合の拡張になっている)。
ランダムな追跡不能例がある場合やもっと複雑な場合についても、適宜設定を変更して関数 $S$ や $\Pr(U \gt t)$な どに仮定を置いて定めることで、数値積分さえできれば同様に計算可能である。
シミュレーション
論文中では第一種の過誤の確率のシミュレーションと、検出力のシミュレーションを行っている。 第一種の過誤の確率のシミュレーションは既存の研究と同様で、恒等変換と対数変換は$n=100$であっても第一種の過誤の確率を制御できないことが確認できる (解析に使用しない方よい)。
ここでは、検出力のシミュレーション結果の一部だけを簡単に紹介する。 シミュレーション条件は
- 生存時間の確率変数が$\exp(\lambda), \mathrm{Weibull}(\lambda, k = 0.5), \mathrm{Weibull}(\lambda, k = 2)$の3パターン ($\exp(\lambda)=\mathrm{Weibull}(\lambda, k = 1)$)
- 各仮説のもとでのハザードは$\lambda=\{-\log S(t)\}^{1/k}/t$で計算できる
- 解析時点$t=12$、組入期間$a=24$で組入れられる時間は一様分布に従う、追跡期間$b=6, 12$
- ランダム打切りがない場合とある場合。打切時間の確率変数はそれぞれ$\exp(4\lambda), \mathrm{Weibull}(4^{1/k}\lambda, k = 0.5), \mathrm{Weibull}(4^{1/k}\lambda, k = 2)$で約20%打切られる状況とした)
- $S_0(t)=0.1, 0.4, 0.7$、$S_1(t)-S_0(t)=0.1$
- $\alpha=0.1, 0.05$、$1-\beta=0.8, 0.9$
とした。 各組み合わせの条件下で提案法とSWOGの方法でサンプルサイズ設計し、各サンプルサイズの乱数を100万セット生成し、それぞれ対応するKaplan–Meier推定量の変換推定量に基づく検定を実施してempirical power \[ \hat{P}=\frac{\{\# \text{ of replications with confidence intervals not covering the null hypothesis } S_0(t)\}}{\{\# \text{ of replications } (1,000,000)\}}. \] を計算した。 生存時間の確率変数が$\exp(\lambda)$, $\alpha=0.05$, $1-\beta=0.8$の場合の結果を以下に示した。
打切り | $b$ | $S_0(t)$ | $S_1(t)$ | サンプルサイズ | 検出力 | ||||||||||
SWOG | 恒等 | 対数 | 二重対数 | ロジット | 逆正弦 | SWOG | 恒等 | 対数 | 二重対数 | ロジット | 逆正弦 | ||||
なし | 12 | 0.1 | 0.2 | 71 | 99 | 52 | 75 | 59 | 77 | 0.786 | 0.861 | 0.739 | 0.761 | 0.769 | 0.794 |
0.4 | 0.5 | 144 | 155 | 125 | 166 | 151 | 153 | 0.820 | 0.789 | 0.762 | 0.803 | 0.792 | 0.791 | ||
0.7 | 0.8 | 106 | 99 | 87 | 142 | 137 | 115 | 0.791 | 0.755 | 0.719 | 0.857 | 0.845 | 0.795 | ||
6 | 0.1 | 0.2 | 80 | 111 | 58 | 84 | 66 | 86 | 0.814 | 0.832 | 0.716 | 0.784 | 0.739 | 0.785 | |
0.4 | 0.5 | 158 | 170 | 136 | 181 | 165 | 167 | 0.808 | 0.801 | 0.760 | 0.815 | 0.798 | 0.799 | ||
0.7 | 0.8 | 115 | 107 | 94 | 153 | 144 | 125 | 0.818 | 0.777 | 0.755 | 0.850 | 0.838 | 0.809 | ||
あり | 12 | 0.1 | 0.2 | 98 | 129 | 67 | 97 | 77 | 100 | 0.829 | 0.832 | 0.713 | 0.782 | 0.740 | 0.785 |
0.4 | 0.5 | 161 | 171 | 137 | 183 | 166 | 169 | 0.813 | 0.802 | 0.761 | 0.818 | 0.798 | 0.801 | ||
0.7 | 0.8 | 110 | 102 | 90 | 146 | 137 | 119 | 0.822 | 0.779 | 0.762 | 0.851 | 0.839 | 0.811 | ||
6 | 0.1 | 0.2 | 111 | 145 | 76 | 109 | 87 | 113 | 0.856 | 0.861 | 0.747 | 0.812 | 0.773 | 0.817 | |
0.4 | 0.5 | 178 | 188 | 151 | 201 | 183 | 185 | 0.852 | 0.841 | 0.801 | 0.855 | 0.839 | 0.839 | ||
0.7 | 0.8 | 119 | 111 | 97 | 158 | 149 | 129 | 0.843 | 0.804 | 0.781 | 0.874 | 0.864 | 0.835 |
検出力の観点のみから結果を要約する。
- SWOG の方法は目的の検出力に概ね近い値に収めることができる
- 提案法で恒等変換を使った場合は多少目的の検出力とはずれる程度
- 対数変換を使った場合は上手くコントロールができない
- 概ねロジット変換の検出力 $\lt$ 二重対数変換の検出力の関係があり、ロジット変換はたまに検出力が過小な傾向があるため、二重対数変換の方が望ましい
- 逆正弦平方根変換は目的の検出力に一番近くコントロールできる
第一種の過誤の確率も考慮し、解析手法としての観点で結果を要約する。
- 恒等変換と対数変換を用いた方法は第一種の過誤の確率が制御できないことから、恒等変換と対数変換に基づいた方法を用いることは望ましくないといえる
- 逆正弦平方根変換を用いるか、やや保守的な二重対数変換を用いることが望ましいと考えられる
ただし、二重対数変換の方が必要サンプルサイズがやや多い傾向にあるので、この点が気になる場合もあるだろう。
※興味深いことに SWOG の方法は設定した検出力に近い値にコントロールできる。なぜ???
まとめ
本研究では Kaplan–Meier 推定量をベースにした単群の生存関数の検定の妥当なサンプルサイズ計算式を提案し、どの変換が適切なのかを検討した。
検討の結果、逆正弦平方根変換かそれよりも保守的な二重対数変換を使った方法によるサンプルサイズ設計を推奨できると考えられた。
なお、解析時はサンプルサイズ設計時に使った変換と同じ変換を用いるべきである点には注意が必要である。
そして、統計ソフトウェアによってデフォルトで使用される変換の種類は異なっており、R と S-PLUS は対数変換、SASとStataは二重対数変換が使用されるのでこの点にも注意が必要である
(この文書は2021年に書いたものなので変更されている可能性はある)。
謝辞
本研究の初期の検討を行っていた際に、浅野淳一博士および佐藤宏征博士から有益な提案と助言をいただきました。深く感謝申し上げます。 本研究はJSPS科研費 16K16014 の助成を受けたものです。
参考文献
- Nagashima K, Noma H, Sato Y, Gosho M. Sample size calculations for single-arm survival studies using transformations of the Kaplan–Meier estimator. Pharmaceutical Statistics 2021; 20(3): 499–511. DOI: 10.1002/pst.2090. [arXiv:2012.03355]
- Nagashima K. A sample size determination tool for one sample non-parametric tests for a survival proportion [Internet]. 2016 Mar 21 [cited 20XX YYY ZZ]; Available from: https://nshi.jp/en/js/onesurvyr/, https://nshi.jp/en/js/onesurvmst/.
- Breslow N, Crowley J. A large sample study of the life table and product-limit estimates under random censorship. Annals of Statistics 1974; 2(3): 437–453.
- Bie O, Borgan Ø, Liestøl K. Confidence intervals and confidence bands for the cumulative hazard rate function and their small sample properties. Scandinavian Journal of Statistics 1987; 14(3): 221–233.
- Borgan Ø, Liestøl K. A note on confidence intervals and bands for the survival function based on transformations. Scandinavian Journal of Statistics 1990; 17(1): 35–41.
- Chang MN. Exact distribution of the Kaplan–Meier estimator under the proportional hazards model. Statistics & Probability Letters 1996; 28(2): 153–157.
- Wu J. Sample size calculation for the one-sample log-rank test. Pharmaceutical Statistics 2015; 14(1): 26–33.
- Kaplan E, Meier P. Nonparametric estimation from incomplete observations. Journal of the American Statistical Association 1958; 53(282): 457–481.
- Nelson W. Hazard plotting for incomplete failure data. Journal of Quality Technology 1969; 1(1): 27–52.
- Aalen O. Nonparametric inference for a family of counting processes. Annals of Statistics 1978; 6(4): 701–726.
- Fleming TR, Harrington DP. Counting Processes and Survival Analysis. New York: Wiley, 1991.
- Gill R, Johansen S. A survey of product-integration with a view toward application in survival analysis. Annals of Statistics 1990; 18(4): 1501–1555.
- Andersen PK, Borgan Ø, Gill RD, Keiding N. Statistical Models Based on Counting Processes. New York: Springer-Verlag, 1993, 176–287, Section IV.1–3.
- Kalbfleisch J, Prentice R. The Statistical Analysis of Failure Time Data. New York: Wiley, 1980.
- Meeker W, Escobar L. Statistical Methods for Reliability Data. New York: Wiley, 1998.
- Thomas D, Grunkemeier G. Confidence interval estimation of survival probabilities for censored data. Journal of the American Statistical Association 1975; 70(352): 865–871.
- Aalen O, Borgan Ø, Gjessing H. Survival and Event History Analysis. New York: Springer, 2008.
- Cancer Research and Biostatistics. One arm survival power/sample size calculator; 2002. Accessed: March 3, 2020. https://stattools.crab.org/Calculators/oneNonParametricSurvivalColored.htm.
履歴
- 2021/8/23 公開しました