No.21134　Cox比例ハザードモデルにおける内生性の問題と社会的な共通因子の採用　　【MWakaba】　2014/07/04(Fri) 09:35

統計手法にあまり詳しくなく，にわか勉強でCox比例ハザードモデルを利用している者です．非常に初歩的な質問かもしれませんが，ご容赦ください．
質問は以下の2点です．

1）内生性について
Cox 比例ハザードモデルは，医療分野でよく使われていますが，これを社会科学分野で応用したいと考えています．具体的には，会社や個人がある行動をとるまでの 時間をイベント（ハザード）の発生確率としてモデル化しています．分析結果を論文として投稿したところ，レビューアーから説明変数（共変数）と被説明変数 （イベント）の関係が原因と結果，すなわち因果関係にあるということが明確ではない，との指摘を受けました．社会学では「内生性」として知られる問題で， 被説明変数である企業あるいは個人の行動が，説明変数としている共変数に影響を与えている可能性があるため，得られたパラメータにバイアスが生じているの ではないか，というものです．

こうした問題は，もともとの応用分野である医療の分析では，考慮する必要がないのでしょうか？もし，Cox比例ハザードモデルの中で内生性を扱った研究があれば，参考にさせていただきたく，ご教示いただければ大変有り難いです．

2）社会的共通因子について
通 常，共通因子としては個々の変数，例えば患者の特性や薬の投与の有無などを使いますが，個々の患者に共通の変数を採用することはありますでしょうか？　私 が考えているモデルでは，社会的な変動要因として，全ての固体に共通だが時間とともに変化する変数を使いたいと思っています．　医療分野でも，例えば特定 地域の患者の疾患の分析で，その地域の大気汚染の状況など，個別の患者間では変化しない変数が考えられるかと思います．
今回，STATAというソ フトを使って，こうした変数を取り入れてCox比例ハザードモデルを当てはめたところ，パラメータは求められるのですが，標準偏差が計算できないという結 果が得られました（下記がその出力結果で，eclcvrgという変数が問題としている共変数です）

Cox regression -- Breslow method for ties

No. of subjects = 884
No. of failures = 597
Time at risk = 14989
Number of obs = 14989
Log likelihood = -3629.8262
LR chi2(28) = 332.93
Prob > chi2 = 0.0000

------------------------------------------------------------------------------
　　　　　　|　Haz. Ratio　　Std. Err.　　z　　P>|z|
-------------+----------------------------------------------------------------
1.grtarget |　1.653434　.2807315　2.96　0.003
eclcvrg |　1.490999　　.　　　.　　　.
peerprs |　1.254075　.4252425　0.67　0.504

これはどういう理由によるのか，そもそもCox比例ハザードモデルでこうした共通因子を用いることは不適切なのか，ご教示いただければ有り難いです．

以上，どうぞよろしくお願いいたします
No.21136　Re: Cox比例ハザードモデルにおける内生性の問題と社会的な共通因子の採用　　【青木繁伸】　2014/07/04(Fri) 11:49

eclcvrg は，データファイル中の全てのケースにおいて，同じ値を取るということですか？
そういう変数は，分析から省くべきでしょう。省かなくても，解析プログラムが除外してくれると思いますが。
以下のテストデータでは，y が全て同じ値を取るので，解析プログラムが文句を言っています。
> library(survival)
> n <- 10
> d <- data.frame(time=rexp(n), event=sample(0:1, n, prob=c(0.6, 0.4), replace=TRUE),
                  x=rnorm(n), y=rep(3.5, n))
> d
        time event           x   y
1  0.5976632     1 -1.82862663 3.5
2  0.5691505     1 -1.81353336 3.5
3  4.3659109     0  1.37261371 3.5
4  0.5581318     0 -0.56426954 3.5
5  0.4587813     1  0.97031123 3.5
6  0.2602637     1 -0.01863398 3.5
7  0.6988417     0  0.36237035 3.5
8  0.9631884     0  2.01130559 3.5
9  0.7314796     1 -1.17796750 3.5
10 1.1122796     0 -0.75517048 3.5
> coxph(Surv(time, event) ~ x + y, d)
Call:
coxph(formula = Surv(time, event) ~ x + y, data = d)

    coef exp(coef) se(coef)     z    p
