本日の内容 前回と同じ本 G.Verbeke, and G.Molenberghs. 著 Linear Mixed Models for Longitudinal Data (Springer, 以下 テキスト と呼びます ) の 10 章をまとめます 誤差が独立でない場合を扱います 参考文献は 今

Transcription

1 経時データのモデリング (2) 第 11 回 BioS 継続勉強会 土居正明 1

2 本日の内容 前回と同じ本 G.Verbeke, and G.Molenberghs. 著 Linear Mixed Models for Longitudinal Data (Springer, 以下 テキスト と呼びます ) の 10 章をまとめます 誤差が独立でない場合を扱います 参考文献は 今回と次回の内容を合わせたものです 2

3 本日の内容 はじめに ( データの紹介 モデリングの基礎 ) 10.1 Introduction 10.2 An Informal Check for Serial Correlation 10.3 Flexible Models for Serial Correlation 10.4 The Semi-Variogram 10.5 Some Remarks まとめと参考文献 3

4 注意 テキストは 記号の使い方等が結構いい加減です ( 確率変数とその実現値の使い分け など ) 数理的な部分で たまに間違い ( 誤植と呼ぶには大きいもの ) もあります 本資料ではテキストより記号等を変更しています 大抵は意図的に 修正 していますので 大体はこっちの方が正しいです ( が 誤植は結構あるはずです ) LMM は線形混合効果モデル (Linear Mixed Model) の略です 分散成分の推定には 特に断りがなければ REML 法を用いています 4

5 はじめに ( データの紹介とモデリングの基礎 ) 5

6 はじめに 症例の個の時点のデータが 1 症例の経時変化 である 10 章のテーマは 1. 時点ごとのデータの相関 ( 系列相関 ) の有無の判断 2. 系列相関がある場合はそのモデリングである 6

7 本書を通して大事なこと 目的は探索解析 本当のモデルは分からないので 探索しましょう 使える モデルをみつけるための 試行錯誤 の方法論を考えましょう Essentially, all models are wrong, but some are useful. (by G.E.P. Box) 7

8 本質的問題と暫定的解決 時点ごとの相関が存在するとして どういう関数で表されるか を正確に知るのは一般に不可能 複雑な関数になることも想定されるが 通常 指数関数を用いたモデリングで十分なことが多い 指数関数は 本質的に時点ごとの相関関係を表現している 訳ではなく ある種の微分方程式の解として自然に現れる 現実的に いくつかのデータに対する Fitting が悪くないことなどから 便宜的に使用される True Model ではなくて Useful Model の探索 8

9 統計で最も重要なこと ( 私見 ) 全ての tool は不完全である プロットして目視 尤度比検定 情報量規準など どれも 決定的 な知識はくれない これさえしておけば絶対大丈夫 という人は信用してはいけません しかし 不完全だからいい加減でよい わけではない むしろ 不完全なので 一層注意深く観察 検討しよう という態度が重要 最終モデルに対しても もしかしたら不完全かも と思うこと ただし 現在の知識 情報 技術ではここが限界 というまで検討することが大事 9

10 モデリングとは ( 私見 ) 全く分からない から入って 少しは分かったかも で終わるもの 完全に分かる = 真のモデルが分かる は ( 人間の力では ) あり得ない 情報がさらに蓄積すれば 構築されるモデルは変わりうる あまり自信のない判断をせざるを得なくなるときもあるが その時は ここで微妙な判断をした ということを忘れないことが重要 自分のモデルを冷静に評価すること 10

11 経時データのモデリング 決めなければいけないものがたくさんある ( は 9 章までのテーマ ) 平均構造 どの要因を選ぶか? 変量効果の分散構造 誤差の相関構造 10 章のテーマ (3.8) 11

12 今回のモデルの 固定効果 変量効果 誤差 固定効果 : に依存しない 全症例 投与群など 複数の個人に共通の部分 変量効果 : に依存する 被験者特有の部分 被験者内変動 誤差 : 固定効果 変量効果で説明できない データのばらつき 被験者間変動 12

13 症例は固定効果か変量効果か? どちらでも解析することは可能 症例内のデータの相関を考慮した解析を行うため 変量効果を用いることが多い テキストでも 変量効果として扱う 13

