記事の目次

この記事でわかること

💡 ポイント

  • 同じデータから3つの数字が出る理由:発現割合49.6%、100人年あたり111.2件、Kaplan-Meierによる1年時点の累積発現65.7%。どれも正しく、分母が違うだけだということ
  • 発現割合の信頼区間の選び方:Wald・Wilson・Clopper-Pearsonの違いと、発現1〜2例のPTが並ぶAE集計でWaldを使ってはいけない理由
  • リスク差の区間推定:Newcombeのハイブリッドスコア法をRの自作関数で実装し、14.4%ポイント[5.7, 22.8]という形で差を語る方法
  • 曝露調整発現率(EAIR)が効く場面:実薬群のほうが早く中止して曝露が短いとき、単純な発現割合が1.41倍でも、曝露で割ると2.02倍になること
  • 再発イベントと過分散offset(log(exposure)) で率のモデルにする方法と、過分散を考慮すると標準誤差が約9%大きくなり、無視すれば差を過信すること

はじめに

臨床試験の統計解析担当者にとって、有効性の解析と安全性の解析はまったく性格の違う仕事です。有効性には事前に決めた主要評価項目があり、検定の多重性を統制し、有意水準を守って結論を出します。ところが安全性には、そのような単一の主要評価項目がありません。数百のPT(基本語)にわたる集計表が並び、そのどれもが「効果があった」と言うための検定ではなく、「何かが起きていないか」を見つけるための記述です。

だからこそ、安全性の数字は作り方をひとつ間違えると、そのまま誤った印象を与えます。たとえば「実薬群のAE発現割合は49.6%、プラセボ群は35.2%でした」という一文は、それ自体は正しくても、実薬群の平均曝露が39.8週、プラセボ群が45.0週だったという事実を添えなければ、リスクの大きさを過小に見せてしまいます。観察されていた時間が短い群のほうが、単純な割合の上では安全に見えるからです。

本記事では、第III相・プラセボ対照・計画曝露52週の試験を模した仮想データ(各群250例)を使って、有害事象データの集計と区間推定をRで一通り実装します。扱うのは、発現割合とその信頼区間(Wald・Wilson・Clopper-Pearson)、リスク差の信頼区間(WaldとNewcombe)、1件も観察されなかったときの上限、曝露調整発現率(EAIR)とポアソン率比、初回発現までの時間のKaplan-Meier推定、そして再発イベントに対する負の二項回帰です。すべてのコードはそのまま実行でき、掲載している出力はローカルのR 4.5.1で実際に得た値です。

想定読者は、製薬企業やCROで安全性の集計・解析を担当されている生物統計担当者と統計解析プログラマ、そして安全性情報に関わる方です。統計検定準1級レベルの学習者にとっても、二項割合の区間推定・ポアソン過程・過分散といったテーマが実務の文脈でどう使われるかを見る材料になります。

前提として、有効性の解析対象集団の考え方や欠測データの扱いについては【完全理解】ICH E9「臨床試験の統計的原則」と補遺(Estimand)を徹底解説を、治験中の安全性情報の報告義務についてはICH E2Aを読み解く:治験中に得られる安全性情報の取り扱いとは?をあわせてご覧いただくと、本記事の位置づけがはっきりします。

有害事象データの構造と発現割合の基本

有害事象(adverse event、AE)の解析は、有効性の解析とはデータの形からして違います。有効性の主要変数が「1被験者につき1つの値」であるのに対し、AEは1人の被験者が試験期間中に何度でも、何種類でも起こしうるイベントです。しかも被験者ごとに観察されていた期間(曝露期間)が違います。ここでは、実務で手元に来るAEデータの構造を確認したうえで、本記事共通の仮想データをシミュレーションで作り、最初の集計表を読みます。

AEデータは「1AE1行」で来て、「1被験者1行」で集計される

CDISCのADaMでAE解析に使うデータセットはADAEです。ADAEはOccurrence Data Structure(OCCDS)に基づく構造で、1被験者・1有害事象につき1行を持ちます。5件のAEを報告した被験者は5行占めます。ここにMedDRAでコーディングされたPT・SOC、重篤度、発現日、因果関係の判定といった変数が並びます。

一方、安全性集計表の1行目に載る「AEが1件以上発現した被験者数と割合」は、被験者を単位とした数字です。したがって解析の実体は、ADAEを被験者単位に集約する作業から始まります。集約後のデータセットが持つべき最小限の情報は次の3つです。

  • 曝露期間(exposure):その被験者が治験薬に曝露されていた時間。中止した被験者は中止時点まで。
  • 発現件数(n_ae):その被験者が起こしたAEの延べ件数。0件もあり得る。
  • 発現フラグ(any_ae):1件以上発現したか否かの0/1。

本記事のシミュレーションデータは、この「集約後」の姿を最初から作ります。ADaMの初回発現フラグ(AOCC系の変数)は、集計表の被験者数がADAEのどの行に対応するかをたどれるようにするためのものです。

💡 ポイント
AE集計の分子は原則として「被験者数」であり「件数」ではありません。同じ被験者が同じPTのAEを3回起こしても、そのPTの発現例数は1です。SASでもRでも、被験者IDのユニーク数で数える処理(length(unique(subjid)) 相当)が入っているかどうかは、レビューで最初に見るべき点です。

TEAEと安全性解析対象集団

集計の対象を決める定義が2つあります。ひとつはイベント側の定義、もうひとつは被験者側の定義です。

イベント側の定義がTEAE(treatment-emergent adverse event、治験薬投与下で発現した有害事象)です。ICH E9の用語集は treatment emergent を「投与前には存在せず投与中に発現した事象、または投与前の状態と比べて悪化した事象」と定義しています。実装上は「AE発現日が初回投与日以降、かつ最終投与日+一定の追跡期間以内」という日付ルールになり、投与前から続く事象は重症度が投与開始後に上がった時点で拾います。この境界の切り方(最終投与後何日まで拾うか、発現日が欠測のAEをどう扱うか)はSAPで明示すべき事項です。

被験者側の定義が安全性解析対象集団です。ICH E9の6.3節は、全般的な安全性・忍容性の評価対象は通常「治験薬を少なくとも1回投与された被験者」の集合として定義される、としています。有効性のITT/FASが無作為化された被験者を割り付けどおりに扱うのに対し、安全性の集計は実際に受けた治療で被験者を分類するのが実務の標準です。割り付けは実薬だが誤って対照薬しか投与されなかった被験者は、安全性では対照群に数えます。AEを起こしうるのは実際に体内に入った薬だからです。解析対象集団の考え方はICH E9「臨床試験の統計的原則」と補遺の解説記事で扱っています。

なおICH E9は、因果関係の判定によらずすべてのAEを報告することも求めています。個別症例の報告義務や重篤性の判断基準はICH E2A(治験中に得られる安全性情報の取扱い)の記事、長期の後期開発試験で収集する安全性データを合理的に絞る考え方はICH E19(選択的な安全性データ収集)の記事にまとめてあります。

集計の単位はPTとSOC

AEの用語はMedDRA(ICH M1)でコーディングされます。MedDRAは上からSOC(器官別大分類)、HLGT(高位グループ語)、HLT(高位語)、PT(基本語)、LLT(下層語)の5階層で、SOCは27あります。日本語版のMedDRA/Jは英語版と同期して維持されています。

ICH E9の6.2節が述べているとおり、AEは通常PTの水準で集計し、同じSOCに属するPTをまとめて記述的に提示します。安全性集計表が「SOCごとにPTがぶら下がる」形になっているのはこのためです。ここではデータ全体を1つのまとまりとして扱います。

共通データを作る

以下のコードで、第III相・2群並行群間・計画曝露52週の試験を模した仮想データを作ります。各群250例、合計500例です。

set.seed(4869)

n_per_arm <- 250
max_follow <- 1.0          # planned exposure: 1 year (52 weeks)

sim_arm <- function(n, rate_ae, rate_drop, arm) {
  frailty  <- rgamma(n, shape = 2, rate = 2)          # between-subject variability
  drop     <- rexp(n, rate = rate_drop)
  exposure <- pmin(drop, max_follow)
  lambda   <- rate_ae * frailty
  n_ae     <- rpois(n, lambda * exposure)
  t_first  <- ifelse(n_ae > 0, exposure * runif(n)^(1 / pmax(n_ae, 1)), NA)
  data.frame(
    subjid   = paste0(arm, seq_len(n)),
    arm      = arm,
    exposure = exposure,
    n_ae     = n_ae,
    any_ae   = as.integer(n_ae > 0),
    t_first  = ifelse(n_ae > 0, t_first, exposure),
    event    = as.integer(n_ae > 0)
  )
}