x -0.471     0.624    0.406 -1.16 0.25
y     NA        NA    0.000    NA   NA

Likelihood ratio test=1.61  on 1 df, p=0.204  n= 10, number of events= 5 
 警告メッセージ: 
In coxph(Surv(time, event) ~ x + y, d) :
  X matrix deemed to be singular; variable 2
No.21137　Re: Cox比例ハザードモデルにおける内生性の問題と社会的な共通因子の採用　　【MWakaba】　2014/07/04(Fri) 13:06

青木先生

恐れ入ります．表現が適切ではなく，申し訳ありません
時間とともに変化する変数ですが，固体による差はありません

下記にデータサンプルをつけました　いずれのサンプルにおいても，時間ごとに同じ値が繰り返されています

ご指導，よろしくお願いいたします

. list samplenumber grtarget eclcvrg

+--------------------------------+
| samplenumber　grtarget　eclcvrg |
|--------------------------------|
1. | 1 0 1 |
2. | 1 0 1.014673 |
3. | 1 0 1.019139 |
4. | 1 0 1.013716 |
5. | 1 0 1.102392 |
6. | 1 0 1.119936 |
7. | 1 0 1.12185 |
8. | 1 0 1.130782 |
9. | 1 1 1.162998 |
10. | 1 1 3.130144 |
11. | 1 1 3.20925 |
12. | 1 1 3.255183 |
|--------------------------------|
21. | 2 0 1 |
22. | 2 0 1.014673 |
23. | 2 0 1.019139 |
|--------------------------------|
44. | 4 0 1 |
45. | 4 0 1.014673 |
46. | 4 0 1.019139 |
No.21140　Re: Cox比例ハザードモデルにおける内生性の問題と社会的な共通因子の採用　　【青木繁伸】　2014/07/04(Fri) 14:58

> いずれのサンプルにおいても，時間ごとに同じ値が繰り返されています

まだよくわからないんですが。
例に挙げられた 1〜12のデータは，同じサンプルのデータなんですか？
独立なデータでないにもかかわらず，独立なデータとして扱ったんですか？
そうでないなら，どうやったのですか？
No.21141　Re: Cox比例ハザードモデルにおける内生性の問題と社会的な共通因子の採用　　【MWakaba】　2014/07/04(Fri) 15:25

適切な説明ができず，申し訳ありません

sample　の数字が同じものは，同じサンプルのデータです．時間tによって値が変化していて，targetが0->1になるところでイベントが発生しています．
tは0から最大23まであり，targetが1になった時点で打ち切られます（sample=1のデータはt=20で打ち切りとなっています）

変数eclcvrgというは，時間によって逓増していますが，社会的な変数であり，それぞれのサンプルのtが同じであれば，同じ値が入っています．