14 モデリングで決めなければいけないこと どの変数が入るか? 時間 性別 etc どういう風に入るか? 時間の 1 次関数?2 次関数? 変数はそのままで OK? 変換してから入れる? 誤差は独立? 相関あり? 変量効果同士の相関は? 14

15 本章であつかう実データの紹介 (Prostate Data: 前立腺癌のデータ ) 詳しくは 節参照 前立腺癌は アメリカでは男性の癌による死亡の 2 番目の原因 治療にお金もかかる 早期発見が重要 PSA (prostate-specific antigen) がマーカーになる PSA は正常細胞にも癌性前立腺細胞にも含まれる酵素 前立腺組織の体積と関係がある 15

16 本章であつかう実データの紹介 (Prostate Data: 前立腺癌のデータ ) BPH(benign prostatic hyperplasia, 良性前立腺過形成 ) でも PSA が大きくなる Pearson et al. (1991) によると PSA の値だけで判断した場合 最大 60% の BPH 患者が前立腺癌と誤診される 現状では PSA は前立腺癌のマーカーとして不十分 前立腺癌だけを検出する判断基準を導けるようなモデルを作りたい ( 目的 ) 16

17 Prostate Data 個人ごとのデータの推移 ( 図 2.3) ( 横軸は診断前の時間 : 右が過去 ) これを区別するモデリングをしたい テキスト 13P より引用 17

18 本章であつかう実データの紹介 (Prostate Data: 前立腺癌のデータ ) 単位時間あたりの変化 : 正しい意味での 率 仮説 PSA の変化率を見れば 前立腺癌の早期発見ができるのではないか? PSA の数値だけではなくて経時変化の仕方までよく見れば よりよい前立腺癌のマーカーになるのでは? 時間に依存した部分に注目 18

19 本章であつかう実データの紹介 (Prostate Data: 前立腺癌のデータ ) BLSA (Baltimore Longitudinal Study of Aging) のデータ Pearson et al. (1994) 参照 デザイン 後ろ向き Case-Control 研究 凍結させた血清サンプルを使用 被験者の内訳は 前立腺癌 :18 例 局所浸潤性癌 (L/R Cancer):14 例 転移性癌 (Metastaic Cancer):4 例 良性前立腺過形成 (BPH):20 例 対照群 (Control, 前立腺癌の兆候なし ):16 例 19

20 本章であつかう実データの紹介 (Prostate Data: 前立腺癌のデータ ) 選択基準 1. 泌尿器科医によって 前立腺癌 BPH による単純前立腺摘出術 前立腺の病気はない と診断されるまでに 7 年以上の追跡調査のデータがある 2. 病理学的診断により確認されている 3. 診断の前に前立腺の手術がない 20

21 デザインの詳細 診断時の年齢 追跡期間は対照群 BPH 群 前立腺癌群でマッチングした 50 歳以上では BPH の罹患率が高すぎるため 対照群を見つけるのが難しかった 対照群は BPH 群に比べて 初回来院や診断時の年齢がだいぶ若い 局所浸潤性癌と転移性癌を分けて考えた PSA は指数関数的に増加するので 対数を考え 0 に近い値であることも考慮して をプロットした 21

22 Prostate Data 人口統計学的データなど ( 表 2.3) テキスト 12P より引用 22

