コホート研究で発生頻度を計算するとき、分母を「人数」にするか「観察時間」にするかで、得られる量の意味が変わります。登録が2年、3年と続く研究では、最初に登録された患者を5年近く追跡できても、最後に登録された患者は2年ほどで研究終了を迎えることがあります。このデータでイベント数を単純に登録人数で割ると、イベントの起こりやすさだけでなく、各群がどれだけ長く観察されたかまで数字に入り込みます。
そこで使われるのが、人時間（person-time）を分母にした発生率です。ただし、3年以内にイベントを経験する確率を知りたい研究と、1人年あたりどの程度の速さでイベントが発生するかを知りたい研究では、そもそも推定したい量が違うため、追跡期間がそろっていない＝人時間を使うということではない点だけご留意ください。
ここでは、対照群600例、治療群600例のデータセットをランダムに生成し、人数を分母にした割合と、人時間を分母にした率を比べます。前提となる打ち切りと追跡時間の考え方は、生存時間解析入門：打ち切りがあるデータで何を推定するのかで扱っています。対象者$i$の観察期間を$T_i$、その期間に初回イベントが観察されたかどうかを$D_i$、群を$g$と書きます。
人数の割合には追跡期間も入り込む
今回のデータでは、心血管イベントの初回発生をアウトカムとしました。対照群は研究開始から0〜2年の間に登録し、治療群は1〜3年の間に登録する設定です。研究は開始から5年で終了するため、治療群の方が構造的に短い期間しか観察できません。生成されたデータでは、追跡期間の中央値が対照群3.27年、治療群2.66年でした。
イベントは対照群195例、治療群90例でした。これを各群600例で割ると、イベントを経験した割合は0.325と0.150で、その比は0.462です。一方、データ生成時の1人年あたり発生率は対照群0.120、治療群0.072なので、真の率比は0.600です。
0.462が0.600より小さくなった理由には、治療効果だけでなく観察機会の違いも入っています。今回の設定では治療群の登録が遅いため、イベントを経験するまでの時間そのものが短くなりました。人数を分母にした0.150は、その短い追跡まで含めた割合です。
ここで「追跡期間が短い群では、イベント割合が必ず小さくなる」と一般化するのは避け、群間で本来の発生率が大きく異なれば、その影響が追跡期間の差を上回ることもあります。提示できることは、同じような発生過程なら観察時間が短いほどイベントを経験する機会が減り、単純なイベント割合を小さくする方向へ働くということです。
人時間は観察された時間を足し合わせる
追跡時間の違いを分母へ反映するには、イベント数を観察時間の合計で割ります。群$g$の発生率は次のように書けます。
$$\hat{\lambda}_g=\frac{\sum_{i\in g}D_i}{\sum_{i\in g}T_i}\tag{1}$$
式(1)の分子は群$g$で観察された初回イベント数、分母はイベントまたは打ち切りまでに各対象者が観察された時間の合計です。時間を年で測れば、分母の単位は人年になります。
人時間には足し算ができます。ある患者を3年間観察した場合、その3年を1年ずつ3区間に分けても、合計すれば3人年です。この性質があるため、年齢、暦年、治療開始後の期間などで追跡を分割し、それぞれの区間にイベント数と人時間を割り振れます。Breslow（1984）は、こうしたincidence density dataから共通率比を推定し、層ごとの率比の異質性や傾向を調べる方法を扱っています。
一方、600人を1年間観察して得た600人年と、120人を5年間観察して得た600人年は、同じ情報を持っているとは限りません。イベント発生の速さが追跡開始後の時間によって変わるなら、前者は主に早期、後者はより長期の発生過程を反映します。私は人時間だけを単独では報告せず、対象者数と追跡期間の分布も併記します。
show <- c(which(cohort$group == "control")[1:40],
which(cohort$group == "treatment")[1:40])
y <- seq_along(show)
plot(cohort$entry[show], y, type = "n", xlim = c(0, 5),
xlab = "Calendar time (years)", ylab = "Selected participants")
segments(cohort$entry[show], y, cohort$entry[show] + cohort$followup[show], y,
lty = ifelse(cohort$group[show] == "control", 1, 2))
abline(v = 5, lty = 3)
率比と一定時点のリスク比は別の量
人時間を分母にした2群の発生率を比べたものが発生率比（incidence rate ratio）です。リスク比と区別するため、ここでは$\mathrm{IRR}$と書きます。
$$\mathrm{IRR}=\frac{\hat{\lambda}_1}{\hat{\lambda}_0}
=\frac{D_1/T_1}{D_0/T_0}\tag{2}$$
式(2)の$D_g$は群$g$のイベント数、$T_g$は人時間の合計です。治療群の1人年あたり発生率を対照群の1人年あたり発生率で割った値になります。
今回のデータでは、対照群は1,761.3人年で195イベントなので100人年あたり11.07、治療群は1,526.2人年で90イベントなので100人年あたり5.90でした。率比は0.533で、95%信頼区間は0.415から0.684です。この1回のランダム標本では、人数を分母にした比0.462より、データ生成時に置いた率比0.600に近い値になりました。
ただし、率比0.600と「3年以内にイベントを経験する確率の比」は同じ数字にはなりません。今回のように発生率を時間によらず一定とした指数分布なら、時点$t$までのイベント発生確率は$1-\exp(-\lambda t)$です。真の発生率0.120と0.072から計算すると、リスク比は1年時点0.614、2年時点0.629、3年時点0.643となります。
率比は0.600のままなのに、リスク比は時点とともに変わっています。
この違いを考えると、追跡期間が患者ごとに異なるという理由だけで、率比を主要指標に切り替えるのは適切ではありません。研究目的が「1人年あたりの発生速度を比べる」なら率比が対応しますが、一方、「治療開始後3年以内に心血管イベントを経験する確率を比べる」のであれば、3年時点のリスク差やリスク比の方が研究目的に直接対応します。VandenbrouckeとPearce（2012）が区別しているincidence rateとcumulative incidenceも、この違いです。
proportion <- tapply(cohort$event, cohort$group, mean)
proportion_ratio <- proportion["treatment"] / proportion["control"]
plot(c(proportion_ratio, rate_ratio), c(2, 1), log = "x",
xlim = c(0.35, 0.8), yaxt = "n", xlab = "Effect measure", ylab = "")
axis(2, at = c(2, 1), labels = c("Event proportion ratio", "Incidence rate ratio"))
segments(ci_irr[1], 1, ci_irr[2], 1)
abline(v = 0.6, lty = 2)
全期間の率にまとめる前に時間変化を見る
人時間を使えば追跡期間の違いが全て解消される、という理解にも注意が必要です。式(1)は全期間のイベント数と人時間をまとめて割っているため、時間とともに発生率が変化していても、最後には1つの値が返ってきます。その値は、各期間の率を人時間で重み付けした要約と考える方が近くなります。
ここで群ごとの追跡時間の分布が違うと、問題が見えやすくなります。対照群には4〜5年目の人時間がある一方、今回の治療群は登録開始が1年遅いため、4年以上追跡された患者がいません。もし4〜5年目に発生率が急に上昇する疾患なら、対照群の全期間率にはその高い率が入り、治療群には入りません。この状態で2つの全期間率を割ると、治療だけでなく比較している追跡時期の違いまで率比へ入ります。
今回のデータでは、1年ごとの率を並べると、対照群は100人年あたり10.73、9.52、11.93、12.98、12.29、治療群は5.33、7.31、5.19、4.54でした。データ生成時には発生率を対照群12.0、治療群7.2で一定にしているため、これらの上下は標本変動によるものです。少なくとも今回の値から、時間とともに一方向へ大きく動く様子は出ていません。
そして、治療群の4〜5年区間だけは人時間が0なので、率そのものを計算できません。この欠け方がある以上、時間によって率が変わる研究では全期間の0.533だけを見ない方がよいでしょう。私なら、まず期間別のイベント数と人時間を確認し、率比が時期によって大きく変わっていないかを見ます。Breslow（1984）が共通率比とともに層別率比の異質性や傾向を扱っている理由もここにあります。
split_cohort <- survSplit(Surv(followup, event) ~ ., data = cohort,
cut = 1:4, episode = "period")
py_period <- tapply(split_cohort$followup - split_cohort$tstart,
list(split_cohort$group, split_cohort$period), sum)
ev_period <- tapply(split_cohort$event,
list(split_cohort$group, split_cohort$period), sum)
period_rate <- 100 * ev_period / py_period
plot(1:5, period_rate["control", ], type = "b", ylim = c(0, 15),
xlab = "Time since entry (1-year interval)", ylab = "Rate per 100 person-years")
lines(1:5, period_rate["treatment", ], type = "b", pch = 17)
abline(h = c(12, 7.2), lty = 2)
人時間の割り付けと信頼区間にも前提がある
人時間を使う解析では、「何年観察したか」だけでなく、その時間をどの群へ割り付けるかが結果を左右します。薬剤疫学で典型的なのがimmortal time biasです。ある薬を処方された患者を最初から治療群として扱い、実際の初回処方より前の期間まで治療群の人時間へ入れる状況を考えます。
最終的に治療群へ分類されるには、患者は処方時点までイベントを起こさず追跡下に残っていなければなりません。その処方前期間を治療曝露後の人時間として数えると、治療群にイベントのない時間が余計に積み上がり、治療が有利な方向へ率が動きます。Suissa（2008）は、この時間の誤分類によって生じるimmortal time biasを薬剤疫学の例から整理しています。
Levesqueら（2010）は、スタチンと糖尿病進行を題材に、このバイアスが実際のコホート解析でどのように生じるかを示しました。人時間の集計は最後に行う単純な足し算ではなく、追跡開始、曝露開始、イベント、打ち切りを時間軸上でどう定義したかの結果です。
脱落も同じです。患者が脱落した時点までの人時間を数えること自体はできますが、脱落の起こり方がその後のイベントリスクと結びついていると、観察された人時間だけから対象集団の発生率を推定する根拠が弱くなります。この点は情報的打ち切り：脱落した患者は他の患者と同じだったかで扱っています。
率比の信頼区間についても、条件を付けて考えます。独立なPoissonイベント数を仮定し、観察された人時間を所与として扱う単純な2群比較では、対数率比の標準誤差を次の形で近似できます。
$$\mathrm{SE}\{\log(\widehat{\mathrm{IRR}})\}
=\sqrt{\frac{1}{D_1}+\frac{1}{D_0}}\tag{3}$$
式(3)では、観察されたイベント数$D_1$と$D_0$が信頼区間の幅を直接決めます。今回の90イベントと195イベントを入れると標準誤差は0.127となり、率比0.533の95%信頼区間は0.415から0.684です。
式(3)が表しているのは、実際に観察されたイベント数を条件にしたときの近似標準誤差です。そのため、「3,000人年以上を追跡したから精度が高い」と人時間だけを根拠にする書き方は勧めず、何件のイベントが観察されたのかを一緒に見ます。また、再発イベント、個人内相関、施設内相関、過分散がある場合には、この単純なポアソン標準誤差をそのまま使えないことがあります。
Rで率比を確認し、論文の数字へつなげる
ここまでのデータは、base R、stats、survivalパッケージだけで再現できます。対照群は研究開始から0〜2年、治療群は1〜3年の間に登録し、研究終了を5年としました。イベント時間は群ごとに一定の率を持つ指数分布、脱落時間も指数分布から生成しています。
library(survival)
set.seed(42)
n_per_group <- 600
entry_time <- c(runif(n_per_group, 0, 2), runif(n_per_group, 1, 3))
group <- factor(rep(c("control", "treatment"), each = n_per_group))
true_rate <- ifelse(group == "control", 0.12, 0.072)
event_time <- rexp(2 * n_per_group, rate = true_rate)
dropout_time <- rexp(2 * n_per_group, rate = 0.05)
admin_time <- 5 - entry_time
cohort <- data.frame(group, entry = entry_time,
followup = pmin(event_time, dropout_time, admin_time),
event = as.integer(event_time <= pmin(dropout_time, admin_time)))群ごとのイベント数と人年を合計すると、率と率比をそのまま計算できます。あわせて、同じ値が人時間の対数をオフセットにしたポアソン回帰から得られることも確認します。
n_event <- tapply(cohort$event, cohort$group, sum)
person_years <- tapply(cohort$followup, cohort$group, sum)
incidence_rate <- n_event / person_years
rate_ratio <- incidence_rate["treatment"] / incidence_rate["control"]
se_log_irr <- sqrt(1 / n_event["treatment"] + 1 / n_event["control"])
ci_irr <- exp(log(rate_ratio) + c(-1, 1) * qnorm(0.975) * se_log_irr)
risk_ratio_theory <- (1 - exp(-0.072 * (1:3))) / (1 - exp(-0.12 * (1:3)))
rate_difference <- 100 * (incidence_rate["treatment"] - incidence_rate["control"])
agg <- data.frame(group = names(n_event), events = as.numeric(n_event),
py = as.numeric(person_years))
poisson_fit <- glm(events ~ group, family = poisson(link = "log"),
offset = log(py), data = agg)
round(c(rate_ratio, ci_irr, rate_difference, risk_ratio_theory,
exp(coef(poisson_fit)[2])), 3)#> Incidence rate ratio 0.533
#> 95% CI lower 0.415
#> 95% CI upper 0.684
#> Rate difference -5.171
#> Risk ratio at 1 year 0.614
#> Risk ratio at 2 years 0.629
#> Risk ratio at 3 years 0.643
#> Poisson rate ratio 0.533ポアソン回帰の群係数を指数化した値は0.533で、イベント数を人年で割って求めた率比と一致します。Frome（1983）が扱った率のポアソン回帰では、オフセットに人時間の対数を置くことで、イベント数ではなく単位時間あたりの率をモデル化します。
もう一つ見ておきたいのが率差です。治療群と対照群の差は100人年あたり-5.17でした。これは「治療によって100人年あたり5.17件防げた」という因果効果を自動的に意味する数字ではなく、今回の2群で観察された発生率の差です。因果効果として解釈するには、群の比較可能性や研究デザインについて別の根拠が必要になります。
論文では、私は「100人年あたり11.07対5.90、率比0.53、95%信頼区間0.41〜0.68」と率だけを書くより、イベント数195対90、人年1,761.3対1,526.2も同じ表に置きます。これなら読者は分子と分母を確認できます。追跡期間の中央値も併記すれば、同じ人年がどのような追跡から積み上がったのかまで見えます。
最終的に率比を主要な効果指標にするかどうかは、追跡期間がそろっているかだけでは決まりません。1人年あたりの発生速度を比較したいのか、3年時点の発生確率を比較したいのか、それとも時間とともに変わるハザードを扱いたいのかで選ぶ量が変わります。率と瞬間的なハザードの違いは、生存関数とハザード関数：ハザード比の正しい読み方で扱っています。
実際の研究では、登録時期、曝露開始時点、脱落、時間依存曝露、期間ごとの発生率まで見たうえで、人時間をどう区切るかを決めます。Dr.データサイエンスでは、こうした追跡データの構造と解析方針についてご相談を承っています。
参考文献
- Frome EL. The analysis of rates using Poisson regression models. Biometrics. 1983;39(3):665-674.
- Breslow NE. Elementary methods of cohort analysis. International Journal of Epidemiology. 1984;13(1):112-115.
- Suissa S. Immortal time bias in pharmacoepidemiology. American Journal of Epidemiology. 2008;167(4):492-499.
- Levesque LE, Hanley JA, Kezouh A, Suissa S. Problem of immortal time bias in cohort studies: example using statins for preventing progression of diabetes. BMJ. 2010;340:b5087.
- Vandenbroucke JP, Pearce N. Incidence rates in dynamic populations. International Journal of Epidemiology. 2012;41(5):1472-1479.