下記にデータリストの一部をつけます
       +------------------------------------------------+
       | sample　t　target　grtarget　eclcvrg 
       |------------------------------------------------|
    1. |        1    1          0          0          1 |
    2. |        1    2          0          0   1.014673 |
    3. |        1    3          0          0   1.019139 |
    4. |        1    4          0          0   1.013716 |
    5. |        1    5          0          0   1.102392 |
    6. |        1    6          0          0   1.119936 |
    7. |        1    7          0          0    1.12185 |
    8. |        1    8          0          0   1.130782 |
    9. |        1    9          0          1   1.162998 |
   10. |        1   10          0          1   3.130144 |
   11. |        1   11          0          1    3.20925 |
   12. |        1   12          0          1   3.255183 |
   13. |        1   13          0          1   3.257416 |
   14. |        1   14          0          1   3.348644 |
   15. |        1   15          0          1   3.434769 |
   16. |        1   16          0          1   3.569059 |
   17. |        1   17          0          1   4.322488 |
   18. |        1   18          0          1   4.502711 |
   19. |        1   19          0          1   4.689952 |
   20. |        1   20          1          1   4.699203 |
       |------------------------------------------------|
   21. |        2    1          0          0          1 |
   22. |        2    2          0          0   1.014673 |
   23. |        2    3          0          0   1.019139 |
   24. |        2    4          0          0   1.013716 |
   25. |        2    5          0          0   1.102392 |
   26. |        2    6          0          0   1.119936 |
   27. |        2    7          0          0    1.12185 |
   28. |        2    8          0          0   1.130782 |
   29. |        2    9          0          1   1.162998 |
   30. |        2   10          0          1   3.130144 |
   31. |        2   11          0          1    3.20925 |
   32. |        2   12          0          1   3.255183 |
   33. |        2   13          0          1   3.257416 |
   34. |        2   14          0          1   3.348644 |
   35. |        2   15          0          1   3.434769 |
   36. |        2   16          0          1   3.569059 |
   37. |        2   17          0          1   4.322488 |
   38. |        2   18          0          1   4.502711 |
   39. |        2   19          0          1   4.689952 |
   40. |        2   20          0          1   4.699203 |
   41. |        2   21          0          1   4.671452 |
   42. |        2   22          0          1   4.739394 |
   43. |        2   23          0          1   4.745773 |
       |------------------------------------------------|
   44. |        4    1          0          0          1 |
   45. |        4    2          0          0   1.014673 |
   46. |        4    3          0          0   1.019139 |
   47. |        4    4          0          0   1.013716 |
   48. |        4    5          0          0   1.102392 |
   49. |        4    6          0          0   1.119936 |
   50. |        4    7          0          0    1.12185 |
   51. |        4    8          0          1   1.130782 |
   52. |        4    9          0          1   1.162998 |
   53. |        4   10          0          1   3.130144 |
   54. |        4   11          0          1    3.20925 |
   55. |        4   12          0          1   3.255183 |
   56. |        4   13          0          1   3.257416 |
   57. |        4   14          0          1   3.348644 |
   58. |        4   15          0          1   3.434769 |
   59. |        4   16          0          1   3.569059 |
No.21142　Re: Cox比例ハザードモデルにおける内生性の問題と社会的な共通因子の採用　　【青木繁伸】　2014/07/04(Fri) 15:42

STATA での分析のためにどのようにデータを準備する必要があるか分からないのですが，普通のソフトでは，

> sample　の数字が同じものは，同じサンプルのデータです．時間tによって値が変化していて，targetが0->1になるところでイベントが発生しています．
tは0から最大23まであり，targetが1になった時点で打ち切られます（sample=1のデータはt=20で打ち切りとなっています）

の ようなデータは，以下のようになりますよ（No. 21136 で例示したデータも同じ形式）。つまり，1サンプル1行。時間 time は，あなたの言うようにイベントの発生時間，または，sample 2，4 のような，イベントが起きないまま追跡が完了する「打ち切り」時間（あなたのいう『sample=1のデータはt=20で打ち切りとなっています』は用法 が間違えています）。event は，事象が起きれば1，打ち切りなら 0。
id time  event covariate
1    20      1   ??
2    23      0   ??
4    16      0   ??
STATA の分析プログラムが，あなたが挙げた1サンプル複数行のデータから，上のようなデータにまとめてくれるのならそれでよいのですが。
しかし，そうであっても，covariate をどうするかは解決されませんね。いつの時点のeclcvrg を使うのか，最初？最後？平均？それともそれらを含めていくつか。
No.21143　Re: Cox比例ハザードモデルにおける内生性の問題と社会的な共通因子の採用　　【MWakaba】　2014/07/04(Fri) 17:39

度々申し訳ありません
基本的なところで理解できていない可能性があるのですが，1サンプル1行ではなく，1サンプルで複数（観測時間数）のデータを想定していました．

異なる例になって恐縮ですが，簡単化のため，以下のようなデータで質問を続けさせてください．

id=1 t=0から5までのcovariateはx11，
t=5から8までのcovariateはx12であり，
t=8 でイベントが発生する

先生の示してくださったデータ形式で書くと，次のようなになると思います

id t1 t2 event covariate
1 0 5 0 x11
1 5 8 1 x12
2 0 23 1 x2
4 0 16 0 x4

この2つのレコード（x11とx12）は，いずれも推定に用いられると考えていたのですが，間違っていますでしょうか？
No.21144　Re: Cox比例ハザードモデルにおける内生性の問題と社会的な共通因子の採用　　【青木繁伸】　2014/07/04(Fri) 19:10