23 国立がんセンターの Web ページ ( より タンデム R 法の グレーゾーン の ~ の値は 23

24 Prostate Data 個人ごとのデータの推移 ( 図 2.3) ( 横軸は診断前の時間 : 右が過去 )( 再掲 ) これを区別するモデリングをしたい テキスト 13P より引用 24

25 10.1 INTRODUCTION 25

26 誤差とは? 誤差というからには 独立同分布が基本では? 少なくとも独立性は欲しい 誤差が相関する というのは不思議な表現 モデリングが不十分 と解釈するべきでは? 誤差が相関がある場合 その要因を取り出して 変量効果によるモデリング + ( 独立同分布の ) 誤差 となるまでモデリングを続けたい ( 理想 ) 26

27 探索 : モデルを段々複雑にしていく (1) 固定効果のみのモデル ( 固定効果以外は全て誤差 ) (2) 変量効果をいくつか入れたモデル 相関構造が複雑 別の要因の影響では? を分解 この相関構造は? 複雑なら 別の要因を考える 27

28 系列相関 (Serial Correlation) いくつかの変量効果を入れた後 個人のデータの経時推移による相関がに残っているか? 残っているなら その部分を 時間依存する要因 として取り出したい 逆に とりあえず取り出して見て 相関が無視できるか を考えてみては? この部分の相関が十分大きいかどうかをみる 誤差時間に依存する部分 ( 変量効果 ) 28

29 本章で扱うモデル ( 一般形 ) 誤差 : 独立同分布 独立 時点ごとの相関 ここが本題 誤差 系列相関 (10.1) 29

30 10.2 An Informal Checks for Serial Correlation 30

31 モチベーション 系列相関があるかどうかを検討したい 誤差が独立同分布かどうかが気になる ( 注意 ) 前節と少し記号使いが異なり ます 31

32 本章で何度か出てくる重要なこと ( 誤差と残差の関係 ) 真のモデル推定したモデル y 9 y 誤差 3 2 残差 x 誤差は未知 データは同じ x 残差は既知 32

33 本章で何度か出てくる重要なこと ( 誤差と残差の関係 ) 誤差 と 残差 は違います 誤差 : 真のモデルに入っているもの 未知 残差 : モデルを当てはめた後 データと推定値 ( や予測値 ) とのズレ 既知 残差は誤差の予測値 と考えることができます 誤差の性質を検討して 誤差は未知だから 代わりに残差でその性質を満たすかどうか検討する 大体性質同じでしょ? という論法が 本章でよく用いられます 33

34 系列相関 が必要かどうかの検討 誤差 (3.11) 時点ごとの相関 ここが邪魔 ここがなければ 時点ごとで独立等分散ならと判断できそう はいらない 34

35 を 穏便に ( 他の要素に影響のないように ) 消したい 固定効果を除いた部分 の各列と直交するベクトル空間へ射影 35

36 直交補空間 とおく の 1 次独立と仮定 における直交補空間 次元ベクトル空間 をとる はの正規直交基底 ( の構成方法 ) を含むの基底を1つとり Schmidt の直交化法で正規直交基底を作る はから具体的に構成できる 36

37 行列 の性質 とおくと 正規直交基底の性質より となる よって である 37

38 また 直交補空間の性質から となる 従って となる 38

39 の を に射影する作用素は だが 簡単のため 最初の を除いて を考える 39

40 となる これより となる 独立 等分散性からのズレ 40

41 よって が独立同分布に従うかどうかを調べれば が 0 かどうか ( 系列相関があるかどうか ) 分かる しかし 誤差 予測 から 残差 は未知 そこで 最小 2 乗推定量は平均構造が正しく特定されていれば不偏推定量 残差は誤差の予測値 を の代わりにする が独立同分布か 調べる 41

42 系列相関が必要かどうかの 手順のまとめ 1. からを求める 2. 外れの少なそうな平均モデルを仮定して を求める 3. が独立同分布か調べる 独立同分布っぽい : 系列相関なさそう独立同分布っぽくない : 系列相関ありそう 42

43 p136 の 1 行目で テキストの修正 1 ( も同じ ) 残差 とし 同ページ 10.2 の 8 行目で として この分散を を考える ではなく としているが これは間違い 残差の部分を誤差にする ( をにする ) と正しい 43

44 10.3 Flexible Models for Serial Correlation 44

45 モチベーション 系列相関がある場合 パラメトリックなモデリングを行い LMM の枠組みで推定を行いたい 45

46 系列相関をモデル化したい たとえば 時間が離れるほど ( 通常は ) 時点間の相関は減少する 時点 1 時点 2 時点 3 時点 4 時点 1 時点 2 時点 3 時点 4 相関の減少の仕方 をモデル化したい いくつか関数を作ってみる 多項式よりは柔軟な関数にしたい 46

47 仮定と関数の導入 何か仮定をおかないと漠然としすぎるので 以下の仮定をおく 仮定 2 時点間の相関の強さは どれだけ離れているか のみに依存する 時点と時点の間の相関が 関数を用いて と表せるとする 自身には依存しない 47

48 たとえば のとき の関数形を決めたい 48

49 Royston and Altman(1994) の関数 (fractional polynomial) に対して とおく ただし 0 も負の数も無理数も OK 多項式より柔軟 (10.2) とする 49

50 Lesaffre, Asefan, and Verbeke (1999) の関数 Royston and Altmanの関数は 以下のように変更する 初回観察日 のとき困る そして を (10.3) とモデル化する ( なお ) 50

51 P 行目 テキストの修正 2 は の誤植 51

52 SAS の現状 (ver 6.2) など 系列相関に式 (10.3) を用いて ML REML 推定を行うことは SAS ver 6.2 では無理 他の数値計算ソフトなら可能らしい ( 確認してません ) 最尤法の数値計算は収束しないことがある SAS で扱えるモデルに対しては PARM オプションで初期値を変更してみる 著者らの経験では 一般にが大きくなると収束しない場合が増える を小さく固定して を変化させる のがよいのでは? 52

53 具体例 :Prostate Data これらは 3 章と同じ 固定効果 ( 平均構造 ) 切片 時点 時点の 2 乗 年齢 変量効果 切片 時点 時点の 2 乗 53

54 系列相関の構造 に対して やを変更した色々なモデルを当てはめる モデルによっては 計算が収束しないこともある ( 複雑すぎる ) 54

55 ( いくつか試した中で ) 尤度関数が最大なモデル 表 10.1 (10.4) 含まれる 指数関数 Gaussian テキスト 139P より引用 大 系列相関の推定は正しくできていなさそう 55

56 重要なポイント 尤度関数が最大のモデル = 最もよいモデル と条件反射するのは 思考停止 尤度関数の値はあくまで 目安 ここから先で モデル構築 の腕前が試される 変数の数が多い場合 fitting はよくて当たり前 もっと変数が少なくて fitting がほとんど変わらない かつ 系列相関の推定がもっと妥当 なモデルはないか? 他のモデルも検討してみる 指数関数 Gaussian 56

57 例 :Prostate Data( 図 10.1) モデルによる推定曲線の違い モデルによって結構違うように見える 実線 : 尤度最大のモデル点線 ( 短い方 ):Gaussian 点線 ( 長い方 ): 指数関数 3 つとも大体近いのでは? ( 指数関数の方が 少し尤度最大のモデルに近い ) テキスト 140P より引用 57

58 表 10.2 SAS の出力 : 系列相関は指数関数 少しは改善 SAS ではは出ない テキスト 141P より引用 58

59 推定量の相関 59

60 10.4 The Semi-Variogram 60

61 モチベーション 前節では 系列相関の関数に対してパラメトリックなモデリングを行った しかし パラメトリックモデルでは扱いきれない パラメトリックモデルの妥当性を保障したい 場合なども考慮して ノンパラメトリックな推定方法を考えたい 61

62 Semi-Variogram とは 10.3 節のアプローチ 線形混合モデルでパラメトリックに考えた Semi-Variogram はノンパラメトリックな方法 ( 論文 ) Diggle (1988) : random-intercept の場合のみ Verbeke, Lesaffre and Brand (1998) : 変量効果が intercept 以外にもある場合に拡張 62

63 まずは簡単なモデルで考える 変量効果を切片項にのみ含むモデルを考える 変量切片 とおくと (10.5) となる 63

64 例 ) 4 時点の場合 64