ae <- rbind(
  sim_arm(n_per_arm, rate_ae = 1.10, rate_drop = 0.55, arm = "Active"),
  sim_arm(n_per_arm, rate_ae = 0.60, rate_drop = 0.28, arm = "Placebo")
)
ae$arm <- factor(ae$arm, levels = c("Placebo", "Active"))

summ <- data.frame(
  arm        = levels(ae$arm),
  N          = as.vector(table(ae$arm)),
  n_with_ae  = as.vector(tapply(ae$any_ae, ae$arm, sum)),
  events     = as.vector(tapply(ae$n_ae,   ae$arm, sum)),
  PY         = round(as.vector(tapply(ae$exposure, ae$arm, sum)), 1)
)
summ$pct  <- round(100 * summ$n_with_ae / summ$N, 1)
summ$EAIR <- round(100 * summ$events / as.vector(tapply(ae$exposure, ae$arm, sum)), 1)
print(summ)

生成の仕組みは4段階です。

第1に、rgamma(n, shape = 2, rate = 2) で被験者ごとの「起こしやすさ」を1つ引きます。このガンマ分布は平均1・分散0.5で、AEを起こしやすい人と起こしにくい人が混在する状況を作ります。これはfrailty(脆弱性)と呼ばれる被験者間のばらつきで、件数データに過分散をもたらす源です。

第2に、rexp(n, rate = rate_drop) で中止までの時間を指数分布から引きます。中止率は実薬0.55/プラセボ0.28で、実薬群のほうが早く中止する人が多い設定です。第3に、pmin(drop, max_follow) で計画曝露の1年を上限に打ち切ります。1年より前に中止すればその時点まで、しなければちょうど1年が曝露期間です。exposure の単位は年で、1.0が52週にあたります。

第4に、AE件数を平均 \( \lambda_i\times\text{exposure}_i \) のポアソン分布から発生させます(\( \lambda_i \) は群ごとの基準率にfrailtyを掛けた値)。曝露が長い被験者ほど件数の期待値が大きくなるという、実際のAEデータが持つ構造をここで作り込んでいます。1件以上発現した被験者には初回発現時刻 t_first を、発現しなかった被験者には曝露終了時点を打ち切りとして入れています。件数データとポアソン分布の関係はポアソン回帰の記事で詳しく扱っています。

出力は次のとおりです。

      arm   N n_with_ae events    PY  pct  EAIR
1 Placebo 250        88    119 216.2 35.2  55.1
2  Active 250       124    213 191.5 49.6 111.2

3つの数字が示していること

この6列の表には、性質の違う3種類の数字が並んでいます。

1つ目は発現割合pct)です。1件以上のAEが発現した被験者は実薬124例で49.6%、プラセボ88例で35.2%。分母はどちらも被験者数250で、実薬群のほうが高い結果です。

2つ目は延べ件数events)です。実薬213件、プラセボ119件。発現した被験者数の差(124対88)に比べて件数の差はさらに開いており、実薬群では同じ被験者が複数回発現しているケースが多いことがわかります。

3つ目は曝露人年PY)です。実薬191.5人年、プラセボ216.2人年。ここで表の様子が変わります。実薬群のほうがイベントは多いのに、観察されていた時間の総量は少ないのです。原因は中止で、実薬群では早く中止する被験者が多く、その分だけAEを起こす機会のある時間が短くなっています。

実薬群は「より短い観察時間の中で、より多くのイベントを起こした」ことになります。にもかかわらず、被験者数を分母にした49.6%と35.2%という発現割合は、この観察時間の差をまったく考慮していません。

⚠️ 注意
延べ件数を被験者数で割った値を「発現割合」と呼んではいけません。実薬群なら213/250=85.2%になりますが、これは割合ではありません。同じ被験者の複数回発現を数えているため、値は容易に100%を超えます。AEの数字を出すときは、分子が被験者数なのか件数なのか、分母が人数なのか人時間なのかを、必ずセットで言えるようにしてください。

用語を分母で区別する

混同されやすい3つの用語を、分子と分母の観点で整理します。

用語分子分母値域と読み方
発現割合
incidence proportion
1件以上発現した被験者数(124 / 88)被験者数(250 / 250)0〜1(0〜100%)。「試験期間中に1回でも起こした人の割合」。観察期間の長短は反映されない
発現率
incidence rate
イベント数または初回発現の被験者数(213 / 119)曝露人時間(191.5 / 216.2人年)0以上、上限なし。単位は「100人年あたり何件」。時間あたりの起こりやすさ
延べ件数
number of events
―(割り算をしない生の数)0以上の整数。単独では群間比較に使えない。分母とセットで初めて意味を持つ

発現割合は、被験者 \( n \) 人のうち \( x \) 人が1件以上発現したとき

\[ \hat{p}=\frac{x}{n} \]

で推定します。実薬群なら \( \hat{p}=124/250=0.496 \)。分子の \( x \) はあくまで人数であって件数ではありません。一方の発現率は、総イベント数 \( d \) と総曝露時間 \( T \) から

\[ \hat{\lambda}=\frac{d}{T} \]

で推定します。実薬群は総イベント数213件を総曝露191.5人年で割って、100人年あたり111.2件(出力の EAIR 列)です。発現割合は無次元、発現率は「時間あたり」の量で、単位からして別物です。

ICH E9の6.3節は、この点を明示的に扱っています。AEの発現状況は通常、リスクに曝された被験者数に対する発現被験者数の割合として表されます。ただし分母に曝露被験者数を取るか曝露量(人年)を取るかは状況によるため、その定義をプロトコルに書いておくことが重要であり、長期投与が計画され中止や死亡が相当数見込まれる場合には特に重要だ、というのが同節の指摘です。そのような状況では過小評価を避けるために生存時間解析の手法を検討し、累積の発現率を算出すべきだとも述べられています。本記事の設定はまさにこの「相当数の中止が見込まれる長期投与試験」であり、手元にあるのは被験者数250という同じ分母を共有しながら曝露人年は191.5と216.2で1割以上違う2つの群です。

発現割合とリスク差の信頼区間をRで求める

AE集計表に並ぶ「49.6%」「35.2%」という数字は、あくまで500例という有限のサンプルから得られた推定値です。同じ試験をもう一度やれば違う数字が出ます。その揺らぎの幅を示さないまま割合だけを並べると、読み手は点推定を確定値のように受け取ってしまいます。ここでは二項割合の信頼区間を3種類の方法で計算し、続いて群間差(リスク差)の区間、最後に「1件も発現しなかった」ときに何が言えるかを、Rの実装とともに整理します。

二項割合の信頼区間 ― Wald・Wilson・Clopper-Pearson

\( n \) 例中 \( x \) 例でAEが発現したとき、発現割合の推定値は \( \hat{p}=x/n \) です。この \( \hat{p} \) の95%信頼区間として実務で使われる方法は、大きく3つに分かれます。

もっとも素朴なのがWald区間で、二項分布を正規分布で近似し、\( \hat{p} \) を中心に標準誤差の1.96倍を足し引きします。

\[ \hat{p}\pm z_{1-\alpha/2}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}} \]

ここで \( z_{1-\alpha/2} \) は標準正規分布の上側 \( \alpha/2 \) 点(95%なら1.96)です。式が単純で手計算もできますが、標準誤差の中の \( \hat{p} \) を真値の代わりに使っている点が弱点で、\( \hat{p} \) が0や1に近づくと標準誤差が過小になり、区間が短くなりすぎます。

Wilsonのスコア区間は、標準誤差に推定値ではなく仮説上の値を使い、\( |\hat{p}-p|/\sqrt{p(1-p)/n}=z \) という方程式を \( p \) について解いて得られます。

\[ \frac{\hat{p}+\dfrac{z^{2}}{2n}\pm z\sqrt{\dfrac{\hat{p}(1-\hat{p})}{n}+\dfrac{z^{2}}{4n^{2}}}}{1+\dfrac{z^{2}}{n}} \]

分子の \( \hat{p}+z^{2}/(2n) \) は、点推定を0.5の方向へわずかに引き寄せた「連続修正つきの中心」にあたります。分母の \( 1+z^{2}/n \) で全体を縮めることで、区間は必ず0から1の内側に収まります。\( n \) が大きくなれば \( z^{2}/n \) は0に近づき、Wald区間と一致していきます。

Clopper-Pearson区間は正規近似を使わず、二項分布の裾確率をそのまま釣り合わせる厳密(exact)な区間です。ベータ分布の分位点を使って

\[ L=B(\alpha/2;\,x,\,n-x+1),\qquad U=B(1-\alpha/2;\,x+1,\,n-x) \]

