2型糖尿病の治療比較のように同じ患者のHbA1cを0週、12週、24週、52週と繰り返し測る研究では、4つの値が独立しないため相関をモデルで表す必要があります。この表し方を指定するのが共分散構造で、解析計画にはほぼ必ず構造名が書かれます。
ところが、構造を選び替えると何が変わるかは解析計画の段階では見えにくく、点推定値がほとんど動かない一方で、標準誤差と信頼区間の幅、検定の水準は構造の指定に左右されます。
そこで、この3つを同じ反復測定データへ当てはめ、当てはめられた相関、情報量規準、群間差の標準誤差がどう変わるかを数値で並べます。
反復測定データで共分散行列が担う役割
この試験では各被験者に4つの観測値が並び、平均構造は来院ごとの群平均で表せますが、残差の相関をどう置くかはそれとは別に決めることになります。
ここで、被験者$i$の観測値ベクトルを$\mathbf{y}_i$、計画行列を$\mathbf{X}_i$、固定効果を$\boldsymbol{\beta}$と書くと、モデルは次の形です。
$$\mathbf{y}_i = \mathbf{X}_i\boldsymbol{\beta} + \boldsymbol{\varepsilon}_i, \qquad \boldsymbol{\varepsilon}_i \sim N(\mathbf{0}, \mathbf{V}_i) \tag{1}$$
式(1)は、平均部分を$\mathbf{X}_i\boldsymbol{\beta}$、被験者内のばらつきと相関を$\mathbf{V}_i$で表します。$\mathbf{V}_i$は4行4列の対称行列で、対角に各来院の分散、非対角に来院間の共分散が入ります。
この$\mathbf{V}_i$へどんな形を課すかが共分散構造です。相関を無視して対角行列と置くと、繰り返し得た情報を独立として数え、群間差の標準誤差が実際より小さく出ます。
そのため一般化最小二乗法では$\mathbf{V}_i$の逆行列で重みを付けて$\boldsymbol{\beta}$を推定します。Laird and Ware(1982)は、ランダム効果が相関を生む仕組みを定式化し、$\mathbf{V}_i$を平均構造と切り離す術を検討しました。
なお、$\mathbf{V}_i$の指定は生存アウトカムと結びつける場面では生存モデル側にも影響しますが、そちらは縦断データと生存データの結合モデルで扱っているので、本記事は線形モデルに絞ります。
compound symmetryが置く一様な相関
最も単純な構造がcompound symmetryで、どの来院の分散も等しく、どの2つの来院の相関も同じ値だと仮定します。
$$\mathbf{V}_i = \sigma^2\begin{pmatrix} 1 & \rho & \rho & \rho \\ \rho & 1 & \rho & \rho \\ \rho & \rho & 1 & \rho \\ \rho & \rho & \rho & 1 \end{pmatrix} \tag{2}$$
式(2)は、4来院とも残差分散が$\sigma^2$、組み合わせによらず相関が$\rho$と置いた形で、推定するのはこの2つだけです。
この単純さには由来があり、ランダム切片だけの線形混合モデルは周辺で見ると式(2)と同じ共分散を与えます。Laird and Ware(1982)の枠組みでは、被験者間のばらつきを1つの切片に集約した場合にあたります。
ランダム切片を入れれば相関は扱えている、と受け取られがちですが、それは0週と52週の相関が12週と24週の相関に等しいと宣言するのと同じで、時間が離れるほど関連が弱まる性質を表現できません。
実際、今回のデータではcompound symmetryが返す相関が組み合わせによらず0.393で一定でした。観測された相関は0.562、0.306、0.127と減っており、近い来院を低く、遠い来院を高く見積もる形です。
AR(1)構造が表す時間による減衰
時間が離れるほど相関が弱まる性質を2つのパラメータで表すのがAR(1)です。Diggle(1988)は反復測定の共分散を被験者間の変動、系列成分、測定誤差へ分解し、系列成分にこの形を置く定式化を提案しました。
$$\operatorname{Cov}(\varepsilon_{ij}, \varepsilon_{ik}) = \sigma^2 \rho^{|j-k|} \tag{3}$$
式(3)は、共分散がズレ$|j-k|$の1増ごとに$\rho$倍ずつ小さくなることを表します。$\varepsilon_{ij}$は残差、$\rho$は隣り合う来院の相関で、$|\rho| < 1$の範囲を動きます。
今回は$\rho$が0.580と推定され、当てはめ相関は0.580、0.336、0.195と減衰して観測値0.562、0.306、0.127に近い値でした。赤池情報量規準1004.341、ベイズ情報量規準1045.911は、compound symmetryの1052.707と1094.277を下回ります。
ただし、見落としやすいのが来院間隔の扱いで、今回の0週、12週、24週、52週は等間隔ではありません。来院番号をそのまま$j$とすると、24週から52週の28週間を12週から24週の12週間と同じ1期分のズレとして扱うことになります。
この点については、実測の週数を距離に使う連続時間版のAR(1)もあり、来院間隔が不揃いな試験ではそちらが設計に合います。Littell et al.(2000)も、時間の間隔をどう構造へ反映するかが選択の要点だと述べています。
unstructuredが仮定を置かない代償
unstructuredは相関の形を何も仮定せず、4来院なら分散4つと相関6つの計10個を自由に推定します。当てはめた相関は標本相関とほぼ重なり、今回も0.594、0.549、0.544となりました。
自由に動ける分だけ当てはまりは良く、対数尤度は-488.314と最も高い値です。もっとも赤池情報量規準は1012.628とAR(1)の1004.341を上回り、ベイズ情報量規準では1087.454まで離れました。
しかも代償はパラメータ数だけではありません。Kenward and Roger(1997)は、共分散パラメータの推定の不確実性が固定効果の分散へ反映されず、標本が小さいほど標準誤差が小さめに出ることを扱い、調整法を提案しています。
実際、パラメータ数の多さは標準誤差の振れに表れ、unstructuredの群間差の標準誤差は0.1518、0.1286、0.1324、0.1436と来院ごとに違いました。全来院で0.1394に固定されるcompound symmetryとは対照的です。
この振れ幅は、来院ごとに分散が本当に違う研究では拾うべき情報ですが、120例で10個を推定する今回は、差の一部が推定の揺らぎにすぎません。
同じデータへ3つの構造を当てはめた結果
今回のデータは各群60例の計120例、4回の来院で480観測、残差標準偏差0.80%、隣接来院の相関0.60のAR(1)から作りました。群間差は0週で0、12週で-0.15、24週で-0.35、52週で-0.50%と置いてあります。
3つを比べると赤池情報量規準はcompound symmetry 1052.707、AR(1) 1004.341、unstructured 1012.628、ベイズ情報量規準でもAR(1)が最小でした。生成に使った構造がAR(1)なので期待どおりです。
どちらもunstructuredの特別な場合なので尤度比検定でも比較でき、対立側をunstructuredとすると、compound symmetryとの統計量は自由度8で56.079、p値0.001未満、AR(1)とは7.713、p値0.462でした。
後者のp値が大きいことをもって、AR(1)で十分だと確認できたと書くのは行き過ぎです。この検定が見ているのは追加の8パラメータが説明を増やすかどうかであって、有意でないという結果は、今回のデータでAR(1)からの明確な逸脱が確認されなかったという以上のことを意味しません。
一方、52週の群間差の点推定値は3つとも-0.511%で同じでした。標準誤差だけが0.1394、0.1411、0.1436と広がり、p値は0.000272、0.000322、0.000407と並びます。
この点推定値が動かなかったのは、欠測がなく平均構造を飽和させたためで、この条件ではどの重み付けでも群平均は標本平均に一致します。観測数が揃わない研究では、$\mathbf{V}_i$が欠測した来院の情報を他からどう借りるかを決めるため、点推定値自体が構造で動きます。
構造の指定が標準誤差をどちらへずらすか
構造の誤指定が標準誤差をどちらへ外すかは、同じ設定でデータを作り直して確かめられます。真の構造をAR(1)に固定したまま300回生成し、52週の群間差について、各構造が返すモデル由来の標準誤差と、300個の推定値が実際に散らばった幅を比べました。
実際、推定値の平均は3つとも-0.487%、実際のばらつきは0.1483でした。モデルが返した標準誤差の平均は0.1463と0.1464で、どちらも実際のばらつきの0.987倍にとどまります。
さらに、95%信頼区間が真値-0.50%を含んだ割合も0.940と0.933で、名目の0.95を下回りました。unstructuredの被覆が最も低い点は、Kenward and Roger(1997)の指摘と整合します。
ただし、構造の誤指定は実害が小さいという一般則にはなりません。今回は欠測がなく例数が等しく来院4回という条件が揃っており、Littell et al.(2000)は、来院数が多い場合には水準がより大きく歪みうると報告しています。
そのため私は構造を決める前に、来院ごとの残差標準偏差とズレごとの標本相関を並べます。今回は残差標準偏差が0.829から0.783と大きく違わず、相関は0.562、0.306、0.127と落ちていたので、時間による減衰の方に注目しました。
臨床試験の解析計画にどう書くか
連続量を主要評価項目とする試験では、Mallinckrodt et al.(2008)が既定の解析にunstructuredを置き、収束しない場合の切り替え順序を解析計画へ事前に書く進め方を勧めています。順序を固定すれば、結果を見てから選び直す余地がなくなります。
この方針に沿って推奨を条件付きで書きます。来院が5回程度までで、欠測が中程度にとどまり、各群の例数が数十以上ある試験なら、unstructuredを既定の解析に置き、Kenward and Roger(1997)の自由度調整を併用する形を勧めます。
一方、来院が10回あるとunstructuredのパラメータは55個に達し、収束しないか標準誤差が信頼しにくくなります。この規模では、Toeplitz構造やAR(1)のようにパラメータ数が来院数に比例する構造が現実的です。
逆に勧めないのは、複数の構造を当てはめて情報量規準が最小の結果を主要な解析として報告する順序です。今回も赤池情報量規準ならAR(1)、事前規定ならunstructuredとなり、52週のp値は0.000322と0.000407で結論は変わりませんが、境界付近なら選び方が有意性を左右します。
また、構造の名前、代替の順序、自由度の計算方法まで書いておくと、査読で妥当性を問われた時に答えられます。Littell et al.(2000)の指摘どおり構造の指定は検定結果を左右するので、構造名のない縦断解析は報告として不十分です。
Rで3つの構造を当てはめる
ここまでの数値は、base Rとnlmeパッケージで出したものです。乱数から作った反復測定データへ3つの構造を当てはめる部分だけを載せます。
library(nlme)
set.seed(42)
lag_matrix <- abs(outer(1:4, 1:4, "-"))
true_cov <- 0.80^2 * 0.60^lag_matrix # 真の共分散はAR(1)
errors <- matrix(rnorm(120 * 4), nrow = 120) %*% chol(true_cov)
mean_matrix <- matrix(8.00, nrow = 120, ncol = 4)
mean_matrix[61:120, ] <- mean_matrix[61:120, ] +
matrix(c(0, -0.15, -0.35, -0.50), nrow = 60, ncol = 4, byrow = TRUE)
hba1c_matrix <- mean_matrix + errorsここで用意したのは、両群60例ずつ4来院分のHbA1cです。治療群の平均だけを0、-0.15、-0.35、-0.50%ずらし、被験者内の相関は隣接来院0.60のAR(1)から作って、縦長データをrepeated_dataとしました。
# repeated_data は被験者120人×4来院を縦に並べ、visit_numに来院番号を持たせたもの
mean_formula <- hba1c ~ visit + visit:arm - 1
fit_cs <- gls(mean_formula, data = repeated_data,
correlation = corCompSymm(form = ~ visit_num | subject), method = "REML")
fit_ar1 <- gls(mean_formula, data = repeated_data,
correlation = corAR1(form = ~ visit_num | subject), method = "REML")
fit_un <- gls(mean_formula, data = repeated_data,
correlation = corSymm(form = ~ visit_num | subject),
weights = varIdent(form = ~ 1 | visit), method = "REML")
c(cs = AIC(fit_cs), ar1 = AIC(fit_ar1), un = AIC(fit_un))#> cs ar1 un
#> 1052.707 1004.341 1012.628これら3つの当てはめでは平均構造を同じ飽和した形に揃えたので、赤池情報量規準の差は共分散構造だけの差になります。数値が小さいAR(1)が良いという結論はこのデータについての話で、別の試験で同じ順序になる保証はありません。
round(sapply(list(cs = fit_cs, ar1 = fit_ar1, un = fit_un),
function(fit) summary(fit)$tTable["visit4:armtreatment", 1:2]), 5)#> cs ar1 un
#> Value -0.51146 -0.51146 -0.51146
#> Std.Error 0.13941 0.14114 0.14364この出力では、52週の群間差が3つとも-0.51146%で一致し、標準誤差だけが0.13941から0.14364まで動きました。点推定値の一致を構造の選択が結果に影響しない根拠と受け取ると読み違いになり、影響は標準誤差と信頼区間の幅に出ています。
by_lag <- function(cor_matrix) sapply(0:3, function(k) mean(cor_matrix[lag_matrix == k]))
fitted_cor <- function(fit)
matrix(as.numeric(corMatrix(fit$modelStruct$corStruct)[[1]]), nrow = 4)
observed_cor <- cor(hba1c_matrix - mean_matrix)
plot(0:3, by_lag(observed_cor), type = "b", pch = 16, ylim = c(0, 1),
xlab = "Lag between visits", ylab = "Correlation")
lines(0:3, by_lag(fitted_cor(fit_cs)), type = "b", pch = 1, lty = 2)
lines(0:3, by_lag(fitted_cor(fit_ar1)), type = "b", pch = 2, lty = 3)
lines(0:3, by_lag(fitted_cor(fit_un)), type = "b", pch = 5, lty = 4)最後に、ズレごとに相関を平均し観測値と3つの当てはめを重ねたものが次の図の左側、右側は300回生成し直したときのモデル標準誤差と推定値の実際のばらつきの比です。
左の図で見てほしいのは、compound symmetryの線だけが水平な点です。AR(1)とunstructuredは減衰を追えていますが、ズレ3ではAR(1)が0.195、観測値が0.127と離れ、べき乗の形は遠い来院で減衰を緩めに見積もります。
右の図では3つの棒がどれも1をわずかに下回って揃い、この設定では構造の違いより、共分散を推定したこと自体による標準誤差の過小評価の方が大きく効いていたと読めます。
反復測定の解析計画で先に決めること
共分散構造は平均構造の飾りではなく、群間差の標準誤差と検定の水準を決める部分です。今回のデータでは点推定値が-0.511%で動かない一方、標準誤差は0.1394から0.1436まで広がり、欠測がある研究なら点推定値そのものも構造で動きます。
そのうえで判断を分けるのは、来院回数と例数の比、来院間隔、来院ごとの分散の違いです。来院が少なく例数に余裕があればunstructuredを既定に、来院が多い設計では減衰を仮定する構造へ移し、切り替え順序を事前に解析計画へ書く形を勧めます。
避けたいのは、当てはめた複数の構造の中から、結果を見て都合のよいものを選ぶ順序です。
ここまで読んで、自分の試験では来院回数と欠測の程度からどの構造を既定にすべきか、と考えられた方もいるかもしれません。実際の判断では、残差のばらつき、欠測の起き方、解析計画の記載、自由度調整の有無まで含めて決めることになり、Dr.データサイエンスでは、こうした解析方針のご相談を承っています。
参考文献
- Laird NM, Ware JH. Random-effects models for longitudinal data. Biometrics. 1982;38(4):963-974.
- Diggle PJ. An approach to the analysis of repeated measurements. Biometrics. 1988;44(4):959-971.
- Littell RC, Pendergast J, Natarajan R. Modelling covariance structure in the analysis of repeated measures data. Statistics in Medicine. 2000;19(13):1793-1819.
- Kenward MG, Roger JH. Small sample inference for fixed effects from restricted maximum likelihood. Biometrics. 1997;53(3):983-997.
- Mallinckrodt CH, Lane PW, Schnell D, Peng Y, Mancuso JP. Recommendations for the primary analysis of continuous endpoints in longitudinal clinical trials. Drug Information Journal. 2008;42(4):303-319.