65 続き ( ただし 4 時点に限らない一般論 ) 分散は時点によらず 対角成分 時点との相関は 非対角成分 65

66 固定効果 とおく 推定 66

67 67

68 以上をまとめると 同一症例の時点間差 から 右辺のをとした を Semi-Variogram と呼ぶ なお を仮定 遠く離れた 2 時点間の相関はほぼ 0 68

69 図 10.2 指数関数と Gaussian の例 テキスト 143P より引用 個人ごとの変量効果 個人ごとの変量効果 誤差 誤差 経過時間 系列相関がなければ ここがつぶれる 69

70 Semi-Variogram の推定方法 に 前ページのように にモデルを当てはめない から 縦軸 : 横軸 : を全 に対してプロットして Smoothing する ノンパラメトリック 70

71 の推定方法 異なる個人間の場合 全て独立となるので より 71

72 これより は で足されたものの数 なお である 72

73 忘れた頃だと思うので 復習 変量切片 とおくと となる (10.5) 変量効果 誤差 系列相関 相対的に小さいと なくてもよいのでは? 73

74 縦方向のつぶれ具合で 系列相関の有無の Informal Check ができる と比べてが大きい 系列相関必要では? テキスト 143P より引用 74

75 図 10.3:Vorozole Study の例 テキスト 144P より引用 75