と書けます。\( B(q;a,b) \) はパラメータ \( a,b \) のベータ分布の \( q \) 分位点です。Rの binom.test() が返す区間がこれで、被覆確率が名目の95%を下回らないことが保証される代わりに、やや広めになります。

3つをRで実装して、実薬群・プラセボ群それぞれに当てはめます。

ci_wald <- function(x, n) {
  p  <- x / n
  se <- sqrt(p * (1 - p) / n)
  c(est = p, lower = p - 1.96 * se, upper = p + 1.96 * se)
}
ci_wilson <- function(x, n, conf = 0.95) {
  z <- qnorm(1 - (1 - conf) / 2)
  p <- x / n
  d <- 1 + z^2 / n
  c(est   = p,
    lower = (p + z^2 / (2 * n) - z * sqrt(p * (1 - p) / n + z^2 / (4 * n^2))) / d,
    upper = (p + z^2 / (2 * n) + z * sqrt(p * (1 - p) / n + z^2 / (4 * n^2))) / d)
}
ci_cp <- function(x, n) {
  bt <- binom.test(x, n)
  c(est = x / n, lower = bt$conf.int[1], upper = bt$conf.int[2])
}

x1 <- summ$n_with_ae[summ$arm == "Active"];  n1 <- summ$N[summ$arm == "Active"]
x0 <- summ$n_with_ae[summ$arm == "Placebo"]; n0 <- summ$N[summ$arm == "Placebo"]

tab_ci <- rbind(
  Active_Wald    = ci_wald(x1, n1),
  Active_Wilson  = ci_wilson(x1, n1),
  Active_CP      = ci_cp(x1, n1),
  Placebo_Wald   = ci_wald(x0, n0),
  Placebo_Wilson = ci_wilson(x0, n0),
  Placebo_CP     = ci_cp(x0, n0)
)
print(round(100 * tab_ci, 1))
                est lower upper
Active_Wald    49.6  43.4  55.8
Active_Wilson  49.6  43.5  55.8
Active_CP      49.6  43.2  56.0
Placebo_Wald   35.2  29.3  41.1
Placebo_Wilson 35.2  29.5  41.3
Placebo_CP     35.2  29.3  41.5

実薬群の下限はWald 43.4%、Wilson 43.5%、Clopper-Pearson 43.2%で、幅にして0.3ポイントしか違いません。上限も55.8%・55.8%・56.0%で、差は0.2ポイントです。プラセボ群も同様で、3手法の食い違いはどこを取っても0.5ポイント以内に収まっています。理由は明快で、発現割合が35〜50%と0.5に近く、しかも各群250例と分母が大きいからです。この条件では正規近似が十分よく効き、どの方法を選んでも結論は変わりません。

問題は、AE集計表の大半がこの条件を満たさないことです。安全性の要約表は「全体のAE発現割合」1行で終わらず、その下にMedDRAのSOC・PT別の行が数百行並びます。そこに現れるのは250例中2例、3例、1例といった低頻度の数字です。この領域でWaldがどうなるか、式に代入して確かめてみます。250例中2例(\( \hat{p}=0.008 \))なら標準誤差は \( \sqrt{0.008\times0.992/250}=0.00563 \) で、その1.96倍は0.01104。したがってWald区間は−0.30%から1.90%となり、下限が負の値に飛び出します。割合が負になることはあり得ないので、この区間は定義からして破綻しています。同じデータにWilsonを当てれば0.22%から2.87%、Clopper-Pearsonなら0.10%から2.86%で、いずれも0以上に収まり、上限も正しく2.9%近くまで伸びています。Waldは上限まで1.90%と過小に見積もっており、「めったに起きない事象」のリスクを実際より小さく報告してしまいます。

⚠️ 注意
「今回の例では3手法が一致したのでWaldでよい」と読まないでください。一致したのは発現割合が0.5に近く分母が大きいという恵まれた条件だったからです。Wald区間は \( \hat{p} \) が0や1に近い領域で被覆確率が名目の95%を大きく下回ることが知られており、区間が0や1をはみ出すことも珍しくありません。AE集計の主戦場はまさにその低頻度領域です。SAPには WilsonまたはClopper-Pearsonを既定 と明記し、PT別集計まで同じ方法で通すのが安全です。
方法考え方特徴推奨場面
Wald正規近似(標準誤差に推定値を代入)計算が簡単。低頻度で区間が0を割り込み、被覆確率が名目を下回る割合が0.5付近かつ分母が大きい場合の参考値。既定にはしない
Wilson(スコア)スコア統計量を \( p \) について解く必ず0〜1に収まる。平均的な被覆が名目に近く、幅も無駄に広くないAE発現割合の既定。PT別の低頻度集計を含め全表で統一しやすい
Clopper-Pearson(厳密)二項分布の裾確率をベータ分位点で釣り合わせる被覆確率が95%を下回らない保証。やや保守的で広め重篤なAEや0件・極少数例の上限提示など、保守的に述べたい場面

リスク差はNewcombeのハイブリッドスコア法で

群ごとの区間が出たら、次は群間差です。リスク差は \( \hat{p}_1-\hat{p}_0 \) で、単位は「%ポイント」です(比ではないので「%」と書かないよう、CSRの表記規約でも区別されます)。Wald型の区間は2群の分散を足し合わせて

\[ (\hat{p}_1-\hat{p}_0)\pm z\sqrt{\frac{\hat{p}_1(1-\hat{p}_1)}{n_1}+\frac{\hat{p}_0(1-\hat{p}_0)}{n_0}} \]

とします。単群のWaldと同じ弱点をそのまま引き継ぐため、低頻度では−100%から100%の範囲を平気ではみ出します。

これを避けるのがNewcombeのハイブリッドスコア法です(Newcombe RG. Statistics in Medicine. 1998;17(8):873-890)。まず各群のWilson区間 \( (l_1,u_1) \) と \( (l_0,u_0) \) を求め、リスク差の下限には「実薬が下限まで低く、プラセボが上限まで高い」という最も不利な組み合わせから来る変動を、上限にはその逆を割り当てます。

\[ L=(\hat{p}_1-\hat{p}_0)-\sqrt{(\hat{p}_1-l_1)^{2}+(u_0-\hat{p}_0)^{2}} \]

\[ U=(\hat{p}_1-\hat{p}_0)+\sqrt{(u_1-\hat{p}_1)^{2}+(\hat{p}_0-l_0)^{2}} \]

各群のWilson区間が非対称なぶん、リスク差の区間も点推定に対して非対称になります。Newcombeの論文は11の方法を比較したもので、この方法は区間が−100〜100%ポイントの外へ飛び出すといった破綻を起こさず、しかも専用のプログラムを必要としない簡便さを併せ持つ点が評価されています。Rでの実装は次のとおりで、ci_wilson() をそのまま部品として呼び出しています。

rd_wald <- function(x1, n1, x0, n0) {
  p1 <- x1 / n1; p0 <- x0 / n0
  se <- sqrt(p1 * (1 - p1) / n1 + p0 * (1 - p0) / n0)
  c(diff = p1 - p0, lower = p1 - p0 - 1.96 * se, upper = p1 - p0 + 1.96 * se)
}
rd_newcombe <- function(x1, n1, x0, n0) {
  w1 <- unname(ci_wilson(x1, n1)); w0 <- unname(ci_wilson(x0, n0))
  p1 <- x1 / n1; p0 <- x0 / n0
  c(diff  = p1 - p0,
    lower = (p1 - p0) - sqrt((p1 - w1[2])^2 + (w0[3] - p0)^2),
    upper = (p1 - p0) + sqrt((w1[3] - p1)^2 + (p0 - w0[2])^2))
}
print(round(100 * rbind(Wald = rd_wald(x1, n1, x0, n0),
                        Newcombe = rd_newcombe(x1, n1, x0, n0)), 1))

rr <- (x1 / n1) / (x0 / n0)
tab2 <- matrix(c(x1, n1 - x1, x0, n0 - x0), nrow = 2, byrow = TRUE)
ft <- fisher.test(tab2)
cat("risk ratio =", round(rr, 3), "\n")
cat("Fisher p =", round(ft$p.value, 4), " OR =", round(ft$estimate, 3), "\n")
         diff lower upper
Wald     14.4   5.8  23.0
Newcombe 14.4   5.7  22.8
risk ratio = 1.409 
Fisher p = 0.0015  OR = 1.81 

リスク差は14.4%ポイント、Newcombeの95%信頼区間は5.7から22.8%ポイントです。Wald型は5.8から23.0%ポイントで、ここでも両者はほぼ重なりますが、Newcombeの区間はわずかに下側へずれています。各群のWilson区間の非対称性が反映された結果です。リスク比は1.409、オッズ比は1.81で、割合が3〜5割と高いためオッズ比がリスク比より大きく出ています。オッズ比は発現割合が高い事象で効果を誇張して見せるので、AE集計ではリスク差とリスク比を主役に据えます。