> 1サンプル1行ではなく，1サンプルで複数（観測時間数）のデータを想定していました．

あなたの使った STATA の Cox 比例ハザードモデルのプログラムが，どのようなデータを期待しているか，再確認してください。

Cox の比例ハザードモデルに限りませんが，普通の統計解析プログラムは，1サンプル1行のデータを仮定していると思います。repeated measure の場合には，1サンプル複数行で表現されることもありますが，それは，統計解析プログラムがそのようなデータを扱うように書かれているからです。

> この2つのレコード（x11とx12）は，いずれも推定に用いられると考えていたのですが，間違っていますでしょうか？

間違っていると思います。統計で使われるデータ行列がどういうものか，いくつか例を見れば，あなたにも分かってもらえると思います。

> レビューアーから説明変数（共変数）と被説明変数（イベント）の関係が原因と結果，すなわち因果関係にあるということが明確ではない，との指摘を受けました

このレビューアーは，あなたが分析に使ったデータ行列がおかしいということに全く気づかなかったのでしょうか。それとも，データ行列がどういう形になっているかの説明が十分じゃなかったのでしょうか。
No.21145　Re: Cox比例ハザードモデルにおける内生性の問題と社会的な共通因子の採用　　【MWakaba】　2014/07/04(Fri) 19:59

非常に基本的なところで大きな間違いを犯していたようです
ありがとうございます

Stataで用いられるデータ構造を再確認の上，改めて不明な点があれば質問させていただきます．本当に何度もありがとうございました

論文では，データ行列について細かに説明しておりませんので，後者（私の説明が不十分であったため）と思われます．．．
No.21146　Re: Cox比例ハザードモデルにおける内生性の問題と社会的な共通因子の採用　　【青木繁伸】　2014/07/04(Fri) 20:52

データ行列について説明しておきましょう。
ケース（観察単位）は，1行で表します。列は，変数（観察対象）を表します。
行も列も，同じものはありません。
   data-1 data-2 … data-p
A
B
C
:
N
もし，各ケースについて何回かずつ測定されたデータがあるならどうするか。
   data-1 x1 x2 … xm … data-p
A
B
C
:
N
の ようにします。x1, x2, …, xm はたぶん何らかの共通な規則の下に測定されていないとマズイでしょう。たとえば 3 ヶ月ごととか。厳密に3ヶ月でなくてもよいかも知れないが，大きく外れればまずい。場合によってはある時点での測定は無いかも知れない（欠損値）
このデータと同じだからと言って，
   data-1 x … data-p
A1        x1
A2        x2
:         :
Am        xm
B1        x1
:
C1        x2
と すると，前述のようなデータだと思っている統計解析プログラムにとってみれば，「いきなり m 倍のデータかよ!」てなことになる。気配りのある統計解析プログラムなら，「1サンプル当たり m 行だから気をつけてね」と言っておけばよいかも知れないが，data-1, data-2 …, data-p は各行に同じ値を複製しないといけない。だれもそんなことはしたくないから，各サンプルの2行目以降はx2, …, xm の1こずつ（そのほかのものもあって良いけど）記述なんて，おかしなことが。
更におかしなことは，1サンプルあたりの行数が同じならよいけど，場 合によっては行数がことなる。そんなことがあると益々おかしなことになる。そこまで，細やかに対応する統計解析プログラムがあるかどうか（だれも，そんな 変なデータを扱えるようなプログラムなんて，書く気にならない）。
こんなことがあると，次のような極端なデータの例も考えられる。
A さんの体重1
A さんの体重2
　:
A さんの体重na
B さんの体重1
B さんの体重2
　:
B さんの体重nb
  :
