morの解析ブログ

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

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

道草 Rで 全因子26データからMHRDを作る記述

保存したデータ
 csv ■◇ ただしい症例定義・・ .csv ;   全メニュー
 cludeなRDデータ名:crd
------------ 
 data = kanzi   記述  
dd<- NULL  
dd<- kanzi
     #  dd= read.csv(file.choose())  # dd:csv 新データ
   # 空白削除
  d26 <-complete.cases(dd)
  d <- d[d26,]         # d :Rに保存された全因子データ


【cludeなRD】曝露数、部分リスクをも計算


   crd<- NULL
 t_ec<- NULL
 hidari<-NULL
     migi<-NULL
  for (x in 1:27) {          # x   全因子を調整対象とし、、、、       
  a1 <- sum( d[,1] *d[x] )    # 発生有
  b1 <- sum( d[,1] *(1-d[x]) )   # 
    c1 <- sum( (1-d[,1]) *d[x] )     # 発生なし
  d1 <- sum( (1-d[,1]) *(1-d[x]) )
                                e_count <- a1+c1   # 曝露数
  n1 <-a1+b1+c1+d1 
    k11<-a1+c1 
      k10<- b1+d1
     rd1 <- a1/k11-b1/k10
                            hidari<- c(hidari,a1/k11)   # リスク差 左
                              migi <- c(migi,b0/k10)
crd<-c(crd,rd1) 
                           t_ec<-c(t_ec,e_count)    # 曝露数ベクトル
}
    ----------
 書き出し 
 例 write.csv( データフレーム名, "ファイル名.csv")
 例 write.csv(data4, "data4.csv")
  write.csv(hidari , "hidari.csv")
  write.csv(migi , "migi.csv")
    
 【生起因子で26因子を調整してMHRD】
    # 26因子データから tで調整したMHOR を計算
   y <-13  # 固定する * yについて全因子を指定すると全因子×全因子のMHRD
mhrd <-NULL
tmhrd<-NULL
for (x in 1:27) {       # x   全因子を調整対象とし、、、、       
# for (y in 1:27) {      # y 13 ; t 因子固定* してt 有無についてtableを  
 a1 <- sum( d[,1] *d[x]*d[y] )   # y t あり
 b1 <- sum( d[,1] *(1-d[x]) *d[y])  # xの2×2を 固定yの tで層化
 c1 <- sum( (1-d[,1]) *d[x] *d[y])  
 d1 <- sum( (1-d[,1]) *(1-d[x])*d[y] )
  a0<- sum( d[,1] *d[x]*(1-d[y]) )  # y t なし 
  b0 <- sum( d[,1] *(1-d[x]) *(1-d[y] ) )
  c0 <- sum( (1-d[,1]) *d[x] *(1-d[y]) )  
  d0 <- sum( (1-d[,1]) *(1-d[x])*(1-d[y]) )
# 式 ;       ありなしによる、 の mhrd;;
  n1 <-a1+b1+c1+d1
   n0 <-a0+b0+c0+d0 
   k11<-a1+c1
k10<- b1+d1
k01<-a0+c0
k00<- b0+d0
rd1 <- a1/k11-b1/k10
rd0 <-a0/k01 -b0/k00
wei0 <- k01*k00/(k01+k00)
wei1 <- k11*k10/(k11+k10)
 mhrd<- (wei0*rd1+wei1*rd0)/(wei0+wei1)
mhrd
tmhrd<-c(tmhrd,mhrd)
 # print( tmhrd ) 
}
# * }

《前座》 交互作用モデルの方がよいのか 4,5因子モデル 全26データ による

■ 元データ
・4因子は客観的に選択した steなど以外を加えるモデルを作るため、26因子全部をデータとした
・症例定義2 ただし症例定義1は考慮せず発生としたもの
■ 4因子モデルと”+適当な因子”で5因子モデル さらに独立と交互作用
 切片とRsqで評価する
   
 4因子モデル
・ 因子   cRDの絶対値の大きな、または曝露の大なもの;4データ
・モデル  独立モデル、交互作用項*モデル
  *交互作用は、抑制因子のある3因子以上のモデルでは必ずあるはずだと想定


4データ独立m
 ko6 lm(formula = y ~ wat + mesi + tam + potesara, data = d)
  (Intercept) 0.10317 Multiple R-squared: 0.08622, Adjusted R-squared: 0.07039
 
4データ交互
 ko6k lm(formula = y ~ wat + mesi + tam + potesara + mesi:tam + mesi:tam:wat,
 data = d)
  (Intercept) -0.003488  Multiple R-squared: 0.09202, Adjusted R-squared: 0.06823
 →交互モデルでIC改善


■ 5因子モデル
・因子  4データに適当な因子を加えた ;+?? 5因子
・独立、交互


+steakで5因子独立
 ko6l<-lm(formula = y ~ wat+mesi+tam+potesara+steak , data = d)
  (Intercept) 0.09585  Multiple R-squared: 0.08684, Adjusted R-squared: 0.06698
+steakで5因子交互 
 ko6st lm(formula = y ~ wat + mesi + tam + potesara + steak + steak:mesi:tam +
 mesi:tam, data = d)
  (Intercept) 0.009791 Multiple R-squared: 0.09739, Adjusted R-squared: 0.06968 

+sake交互
 ko6ss lm(formula = y ~ wat + mesi + tam + potesara + sake + sake:mesi:tam +
 mesi:tam, data = d)
  (Intercept) 0.01796  Multiple R-squared: 0.09732, Adjusted R-squared: 0.0696
 → +steakまたは⁺sakeでIC,Rsq 改善  


■ 結果
・独立より交互作用モデルは有利
  ICは小さめで好ましく、Rsq は改善されている
  うち、steakとsakeでの違いは鮮明でない
・4因子より多い5因子モデルは有利
  思い付きで*steakなどを追加しただけで交互作用モデルは、さらに好ましくなる
  これはsteakに限らない
■ 課題
・症例定義による、ID削除
・lmでのモデルのよさをRsqで評価することの意味
 ”適当な因子選択”は、別記事でより検討