76 図 10.3 Vorozole Study 固定効果は以下のものを考えた 時間 時間 ベースライン 時間の 2 乗 時間の 2 乗 ベースライン 変量効果は切片のみ 系列相関は Gaussian + 独立同分布 ( 正規分布 ) に従う誤差 76

77 補足 より詳しい話は Diggle, Liang, and Zager (1994) Analysis of Longitudinal Data の Chapter 5 参照 77

78 変量効果を追加したモデル 前節では変量効果が切片項だけに含まれるモデルを考えた 本節では 他の変量効果も追加したモデルを考えたい 78

79 一般のモデルの場合 の部分を消したい 考えられる方法 1. を考える (10.2 節参照 ) 2. を考える 何か推定量 経験 Bayes 推定量 79

80 どっちがよいか? (2. の欠点 ) はの正規性に強く依存する (7.8.2 節参照 ) 分散共分散構造の特定に大きく依存する 本書では 1. のアプローチを支持 80

81 記号の整備 81

82 ここで 以下のようにおく 成分を とおく このとき 82

83 同一症例の時点間差の射影 83

84 (i) (ii) いま から より (i) = 84

85 (ii) の計算 対角成分 非対角成分 2 倍は が対称行列だから 85

86 まとめ { (i) + (ii) } より 既知定数 これの数はデータ依存 (10.6) 未知パラメータ 期待値が 未知パラメータの線形結合 ( データの数 ) > ( 未知パラメータ数 ) なら線形回帰できる 86

87 パラメータ推定 ( 理想 ) の数が少ない場合 ( 規定来院日からのズレがほとんどない場合 ) のうち 値の異なるものを小さい順に並べて とおく とし を応答変数 パラメータを とした線形回帰分析を行う 87

88 を 曲線として推定 するのではなくて 離散的な点を推定する (10.6) より メータにする方が便利 を 1 パラ 88

89 パラメータ推定 ( 現実 : 方法 ) 実際は のバリエーションが多すぎる 全体からいくつか取り出して 小さい順に並べ とする 離散的な点しか推定されない の最大値 残りのに対応する部分は 線形補間する 要は 折れ線で回帰分析を行っている 89

90 パラメータ推定 ( 現実 : 論点 1) をどうやって決めるか? 増やすと補間の精度が上がるが パラメータが増えて パラメータ推定の精度は下がる を決めた後 どの点をとするか ( 著者らの推薦する方法 ) を全て計算し その % 点をとる 90

91 パラメータ推定 ( 現実 : 論点 2) Verbeke (1995) によると シミュレーションの結果 データごとに相関が高すぎて 多重共線性が起きる Ridge 回帰で回避できる の選び方は一意的ではない 系列相関がない場合 式 (10.6) はしない ( 対偶 ) に依存 を変更したとき 推定された様子が変わるようなら 系列相関がある 91

92 パラメータ推定 ( 現実 : 論点 3) Verbeke et al(1995) によると 経験的にの変更にあまり影響を受けないの推定方法がある ただし 系列相関の存在の有無 タイプについての結論は上の方法と変わらないことが多い 92

93 例 :Prostate Data ( 図 10.4) ほぼ一致 実線 :Semi-Variogram( ノンパラ ) 点線 :Gaussian を用いたパラメトリックな推定 Gaussian で十分よさそう テキスト 148P より引用 ノンパラを使ってパラメトリックな解析の妥当性を評価 93

94 10.5 Some Remarks 94

95 Prostate Data の解析結果のまとめ パラメトリックな検討結果 指数関数がよさそう ノンパラメトリックな検討結果 Gaussian がよさそう 変量効果がある場合 系列相関のモデルを 1 つに特定するのは難しい (True Model の特定はやっぱり難しい ) 著者らの経験によると 重要なことは 系列相関をモデルに入れること ( 最重要 ) 正しく系列相関をモデル化すること ( あくまでその次 ) とにかく 系列相関の有無 が最重点の検討課題 これが正しくできれば そこそこ Useful Model では? 95

