欠測を含む順序カテゴリカル経時データの解析 -GEE プロシジャの有用性 - 駒嵜弘 藤原正和 2 ( マルホ株式会社 2 塩野義製薬株式会社 ) Ordinal longitudinal data analysis with missing data -Usefulness of Proc GEE- Hiroshi Komazaki,Masakazu Fujiwara 2 Maruho Co, Ltd., 2 Shionogi & Co., Ltd.
2 発表の構成 カテゴリカル経時データの発生 前発表 欠測データの発生 多重補完 (REG MCMC FCS) 本発表 解析 (CATMOD GEE) 2
要旨 : 順序カテゴリカル経時データの解析方法として GEE プロシジャによる方法及び CATMOD プロシジャによる方法をそれぞれ紹介するとともに シミュレーションデータを用いて性能評価を行う キーワード : CATMOD プロシジャ GEE プロシジャ Generalized Estimating Equations (GEE) alternating logistic regression(alr) 3
Outline 順序カテゴリカル経時データの解析仕様 ( CATMODプロシジャ GEEプロシジャ ) シミュレーション結果 考察 4
Background 及び Motivation 順序カテゴリカル経時データに対する解析方法 PROC 特徴 CATMOD 重み付き最小二乗法 ( スライド の式参照 ) GENMOD GEE GEE による推定相関構造は IND しか使用できず ( ただしロバスト分散の使用が可 ) GEE による推定相関構造は EXCH が指定可能 ( スライド 8 4 の式参照 ) SAS 9.4(stat 4.) より PROC GEE は応答変数が順序カテゴリカルデータに対して も実行できるよう DIST=MULTINOMIAL が さらに相関構造も指定できるよう LOGOR= EXCH が追加された PROC CATMOD と PROC GEE を用いてカテゴリカル経時データの解析方法を 紹介し さらに性能評価を行う 5
順序カテゴリカル経時データにおいて CATMOD プロシジャでの解析方法と GEE プロシジャでの解析方法は Liang and Zeger (986) 及び Prentice(988) の GEE 法に類似 Liang and Zeger (986) 及び Prentice(988) の GEE 法との相違点 CATMOD プロシジャ 重み付き最小二乗法分散構造は saturated correlation structure GEE プロシジャ 相関パラメータ α の推定方程式を改良 Liang and Zeger (986) 及び Prentice(988) の GEE 法を提示後 それぞれのプロシジャでの推定方程式を紹介する 6
* 順序カテゴリカルデータを表記する上での留意点 本発表にて 順序カテゴリカルデータは指示変数 (indicator variable) を用いて表記 例 )5 カテゴリの応答変数の場合 Yのカテゴリ指示変数 2 3 4 5 y ( O c) [,2,3,4 ] = I c ( y ) y 2 y 3 y 4 ( 0 0 0) ( 0 0) ( 0) ( ) ( 0 0 0 0) P : 群 i : 被験者 t : 時点 7
順序カテゴリカルデータのGEE( Liang and Zeger (986) Prentice(988) ) ( 主効果 β の推定方程式 ) U ( β ) D V ( α) ( Y ) = 0 = = 0 0k β X : カテゴリ k の切片 : デザイン行列 ( 群 時点 群 時点 ) Y : y k D 各群 (=0,) 各症例 各時点 (t=,2,3,4) の指示変数 例 : 症例 i の指示変数 Y の期待値 β φ y = { y y y y } 2 3 4 ) = ln ( k ) = β 2 k k = L( k 0k 5 累積ロジットモデル X β V (α ) 各群 (=0,) の各時点 各カテゴリ間の分散共分散行列 8
9 5 4 3 2 ln β β X = 0 ln ロジットモデル累積ロジットモデルカテゴリ数 :5 の場合 β β X k k k = 0 5 ) ( 2 ln 5 4 3 2 ln 5 4 3 2 ln 5 4 3 2 ln 各カテゴリにて 共変量は同じと仮定 比例オッズモデル β
順序カテゴリカルデータのGEE ( 相関パラメータ α の推定方程式 ) U ( α ) E W ( Z η ( α)) = 0 = = 0 Z η E 相関係数 η α ( ) Z Z の期待値 ( tk )( t k ) α η = ( y k k )( y k k [ ( ) ( )] 2 k k k k ex( f α ) i( tk )( t k ) ( α) = i( tk ) ( t k ) ex( f α ) i( tk ) f ( t k ) ( i( tk t k X と X の関数 ) )( ) : ) Miller(993) W Z i( tk )( t k ) の分散共分散行列 0
順序カテゴリカルデータの GEE 推定量 ( ) = = = 0 ) ( 0 ) ( t t t t t t X G X X G X φ β β CATMOD プロシジャの推定量分散共分散構造が saturated correlation structure* のとき CATMOD の推定量は 上式の推定量の赤線部分に対応している (Miller,993) *saturated correlation structure 群間 時点間 カテゴリ間で異なる相関パラメータ α を指定 β ( ) = t t t t V n G φ α φ ) ( β
PROC GEEではHeagerty and Zeger(996) の方法を使用 Liang and Zeger (986) Prentice(988) のGEE β の推定方程式 ( スライド 8) 相関パラメータ α の推定方程式 ( スライド0) Heagerty and Zeger(996) の GEE β の推定方程式 ( スライド 8) 相関パラメータ α の推定方程式 ( スライド4) ξ, c2 ) ( y y ) = E i( t, t )( c itc it c 2 2
GEE プロシジャの順序カテゴリカルデータの GEE 法は Heagerty and Zeger(996) の手法に基づいている 相関パラメータ算出の際 alternating logistic regression(alr) を使用 β 主効果の推定方程式は Liang and Zeger (986) Prentice(988) の GEE と同じ α 相関パラメータは 時点間 カテゴリ間の条件付き期待値を利用して推定 条件付き期待値の回帰パラメータが対数オッズ比となる ( スライド 7 参照 ) ξ, c2 ) ( ) y y = E i( t, t )( c itc it c 2 3
ALRによる相関パラメータ α の推定方程式 U ξ = = 0 α t ( ) * α V ( Y ξ ) = 0 y * ( y,, y ) = 4 ( ) n i 4 各群 (=0,) の各時点 (t=,2,3,4) の指示変数 各データを ( 全時点 - 当該時点 ) 個分 ( カテゴリ数 -) 個分 重複生成している ( 詳細は次ページ参照 ) ξ, c2 ) ( y y ) = E ( t, t )( c tc t c 2 時点 t でカテゴリ c2 が得られた条件のもとで 時点 t でカテゴリ c が得られる期待値 (t >t) V i = diag [ ξ ξ )] i ( i 4
例 : 説明の簡略化のため時点数 t:3 カテゴリ数 k:3 を想定 以下の 2 時点のデータが得られたとする ( 時点 3 のデータは説明変数で使用し 応答変数には用いないため省略 ) 時点 時点 2 yi y i 2 カテゴリ 2 指示変数 (,0) (,) n i : 被験者 iの時点数 y 各時点 各カテゴリごとの条件付き期待値を求めるにあたり 応答変数は以下の通りデータを重複発生する必要がある * i n n 2 i ( y,, y, y,, y,, y ) = i ( k ) i ( k ) i2 ( k ) i2 ( k ) i( n ) ( k ) y = = [ yi 2 yi 2 yi2 2 ] [(,,0,0) (,,0,0) (,,,) ] * i i i 5
例 : 各時点 各カテゴリごとの条件付き期待値の構造は以下の通り ξ = = i [ ξ ξ ξ ] 2 [ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ] 2() 3 2(2) 23 2(2) 従って本例の場合 残差 y = 2(22) 3() 3(2) = [ yi 2 yi 2 yi2 2 ] [(,,0,0) (,,0,0) (,,, )] * i ξ = = [ ξ ξ ξ ] 3(2) 3(22) 時点 2 で が得られた条件のもと 時点 で が得られる期待値 i 2 23() は 2 個のデータで構成される [ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ] 2() 3 2(2) 23 2(2) y * 2(22) ξ 3() 3(2) 3(2) 3(22) 23() 23(2) 23(2) 23(2) 23(2) 23(22) 23(22) 6
ξ ( y y ) = E ( t, t )( c tc t c, c2 ) 2 のモデル [ E( Y Y )] log it = itc it c γ i( tt )( c c ) Y 対数オッズ比二項目は切片 it c log 2 2 2 µ µ itc itc E Y µ it c ( ) itc Y it c2 E( Y Y ) 2 itc it c 2 γ i( tt )( c c ) = log ρ ( ( ) t i( tt )( cc 2 ) OR Y = = itc, Y it c Z 2 i( tt )( cc 2 ) α log ρi( tt )( c ) c2 2 時点間 カテゴリ間の対数オッズ比より相関パラメータ ρ i t Z i = を推定 相関係数 ( tt )( c c2 ) が求まる ( のとき相関構造はexchangeable) α 数式の詳細は Heagerty and Zeger(996) Lisitz(99) SAS HELP 参照 7
シミュレーションの流れ カテゴリカル経時データの発生 欠測データの発生 多重補完 (REG MCMC FCS) 解析 (CATMOD GEE) 8
シミュレーションデータ Donneau(205) 及び Ibrahim and Suliadi(20) を参考にシミュレーションデータを作成 [ Pr( Y k x, t) ] logit = β0 β x β I( t = j) β x I( t j) i y x t ij i k x i t tx i = j : 被験者 : 時点 (,2,3,4) : 応答変数のカテゴリ (,2,3,4) : 応答変数 : 投与群 (0: プラセボ : 実薬 ) : 時点 (,2,3,4) シミュレーション回数 :000 回 k β 0 β 02 β 03 β 04 β x β t β xt パラメータ真値 -.39-0.4 0.4.39 0.5-0.4-0.2-0. 0-0.4-0.2-0. 0 9
シミュレーションデータ各時点 各カテゴリの周辺分布 プラセボ群 VISIT K 2 3 4 4 7 8 20 2 6 8 9 20 3 9 20 20 20 4 23 2 2 20 5 27 23 22 20 % 実薬群 % VISIT K 2 3 4 6 22 25 29 2 7 2 22 23 3 20 20 20 9 4 22 9 7 6 5 25 8 6 3 20
シミュレーションデータ 観測確率モデル logit( ) Y t = 2 0. t 時点前のデータが欠測の有無に影響 各群 最終時点までに 40% 脱落 2
CATMOD プロシジャの指定方法 データセットの構造 時点 ~ 時点 4 のカテゴリのパターン sim Imutation Number arm 2 _3 _4 count 0 0.0 0 2 0 3 0.0 0 4 0.0 0 5 : : : : : : : : 000 20 5 5 5 5 0.0 各時点の応答変数パターンを度数で格納度数 0 のパターンは代わりに 0.0 を格納 ( 値の妥当性は未検証 ) 分散共分散行列のランク落ち 累積ロジスティックモデルの準完全分離を防ぐため ERROR: Logits cannot be comuted. One of the marginal robabilities is equal to zero. ERROR: The covariance matrix of the resonse functions is singular in oulation. 22
CATMOD プロシジャの指定方法 roc catmod data=in_data ; resonse clogits; oulation arm; weight count ; model _*_2*_3*_4 = ( 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0, ~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~ 0 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0, 0 0 0 0 0 0 0, 0 0 0 0 0 0 0, 0 0 0 0 0 0 0, ~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~ 累積ロジットモデルの指定投与群の指定度数変数の指定 列 ~4 : ( 応答変数のカテゴリ数 -) 個の切片 列 5 : 投与群 列 6~8 : 時点 列 9~ : 投与群 時点の交互作用 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0, 0 0 0 0 0 0 0 0 0 ) ( 2 3 4 = 'intercet', 5='arm',6 7 8='VISIT', 9 0 ='arm*visit') ; run; 23
GEE プロシジャの指定方法 データセットの構造 sim Imutation Number arm ID VISIT AVAL 0 0 2 2 0 3 2 0 4. : : : : : : 000 20 00 4 4 応答変数の順序カテゴリカルデータ (AVAL:,2,3,4,5) 被験者ごと VISIT ごとに応答変数データを格納 24
GEE プロシジャの指定方法 roc gee data=in_data; class arm ID; model AVAL_M =visit arm visit*arm / dist=multinomial; reeated subject=id / logor=exch; run; 相関構造に exchangeable を指定 累積ロジットモデルの指定 arm : 投与群 (0 or ) visit : 時点 (,2,3,4) Visit*arm : 投与群 時点の交互作用 25
結果 : モデルを正しく特定 ver CATMOD 真値推定値 (MSE) BIAS RB αエラー MCMC 0.5 0.48 (0.078) 0.082 83.6 2.80 REG 0.5 0.436 (0.090) 0.064 87. 2.50 FCS 0.5 0.433 (0.09) 0.067 86.6 2.50 補完なし 0.5 0.420 (0.092) 0.08 84.0 2.70 完全データ 0.5 0.450 (0.059) 0.05 90. 2.50 GEE 000 ( µ ) 真値推定値 (MSE) BIAS RB αエラー MCMC 0.5 0.460 (0.086) 0.040 92.0 2.20 REG 0.5 0.479 (0.095) 0.02 95.8 2.60 FCS 0.5 0.477 (0.095) 0.023 95.4 2.00 補完なし 0.5 0.490 (0.094) 0.00 98.0 2.40 完全データ 0.5 0.497 (0.059) 0.003 99.3 2.30 000 i 000 µ µ 000 i µ 26
結果 2: 相関構造誤特定 ver CATMOD 真の相関構造 真値推定値 (MSE) BIAS RB αエラー MCMC 0.5 0.448 (0.080) -0.052 89.6 2.0 REG 0.5 0.462 (0.085) -0.038 92.5.50 FCS 0.5 0.45 (0.086) -0.049 90.2.80 補完なし 0.5 0.452 (0.090) -0.048 90.3.80 完全データ 0.5 0.47 (0.058) -0.029 94.3 3.00 0.6 0.4 0.2 0.6 0.6 0.4 0.4 0.6 0.6 0.2 0.4 0.6 GEE 真値推定値 (MSE) BIAS RB αエラー MCMC 0.5 0.493 (0.094) 0.007 98.5 2.60 REG 0.5 0.509 (0.098) 0.009 0.8.70 FCS 0.5 0.497 (0.097) 0.003 99.3.40 補完なし 0.5 0.524 (0.08) 0.024 04.7 2.60 完全データ 0.5 0.58 (0.065) 0.08 03.6 2.70 27
結果 3: 欠測メカニズムがMNAR ver CATMOD 真値推定値 (MSE) BIAS RB αエラー MCMC 0.5 0.45 (0.20) -0.085 82.9 2.00 REG 0.5 0.426 (0.46) -0.074 85. 2.60 FCS 0.5 0.422 (0.47) -0.078 84.5 2.50 補完なし 0.5 0.396 (0.37) -0.04 79.3 2.50 GEE 真値推定値 (MSE) BIAS RB αエラー MCMC 0.5 0.463 (0.37) -0.037 92.6 2.60 REG 0.5 0.475 (0.58) -0.025 94.9 3.00 FCS 0.5 0.473 (0.57) -0.027 94.6 2.30 補完なし 0.5 0.498 (0.58) -0.002 99.6 2.50 28
結果のまとめ モデルを正しく特定できたとき 性能 : GEE CATMOD GEEでは 補完によりバイアス増加 相関構造を誤特定したとき 性能 : GEE > CATMOD GEEの頑健性を確認できた GEEでは 補完によりバイアス減少 欠測メカニズムを誤特定したとき 性能はGEE> CATMOD 相関構造の誤特定よりはバイアス小 MSEはいずれも増加 GEEでは 補完によりバイアス増加 α エラーの増大 (>>2.5) は特に確認できなかった 29
発表のまとめ カテゴリカル経時データの解析方法の紹介 (CATMOD GEE) GEEプロシジャでは相関パラメータαを対数オッズ比より推定 性能 : GEE>CATMOD GEEでは補完した方がバイアスが増大することもあった 累積ロジットモデルのとき Weighted-GEEによる解析はSAS のPROC GEEでは実装されていないため 今後の導入が期待される 30
参考文献 Donneau, A. F., Mauer, M., Molenberghs, G., & Albert, A. (205). A simulation study comaring multile imutation methods for incomlete longitudinal ordinal data. Communications in Statistics-Simulation and Comutation, 44(5), 3-338. Heagerty, P. J., & Zeger, S. L. (996). Marginal regression models for clustered ordinal measurements. Journal of the American Statistical Association, 9(435), 024-036. Ibrahim, N. A., & Suliadi, S. (20). Generating correlated discrete ordinal data using R and SAS IML. Comuter methods and rograms in biomedicine, 04(3), e22-e32. Liang, K. Y., & Zeger, S. L. (986). Longitudinal data analysis using generalized linear models. Biometrika, 73(), 3-22. Lisitz, S. R., Laird, N. M., & Harrington, D. P. (99). Generalized estimating equations for correlated binary data: using the odds ratio as a measure of association. Biometrika, 78(), 53-60. Miller, M. E., Davis, C. S., & Landis, J. R. (993). The analysis of longitudinal olytomous data: Generalized estimating equations and connections with weighted least squares. Biometrics, 033-044. Prentice, R. L. (988). Correlated binary regression with covariates secific to each binary observation. Biometrics, 033-048. 3