読み方は「実薬群の発現割合はプラセボ群より5.7ポイントから22.8ポイント高い範囲にある」です。250例ずつの設計では、リスク差の精度はせいぜい±8ポイント程度しかありません。この幅の広さこそが、安全性データを検定でなく区間で語るべき理由です。

💡 ポイント
Fisherの正確検定のp値は0.0015でしたが、この記事でp値に触れるのはここ一度だけにします。ICH E9は第6.4節で「多くの試験において安全性および忍容性の含意は、記述的な統計手法をデータに適用し、解釈の助けとなる場合に信頼区間の計算で補うことによって最もよく扱われる」と述べています。さらに、群間に臨床的に許容できない差が無いことを確かめようとする状況では、仮説検定より信頼区間の使用が好ましいとも明言しています。区間の幅を示せば、発現頻度が低いときに生じる推定精度の乏しさがそのまま見えるからです。ICH E9「臨床試験の統計的原則」と補遺の全体像はこちらで整理しています。

比率の差そのものの理論的な扱いは比率の差の検定と信頼区間について徹底解説で、fisher.test() の背後にある超幾何分布の考え方はFisherの正確検定で扱っています。区間とp値が同じ情報の別の見せ方であることは信頼区間とp値の関係を図解で理解するにまとめました。

1件も観察されなかったときに何が言えるか

安全性で最も誤読されるのが、AE集計表の「0(0.0%)」という行です。0%という表示は「起こらない」ではなく「今回の標本では観察されなかった」に過ぎません。ではその上限はどこにあるのか。分子が0のときのClopper-Pearson上限と、その近似である3/nを並べてみます。

zero_event <- t(sapply(c(50, 100, 250, 1000), function(n) {
  c(n = n,
    exact_upper_pct = round(100 * binom.test(0, n)$conf.int[2], 2),
    rule_of_three_pct = round(100 * 3 / n, 2))
}))
print(zero_event)
        n exact_upper_pct rule_of_three_pct
[1,]   50            7.11               6.0
[2,]  100            3.62               3.0
[3,]  250            1.46               1.2
[4,] 1000            0.37               0.3

今回の試験と同じ250例で1件も発現がなくても、真の発現割合の95%信頼区間の上限は1.46%です。これは約68人に1人という水準で、決して「まれだから無視してよい」と言える大きさではありません。50例なら上限は7.11%、1,000例まで増やしてようやく0.37%です。第III相で0件だったからといって、市販後に数万人規模で使われれば相当数の症例が出うる、という直感がこの表から得られます。

3/nという近似はrule of threeと呼ばれ、Hanley JA, Lippman-Hand A. JAMA. 1983;249(13):1743-1745で広く知られるようになった経験則です。発現確率 \( p \) で \( n \) 例すべてに事象が起きない確率は \( (1-p)^{n} \) で、これが0.05になる \( p \) を \( \exp(-np)=0.05 \) の近似から解くと \( p\approx3/n \) が得られます。上の表では6.0・3.0・1.2・0.3と厳密上限(7.11・3.62・1.46・0.37)よりやや小さめに出ますが、「n=250で0件なら上限はおよそ1.2%、実際はもう少し上」と暗算で当たりを付ける道具としては十分です。

この「0件の上限」は、稀だが重篤なリスクの議論、CTD 2.7.4での曝露規模の記述、DSURでの累積曝露の提示で必ず要求される観点です。観察されなかった事象について語れるのは「上限がどこにあるか」だけであり、それは分母の大きさで決まります。AEの数字を報告するときは、割合と区間に加えて分母である症例数と曝露規模を必ずセットで示してください。

曝露調整発現率(EAIR)で観察期間の違いを吸収する

発現割合の分母は「人」です。1週間で中止した被験者も、52週を完走した被験者も、同じ「1人」として数えられます。短期の試験であれば実害はほとんどありません。しかし今回のように計画曝露が52週あり、しかも中止の起こり方が群間で違う試験では、この分母がそのまま結論を歪めます。有害事象は、薬を飲んでいる時間が長いほど起こる機会が増えるからです。ここで登場するのが、分母を「人」から「人時間」に置き換える曝露調整発現率(EAIR: exposure-adjusted incidence rate)です。

分母を人から人時間に変える

EAIRは、分子にイベント(または発現した被験者)の数、分母に総曝露人時間(person-time)を置いた指標です。被験者 \( i \) の曝露期間を \( T_i \)(年)、対象集団の人数を \( N \)、イベント数を \( e \) とすると、100人年あたりのEAIRは次のように書けます。

\[ \mathrm{EAIR}=\frac{e}{\sum_{i=1}^{N}T_i}\times 100 \]

分母の \( \sum T_i \) が総曝露人年(person-years, PY)です。ここで単位の明示は必須です。「100人年あたり」とは、100人を1年ずつ観察したときに何件のイベントが期待されるか、という意味です。同じ人時間は「1人を100年観察する」「200人を半年ずつ観察する」でも作れますから、人時間は人数と期間を掛け合わせて均したものだと理解してください。ジャーナルやCSRでは「per 100 patient-years」「患者年あたり」といった表記が混在しますが、100倍しているのか、1人年あたりの生の値なのかで数字が2桁変わります。表のヘッダーに単位を書き忘れた集計表は、それだけで読めません。

人時間(person-time)の作り方には二つの流儀がある

人時間は、被験者ごとの曝露期間を足し上げて作ります。難しいのは、どこまでを「曝露期間」と呼ぶかです。ここには実務上はっきり分かれた二つの流儀があります。

ひとつは、分子を「初回発現した被験者の数」に限り、分母を「初回発現までの時間」(発現しなかった被験者は曝露終了まで)とする定義です。同じ被験者が同じ事象を3回起こしても分子は1です。この定義は「まだその事象を起こしていない人だけがリスク集団である」という発生率(incidence rate)の考え方に忠実で、PHUSEの標準解析でもEAIRはこちらを指します。

もうひとつは、分子に全イベント数を数え上げ、分母に総曝露期間をそのまま使う定義です。こちらは曝露調整イベント率(EAER: exposure-adjusted event rate)と呼び分けられることが多く、再発を含めた負担量を表します。

同じ「曝露調整」という言葉で、値も解釈も違う二つの指標が流通しているわけです。したがって、SAPには「分子は初回発現のみか全件か」「分母は初回発現までか総曝露か」「単位は100人年か」を必ず書き切ってください。CTDのモジュール2.7.4では、2.7.4.1.2(曝露の全体像)で投与期間別・用量別の曝露量を示したうえで有害事象の集計を提示することが求められており、曝露量とAE集計は本来セットで読ませる建て付けになっています。以下のRコードは後者(全イベント数 ÷ 総曝露人年)で計算しています。

RでEAIRとポアソン率比を求める

まず各群のEAIRと、その比(率比)を求めます。

e1 <- summ$events[summ$arm == "Active"];  T1 <- sum(ae$exposure[ae$arm == "Active"])
e0 <- summ$events[summ$arm == "Placebo"]; T0 <- sum(ae$exposure[ae$arm == "Placebo"])
cat(sprintf("Active : events = %d, PY = %.1f\n", e1, T1))
cat(sprintf("Placebo: events = %d, PY = %.1f\n", e0, T0))

pt1 <- poisson.test(e1, T1)
pt0 <- poisson.test(e0, T0)
eair <- rbind(
  Active  = c(EAIR = 100 * e1 / T1, lower = 100 * pt1$conf.int[1], upper = 100 * pt1$conf.int[2]),
  Placebo = c(EAIR = 100 * e0 / T0, lower = 100 * pt0$conf.int[1], upper = 100 * pt0$conf.int[2])
)
print(round(eair, 1))

rt <- poisson.test(c(e1, e0), c(T1, T0))
cat(sprintf("rate ratio = %.3f  95%%CI [%.3f, %.3f]  p = %.4g\n",
            rt$estimate, rt$conf.int[1], rt$conf.int[2], rt$p.value))
Active : events = 213, PY = 191.5
Placebo: events = 119, PY = 216.2
         EAIR lower upper
Active  111.2  96.8 127.2
Placebo  55.1  45.6  65.9
rate ratio = 2.021  95%CI [1.607, 2.551]  p = 3.596e-10

実薬群は213件を191.5人年で観察したので、100人年あたり111.2件、95%信頼区間は [96.8, 127.2] です。プラセボ群は119件を216.2人年で観察して55.1件、区間は [45.6, 65.9] です。率比は2.021、95%信頼区間は [1.607, 2.551] でした。区間の下限が1.607と1をはるかに上回っており、曝露時間あたりのイベント発生が実薬群で明確に多いことが読み取れます。

