morの解析ブログ

解析疫学、リスクにまつわるメモや計算

「推定」のまわりをさぐる.教科書では「解析はMHにより行う、因子が多ければ重回帰を用いる」という風で詳しい例は少ない.独自(のつもり)な思いつきで具体に試行.
 数理を用いるべきアセスメントにも切り込む.

余興 説明変数と目的変数を入れ替える

■ 結果yが、原因と密接な関連を持つxにより説明されるかを調べるのは、”ふつう”のモデリングである.だが xが、むしろyの原因としても不自然でない場合があり、因子と結果をモデル内で入れ替えた際に何が起こるか、興味が湧いた.
 推定係数とORについて復習しながら手元のデータでこれを調べた.
 結果、yとx:原因/結果を入れ替えた場合、入れ替わって因子としたyの推定値とほぼ一致した.これは生起因子のORの性質から自明と考えられる.xを結果としたモデルで推定される他の因子の係数は、xとの暴露程度を示すから、生起因子が他の因子におよぼすみかけの効果を調整する数値として扱えそうだ.
■ まず、層化によるORは入れ替えの影響でどうなるか.
 2×2表では、
      x暴露あり なし
y発生あり   a    b
なし       c      d     OR1 = ad/bc


 なのだが、真の因子では、


      y暴露あり なし
x発生あり   a    c
なし         b      d     OR2 = ad/cb =OR1


 のように、単なる転置となり、どの因子を交換してもORは変化しない.これに呼応してモデルからの推定値は交換によって変化しないことがわかる.


・副因子x2がモデルに含まれる場合を考える.x2は、結果に対するある程度のリスク値を示し、”交絡因子”の候補ともなる.
 yとxの入れ替えによってx2の曝露如何による、x「発生」を表す2×2表へと替わる.結果、当然ながら新たなORは大きく変化し、xとx2との曝露パターンを表すことになる.


モデリングによる推定係数を調べてみる.
■ 実例

 公開された実例データで、因子交換してみる.
・観光事例
元モデル 
glm(formula = s ~ wat + tya + tuk + tam + pot, family = binomial(link = logit), data =kanko)
Coefficients:                 ーーーー  s:発症
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -2.60178    0.71871  -3.620 0.000295 ***
wat         -0.58586    0.48841  -1.200 0.230322       水 
tya         -0.03148    0.32043  -0.098 0.921731       お茶
tuk          0.63870    0.35570   1.796 0.072560 .     佃煮
tam          2.36972    0.63520   3.731 0.000191 ***   卵焼き;生起因子
pot         -0.03996    0.43921  -0.091 0.927506       フライドポテト
---   Null deviance: 296.77  on 215  degrees of freedom
Residual deviance: 259.81  on 210  degrees of freedom
  (50 observations deleted due to missingness)
AIC: 271.81


交換モデル
> kank_r <- glm(tam~ s+wat+tya+tuk+tam+pot,family=binomial(link=logit,data=kanko)
Coefficients:          ーーーーこちらの y; tam 卵焼き
            Estimate Std. Error z value Pr(>|z|)   
(Intercept) -0.20871    0.57938  -0.360 0.718677   
s            2.38721    0.63644   3.751 0.000176 *** s;卵焼きを結果とした、”発症”の係数
wat          0.05779    0.66729   0.087 0.930991   
tya         -0.10029    0.47880  -0.209 0.834089   
tuk          1.70670    0.44607   3.826 0.000130 ***
pot          0.42120    0.47618   0.885 0.376406   


* 元モデルの卵焼きtam:2.37,交換モデルの”発症” s:2.39 はほぼ一致する.  
 ちなみに、元モデルデータから 
               卵焼き    佃煮
  cludeなOR        14.7     2.91
  各因子に対するMHOR  7.48-12.04  2.25-6.66 
   であった.
* このモデルで、水、お茶、佃煮、フライドポテト の係数は、卵焼きとの暴露重複状況を示し、生起因子の曝露効果(食い合わせ)による、見かけの発症リスク”交絡”を演出する強さとみなせる.これをふまえれば、佃煮の独自なリスクは小さいと予想できる.
   


・納豆オクラ事例
元モデル
 glm(formula = y ~ natto + mesi + pan + milk, family = binomial(link = "logit"), data = natto)
 Deviance Residuals:                y:発症   
    Min       1Q   Median       3Q      Max 
-2.5017  -0.4014   0.2991   0.2991   2.2621 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -2.8765     1.3156  -2.186   0.0288 * 
natto         3.9021     0.7173   5.440 5.34e-08 ***    納豆オクラ:生起因子
mesi          2.0588     1.2461   1.652   0.0985 .      
pan           0.3985     1.2045   0.331   0.7408   
milk         -0.8564     0.7393  -1.158   0.2467  
 ・・・・
 * ”交絡”とされることのある、飯の係数は、納豆オクラに次いで大きいままである.


交換モデル
 glm(formula = natto ~ y + mesi + pan + milk, family = binomial(link = "logit"), data = natto)
 Deviance Residuals:       ---納豆オクラ を”発症”と他の食材で説明させる
    Min       1Q   Median       3Q      Max 
-2.3323  -0.1968   0.3692   0.4458   2.8121 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -2.1932     1.4218  -1.543   0.1229   
y             3.8621     0.7056   5.473 4.41e-08 ***
mesi          0.5899     1.4102   0.418   0.6757    ご飯:納豆オクラとの暴露重複は大
pan          -2.1343     1.1614  -1.838   0.0661 . 
milk          0.3928     0.7802   0.503   0.6146   
・・・
* 納豆オクラ:3.90 , ”発症”:3.86 となっており、2者はほぼ一致する.
 なお、交絡因子とされる(ことがある)「ご飯」の係数の意味は異なる.納豆オクラとの暴露が重複することによる、みかけののリスク;”交絡”の強さを示すことになる.


■ 
・生起因子の推定値は、交換モデルで因子とした”結果”の推定値とほぼ一致する.
・交換されたモデルでは、他の因子は、生起因子との暴露重複程度を表す係数となり、”交絡”を演出する強さがわかる.

×

非ログインユーザーとして返信する