二項定理を使って「発生頻度(連続変数)を二値変数に変換する」
線形予測子の和は連続変数であり、それがなぜbinary(二値変数)にできるのかはなかなかイメージがつきにくい。しかし自分でやってみるとイメージが湧くようになる。
発生頻度を1%から100%まで動かして、その発生頻度に従って二項定理を使用して乱数を発生させる。
bino <- 100 plot(1:bino,rbinom(bino,1,1:bino/bino))
線形予測子を求める(-∞~∞)
-> probability(発生割合)の値を推測する(0-1)
-> 二項定理で{0,1}に変更する
事で自作データベースを作成できる。
二次関数に直線を引くとどうなる？
データベース(DB)作成
まずは年齢だけで説明できる血圧が入ったDB(random errorなし)を作成
num <- 40 age <- sample(65:85, num, replace=TRUE) XX <- data.frame(age) XX$sbp <- 130 + (XX$age-65)*1.2 XX$sbp2 <- 100 + (XX$age-72)^2*0.5 head(XX) ## age sbp sbp2 ## 1 80 148.0 132.0 ## 2 83 151.6 160.5 ## 3 80 148.0 132.0 ## 4 83 151.6 160.5 ## 5 74 140.8 102.0 ## 6 68 133.6 108.0
サンプル数は気軽に集めれる40程度にしておく
直線の関係に直線を書いた
res1 <- lm(sbp ~ age, data = XX) summary(res1) ## Call: ## lm(formula = sbp ~ age, data = XX) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.442e-14 -5.835e-15 -4.310e-16 5.494e-15 1.349e-14 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 5.200e+01 1.611e-14 3.229e+15 <2e-16 *** ## age 1.200e+00 2.135e-16 5.620e+15 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 7.985e-15 on 38 degrees of freedom ## Multiple R-squared: 1, Adjusted R-squared: 1 ## F-statistic: 3.159e+31 on 1 and 38 DF, p-value: < 2.2e-16
まあそりゃ有意になる。
しかもデータベース作成の際の設定通り、ちゃんと年齢が１上がる度に血圧が1.2 mmHg上昇するというデータを作成した時の通りの結果になっている。
散布図に直線を引くと、
plot(XX$age,XX$sbp) abline(res1,col="red")
次、
二次関数の関係に直線を書いた
二次関数の関係に線形回帰モデルをすると、二次関数の関係に関しても直線を引くというちょっと意味のわからない事になる。
それをあえてやってみる
res2 <- lm(sbp2 ~ age, data = XX) summary(res2) ## ## Call: ## lm(formula = sbp2 ~ age, data = XX) ## ## Residuals: ## Min 1Q Median 3Q Max ## -17.491 -15.402 -5.213 14.154 32.952 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) -106.3341 34.4494 -3.087 0.00377 ** ## age 3.0443 0.4567 6.666 6.99e-08 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 17.08 on 38 degrees of freedom ## Multiple R-squared: 0.539, Adjusted R-squared: 0.5269 ## F-statistic: 44.44 on 1 and 38 DF, p-value: 6.986e-08
この結果は、モデルの仮定などを全く気にしない場合、年齢が１上昇するたびに血圧が3 mmHg上昇する(かなり明らかな有意差を以て上昇する)と解釈されうる。
40人しかサンプル数が無いのにここまでの有意差が出るのは大したものである。
散布図に直線を引く。
plot(XX$age,XX$sbp2) abline(res2,col="red") legend("topleft",legend=paste("coef:",round(coef(res2)[2],3))) legend("top",legend=paste(paste("p value:",round(summary(res2)$coefficients[ 2,"Pr(>|t|)"],5),"(signif digits=5)")))
こういう事は論文で結構行われてしまっている(特に多変量解析にするとわかりにくくなるので、本当によくある。メタアナリシスでもある)。例えば、妊娠期間と出生児生命予後の関係などを見る時に、妊娠期間を連続変数で入れるとこういう解析を行っているという事。つまり妊娠期間が長ければ長いほど、40週を超えて41週、42週となればなるほど、できればもっと長い方が予後はよくなるという仮定をしているという事になる。
じゃあどうすれば良いかというと、妊娠期間のように、短すぎても長すぎても予後が悪くなるだろうな、と思われる時には例えば3群に分ける(早産・正期産・過期産など)方が解釈上も仮定上もうまくfitする印象がある。
他にこういった例がどのような分野にあるだろうか。
table.1がいかに重要か
※ この記事はマメに追記・訂正していつか完成する予定です。
カンファレンスなどで「この人賢いなあ」と思う事がよくあるが、一貫してそう言う人たちに共通する性質を見つけた。
それはResultでは「Figure.1とtable.1を半端なく見ている」という所である(医療英語論文の書き方だが)。sampling bias (selection biasでまとめる人もいるかもだが)を検討するだけではなく、その後の解析・解釈に与える影響もきちんと辿っている印象がある。
賢いと思う人の指摘に感化されながら記録したどうすれば自分もそうなれるかメモから引用すると、
- 代表性があるのか(今までの論文と大きな違いが無いか)
- 解析対象集団に偏りが無いか(臨床的にイメージされるような集団か)
- 使用されている統計モデルを使用するための仮定と矛盾はないか
- table.1以外にどんな素データの見せ方が必要か
というあたりを意識して聞いていると賢い指摘ができるのではないかと思った。最近気をつけている事を自戒の念を込めてつらつら述べていく。
代表性があるのか(今までの論文と大きな違いが無いか)
に関しては、みんなやっていることである。いわゆる「素データの確認」と聞いてまず思い浮かぶ確認事項かと思う。(過去の論文と見比べなければいけないので時間がかかるため、oralの発表でパッと言われても正直わからない。)
discussionの３段落目くらいがだいたいそれに相当する印象で、なぜそのような違いが出たのかを考察する必要がある。
普通やっている事なので省略する。
解析対象集団に偏りが無いか(Methodを読んでイメージされているような集団か)
一旦データをもらって素データと向き合うと年齢に関してだけでも「若い人がメインの病気のはずやけど今回は意外と高齢者ばっかりだな」とか「死亡割合と年齢とがめちゃめちゃ相関してる(ほとんどそれで説明できてしまっている)な」1とか「糖尿病で高齢の人にはこの薬全然使われてないな」とか、気づきが生まれる。そういった気づきを経て(思っていなかったデータと向き合う事で)、まずは全項目の違和感の原因を列挙する。
次に違和感を解きほぐしていく。すると、データの抽出の段階でスクリプトに注意点しなくてはならない点2が実はあったり、ある属性の施設のデータがごっそり抜けていたりとかデータにそもそも問題がある事が本当によくある。
解析をやってる側からすると日常茶飯事であるが、それをそもそも調べていない事も多いし、調べて気づいた所でどのように扱うか(bias覚悟で除外するか、真実のはずと強行突破するかとか)は特に定まっていない。実際問題として案外論文には書いていないし、reviewerもわからないから結局は悪い事してもバレないと思う。
経験豊富なデータサイエンティストはデータの背景を理解した上で解析前に「このデータおかしくないですか？」とか言ってきてくれて、矛盾点(>疑問点)を指摘してくれる。また、素データの表示方法を変える事で様々な知りたい事を教えてくれる。医療系の論文でそこまでやっている人はどのくらいいるのだろうか。
医療系の論文をアブストラクトだけ読むという事は？
最近はAIでこんな事がわかりましたという言葉が一人歩きして、「基になった学習データはどういうデータだったのか」などに対する理解が乏しい印象がある。更にそういったニュースに対する批判も経験論や感情論が先行して、データに基づいて批判していない(そもそも開示されていないからそんな事できないのかもしれないが)印象がある。
医学系の論文でアブストラクトだけ読む場合3は「AIでこんな事がわかりました」と同じようにならないよう注意せねばならない。ロジスティック回帰でオッズ比がこうなりました、は全く同じ枠組みになりうる4と思うので。。
使用されている統計モデルを使用するための仮定と矛盾はないか
これも臨床やってた頃には無かった視点だった。今でこそ素データをきっちり見て同時に意識するようになった。数理モデルの仮定・外れ値の扱いなどがこれに相当する。
数理モデルの仮定に関しては、有名な所では、線型性・正規性・独立性などが挙げられる。例えば下記の記事である。
他にも「外れ値があればノンパラ」とかよく言うが、外れ値があるのかどうかの判断が難しい。外れ値の症例がどのような属性でどのような人かを見て、現在自分が解析しようと思っている集団から流石に離れているだろう、と思えばそれは外れ値として除外することもある。少ない症例数であれば、外れ値一人のために有意差が出る事などもありえるので、散布図を確認して解析対象とするかどうかは案外大事である。
外れ値どうするか問題でよく困るのは、BMIが60を超えていたらどのように扱うか考える事などが挙げられるが、これに相当する5。
table.1以外にどんな素データの見せ方が必要か
これが意外と大事だと思う。そもそもロジスティック回帰さえすべきではないという人もいるし、確かにそう思う事もある。しかし、解析していないと論文にならないときには、ロジスティック回帰の結果を書いて、table.3辺りに必要な素データの見せ方だと思うテーブルを書いたりする。
データに二峰性があるときや、年齢との交互作用がある時の追加解析などもこれに相当するのかと思う。また別の機会に記事にしたいと思う。
結論
データサイエンティストとしての実力はtable.1の見方を聞けばだいたいわかるのではないかと思うくらい重要である。
-
こういった相関が強い時に、年齢を連続変数から適当にカテゴリ分けすると、他のオッズ比などで、それっぽい有意差が出てしまったりする。逆に、年齢と相関が強い因子を入れておけば、ビックデータを使用するとどちらも有意差が出て、ほとんど年齢の影響をみているだけになる可能性がある。↩
-
ミスという単語を連想すると語弊がある。データの定義書を読み込めていなかったため条件付きのデータである事をわかっておらず、欠測に大きな偏りが出現するなどは本当によくある話である。他にも解析ソフトによって欠測値の扱われ方が違うために意図しない解析になっているなどもよくある。↩
-
自分の専門外だとそこまで深く読み込めない。逆に「その論文の問題点としてこういった問題が考えられる」と指摘しているがdiscussion読んでない事をさらけ出しているだけという事もよくある話である。↩
-
AIでしばしば使われる機械学習は魔法ではなく、結局は複雑な統計である。(だから概念的な話でなく具体的な話でAIという言葉を使うと、AIってそもそもどれを言ってるの？となる)ビッグデータの時代になってp値が出やすくなっているので、未測定交絡因子の存在など含めたロジスティック回帰の仮定や基データも意識せずに導き出された結論を使うのは「AIがこう言っている」という事を現実世界に当てはめるのと同じなので、気をつけた方が良いと思う。どうでも良いがスパースは逆方向のベクトルなので個人的にはどんどん使っても問題ないと思う。↩
-
身長160で体重160とか書いてあると、BMI62.5だが、データを記入する人が体重に身長を書いてしまっただけなのか本当にそうなのかなど迷う。答えはない、という結論です。↩
医療統計と疫学の雰囲気の違い
今回はいつも以上に独りよがりな意見ですが。。
公衆衛生大学院に入ると、疫学の教室と医療統計の教室で全く考え方が違う印象を受けた。
| 違い | |
|---|---|
| 医療統計 | 「統計学的な正しさ」にかなり価値がある(主眼を置く)イメージ。「誰にも文句を言われないデザイン」(RCTなど)と数理モデルを目指すイメージ |
| 疫学 | まずそのデータを集めた事に価値があるが、数理モデルはロジスティック回帰や傾向スコアなどで十分「こんなデータ使ってみました、もしかしたらこうかもしれません」という方法を多数扱う。 |
一般に、医療統計の教室は製薬会社の治験などにも関わっている事が多い印象。新薬の開発などはベースとなる統計学的な正しさにスキがあってはいけない、という気迫が伝わってくる。
一方で、疫学はエビデンスレベルが高いような研究がそもそも難しいというのもあるが、モデルの仮定にも結構ルーズで、「こうやったらこうなったよー、前からこう言われているから1まあいいんじゃない？」という人が多い気がする。それよりもそのデータを扱える事や料理の仕方の方に価値が置かれている印象がある。
-
医学系の論文って、introductionやdiscussionなどでは、一つでも論文を引用していればOKみたいな所がある気がします。↩
問診をする時何か忘れる、というのは普通の事らしい。
患者に問診に行く時、「あー、喫煙聞くの忘れてたー」とか、「そういや子供いるんやっけ」などとなる事が多い。
専門領域に関しては特に「そういえばアミロイドーシスの可能性ってないの？」とか指摘されて何も答えられないとなると自分の信頼に関わる問題になるので、最終的にはスマートな医師ほど「毎回同じように聞く自分問診リスト」のようなものができてくる。(専門修練医の指導をする際には必ず最初に作らせる)
しかし研修医の段階では、どのような事を聞かねばならないのかさえわからないので、様々な先人の知恵がある。
主訴に関しては「OPQRST」
- Onset：「いつから？」
- Pallative, Provoke：「どんな時に改善/悪化する？」
- Quality, Quantity：「どんな表現ができる？」
- Region：「どこに症状が？」
- aSsociated symptom：「他にどんな症状を伴う？」
- Time course：「発症してから今までの経時変化は？」
とそれ以外に忘れそうな事を連ねた「ABCDEFG」
- Anamne＆Allergy：既往歴・アレルギー
- Back Ground：職業・住居・家族構成
- Cook：食事・飲酒・喫煙歴
- Drug：服薬歴
- Exposure：曝露歴（ペット・旅行・クーラーとか）
- Family：家族歴(主に父母兄弟)
- Gynecology：婦人科的事項（妊娠、月経など）
などがある。
ちなみに漢方でも張景岳の問診十問歌という有名な問診リストがある。
張景岳の問診十問歌
こういうのを各学会で雛形作れば1、それを電子カルテに乗せてクローズドクエスチョンとして(テキストマイニングなんてしなくても)きれいにデータがとれるのにといつも思う。
マルチレベル分析がなぜ必要かの説明をRで簡単に実装して行う
臨床をやっている時には意識していなかったが、施設差を解析に入れないと解析結果が明らかにmisleadingとなる事があるので、Rで実装しておく
個人的には、「その治療法は施設によってやり方が全然違うんちゃうん？」もしくは「その病気の予後は地域によって全然違うんちゃうん？」と思ったら、マルチレベルはやるべきだと思う。
マルチレベル分析とは？
マルチレベル分析とは施設差を考慮する際によく使われる解析方法である。過去にある人の研究発表を聞いた後に「施設差を考慮しないといけないのでは？」と質問した際に「いや、これだけの症例数あるから大丈夫でしょ」みたいな事を言われた事があるのだが、この記事を読んで(Rで実験して)症例数を増やすだけでは全く大丈夫ではない事が伝わればと思う。
まずはデータベース作成のため、以下をコピペする。
library(mvtnorm) N <- 1000; se1 <- 150; se2 <- -90 sigma = matrix(c(se1, se2, se2, se1), ncol = 2) rand = rmvnorm(n = N, mean = c(130, 360), sigma) x0 <- rep("A",N) x1 = rand[, 1] x2 = rand[, 2] X1 <- cbind(data.frame(x0),x1,x2) rand_2 = rmvnorm(n = N, mean = c(140, 380), sigma) x0 <- rep("B",N) x1 = rand_2[, 1] x2 = rand_2[, 2] X2 <- cbind(data.frame(x0),x1,x2) rand_3 = rmvnorm(n = N, mean = c(120, 340), sigma) x0 <- rep("C",N) x1 = rand_3[, 1] x2 = rand_3[, 2] X3 <- cbind(data.frame(x0),x1,x2) X_total <- rbind(X1,X2,X3) colnames(X_total) <- c("state","sBP","time")
作成したデータベースの説明と解析目的
解析目的：「食事指導にかけた時間」と「血圧」の関係を調べたいとする1
データベースには
str(X_total) #'data.frame': 3000 obs. of 3 variables: # $ state: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ... # $ sBP : num 149 141 111 122 123 ... # $ time : num 341 359 372 378 358 ... head(X_total,3) # state sBP time #1 A 148.8035 341.0151 #2 A 140.5801 358.7447 #3 A 110.9342 371.7364
- “state”は3つの地域の名称が
- “sBP”には収縮期血圧が
- “time”には年あたりの減塩にかけた指導時間が
入っているとする
stateのイメージとしては、Aは普通の地域、Bは医療リテラシーが低い(不健康な)地域、Cはリテラシーが高い地域とする
実際解析してみる
施設を無視して(多施設をまぜこぜにして)血圧と指導時間を単変量解析し、散布図に直線を引く
res_total <- lm(sBP~time,data=X_total) summary(res_total) ## ## Call: ## lm(formula = sBP ~ time, data = X_total) ## ## Residuals: ## Min 1Q Median 3Q Max ## -58.428 -9.577 0.041 10.122 48.306 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 90.94149 4.64591 19.575 <2e-16 *** ## time 0.10935 0.01288 8.487 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 14.44 on 2998 degrees of freedom ## Multiple R-squared: 0.02346, Adjusted R-squared: 0.02314 ## F-statistic: 72.03 on 1 and 2998 DF, p-value: < 2.2e-16 plot(X_total[,"time"],X_total[,"sBP"]) abline(res_total,col=2)
timeの係数が有意差をもって正になったので、 食事指導にかけた時間が長ければ長いほど、血圧が上がってしまうという事が検証されてしまった。(plotでも傾きは正となっている)
これは施設の違いを考慮せずに解析した事が影響している。
施設を考慮すると？
施設毎の性質(の分の分散)は別のパラメータに引き取ってもらうので、きちんと減塩指導時間の効果が解析できるというイメージ。
まず最初に素データを施設毎に色分けして散布図をみると、
plot(X_total[,"time"],X_total[,"sBP"],col=as.numeric(as.factor(X_total[,"state"])))
散布図を見るだけでもなんとなく違うんじゃないかと思えてくるが、ここに施設ごとに回帰式を作って線を引いてみる
d_a <- subset(X_total,X_total$state %in% "A") res_a <- lm(sBP~time,data=d_a) d_b <- subset(X_total,X_total$state %in% "B") res_b <- lm(sBP~time,data=d_b) d_c <- subset(X_total,X_total$state %in% "C") res_c <- lm(sBP~time,data=d_c) plot(X_total[,"time"],X_total[,"sBP"],col=as.numeric(as.factor(X_total[,"state"]))) c="magenta" abline(res_a,col=c) abline(res_b,col=c) abline(res_c,col=c)
施設毎にみると、指導時間が長い方が血圧は下がるという結果になる。
初めのスクリプトで仮にN=1000000としても変わらないので対象患者数を増やせば良いという話ではないという事はわかるだろう。
解析結果として出力してみる
まずは施設を考慮せずに解析した結果は、
summary(res_total) ## ## Call: ## lm(formula = sBP ~ time, data = X_total) ## ## Residuals: ## Min 1Q Median 3Q Max ## -58.428 -9.577 0.041 10.122 48.306 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 90.94149 4.64591 19.575 <2e-16 *** ## time 0.10935 0.01288 8.487 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 14.44 on 2998 degrees of freedom ## Multiple R-squared: 0.02346, Adjusted R-squared: 0.02314 ## F-statistic: 72.03 on 1 and 2998 DF, p-value: < 2.2e-16
となり、有意差をもって指導時間が増えると血圧が高いという結果になる。
次に施設差を加えてマルチレベル分析をすると、
summary(lme4::lmer(sBP~time+(1|state), data=X_total)) ## Linear mixed model fit by REML ['lmerMod'] ## Formula: sBP ~ time + (1 | state) ## Data: X_total ## ## REML criterion at convergence: 22103.1 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -3.8264 -0.6823 0.0100 0.6831 3.6717 ## ## Random effects: ## Groups Name Variance Std.Dev. ## state (Intercept) 490.23 22.141 ## Residual 92.02 9.593 ## Number of obs: 3000, groups: state, 3 ## ## Fixed effects: ## Estimate Std. Error t value ## (Intercept) 346.15689 13.78917 25.10 ## time -0.59956 0.01435 -41.77 ## ## Correlation of Fixed Effects: ## (Intr) ## time -0.375
となり、初めに設定した共分散の値に近づいたのがわかる。2
マルチレベルで考慮するのは、地域変数・病院変数・外科の手術なら術者毎のカテゴリ・内科なら外来毎のカテゴリなどが挙げられるが、臨床やってる頃は知らんかった。これはRCTをやる時もランダム性を保てない要因に関しては同様に必要だと思う。