poisson.test() が返す区間は正規近似ではなく厳密区間です。単群の場合はガンマ分布の裾確率を使ったGarwood型の中心区間で、二項分布の正確区間(Clopper-Pearson)を \( n \) が非常に大きい極限に飛ばしたものに相当します。二群の率比の場合は、総イベント数 \( e_1+e_0 \) を固定したうえで「そのうち何件が実薬群で起きたか」という条件付き二項分布に変換して区間を作ります。件数が数件しかない稀な事象でも破綻しない点が、AE集計でこの関数を既定にする理由です。なお p値(3.596e-10)は参考値であり、安全性の評価は区間と大きさで語るのが原則です(信頼区間とp値の関係を図解で理解する)。

曝露で割ると、比が1.41倍から2.02倍に変わる

ここからがこの記事の山場です。同じデータを、発現割合とEAIRの両方で並べてみます。

cmp <- data.frame(
  arm      = c("Active", "Placebo"),
  pct      = c(100 * x1 / n1, 100 * x0 / n0),
  mean_exp = c(mean(ae$exposure[ae$arm == "Active"]), mean(ae$exposure[ae$arm == "Placebo"])),
  EAIR     = c(100 * e1 / T1, 100 * e0 / T0)
)
cmp$mean_exp_weeks <- cmp$mean_exp * 52
cmp$pct <- round(cmp$pct, 1); cmp$EAIR <- round(cmp$EAIR, 1)
cmp$mean_exp <- round(cmp$mean_exp, 3); cmp$mean_exp_weeks <- round(cmp$mean_exp_weeks, 1)
print(cmp)
cat(sprintf("proportion ratio = %.3f / EAIR ratio = %.3f\n",
            (x1 / n1) / (x0 / n0), (e1 / T1) / (e0 / T0)))
      arm  pct mean_exp  EAIR mean_exp_weeks
1  Active 49.6    0.766 111.2           39.8
2 Placebo 35.2    0.865  55.1           45.0
proportion ratio = 1.409 / EAIR ratio = 2.021

注目すべきは mean_exp_weeks です。実薬群の平均曝露は39.8週、プラセボ群は45.0週で、実薬群のほうが5週以上短く観察されています。にもかかわらず、発現割合は実薬49.6%・プラセボ35.2%で、比は1.409倍でした。曝露時間で割ると比は2.021倍に跳ね上がります。

理屈は単純です。観察している時間が長ければ、それだけイベントに遭遇する機会も増えます。実薬群は早く中止する人が多く、平均して短い時間しか観察されていません。つまり実薬群は「イベントを起こすチャンスを与えられないまま集計を打ち切られた」状態です。分母を人数のまま数えると、この不利(薬にとっては有利)が補正されないので、実薬群の危険性が小さく見積もられます。曝露が短い群のほうが、単純な発現割合では安全に見えてしまうのです。

💡 ポイント
発現割合とEAIRのどちらかが「正しい」わけではありません。分母が「人」か「人時間」かが違うだけです。同じデータから1.409倍と2.021倍という2つの数字が出るという事実こそが、AEの数字を出すときに分母を言えなければならない理由です。
指標実薬群プラセボ群分母
発現割合49.6%35.2%1.409被験者数(各250例)
EAIR(100人年あたり)111.255.12.021総曝露人年(191.5 / 216.2)
平均曝露期間39.8週45.0週0.885被験者数(各250例)

この現象が効いてくるのは、長期投与試験と、群間で中止の起こり方が違う試験です。有効性が不足して早期に脱落する被験者が多い群、忍容性の問題で減量・中止が多い群、実薬群の効果が高くイベント再発による中止が少ない群──いずれも群間で人時間が揃いません。extension試験や非盲検継続投与期を含む長期安全性データベース、あるいは複数試験を統合する統合解析(ISS)では、曝露期間の分布が試験ごとに大きく違うため、曝露調整なしの単純合算は誤読を招きます。開発後期の試験で収集する安全性データの範囲を絞る考え方を示したICH E19(2022年9月にStep 4到達)の下でも、集めたデータをどの分母で提示するかという設計は残ります。

⚠️ 注意
EAIRは「発現ハザードが期間を通じて一定」というポアソン仮定に乗った指標です。イベント発生が期間に比例するからこそ、時間で割ることに意味が生まれます。ところがinfusion reactionや初回投与時の悪心のように、投与初期に集中して以後ほとんど起こらないAEでは、この仮定は成り立ちません。曝露が長い被験者の「何も起こらなかった時間」が分母を膨らませ、EAIRは実態より小さく出ます。平均化しすぎて、危険な最初の数週間を薄めてしまうわけです。また、曝露が極端に短い群やサブグループでは分母の人年が小さくなり、EAIRの区間が一気に広がって不安定になります。分母の人年と件数を必ず併記し、EAIRだけを単独で載せないでください。

発現割合とEAIRをどう使い分けるか

実務では両方を出すのが基本です。そのうえで、主たる指標をどちらに置くかを場面で決めます。

場面主に見る指標理由
曝露期間が群間でほぼ揃う短期試験発現割合解釈が直感的で、規制当局の標準集計表と揃う
長期投与試験・中止が群間で偏る試験EAIR(割合も併記)観察時間の差が割合を歪めるため
投与初期に集中するAE発現割合+時期別の集計一定ハザードの仮定が成り立たない
同一被験者の再発が臨床的に重要なAE全件を分子にした率(EAER)初回のみでは負担量を過小評価する

なお、ここで求めた率比2.021は、群という1つの因子だけで比較した粗い推定値です。年齢や併用薬などの共変量で調整したい場合や、群ごとの人時間をオフセットとして明示的にモデル化したい場合には、対数リンクのポアソン回帰が自然な拡張になります。さらに、同じ被験者が何度も発現する体質差(過分散)が強いデータでは、負の二項回帰のほうが標準誤差を正しく見積もります。ただし、AE集計表の第一の役割はモデル推定値を載せることではなく、件数・人年・割合・率という素の数字を、分母を明示して並べることです。そこを外さないでください。

発現までの時間と再発イベントを扱う

発現割合は「500人のうち何人が起こしたか」に答え、曝露調整発現率(EAIR)は「総曝露人年あたり何件起きたか」に答えます。どちらも有用な要約ですが、共通して捨てている情報がひとつあります。いつ起きたかです。

この一点は実務に直結します。注射部位反応のように投与開始直後に集中するAEと、薬剤性肝障害のように曝露が積み上がるほど出やすくなるAEは、発現割合で見ればどちらも「実薬20%、プラセボ10%」という同じ数字になりえます。しかし前者なら初回投与時の観察体制が、後者なら長期投与時のモニタリング間隔が論点になり、リスク管理の設計はまったく違います。時間軸を持ち込まないと、この区別はできません。加えてAEには「同じ人が何度も起こす」という性質があり、再発を数えるなら被験者ごとの起こしやすさのばらつきをモデルに入れる必要があります。ここでは初回発現までの時間をKaplan-Meier法で、再発イベントの件数を回帰モデルで扱います。

初回発現までの時間をKaplan-Meier法で見る

初回AE発現までの時間 t_first と発現有無 event を使い、群別のKaplan-Meier推定量から3か月(0.25年)・6か月(0.5年)・1年(1.0年)時点の累積発現割合 \( 1-\hat{S}(t) \) を取り出します。

library(survival)
fit <- survfit(Surv(t_first, event) ~ arm, data = ae)
s <- summary(fit, times = c(0.25, 0.5, 1.0))
out7 <- data.frame(
  arm     = sub("arm=", "", s$strata),
  time_yr = s$time,
  n.risk  = s$n.risk,
  n.event = s$n.event,
  cum_inc = round(100 * (1 - s$surv), 1),
  lower   = round(100 * (1 - s$upper), 1),
  upper   = round(100 * (1 - s$lower), 1)
)
print(out7)
print(survdiff(Surv(t_first, event) ~ arm, data = ae))
      arm time_yr n.risk n.event cum_inc lower upper
1 Placebo    0.25    213      21     8.7   5.1  12.2
2 Placebo    0.50    177      19    17.2  12.2  21.9
3 Placebo    1.00    112      48    40.9  33.8  47.2
4  Active    0.25    193      30    12.7   8.3  16.8
5  Active    0.50    151      22    23.2  17.4  28.6
6  Active    1.00     48      72    65.7  57.8  72.2
Call:
survdiff(formula = Surv(t_first, event) ~ arm, data = ae)

              N Observed Expected (O-E)^2/E (O-E)^2/V
