統計家に聞いても、正しい解釈はわからないことが多い

臨床統計家に「この結果はどういう意味なのか？こう解析した方が良いのか？」と聞いても、「(数学的に、)こういうものとしか言えない」と言われてもしっくりはこない事がある。それは解釈の問題だと思う

私見だが、

1. 数学的な理解(線形代数/確率統計の知識を基にモデルの仮定の理解やsimulation researchによるbias判断)

2. 数学的な結果をどう解釈するかの理解

3. 解釈できる範囲内でどちらの方が良いのかの理解

4. その解析をする方が良いのか、しない方が良いのかの理解

が必要だと思うが知識でなんとかなるのは第二段階までな気がする。

ただ、その理解する、というのもどこまで必要かも悩ましい

よくわからないまま統計モデルを使うのは良くない、という点では概ね皆意見が一致するが、

• 詳細な理解が必要？

→　しかし、現実ロジスティック回帰すら理解して使っている人は少ない

• 過去論文の理解？

→　過去の論文でしているからと言っても、全然状況が違うので仮定も違う事が多い

• 結果が変わりそうな事は全て理解？

→　オプションがあっても実際ほとんど変わらない事が多い一方、seedで大きく変わる解析もある。

もう少し勉強して結論を自分なりに出そうと思う

なぜロジスティック回帰分析では検証とか言われないのに機械学習ではとやかく言われるのか (バイアスバリアンストレードオフ

最近、世間的にも人工知能一般に関する認識が変わってきて、AIといった抽象的な言葉から機械学習という言葉に移り変わりつつあるようである。

tjo.hatenablog.com

確かに自分に当てはめて考えると、ある程度専門的な知識を有している人と話す時は「機械学習」と言っているし、そういった統計学的な事に日常触れていなさそうな人には「AI」と表現していることに気がついた。

統計学的な事に日常触れていなさそうな人に機械学習と言うと、なんとなく深層学習deep learningの事を連想されるからである。

じゃあ何を理解している人には機械学習と表現して議論するかというと、bias variance tradeoffとノーフリーランチ定理の二つではないかと考えた。

機械学習は統計学の延長であると解釈しているのだが、その際にバイアスバリアンストレードオフ(bias variance tradeoff)もしくは過学習をどのように避けるかは非常に重要な理解である¹。

scikit-learn.org https://www.jclinepi.com/article/S0895-4356(18)31081-3/abstract#

bias variance tradeoff

自分の理解を簡単に言うと、

• モデルのvarianceを重視すれば、複雑なモデルは作れるがデータの誤差も拾ってしまう。

• モデルのvarianceをある程度犠牲にすれば、モデルは単純になるが大雑把にはあっているモデルになる。

biasvariancetradeoff_1

上図のようなデータがある時に、どのような線を引くかという問題である。もともとのデータは真ん中のような二次関数である。

• 単純なモデル(左下)：例えばロジスティック回帰などは一番左のような直線を引く。これだと細かい事を考えれば間違いだらけだが、大まかな方向性は間違えないだろう。

• 複雑なモデル(右下)：機械学習のような複雑なモデルでは一番右のような曲線を引く事ができる。しかし少しデータを増やしてみると様々な直線の形が変わる事になる。

どちらが良いかはケースバイケースで目的にも依るが、現場感がないままに解析するとだいたい何か間違える、という根拠にもなる気がする。

今回は単変量なので散布図を書けばすぐにわかるが、これが多変量になると散布図をみてもなかなか気づきにくいと思う。

ちなみに筆者は36万以上のデータで説明変数50個くらいを使って検証すると、ロジスティック回帰がほとんどのモデルより(勾配ブースティングを除き)精度が高かったという学会発表をした事がある。

---

1. というか正直解析していると、おそらく未測定交絡因子の影響だろうが、チューニングしてもtestデータではロジスティック回帰の方が機械学習よりよっぽど良い予測をしてたりするのだが、理解してもらなかったのでこのような文章を書こうと思った。。↩

明視域とは、その計算の仕方

今回は全く関係がない、眼の話。自分で実験してみると結構計算通りになるし、老眼鏡が必要になりそうとかを予想してみる。正直、眼科のランドルト環(Cをみて右とか言うやつ)検査よりも自分でやってみる方が実際の見え方に繋がる気がする　(眼科の先生に確認とってないし、私自身専門ではない分野の個人的な理解なので間違ってたら教えて下さい)

あの人は目が良いから、とかいうけれど、①遠くが見えることが可能、という意味と②調整できる幅が長い、とどちらの意味で使用しているのかというお話。

コンタクトとかメガネのレンズの「-2.0D」ってどういう意味か

このDはディオプトリという読み方をする。1m(100cm)からそのディオプトリで割った距離が見える距離である。

例えば裸眼であなたが定規を持って、11cm-33cmくらいは見えたが、それ以上は文字が見えないとしよう。

それならあなたの目は100÷11=9, 100÷33=3 となり、あなたの目は3-9Dの幅が見えていることになる。その幅の事を明視域という

そんなあなたが-2.0Dのメガネをつけたとしよう。

そうすると、引き算してあなたの目は1-7Dが見えるという事になる。

逆に長さを計算すると、100÷7=14, 100÷1=100となり、あなたはメガネをつけると14cmから100cmまでの範囲は見える事になる。(この程度だとパソコン用メガネだが、-3.0だと近視遠視だけであればかなり遠くまで見えるという感じである)

これが明視域である。必要になることはあまりないが、10cm以下を見えるようにしていくにはかなりのDが必要となってくる。

日常生活では30cm-10mくらい見える必要があるので、その範囲になるように調整する。しかし歳をとって40歳くらいになってくると、この調整できる幅が小さくなってくる。

そうなると昔メガネが要らなかった人ほど、2-3Dとかの調整域になれば、近くをみるときに+1.0とかそのくらいの老眼鏡が必要になってくる。

白内障で単焦点レンズを入れても1Dの幅くらいは文字は読めるので、60歳くらいになればまあ手術してもあまり差はないという感じになってくる。(明視域でググれば、年齢によってどのくらいの範囲かはわかるかと。乱視が入っているとそのせいで見えにくいなどはあるが)

傾向スコアと多変量ロジスティック回帰の違い

傾向スコア解析というと、医学分野ではちょくちょく使われているが大学院に来てすぐだと何がどうなっているのかわからず、とりあえずパッケージ使ってやってみた、という感じになってしまうため、lectureを準備することになった。

本質的にはある暴露をみたい時にその暴露になりやすさのprobabilityを求める(統計的にやっている事は傾向スコアを作る所までは本質的には同じ)のだが、ロジスティック回帰を漠然としていると理解が難しいので、この授業をした時の復習から入る方針にした

yiliaojingji.hatenablog.com

• マッチングする際に論文に書く条件を3種類と、その使い分け
• マッチングした後に注意する事(overfittingの仮定と一般化妥当性の確認)
• マッチング以外に傾向スコアの使用方法3種類の方法

とした。1:2マッチングを人数の多さで選択したりなどファシリテートしがいのある間違いが多いので、組み立て方を考えてみようと思う。

測定誤差があると真の関係はわからない

医学系の研究で、検査が毎回変わりうる検査(信頼性が低い)をoutcomeにするとなかなか結果が安定しない事を説明するためにグラフを描いた

num <- 40
age <- sample(65:85, num, replace=TRUE)
XX1 <- data.frame(age)
XX1$sbp <- 130 + (XX1$age-72)^2*0.5
res1 <- lm(sbp ~ age, data = XX1)
XX1$sbp3 <- 130 + (XX1$age-72)^2*0.5 + rnorm(num,mean = 0,sd=5)
res3 <- lm(sbp3 ~ age, data = XX1)
XX1$sbp4 <- 130 + (XX1$age-72)^2*0.5 + rnorm(num,mean = 0,sd=10)
res4 <- lm(sbp4 ~ age, data = XX1)
XX1$sbp5 <- 130 + (XX1$age-72)^2*0.5 + rnorm(num,mean = 0,sd=15)
res5 <- lm(sbp5 ~ age, data = XX1)
par(mfrow=c(2,2)) 
plot(XX$age,XX1$sbp,xlab="",ylab="")
plot(XX$age,XX1$sbp3,xlab="",ylab="")
plot(XX$age,XX1$sbp4,xlab="",ylab="")
plot(XX$age,XX1$sbp5,xlab="",ylab="")

誤差の大きな二次関数

最後の二次関数などはもはや直線の関係と言ってもばれないような印象である。

誤差が大きい(もしくは未測定交絡因子の影響が大きい)と、調べたいものにそれを上回る効果が無ければあまり解析の意味がなくなると思うが、ビッグデータになって有意差があるだけではそのような間違いが起こりうると思う。(N数を大きくしても、systematic errorを拾うだけになる)

一般化線形モデルの仕上げ

後輩向けの最後のlectureをファシリテートした。 linear predictor, link function, probability density functionの話の後半二つの内容だった。

x <- c(1:100)
par(mfrow=c(1,3))
plot(x,type='l',xlab = "x",ylab = "sum of linear predictor")
plot(1/1/(1+exp(-(x-50)/5)),type='l',xlab = "x",ylab = "logit")
plot(x,rbinom(length(x),1,1:length(x)/length(x)),xlab = "x",ylab = "outcome")

logistic_3

乱数発生させて想定通りの結果となるかをこの後実験する

今何をしているのかが理解できるようになる事が重要だという結論に至った

なんでもかんでも交絡因子を調整してはいけない

なんでもかんでもいれてはいけない

医学分野の論文で、「これが代替の説明変数かな、とか、とりあえず多めに調整しておけば良いか、と思って説明変数を多数入れる」事がよく行われている。¹

collider bias （日本語では、合流点バイアス？）というが、あまり知られていない。

実際何も考えずに行うと、どうなるか実験する。²

DAG (Directed Acyclic Graph)

DAG_collider

想定しているモデルの因果ダイアグラムを描いておく

意図としては、交通事故の発生確率(連続変数)と年齢の関係を調べたいが、交絡因子として動脈硬化度(ASOや脳梗塞後などの上流因子として連続変数を想定)を考慮したい。

その際、腎機能に関してのデータがあったので、eGFRと動脈硬化度と年齢と腎機能低下に関するリスクは概ね解っていると仮定する。

腎機能と交通事故にはさすがに直接関係ないけれど、動脈硬化の度合いを他の説明変数で調整できないので、仕方なくeGFRを代替交絡因子として調整しようと考える。というようなイメージ

# サンプルサイズ
n.sample <- 5000
# それぞれの値を標準化したとして連続変数としてランダムに作成
age_v              <- rnorm(n.sample,mean=0,sd=1)
arteriosclerosis_v <- rnorm(n.sample,mean=0,sd=1)
totalrisk_v        <- rnorm(n.sample,mean=0,sd=1)
koutujiko_v <- 0.5*age_v + 1.5*arteriosclerosis_v + rnorm(n.sample,mean=0,sd=1)
egfr_v <- 1*age_v + 0.5*arteriosclerosis_v + 0.5*totalrisk_v + rnorm(n.sample,mean=0,sd=1)
data_v <- data.frame(AGE = age_v, ARTERIOSCLEROSIS = arteriosclerosis_v, TOTALRISK = totalrisk_v
                             ,KOTSUJIKO = koutujiko_v, eGFR = egfr_v)
head(data_v)

seed設定していないので各自確認。

res1 <- lm(KOTSUJIKO ~ AGE, data=data_v)
summary(res1)
## 
## Call:
## lm(formula = KOTSUJIKO ~ AGE, data = data_v)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.496 -1.203 -0.001  1.219  7.214 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.02460    0.02574  -0.955    0.339    
## AGE          0.51232    0.02560  20.010   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.82 on 4998 degrees of freedom
## Multiple R-squared:  0.07417,    Adjusted R-squared:  0.07399 
## F-statistic: 400.4 on 1 and 4998 DF,  p-value: < 2.2e-16

交通事故発生割合と年齢の間の関係は、初めに設定した真実(β =0.5)あたりの結果になる。

eGFRも関係あるかもしれへんやんと思ってしまうと、、

動脈硬化に関して調整しておいた方が良いかと思ったとする。まずは動脈硬化も入れて解析すると、

res2 <- lm(KOTSUJIKO ~ AGE + ARTERIOSCLEROSIS, data=data_v)
summary(res2)
## 
## Call:
## lm(formula = KOTSUJIKO ~ AGE + ARTERIOSCLEROSIS, data = data_v)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4910 -0.6877 -0.0116  0.6815  3.8000 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -0.007593   0.014210  -0.534    0.593    
## AGE               0.519588   0.014131  36.769   <2e-16 ***
## ARTERIOSCLEROSIS  1.513968   0.014173 106.819   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.004 on 4997 degrees of freedom
## Multiple R-squared:  0.718,  Adjusted R-squared:  0.7179 
## F-statistic:  6362 on 2 and 4997 DF,  p-value: < 2.2e-16

どちらも有意差はある。

もしeGFRを代理で用いると、

res3 <- lm(KOTSUJIKO ~ AGE + eGFR, data=data_v)
summary(res3)
## 
## Call:
## lm(formula = KOTSUJIKO ~ AGE + eGFR, data = data_v)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.5586 -1.1474 -0.0204  1.1407  6.1140 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.02056    0.02426  -0.847    0.397    
## AGE          0.01741    0.03117   0.558    0.577    
## eGFR         0.49286    0.01965  25.087   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.715 on 4997 degrees of freedom
## Multiple R-squared:  0.1777, Adjusted R-squared:  0.1774 
## F-statistic:   540 on 2 and 4997 DF,  p-value: < 2.2e-16

年齢の効果が完全に消えてしまう。

参考にeGFRのリスクで調整すると、

res4 <- lm(KOTSUJIKO ~ AGE + TOTALRISK, data=data_v)
summary(res4)
## 
## Call:
## lm(formula = KOTSUJIKO ~ AGE + TOTALRISK, data = data_v)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.5377 -1.2012  0.0051  1.2126  7.2430 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.02476    0.02574  -0.962    0.336    
## AGE          0.51238    0.02560  20.013   <2e-16 ***
## TOTALRISK    0.02504    0.02543   0.985    0.325    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.82 on 4997 degrees of freedom
## Multiple R-squared:  0.07435,    Adjusted R-squared:  0.07398 
## F-statistic: 200.7 on 2 and 4997 DF,  p-value: < 2.2e-16

こっちだとAGEに関しては正しい推定結果が出てくる。

まとめると

解析No	調整内容	交通事故に及ぼす年齢の影響
res1	単変量(年齢のみ)	有意差あり
res2	動脈硬化度で調整	有意差あり
res3	eGFRで調整	有意差なし
res4	eGFRのリスクで調整	有意差あり

個人的な意見としては、

DAG自体にも恣意性が入るので、

• DAGは自分の頭の中で書いて、

• この説明変数を除いた解析も入れておいたら文句を言われまい

という理論武装をしておく(modelを複数作る)のが良いと思う。

---

1. 時としてreviewerに求められたりする↩

2. http://takehiko-i-hayashi.hatenablog.com/entry/20130919/1379560929 を参考にさせて頂いた↩