■　因子をできるだけ客観的に選んで、有利な交互作用をもつモデルで組み合わせを試す
■　因子、モデル　
因子
　　wat mesi tam potesara は　客観的選択
　　steakは　曝露がmesiに次いで大　　
　mesi steakには、阻止を想定した交互
　mesi,wat　は抑制を表す交互作用を期待してみる
　------------
交互作用項モデル
　　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　
の性能を持つ
■　試行
　抑制を正してみる
　　kow<-　lm(formula = y ~ wat + mesi + tam + potesara + steak + wat:mesi:tam +
　wat:tam, data = d) 
だと、
 (Intercept) 0.095232  Multiple R-squared: 0.08684, Adjusted R-squared: 0.06291
となり、すべてが悪化する
　つまり第一抑制だけにとらわれると違うぞと

　なのでmesi:tamを入れ戻す
　　kowm<-　lm(formula = y ~ wat + mesi + tam + potesara + steak + wat:mesi:tam +
　wat:tam+mesi:tam, data = d)
  (Intercept) -0.003739   Multiple R-squared: 0.09215, Adjusted R-squared: 0.06428
 ましになった

　ならば、steak：第二抑制を想定する
　　kowsm<-　lm(formula = y ~ wat + mesi + tam + potesara + steak + wat:mesi:tam +
　steak:tam+mesi:tam, data = d)
  としたら、
　(Intercept) -0.01100　Multiple R-squared: 0.09376, Adjusted R-squared: 0.06182
　もっとましになった
　むしろ、大暴露は抑制を考えとけということか

　で、第二曝露にも阻止を想定すると　
　　 kowsms<-　lm(formula = y ~ wat + mesi + tam + potesara + steak + wat:mesi:tam +
　steak:tam + mesi:tam + mesi:tam:steak, data = d)
   (Intercept) -0.01514   Multiple R-squared: 0.12, Adjusted R-squared: 0.08496
 oh　これはもっともよい
　　　　＊　交互作用項係数からsteak,mesiは、抑制性、かつそれらが互いに阻止性

　　kowst<-　lm(formula = y ~ wat + mesi + tam + potesara + steak　　　 +steak:tam+mesi:tam+mesi:tam:steak, data = d)
　(Intercept) -0.01752　Multiple R-squared: 0.1199, Adjusted R-squared: 0.08888
　ほぼ同程度

・曝露の大きな因子　2個　生起性大きな2つ　抑制因子１つの組み合わせとした
　各因子の組み合わせは、
　　　kowsms　　　　　　　　kowst
　　抑制・曝露大１・生起　　　　　　　　　　
　　曝露大2・生起　　　　　  　〃　　　　 　
　　曝露大1・生起　　　　　　  〃　　　  　　
　　曝露大1・生起・曝露大2　　〃　
　　　　良　　　　　　　　　　次

■　なにがよかったのか　
・交互作用として曝露大の因子を重視したモデルは　よい という傾向だった
・交互作用項4つで最良のモデルができた
　当初、抑制と阻止をあらかじめ想定し交互作用項を設定、試行したとはいえ行き当たりばったりだった　因子をどうカバーしてきたかをベン図；質的な模式でみると次のようであった
   

   　　図　kowsmsの交互作用がカバーする因子

　また、交互作用項で限定したカテゴリーのID数を調べると
　　wat:mesi:tam 　　20
　　steak:tam 　　　176
　　mesi:tam　　　  184
　　mesi:tam:steak　170
　と、比較的IDの大なものとなっていた

　これらのことから、「生起因子と重複する項を増やした」ことが奏功した、と考えた
　だとすると、推定pが0.6あたりに集中し、そこのyが１と０に分かれることを裏付けるものと思える
　しかも、さらにpotesara：tamを追加するとわずかに改善する；結果記載なし　それでいて、t:wを追加するとかえって悪化する；結果記載なし　
　注意は、度数0となるカテゴリーを含む項を入れた予測子だと推定が；それとは別な項が　NAになり異常となる
・交互作用をどのように組んでいたか
　生起因子とできるだけ重複する因子で交互作用を組み、試行するのがよいだろう
　　　（曝露グループを計算して試行するよりは手間が少ない

・曝露が大であれば、他の因子との重複は強く、場合によっては係数の代償性変化を引き起こしているおそれはある　

■　まとめ
　客観的に選んだ5因子で交互作用を試した
　モデルを改善する試行のすえに曝露大なものなど、生起因子に影響する因子組み合わせを重視すべきと気づいた　　

　”曝露の大きな因子と、生起抑制因子をさがし”
　”生起因子に絡んだ交互作用項を過不足なく試行錯誤せよ”　
　　
■　課題
　生起因子と重複する因子を調べるR記述

保存したデータ
　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")
　　　　
　【生起因子で２６因子を調整して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 )　
}
# * }