arm=Placebo 250       88    120.1       8.6        20
arm=Active  250      124     91.9      11.2        20

 Chisq= 20  on 1 degrees of freedom, p= 8e-06 

まず時間の経過に沿って読みます。3か月時点は実薬12.7%・プラセボ8.7%で、差は4ポイントほどしかありません。6か月時点で23.2%対17.2%、1年時点では65.7% [57.8, 72.2] 対40.9% [33.8, 47.2] と、差は時間とともに開いていきます。この試験のAEは投与初期に一気に出るタイプではなく、曝露が続くあいだ一定のペースで積み上がるタイプだ、という読み方ができます。

次に、多くの読者がつまずく点です。1年時点の累積発現割合65.7% / 40.9% は、単純な発現割合49.6% / 35.2% より明らかに高く出ています。同じデータなのに、なぜ増えるのか。理由は、Kaplan-Meier推定量が「打ち切られた被験者も、最後まで試験に残っていれば、残った被験者と同じペースでAEを起こしていたはず」という前提(無情報打ち切り)で補外しているからです。実薬群では多くの被験者が1年を待たずに中止しており、その人たちは「AEを起こさなかった人」として発現割合の分母に丸ごと乗ります。一方Kaplan-Meier法は、その人たちを打ち切り時点でリスク集合から外し、残った被験者の発現ペースで先を埋めます。この差が、49.6%と65.7%の開きの正体です。

したがって「1年時点の累積発現割合」は、観察された事実そのものではなく、打ち切りが無情報であるという仮定の下での推定値です。CSRやCTD 2.7.4に載せるなら、単純な発現割合とKaplan-Meier推定値のどちらを主表に置くかをSAPで先に決め、両方載せる場合はなぜ数字が違うのかを本文で説明してください。Kaplan-Meier推定量そのものの導出は生存時間解析の基礎:カプラン–マイヤー法で扱っています。

⚠️ 注意
1年時点の n.risk は実薬群でわずか48、プラセボ群で112です。実薬群の65.7%という推定値は、残り48人の挙動に強く引きずられており、信頼区間も [57.8, 72.2] と幅広くなっています。KM曲線の裾(右端)は、読者がもっとも誤読しやすい場所です。曲線の見た目は最後まで滑らかに描かれるので、そこで何人が残っているのかは図からは分かりません。リスク集合の人数(number at risk)を必ず曲線の下に併記し、n.riskが十分でない区間の推定値を単独の結論として引用しない。これはAEに限らず、生存時間解析の作図における最低限の作法です。

競合リスクがあるときはKMが過大評価になる

無情報打ち切りの仮定が特に危ういのが、AEと競合する事象が存在する場合です。腫瘍領域の試験で、対象のAEが発現する前に被験者が死亡したとします。この死亡を単なる打ち切りとして扱うと、Kaplan-Meier法は「この人も生きていればいずれAEを起こしたはず」と補外します。しかし現実には、死亡した被験者がその後にそのAEを起こすことはありません。結果として、\( 1-\hat{S}(t) \) は真の累積発現割合を系統的に過大評価します。

この場合に使うのが競合リスクの枠組みで、累積発生関数(cumulative incidence function)を推定します。競合事象が起きた被験者について、その後にイベントが起きたとは数えない形で累積確率を積み上げるため、複数の事象タイプの累積発生関数を足し合わせても1を超えません。上の解析は死亡や競合事象を明示的に扱っていないので、KMの値は「AE以外の理由による中止が無情報である」という前提の上に立っています。長期試験や重篤な疾患を対象とする試験では、この前提を毎回問い直してください。詳しくは競合リスク(competing risks)とはを参照してください。

log-rank検定のp値をどう位置づけるか

survdiff() の出力では、プラセボ群の観測88件に対して期待120.1件、実薬群は観測124件に対して期待91.9件で、\( \chi^2=20 \)(自由度1)、p= 8e-06 となっています。数字だけ見れば非常に小さなp値です。

しかし、このp値を確証的に読んではいけません。ICH E9の第6.4節(安全性・忍容性の統計的評価)は、安全性の評価を記述的な統計手法で行い、解釈の助けになる場合に信頼区間を添えることを基本方針としています。そのなかでp値は、関心のある特定の差を評価する補助や、多数の安全性変数から注目すべき差を拾い上げる「フラグ」という位置づけにとどまります。

この試験でAEの発現時間は主要評価項目でも事前に規定した確証的仮説でもありません。SOC・PT単位まで下ろせば数百のlog-rank検定が並び、いくつかが偶然に小さなp値を示すのは当然です。ここでのp= 8e-06は記述の補助として提示するものであり、有効性の主要解析と同じ重みで扱う数字ではありません。安全性の結論は、区間推定の幅と医学的レビューで組み立てます。

再発イベントを率としてモデル化する

初回発現までの時間だけを見ると、2回目以降のAEは情報として捨てられます。総イベント数が実薬213件・プラセボ119件と、発現者数(124人・88人)を大きく上回るこのデータでは、再発を数え上げるモデルも併せて見る価値があります。ここではポアソン回帰と負の二項回帰を当てます。

library(MASS)
m_pois <- glm(n_ae ~ arm + offset(log(exposure)), family = poisson, data = ae)
m_nb   <- glm.nb(n_ae ~ arm + offset(log(exposure)), data = ae)
cat("Poisson: rate ratio =", round(exp(coef(m_pois)["armActive"]), 3),
    " SE =", round(summary(m_pois)$coefficients["armActive", 2], 4), "\n")
cat("Poisson dispersion (Pearson X2/df) =",
    round(sum(residuals(m_pois, type = "pearson")^2) / m_pois$df.residual, 2), "\n")
cat("NB     : rate ratio =", round(exp(coef(m_nb)["armActive"]), 3),
    " SE =", round(summary(m_nb)$coefficients["armActive", 2], 4), "\n")
print(round(exp(confint.default(m_nb)["armActive", ]), 3))
cat("NB theta =", round(m_nb$theta, 3), " (SE", round(m_nb$SE.theta, 3), ")\n")
cat("AIC: Poisson =", round(AIC(m_pois), 1), " NB =", round(AIC(m_nb), 1), "\n")
Poisson: rate ratio = 2.021  SE = 0.1144 
Poisson dispersion (Pearson X2/df) = 1.21 
NB     : rate ratio = 2.014  SE = 0.1247 
 2.5 % 97.5 % 
 1.577  2.571 
NB theta = 3.698  (SE 1.644 )
AIC: Poisson = 1053  NB = 1047.6 

まず offset(log(exposure)) の役割を確認します。被験者 \( i \) のイベント数を \( Y_i \)、曝露期間を \( T_i \)、期待件数を \( \mu_i=E[Y_i] \) とします。対数リンクの下で、群を表す指示変数 \( x_i \) を使ってモデルを

\[ \log\mu_i=\log T_i+\beta_0+\beta_1 x_i \]

と書くと、両辺から \( \log T_i \) を移項して

\[ \log\frac{\mu_i}{T_i}=\beta_0+\beta_1 x_i \]

となります。左辺は「単位時間あたりの期待件数」、すなわちの対数です。オフセット項を置くことで、モデルの対象が「件数」から「率」へ切り替わります。係数 \( \beta_1 \) は率の対数の差なので、その指数 \( \exp(\beta_1) \) がそのまま率比になります。曝露を共変量として推定するのではなく、係数を1に固定した項として押し込むのがオフセットです。

出力を見ると、ポアソン回帰の率比は2.021、負の二項回帰は2.014で、点推定はほとんど動きません。動くのは精度のほうです。標準誤差は0.1144から0.1247へ、約9%大きくなっています。この差を生んでいるのが過分散です。Pearson \( \chi^2 \) を残差自由度で割った値は1.21で、ポアソン分布が仮定する「分散=平均」より2割ほど散らばりが大きいことを示しています。

負の二項分布は、この余分なばらつきを分散パラメータ \( \theta \) で吸収します。平均 \( \mu \) に対する分散は

\[ \mathrm{Var}(Y)=\mu+\frac{\mu^2}{\theta} \]

で与えられ、右辺第2項が過分散にあたります。\( \theta \) が大きいほど第2項は小さくなり、\( \theta\to\infty \) でポアソン分布に一致します。ここでの推定値は \( \theta=3.698 \)(SE 1.644)でした。たとえば \( \mu=1 \) 件のとき分散は \( 1+1/3.698 \)、およそ1.27となり、先ほどのPearson \( \chi^2 \)/df = 1.21 と整合します。モデル選択の指標であるAICも1053から1047.6へ下がり、負の二項回帰のほうが当てはまりが良いという判断になります。負の二項回帰の率比の95%信頼区間は [1.577, 2.571] です。