96 まとめと参考文献 96

97 本日のまとめ 1. 系列相関が必要かどうかの informal check の方法をみた 変量効果のデザイン行列の列ベクトルから生成されるベクトル空間への射影 ( みたいなもの ) を考えた 2. 系列相関のパラメトリックなモデル化を考えた 指数関数 Gaussian 多項式より柔軟な Lesaffre, Asefan, and Verbeke の関数を用いたモデル化 3. 系列相関のノンパラメトリックな推定方法を考えた パラメトリックモデルよりも柔軟なモデリング パラメトリックモデルが正しいかどうかの検討に使える 97

98 テキストの参考文献 (1) Altham,P.M.E. (1984) Improving the precision of estimation by fitting a model. Journal of the Royal Statistical Society, Series B, 46, Diggle,P.J. (1988) An approach to the analysis of repeated measures. Biometrics, 44, Diggle,P.J., Liang,K.-Y., and Zeger,S.L. (1994) Analysis of Longitudinal Data. Cxford Science Publications. Oxford: Clarendon Press. Lesaffre,E., Asefa,M., and Verbeke,G. (1999) Assessing the goodness-of-fit of the Laird and Ware model: an example: the Jimma Infant Survival Differential Longitudinal Study. Statistics in Medicine, 18, Morrell, C.H., Pearson,J.D., and Brant,L.J. (1997) Linear transformations of linear mixed-effects models. The American Statistician, 51,

99 テキストの参考文献 (2) Pearson,J.D., Kaminski,P., Metter, E.J., Fozard, J.L., Brant,L.J., Morrell,C.H., and Carter,H.B. (1999) Modeling longitudinal rates of change in prostate specific antigen during aging. Proceedings of the Social Statistics Section of the American Statistical Association, Washington, DC, pp Pearson,J.D., Morrell, C.H., Landis,P.K., Carter,H.B., and Brant,L.J. (1994) Mixed-effects regression models for studying the natural history of prostate disease. Statistics in Medicine, 13, Royston,P. and Altman,D.G. (1994) Regression using fractional polynomials of continuous covariates: parsimonious parametric modelling. Applied Statistics, 43, Self,S.G. and Liang,K.Y.(1987) Asymptotic properties of maximum likelihood estimatiors and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82,

100 テキストの参考文献 (3) Stram,D.O and Lee,J.W. (1994) Variance components testing in the longitudinal mixed effects model. Biometrics, 50, Stram,D.O and Lee,J.W. (1995) Correction to: Variance components testing in the longitudinal mixed effects model. Biometrics, 51, Verbeke,G., Lesaffre,E, and Brant,L.J. (1998) The detection of residual serial correlation in linear mixed models. Statistics in Medicine, 17,

101 追加の参考文献 Box,G.E.P (1976). Science and Statistics. Journal of American Statistical Association, 71, Box,G.E.P (1979). Robustness in the Strategy of Scientific Model Building. Robustness in Statistics:Proceedings of a Workshop(1979) edited by R.L.Launer and G.N.Wilkinson Box,G.E.P. and Draper,N.R.(1987). Empirical Model-Building and Response Surfaces. Wiley. 土居正明, 横道洋司, 青山淑子, 五百路徹也, 中村竜児, 吉田和生, 白岩健, 松下勲, 西山毅, 井上永介, 上原秀昭, 山口亨, 酒井美良訳 (2011). 線形モデルとその拡張 - 一般化線形モデル 混合効果モデル 経時データのためのモデル-, 株式会社シーエーシー. (McCulloch,C.E., Searle,S.R, and Neuhaus, J.M. (2008) Generalized, Linear, and Mixed Models 2nd edition. Wiley.) 松山裕, 山口拓洋編訳 (2001). 医学統計のための線形混合モデル-SAS によるアプローチ, サイエンティスト社. (Verbeke,G. and Molenbergh,G.M.ed. (1997). Mixed models in Practice A SAS- Oriented Approach-. Springer) 三中信宏 (2006) 系統樹思考の世界, 講談社