さて，このようなデータにおいて，体重の平均値を求めよという統計解析プログラムはどのような動作が望まれるか？
1 サンプルは1行であると思い込んでいる統計解析プログラムで得られる結果はまるでダメだろう。A さんの体重データが100回分あって，Bさんのデータは2回分しか無いなんてことを考えると明らか。親切に同じ人の体重を識別してくれる統計解析プログラ ムは，1番目に記録されている体重の平均値，2番目に記録されている体重の平均値，を計算してくれるかも知れない。しかし，m番目に記録されている体重が ばらばらの時期でのものならどうする？
さらに，今回のあなたのデータのような場合，時間2での生存データは，時間1での生存データの水増しにすぎない。各時点時点で変化しないデータ（例えば性別なんかその最適なもの）があると，変だということも明らかになる。
Aさん　1年後はまだ生存　性別は男　年齢は35歳
Aさん　2年後はまだ生存　性別は男　年齢は36歳
:
Aさん　n年後はまだ生存　性別は男　年齢は34+n歳
Aさん　n+1年後に死亡　性別は男　年齢は35+n歳
同じようなことがBさん，Cさんと繰り返される。おかしいですよね。例えば A さんのデータは10年分，Bさんのデータは2年分とすると，AさんのデータはBさんのデータの5倍の重みを持つことになってしまいますよ。
A さんのデータは，
A さん n+1 年後に死亡 性別は男　観察対象になったのは 35 歳
この1行だけで十分。Cox の比例ハザードモデルの解析プログラムが期待するのは，このデータ形式。もし，分析対象になってからの年ごとのある変数が必要なら（上の例では各時点での年齢も同じ性質を持っているのだが）
A さん n+1 年後に死亡 性別は男　観察対象になったのは 35 歳　観察対象1年目のある変数の値　観察対象2年目のある変数の値　… 観察対象k年目のある変数の値
という形式です。被検者によって何年まで観察されたかは違うので，ある被検者においては，観察対象 m 年目以降のある変数の値は，欠損値になっているでしょうけど。
あなたのデータにおいては更に制約があって，どの対象者においても，「観察対象k年目のある変数の値」は全て同じ値を取るということ。
そのような場合，観察対象1年目のある変数の値は全て同じになり，私が最初にあげたようなことになり，統計解析プログラムは文句を言うことになる。

以上のことは，あまりにも当たり前すぎて，分かってもらえるように書くことの方が難しかった。
ご理解戴けたかなあと，ちょと不安。

No. of subjects = 884
No. of failures = 597
Time at risk = 14989
Number of obs = 14989

というところも，ちょっと変だなとは思う。
それぞれの数値が，あなたのデータの何を表しているかがはっきりすれば，STATA の Cox比例ハザード分析プログラムが，どのようなデータだと思って分析したのかが分かる，かな？
No.21150　Re: Cox比例ハザードモデルにおける内生性の問題と社会的な共通因子の採用　　【MWakaba】　2014/07/07(Mon) 10:57

青木先生

大変丁寧な解説をありがとうございます！
ご指摘の通り，観測値の数による「水増し」を行っていました

筒井淳也他「Stataで計量経済学入門」ミネルヴァ書房，2011　の9章 サバイバル分析
にその方法が示されているのですが，
expand periodというコマンドでそれぞれのサンプルについて生存時間の長さだけ複製し，eclcvrgなどの時間依存関数についてはt（時間）ごとに値を設定していました

No. of subjects = 884
No. of failures = 597
Time at risk = 14989
Number of obs = 14989

と いう表現は，サンプルは884，うち597でイベントが発生（残りの287は打ち切り），1つのサンプルの時間tの観測値を1と数えた場合に，観測値は t=0からイベントが発生するまで（あるいは，打ち切りデータの場合はt=23まで）あることから，全部で14989のsample period dataがあることを示しています
（以下がプログラムで生存時間分析のデータ宣言をした結果です）

. stset t, failure(target==1) id(sample)

id: sample
failure event: target == 1
obs. time interval: (t[_n-1], t]
exit on or before: failure

------------------------------------------------------------------------------
16388 total obs.
0 exclusions
------------------------------------------------------------------------------
16388 obs. remaining, representing
956 subjects
633 failures in single failure-per-subject data
16388 total analysis time at risk, at risk from t = 0
earliest observed entry t = 0
last observed exit t = 23

上掲書のp.256にも同じような例があります

ご 指摘のように，通常の重回帰モデルのイメージだと，サンプルによって回帰の重み付けが異なる結果になるなぁ．．．とは感じていましたが，Cox回帰モデル ではそのようなデータを用いても結果にバイアスが生じないよう，うまく処理してくれているんだ，と勝手に理解していました．．．

● 「統計学関連なんでもあり」の過去ログ--- 046 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る