💡 ポイント
点推定が2.021と2.014でほぼ同じだからといって「どちらのモデルでも同じ」ではありません。過分散を無視すると、標準誤差が小さく出て信頼区間が狭くなりすぎ、群間差を実際より確からしいものとして提示してしまいます。AEの件数データは、少数の被験者が多数のイベントを起こす偏った分布になりがちで、過分散はむしろ既定と考えるべきです。ポアソン回帰を当てたら、Pearson \( \chi^2 \)/df を必ず確認してください。

なお、ポアソン回帰・負の二項回帰それぞれのモデル診断やリンク関数の詳細はポアソン回帰とはおよび負の二項回帰と再発イベント率解析で扱っています。ここで押さえるべきは、AEの再発を数える以上は個体間のばらつきを明示的にモデル化する必要がある、という一点です。発現のタイミングまで含めて再発を扱うなら、Cox比例ハザードモデルを計数過程の枠組みに拡張したAndersen–Gillモデルや、時点 \( t \) までの平均累積イベント数を表す平均累積関数(mean cumulative function)という選択肢もあります。ベースラインのイベント率の形状を仮定せずに済む点が特徴です。

三つの見方で数字がどう変わるか

同じデータに三つの分母を当てた結果を並べます。どれかが正しくてどれかが間違っているのではなく、答えている問いが違います。

見方(分母)実薬プラセボ群間比較置いている仮定
発現割合(分母=被験者数250)49.6%35.2%比 1.41倍曝露期間の群間差を考慮しない
EAIR(分母=総曝露人年)100人年あたり111.2100人年あたり55.1比 2.02倍発現ハザードが期間を通じて一定
KM 1年累積発現割合(分母=各時点のリスク集合)65.7% [57.8, 72.2]40.9% [33.8, 47.2]差 24.8ポイント打ち切りが無情報・競合リスクなし
負の二項回帰(再発を数える・分母=曝露期間)総イベント213件総イベント119件率比 2.014 [1.577, 2.571]個体間のばらつきを \( \theta \) で吸収

数字は49.6%から65.7%まで、比は1.41倍から2.02倍まで動きます。この幅は解析の失敗ではなく、AEというデータが本来もっている多面性です。安全性の要約表に数字を書くときは、その分母が何で、どの仮定の下で成り立っているのかを常に言えるようにしておいてください。

実務でのポイント

CSRに載る安全性の表は、全体の発現割合が1行あって終わりではありません。MedDRAの器官別大分類(SOC)と基本語(PT)で層別され、行数は数百に達します。ここでは同じデータから5つのPTを取り出し、PTごとに発現割合とリスク差の信頼区間を並べてみます。

set.seed(20260810)
pt_names <- c("Headache", "Nausea", "Diarrhoea", "Dizziness", "Rash")
prob_act <- c(0.30, 0.25, 0.20, 0.15, 0.10)
prob_pla <- c(0.32, 0.22, 0.18, 0.13, 0.15)
long <- do.call(rbind, lapply(which(ae$n_ae > 0), function(i) {
  k <- ae$n_ae[i]
  p <- if (ae$arm[i] == "Active") prob_act else prob_pla
  data.frame(subjid = ae$subjid[i], arm = ae$arm[i],
             pt = sample(pt_names, k, replace = TRUE, prob = p))
}))
inc <- t(sapply(pt_names, function(p) {
  s1 <- length(unique(long$subjid[long$pt == p & long$arm == "Active"]))
  s0 <- length(unique(long$subjid[long$pt == p & long$arm == "Placebo"]))
  rd <- unname(rd_newcombe(s1, n1, s0, n0))
  c(Active_n = s1, Active_pct = 100 * s1 / n1,
    Placebo_n = s0, Placebo_pct = 100 * s0 / n0,
    RD = 100 * rd[1], lower = 100 * rd[2], upper = 100 * rd[3])
}))
print(round(inc, 1))
          Active_n Active_pct Placebo_n Placebo_pct   RD lower upper
Headache        63       25.2        37        14.8 10.4   3.4  17.3
Nausea          36       14.4        25        10.0  4.4  -1.4  10.2
Diarrhoea       32       12.8        25        10.0  2.8  -2.8   8.5
Dizziness       36       14.4        12         4.8  9.6   4.5  14.9
Rash            20        8.0        12         4.8  3.2  -1.2   7.7

分子は件数ではなく人数で数える

最も重要なのは length(unique(subjid)) で数えている点です。同じ被験者が同じPTを2回起こしても、分子には1例として数えます。延べ件数を分子にすると分母の被験者数を超えかねず、割合が1を超えるという事故が起こります。ICH E9も、AEの発現割合は通常「イベントを経験した被験者数を、リスクにさらされた被験者数に関係づけた割合」として表現されると述べています。

数字を読むと、Headacheは実薬25.2%・プラセボ14.8%でリスク差 \( \hat{p}_1-\hat{p}_0 \) は10.4%ポイント、95%信頼区間[3.4, 17.3]。Dizzinessは14.4%対4.8%で9.6[4.5, 14.9]。この2つは区間が0を含みません。一方Nauseaは4.4[-1.4, 10.2]、Diarrhoeaは2.8[-2.8, 8.5]、Rashは3.2[-1.2, 7.7]と0をまたぎます。どのPTでも区間幅が9〜14%ポイントに達することに注目してください。n=250でも、PT単位の解像度はこの程度です。

並べた区間を有意判定に使ってはいけない

この5本には多重性の調整が一切入っていません。実際のSOC/PT集計では、重症度別・因果関係別の内訳まで含めれば比較は数百から千に達します。名目95%の区間を千本並べれば、真の差が皆無でも数十本は0をまたがない区間が現れます。「Headacheの区間が0を含まないから実薬で有意に多い」という読み方は、この表からは成立しません。

では調整すればよいかというと、そこも単純ではありません。ICH E9は、仮説検定を用いるなら第一種の過誤に対する多重性調整は適切だとしつつ、安全性では通常むしろ第二種の過誤のほうが懸念されると述べています。見逃しのほうが怖い領域で強い調整をかければ、真のシグナルまで潰れます。だから実務のシグナル検出は、検定ではなく記述統計と医学的レビューの組み合わせで行います。作用機序から説明がつくか、投与開始からの時間的関係が自然か、関連PTを標準検索式(SMQ)で束ねたときに一貫するか、個別症例のナラティブが納得できるか。この往復で追う・追わないが決まります。

⚠️ 注意
PT別の区間は、差の大きさと精度を伝える記述的な道具です。表に載せるなら「多重性の調整は行っておらず、推定精度を示す目的で提示する」旨をSAPと脚注に明記してください。多重性の考え方は多重比較法の基礎と活用で整理しています。

ICH E9が安全性の解析をどう位置づけているか

ICH E9(1998年)の立場は明確です。多くの試験では安全性・忍容性の示唆は記述統計で最もよく扱われ、解釈の助けになる場合に信頼区間の計算で補う、とされています。群間に臨床的に許容できない差が無いことを示す状況では、仮説検定より信頼区間の使用が好まれるとも書かれています。p値は、特定の差を評価する補助か、多数の安全性変数から注目に値する差を浮かび上がらせる「フラグ付け」の道具という位置づけにとどまります。同じ章は、起こり得る有害作用の範囲は非常に広く新規かつ予見不可能な作用も常にあり得るため、検証的試験から結論的な情報が得られるのは例外だと明言しています。安全性の解析が探索的・仮説生成的だと言われる根拠はここにあります。

ただし「安全性はすべて記述的でよい」は行き過ぎです。ICH E9は、他剤や他用量に対する安全性の優越性・同等性を主張する目的で計画された試験には、有効性と同様の確証的エビデンスを求めています。詳しくはICH E9「臨床試験の統計的原則」と補遺をご覧ください。

SAPに事前に書いておくべきこと

記述的だから事前規定は緩くてよい、というのは誤解です。検定が無いぶん、定義の一つひとつが数字を直接動かします。曝露期間の数え方を変えるだけで、曝露調整発現率は数パーセント動きます。

項目SAPで決めておくこと
安全性解析対象集団治験薬を1回以上投与された被験者と定義する。割付群と実投与群のどちらで分類するかも明記
TEAEの定義発現期間のウィンドウ(初回投与日から最終投与日+何日まで)と、ベースラインから悪化した事象の扱い
曝露期間の算出規則最終投与日までか+1日か、休薬期間を差し引くか、減量期間をどう数えるか
EAIRの分子全発現件数か初回発現のみか。初回発現のみなら分母も初回発現までの時間に揃えるか
重症度・因果関係複数回発現時に最悪の重症度を採用する規則、判定が欠測のときの寄せ方
集計単位と辞書MedDRAのバージョン、SOC/PTの階層、SOC内の並び順
表示閾値いずれかの群で5%以上のPTを表示する、など閾値と適用対象
区間推定と多重性どの手法で区間を作るか、多重性の調整を行わない旨と、区間を有意判定に用いない旨
特に注目すべきAEAESIの定義と対応するPTのリスト、標準検索式を使う場合はその指定

