ムラゴン ホーム
・linear model を2値データに適用したとき summayの Multiple R-squared は、実感とずれがある かなり小さな値がでてくる
・yが0,1のデータがそのまま相関x,yでのyとして計算され Multiple R-squared となっている
・推定結果、得た係数をIDごとの曝露で計算したpと、観察したyとの関係を考え、別な相関を計算してみる
■ lmでのMRsqは、公式通りの単純な相関係数である
plot(jitter(idp),jitter(d$y)) ：変数は後述
図 IDごとのpとy ；とあるモデルから
p：推定係数と曝露からなる期待値
例としたモデルの、推定pと発生yの相関は、
cor((idp),(d$y))
[1] 0.3120816
cor((idp),(d$y))^2
[1] 0.09739494 ；これは Multiple R-squared に一致
つまり、サマリーは２値のyそのまま、相関の計算をしている
また、その数値は、図からも明らかなようにごく弱い相関となってしまい、小さい
■ 別な相関 ・・ましな相関をどう計算するか
相関をみるべきxとy
まずxとなる、p これは推定の結果、coefficients×曝露有無 ；因子vectorの内積；IDごとに計算したvectorである pはもともと、yと、推定された係数によって決まる pの値は、あるID がもつ発生の可能性であって、曝露パターンの一致したサブｇ；ID n個 の pjは 対応するn個 のyのなかで y1/n という発生率に近づけたい性質のものであり、0,1のままでなく、pごとに集めたy1のy0+y1に占める割合として計算する
・記事；《前座》でまあまあよいモデルだった、 steak での計算例
ko6st <- lm(formula = y ~ wat + mesi + tam + potesara + steak + steak:mesi:tam +
mesi:tam, data = d)
summary(ko6st)
・・中略
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.009791 0.141247 0.069 0.94480
wat -0.110899 0.100638 -1.102 0.27164
mesi 0.166222 0.200210 0.830 0.40727
tam 0.620697 0.221751 2.799 0.00556 **
potesara 0.064811 0.095652 0.678 0.49873
steak -0.126427 0.166205 -0.761 0.44764
mesi:tam -0.540110 0.330976 -1.632 0.10409
mesi:tam:steak 0.238870 0.207519 1.151 0.25091
---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4628 on 228 degrees of freedom
Multiple R-squared: 0.09739, Adjusted R-squared: 0.06968
F-statistic: 3.515 on 7 and 228 DF, p-value: 0.00133
kc <- ko6st $coefficients # 推定結果 coef を kcとしておく
■ 手計算で
IDごとの、期待値 p ；
idp = kc[ 1] + kc[2 ]*vw+ kc[3 ]*vm+ kc[4 ] *vt+ kc[5 ]*vp+ kc[6 ]*vste+ kc[7 ]*vm*vt+ kc[8 ]*vm*vt*vste
idp ；期待p の ID分のvec
などとして手作業を織り込んで計算できそうだ
■ 曝露、係数ベクトルから IDごとの期待値pを計算する
pivdf<- data.frame(table( idp))
yp_<- cbind(vy,1-vy,idp ) 236行 3列
table( yp_) 21 1 （項目名 目盛
table(yp_)[1] -0.16・・ 1
-------------------
y1idp<-cbind(idp,vy) # IDごと pとy
y0idp<-cbind(idp,1-vy) # 〃 y0；1 これいらない
y1df<-data.frame(y1idp) # df化
y0df<-data.frame(y0idp) #
tab1 <- table(y1df) # 度数
rdf<- cbind(tab1,tab1[,2]/(tab1[,2]+tab1[,1])) # pごとのyの率
rowr<- as.numeric(rownames(rdf)) # 文字であるpを数値にする
・pの値ごとの、yの率の分布
この段階でpごとのy率について形式的に相関をみる：上の rdf；度数表から直にMRsqが再現される
cor(rowr,rdf[,3]) # pが作るサブgのID 数は相関係数に反映されない
[1] 0.8429513
> cor(rowr,rdf[,3])^2
[1] 0.7105669
図示 plot(rowr,rdf[,3],cex=(log(rdf[,2]+rdf[,3]+1)))
図 pごとのy率の位置 ；プロット大きさは目盛ごとの例数に応じる
* 2×2表でいえば
------------------
・本来的な相関
p.rate<-as.numeric(rownames(rdf)) # 目盛を数値にして取り出し
ratedf<- cbind(p.rate,rdf[,1],rdf[,2],rdf[,3]) # y01、yrateとつなげる
# 項目名をつける
これは、度数表
ID ごとのデータをつくる
n=pごとのy01数；ID 数
rep(ratedf[,1], times = ratedf[,2]+ratedf[,3] ) # p をn個つくる
repp<- rep(ratedf[,1], times = ratedf[,2]+ratedf[,3] ) # p をn繰り返したもの
repy<- rep(ratedf[,4], times = ratedf[,2]+ratedf[,3] ) # yrate を 〃
rept<-rbind(repp,repy) # 2つを合わせ・・
cor(rept[1,],rept[2,]) # 相関をみる これはいわば ｎのウエイトつき相関
[1] 0.8803379
> cor(rept[1,],rept[2,])^2
[1] 0.7749948
図 pによるサブｇごとのy発生率の相関 jitter入り
・この処理はあたかも重みをつけたような相関性を算出できる
■ 考えるに
２値解析はﾛｼﾞｽﾃｨｯｸ回帰が、連続量はlmが連想される
pが比較的大きく推定される一方、観測yは、1に一致せず、ましてかなりがy0になっていることは、交互作用項を含む複数のモデル試行で明らかだが、かつての記事にあるとおり、不思議なことでも不都合なことでもない
よって、その値 2×2表 でいえば c * に直接左右されない指標がほしいところ、発生有無でなくその率とで相関係数風な値を得て視覚化し、実感した pという位置において y1のものの割合が適合しているかの一致具合こそ実は重要なのだろう
ここで試行した”ましな”相関は、１つの思いつきであった
これによる値は、モデルサマリーと比べて、かなり良好な相関を示す
ちょっと残念なのはいちいちモデルごと計算しなければならず、使いづらい 当面、Rsqをみていくことになる
