morの解析ブログ

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

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

Rで 第3の因子をさがす ・・tとs に影響するもの

・生起因子曝露群の中で、抑制因子1つとし、それ以外の因子を調べる.
 粗table、MH指標を一斉に眺める方法はRでも試せた.実発生数率は、効果の大きな因子を探すのに有効だった.
 目立たない因子は、BGとして無関連か、他の効果を持つかを調べる総当たり的方法を模索する.
 まずID数の大きな t曝露下でs有無に対して及ぼされる 任意な i因子の効果をみる.
 
■ 処理
    i はデータ内の全因子とし繰り返し処理でリスク値計算.
 計算不能を増やさないRD風なリスク値を使う.なぜならオッズだとc=0で無限大になり、検討の幅が狭まるから.
 算出するデータは、各値のまま;risk、差をとった値;rd、さらにその差;rddをみる.
 plotもいろいろみる.


■ 計算に用いる各データの考え方 
 まず、t 曝露によるrisk;発生率は、
         rt  y1[7] / dr[7]           #  sum略 
 sありなしで tの発生率の差;rdをとれば、sの効果が残る. 
 s ありに限定した t のriskは、dr[7] ; s 曝露に絞るから分子にその曝露y1[6]をかけ、分母にdr[7]をかけて、
   s1rt = y1[6]* y1[7] / dr[6]*dr[7]        #  sum略  95/158
 また、s なしtのriskは、 
   s0rt = (1-y1[6]) * y1[7] / (1-dr[6])*dr[7]       #  sum略    37/55
となる.
 調べる因子のありなしについて計算し、4通りの「各データ」;riskを観察する.
   それらの差;rd 2通り と最終的にそれらの差;rdd も.;記述下記.


■ 各データ risk
順にplotし、記す.
1)各データをみておく 
 図示              menu
                  1  2   3       4    5   6     7        8
                [1] "y" "wat" "tya" "mesi" "tori" "sake" "tam" "potesara"
                   〇    △    +        ×       ◇      ▽      〼    米

    

            i s1t1   i s0t1         ni s1t1          nis0t1  
 観察:
1   i s1t1  wある場合    sあり群でリスク値は小さめ    w;抑制
3 ni s1t1   p、mなしの場合  sあり群で 〃 小さめ      p、m;阻止
4 nis0t1   mなしの場合    t 群のリスク値大        m;抑制?
 
              ;〇ほか、ここで意味ない因子のマーカーは注目しない
   各データ :   

      

■ 最初からメインなデータまでの記述


    idef <-NULL
   ivec4<-NULL 
   for (i in 1:8) {
# i あり sあり
   i1s1rt <- sum( y1[i]* y1[6]* y1[7]) / sum( dr[i]* dr[6]*dr[7])      # mだと 93/154
# i あり sなし 
   i1s0rt <- sum( y1[i]*(1-y1[6]) * y1[7]) / sum(  dr[i]*(1-dr[6])*dr[7])       # m   31/ 48
#     i あり sの効果:
        iRD <-  i1s1Rt-i1s0Rt  
# i なし sあり
   nis1rt <- sum( (1-y1[i])* y1[6]* y1[7]) / sum((1-dr[i])* dr[6]*dr[7])          # m 2/4
# i なし sなし 
   nis0rt <- sum( (1-y1[i])*(1-y1[6]) * y1[7]) / sum((1-  dr[i])*(1-dr[6])*dr[7] ) # m 6/7
#   i なし sの効果:
       niRD <-  nis1rt-nis0rt  
  idef<- rbind(    idef , c( iRD, niRD ,iRD-niRD  ) )   # ** i についての s効果の差
          # ■ 1,2,3,4を 並べたvec
          ivec <- c( i1s1rt , i1s0rt , nis1rt, nis0rt)
           ivec4<- rbind(ivec4,ivec)  # 4個を 因子数だけ      
        } 
 rownames(idef)<-t(menu)
   colnames(idef) <-c(  "iRD"," niRD" ,"iRD-niRD")
 rownames(ivec4)<-name
   colnames(ivec4) <-c( "ist"," it" ,"i0st","i0t")


* 各データをみる記述   
  for (i in 1:4) { 
  for (j in 1:8) {
  plot( i,ivec4[j,i],xlim=c(0,4),ylim=c(0,1), pch=j )
 par(new = T)
  } } 
 ■ 差データ
 各データの差をとったものをみる.  
2) i ありなしの結果

          

  図示すると、 ;記述、末尾
   

       上段:i 因子ない場合のR、下段: ありのR


・ 因子の効果
 このplotでも個別のデータ観察と同様、特徴ある因子がわかる.
 もともとの、tに対するsの抑制は0.07程だった.
 mについて、これが曝露しない場合;上段、-0.3ほどだが、mあり;下段でsの効果は0近くまで上がっている.wありでは、逆に-0.3 ほど下となる.ということはmは、生起性か抑制性を阻止する可能性、wは抑制か、s抑制を促進する可能性がある、
 抑制因子の定義*、粗なtableでの観察を根拠に、最終的な数値からみれば、mは生起性が弱いから阻止性、wは抑制とみてもよいが、kは小さく、全体として目立たない因子ともいえる.
 potもまた、mと似た挙動で、同様な因子.残る因子は若干の抑制性をもつ、とみる.
                * 「生起因子のもとで効果を減ずる因子」


・結果の散布図
 

      

          ↑  i なし   → i あり    
  散布図記述
   if( dev.cur() > 1 ) dev.off() # 前の図を消す
   plot(idef[,1],idef[,2])
   text( x=0.03+idef[,1]*0.98,y=idef[,2],rownames(idef) )
     散布図は、iありなしでplotすると、軸との離れ具合が sへの類似度;相反度を表す.
   
3) 差の差 結果のみの図示;rdd
   plot(c(1,1,1,1,1,1,1,1),idef[,3])
   text( x=c(1,1,0.8,1,1,1,1,1)+0.1,y=idef[,3],rownames(idef) ,cex=0.75)


        

 注目する i は、生起性、抑制性がはっきりしない因子であるから、i 有無でその差;rddは i ありが + ならsの阻止、i なしが - なら抑制性;弱 などとみることができる. 


■ まとめると
・各データ、差データ、差の差データを試した.plotも試した.
・rddは、
 「i 曝露下でのs効果-i 非曝露下でのs効果」
 であるが、同時に「抑制因子に対する i の効果 - 抑制されない因子への i 効果」でもある.


・risk、rdは、同じ程度の情報を得られる


 大きな効果を示す因子とそれを揺さぶる因子までは、調べられた.Rでそれらを記述することもできた.それらに影響する因子を調べるのは局所にとらわれ難儀したが、1つの方法がみつかった.ただし、tに曝露したIDのみを対象としており、場面によっては異なった方向の検討が必要.
・大きな効果を持つ因子に対する他の因子を検討するという段階を踏む.一方、数個の因子を曝露gとして検討する.このこととどう比較できるか.
-------------------
プロット記述
name<-t(menu)
 plot(c(1,1,1,1,1,1,1,1),idef[1:8,1],xlim=c(0.5,2.5),ylim=c(-0.4,0 ) )
  text( x=1+0.25,y=idef[1:8,1],name,cex=0.75 )
   par(new = T)
 plot(c(2,2,2,2,2,2,2,2),idef[1:8,2],xlim=c(0.5,2.5),ylim=c(-0.4,0) )
    text(x=2+0.25,y=idef[1:8,2],name ,cex=0.75)

×

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