治験中の安全性情報の取り扱いと報告の枠組みはICH E2Aが定めています。SAPの安全性の章は、これと矛盾しない形で書かれている必要があります。

CSR・CTD 2.7.4・DSURでこの数字がどう使われるか

試験単位では、これらの数字はCSRの安全性の章に載ります。曝露量の要約、SOC/PT別の発現被験者数と割合、重篤なAE・中止に至ったAE・死亡の個別記述という構成で、本記事で扱った発現割合と信頼区間、曝露調整発現率はその骨格そのものです。申請段階では複数試験を束ねた要約がCTDのモジュール2.7.4(臨床的安全性の概要)に置かれます。2.7.4.1が治験薬への曝露、そのなかに全体的な曝露の程度を記述する項が独立して立ち、続く2.7.4.2で有害事象を扱います。曝露の程度がAE解析より先に独立項として求められる構成自体が、「AEの数字は曝露量とセットでなければ読めない」という考え方の制度的な表れです。

開発期間中は、これらが年1回のDSUR(開発中の医薬品の安全性最新報告、ICH E2F)に集約されます。データロックポイントを1つ定め、その時点までの全試験の安全性情報を統合し、曝露状況、期間中の重篤な副作用の一覧、新たな知見、ベネフィット・リスクの評価を報告する文書です。個々の試験のAE集計は、そこに至る部品です。安全性データの収集をどこまで絞り込めるかという議論はICH E19で扱われています。

盲検下の安全性レビューとIDMC

試験実施中の安全性の見方には2つの経路があります。1つは盲検下の安全性レビューで、投与群を伏せたままプールしたAE一覧を眺め、コーディングの妥当性、重篤性や因果関係の判定の一貫性、発現期間の運用を点検します。群間比較をしないので試験の完全性は損なわれません。もう1つが群を開けて比較する経路です。ICH E9は、特に公衆衛生上の意義が大きい試験では有効性・安全性の監視責任を外部の独立した集団に委ねるべきだとしており、これが独立データモニタリング委員会(IDMC)やDSMBと呼ばれる組織です。委員会は定められた間隔で進行状況・安全性データ・主要エンドポイントを評価し、継続・変更・中止を依頼者に勧告します。依頼者自身が非盲検情報を扱う場合は、情報共有を管理・限定し、手順書の遵守と議事録の維持を文書で保証することが求められます。

💡 ポイント
PT別の発現割合は必ず被験者数で数え、延べ件数を分子にしないでください。区間は差の精度を示すために付けるものなので、多重性を調整していない数百本の比較を有意判定には使えません。検定に頼れないぶん、安全性解析対象集団・TEAEの定義・曝露期間の算出規則・分子の定義・表示閾値をSAPで決め切ることが、安全性統計の実質的な設計作業になります。

参考書籍

安全性データの解析を独学で固めるには、性格の違う3冊が要ります。臨床試験の設計と解析を通しで扱う本、件数データの過分散をモデルの言葉で理解する本、そして統計解析の外側にある医薬品安全性監視の枠組みを押さえる本です。

『臨床試験ハンドブック ―デザインと統計解析―』丹後俊郎・上坂浩之 編(朝倉書店)
無作為化比較試験のデザインと統計解析の決定版。772ページの大部で、臨床試験の計画・実施・デザイン・統計解析を通しで扱う和書です。本記事との関係で価値が高いのは、有害事象の経時的な発現状況の推測を扱った章があること。発現割合という一時点の要約から、時間軸を持った見方へ踏み出すときの足場になります。安全性解析対象集団や曝露の定義をSAPに落とし込む段階で、用語と考え方を確認する参照本としても使えます。2021年に新装版が出ています。
『データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC』久保拓弥 著(岩波書店)
通称「緑本」。ポアソン分布から始めて、なぜ現実の件数データは分散が平均を超えるのか、その過分散を個体差としてモデルに入れるとどうなるかを、手を動かしながら理解できます。本記事で負の二項回帰の標準誤差がポアソン回帰より大きくなった理由や、分散パラメータ ( \theta ) が何を表すのかが腑に落ちるはずです。オフセット項の扱いも登場するので、曝露時間で割った率のモデルを自分で書けるようになります。
『医薬品安全性監視入門 第2版 ―ファーマコビジランスの基本原理―』Patrick Waller・Mira Harrison-Woolrych 著/久保田潔 監訳(じほう)
英国で長年ファーマコビジランスに携わってきた2人の専門家による入門書(原著第2版・2017年)の邦訳で、通読しやすい分量にまとまっています。本記事で「シグナル検出は記述統計と医学的レビューの組み合わせで行う」と書いた、その医学的レビュー側の作法を学べます。有害事象と副作用の区別、因果関係の評価、リスクの伝達といった統計解析の前後の工程を押さえておくと、集計表の数字を誰に何のために渡しているのかが見えてきます。安全性部門と会話するための共通語彙を得るのに適した一冊です。

まとめ

本記事では、有害事象データの集計と区間推定を、発現割合・曝露調整発現率・時間軸という3つの見方でRを動かしながら追いました。

出てきた数値を並べ直します。各群250例・計画曝露52週の試験で、1件以上のAEが発現した被験者の割合は実薬49.6%・プラセボ35.2%、リスク差は14.4%ポイント、Newcombeの95%信頼区間は[5.7, 22.8]でした。ところが実薬群の平均曝露は39.8週、プラセボ群は45.0週で、実薬群のほうが1割ほど短い期間しか観察されていません。曝露で割った発現率は100人年あたり111.2件と55.1件、率比は2.021[1.607, 2.551]となり、単純な割合の比1.409から大きく動きます。さらに打ち切りを考慮したKaplan-Meierでは、1年時点の累積発現割合が65.7%と40.9%になります。同じ1つのデータセットから、49.6%・111.2件/100人年・65.7%という3つの数字が出てくるわけです。

ここから引き出せる実務的な結論は3つです。第一に、AEの数字を出すときは分母が何かを必ず言えなければならないこと。人を分母にするのか、人時間を分母にするのか、打ち切りを考慮した時間軸で見るのかで、値も解釈も変わります。第二に、区間で語ること。発現1〜2例のPTが数百行並ぶ安全性の集計表では、点推定だけを見て差を論じることに意味はなく、Wilson区間やNewcombe区間で不確実性の幅を示すことが出発点になります。1件も観察されなかったPTでも、n=250なら真の発現割合が1.46%までありうるという事実は、その典型です。第三に、再発を数えるなら過分散を疑うこと。ポアソン回帰から負の二項回帰に切り替えると率比の点推定は2.021から2.014とほとんど動かないのに、標準誤差は約9%大きくなりました。区間の狭さは精度ではなく、モデルの仮定が生んだ見かけであることがあります。

そして、これらすべてに共通する前提が「事前規定」です。TEAEの定義、曝露期間の算出規則、EAIRの分子に何を数えるか、集計に載せる閾値。これらを解析計画書に書き切っておくことが、安全性の数字を後から動かさないための唯一の防御になります。

関連するトピックとしては、イベント件数と発生率のモデリングそのものを扱ったポアソン回帰とは ― 件数・発生率データをRで解析する一般化線形モデル入門 ―、再発イベントの解析を喘息・COPD増悪の文脈で掘り下げた負の二項回帰と再発イベント率解析、打ち切りと競合事象の扱いを整理した競合リスク(Competing Risks)とは ― Fine-Grayモデルと累積発生関数をRで実装 ―、そして開発後期における安全性データ収集の考え方を扱ったICH E19とは ― 開発後期・承認後試験における「選択的な安全性データ収集」をわかりやすく解説 ―もあわせてご覧ください。

まずは担当している試験のAE集計表を1枚開いて、その分母が何になっているかを確認してみてください。人なのか、人時間なのか。そこが言えるようになるだけで、安全性の数字に対する解像度は大きく変わります。

ABOUT ME
tomokichi
外資系製薬会社で生物統計家として働ている1児のパパ。生物統計家とは何か、どのようなスキルが必要か、何を行っているのかを共有していきたいと思っております!生物統計に関する最新情報を皆様にお届けすべく、日々奮闘中です。趣味は筋トレ、温泉巡り、家族と散歩。