医療データ奮闘記

公衆衛生大学院に入った内科系専門医が医師として培った現場感と大学院で培った統計の知識を交えながら、医療や疫学や統計に関する素朴な疑問や本音をつらつら書いています。

Reviewerが10人ついた！！？

論文投稿

この間投稿した医学系の論文だが、最終的にアクセプトされたのだが、はじめにreviewerが10人ついて驚いた。

周りの割と偉い人たちに聞いても誰一人そんな事経験した事がないと言っていた。1

今となってはおそらくEditorが通したかったのだろうという結論に至っているが、当初は意図がわからず結構不安だった2。

初めの４人くらいは否定的だったけど、後半につれてだんだん好意的なコメントが増えて、最後の二人はかなり中立な意見だった。

内容も人によって指摘の仕方も全く違ったし、結局初めの解析の5倍くらい追加解析して通った3。

最終的にはワードファイルで一旦27枚くらい(その後調整して24枚くらい)になったので、いっときは本文より多かった。これで落としたら人としてどうなん、というくらい誠実に、ほとんどの質問に何らかの文献引用か解析結果を用いてきっちり書いた。

major reviseで人数が多くボロボロに書かれていても、諦めなければなんとかなるもんだ、という結論に至った。

  1. 教授も５人くらい聞いたが、みんなびっくりしていた。top journalでもせいぜい４人とか統計家もついて５人だろうと。↩

  2. しかも致命的なスクリプトミスyiliaojingji.hatenablog.comがあって、revise出す時は不安の頂点に達した。↩

  3. 傾向スコア解析もするべき、とか一行で書かれるけど、methodとresultにも書かなあかんし、マッチング前後のSDのテーブルも書かなあかんし。。)↩

ベイズの考え方をRで簡単に実験してみよう。

R 統計

ベイズの勉強をしたいと思ったが、何をみてもいまいち実感が湧かなかったため、Rで実験したくなった。(初めの段階の理解は以下の記事参照)

yiliaojingji.hatenablog.com

頻度論とベイズは同じ式の形をしてても（同じ確率密度関数の形）解釈が違う。例えば以下のHPで説明されている

to-kei.net

しかし私自身が実際やってみないと何を言っているのかわからなかったためRで手を動かして実験してみる1。

まずは頻度論の考え方

目の前にコインがあるとする。一般的にはこのコインが表になる確率は1/2だろうと考える(仮定する)。5回くらいコインを投げて全部表が出た場合に、そんな確率は0.03125(=1/32：3%ちょっと)しかないんだからこのコインはいかさまだ！！と考えるのが頻度論の考えである。

頻度論ではこのコインはまずは「理想的なコインである」と仮定して（考えて）おり、理想的なコインとしてふさわしくない行動を取るとそれはいかさまコインだと考える。

分布を考えて(plotを書いて)、それより低い所の総和(面積)が0.05%以下となっている所なのか否かを検定する事になる。

ベイズ的に考えるとこうなる。

目の前に得体の知れないコインがあるとする。 このコインがどのくらいの確率で表が出るのか(表が出る真の確率であるpの値がいくらと考えるのが妥当か)を実際やってみて(データから)推定したい。

まずは事前情報がないので、事前分布(今回は共役事前分布2としてベータ分布を使う)を一様分布：uniとする

pには「このコインが本当はどのくらいの確率で表が出やすいんだろうか」という確率が入っており、縦軸の高さはわかりやすく言うと「そのp値(表の出やすさ)が妥当な値(真の値)と決めてしまっても良いよーという度合い」が入っている3：自分は理解の上でここが一番難しかった

p <- seq(from=0,to=1,length=100)
uni <- dbeta(p,1,1)
plot(p,uni)

f:id:yiliaojingji:20190509142020p:plain
uni

初めの状態は、何の情報もないので、均一(どのp値も同じくらい真の値でありうる)


次に、このコインを5回回すと、表が3回・裏が2回でたとする。

aft <- dbeta(p,1+3,1+2)
plot(p,aft)

f:id:yiliaojingji:20190509142134p:plain
aft

表が出やすいかもしれないという程度のグラフになる。(細かく言うと、今ある情報だけから考えると、pは0.5よりちょっと大きめなコインなのかもしれないな、でも他の値である確率も十分高いな)

例えばここで、本当はp=0.2じゃないか？(イカサマコインだから表は5回に1回しか出ない,というのがp=0.2の意味)とか言われると、まだその可能性も十分あるよね、もう少しやってみないとわからないよね、と返事する事になる。

さっきは表が3回・裏が2回 (3,2)としたが、自分でこの値を(4,0)とか(4,1)とか(5,1)とか自分で実験していくと、だいぶイメージがついてくると思う4。


そしてそんな数回ではなくコインをどんどん回してみて、最終的には例えば表が1000回、裏が1000回出たとするとこんなグラフに成長する。

matomo <- dbeta(p,1+1000,1+1000)
plot(p,matomo,type="l") #グラフを直線表記にした

f:id:yiliaojingji:20190509142357p:plain
aft1000

だいぶ尖っている。

これをみてさっきのグラフと比較すると、p=0.5なのだろう(表が出る確率が0.5なのだろう)という考えに近づいていく事が可視化される。

つまりこれはp=0.5のまともなコインなのだろう、という結論に近づく。

ここでp=0.2じゃないか？とか言われると、さすがにそれはないんじゃないかな、といえる。

結論

どっちが良いかという結論はないだろうけど、

  • 医学系の論文ではほとんど頻度論の手法がメイン
  • 最近の統計手法(機械学習など)はベイズの手法の方がどんどん使用されている

気がするので、今後置き換わってくるかもしれない。ビジネスに使うなら人間の考え方に近いベイズなのではないだろうか。

  1. 学んだことの備忘録＋自分なりの理解の追記ですので、間違いあれば教えて下さい。↩

  2. データを入れても事後分布と事前分布の属性が変わらないようにする分布を言う。調べてもよく解らんかったら、とりあえずはそんなに気にしなくて良い。↩

  3. ゼロに近いと、まあその値はとらないだろうという解釈になる↩

  4. こういうのはshinyで実装したら面白いので、いつかスクリプトをupする予定。縦軸の絶対値が何かは結構難しいので、今回の記事では割愛↩

公衆衛生大学院後輩指導論文作成解析用資料(疫学)

R 統計 後輩レクチャー

Rで解析するための夏休み強化学習50本ノックというものを作った。

自己学習用の資料で、このままやれば50本ノックした際には理解も深まって論文ができている1というものである。答えもR Markdownで作っているのだが、使用したデータベースは公表できないため、ネットで公開する時にはデータを自作で作る必要がある2。一度解析したら理解を深めるには自分が設定した事が答えとなるデータベースを自作するようになるのが一番の理解の深め方だと思っているので、その辺りも実際後輩に対する授業を進め、反応を見ながら今後は少しづつ記事にしていこうと思う。

どのパッケージ使うかをgoogleとかで検索しながらまあこれだけできたら、医療系の比較疫学研究に関して最低限の解析はできたと言えるのではないだろうか。

  1. 今から行おうとしている研究のデザインを簡潔に述べよ(ADL、退院時後遺症の推定のためにどのような数理モデルを作成するのか)
  2. 練習用データ内のcsvファイルをRに取り込め
  3. 5人分のデータを出して、各説明変数の定義を述べよ
  4. 各説明変数毎のデータの型を確認せよ
  5. 必要に応じて自然数のデータを因子型のデータに変更せよ
  6. 性別を二値変数に変更せよ
  7. 年齢をヒストグラムにせよ
  8. 年齢を4群に分けた説明変数を作成せよ(分け方は根拠があれば何でも良い)
  9. 年齢を10群に各群均等になるように分けた説明変数を作成せよ
  10. 治療法と各説明変数のクロス集計を全て作成せよ
  11. 治療法毎の、年齢のヒストグラムを描け
  12. table.1のデータをまとめよ
  13. table.1をwordファイルで完成させよ
  14. 最小二乗法に関して簡潔に説明せよ
  15. ADLの分布をヒストグラムで表せ
  16. ADLを治療方法のみを使用して予想する線形回帰モデルを作成し、その係数の解釈に関して説明せよ
  17. ADLと年齢(連続変数)の関係を散布図で書き表せ
  18. ADLを年齢(連続変数)のみを使用して予想する線形回帰モデルを作成し、その係数の解釈に関して説明せよ
  19. ADLを年齢(連続変数)のみを使用して予想した線形回帰モデルの当てはまりの良さに関して述べよ
  20. ADLを年齢(9のカテゴリ変数)のみを使用して予想する線形回帰モデルを作成し、その係数の解釈に関して説明せよ
  21. ADLを年齢(連続変数)のみを使用して予想する線形回帰モデルを作成し、その係数の解釈に関して説明せよ
  22. ADLを年齢(連続変数)と性別のみを使用して予想する線形回帰モデルを作成し、その係数の解釈に関して説明せよ
  23. ADLを年齢(連続変数)のみを使用して予想した線形回帰モデルの当てはまりの良さに関して18と比較せよ
  24. ADLを年齢(連続変数)のみを使用して予想した線形回帰モデルの当てはまりの良さに関して18との比較を検定せよ
  25. ADLを年齢(連続変数)、性別、高血圧、糖尿病、脳卒中、心筋梗塞の有無を説明変数とした線形回帰モデルで予想せよ
  26. 一般化線形回帰モデルと一般線形回帰モデルの違いに関して述べよ
  27. 確率分布・リンク関数・線形予測子に関して簡潔に述べよ
  28. ロジスティック回帰モデルに関して簡潔に述べよ
  29. 退院時後遺症の有無をtable.1に追加せよ
  30. 退院時後遺症の有無に関する治療方法の係数をロジスティック単回帰モデルで求めよ
  31. 30の係数をオッズ比に変えよ
  32. 25と同様の説明変数を用いて、31の多変数解析を行い、治療方法によるオッズ比を求めよ
  33. 32のロジスティック回帰分析のROC曲線を描き、c統計量を求めよ
  34. 区間推定を踏まえて、全説明変数のオッズ比をtable.2としてまとめよ
  35. table.2をwordファイルで完成させよ
  36. 傾向スコアの概念を簡潔に述べ、傾向スコアを活用した解析方法を3つ以上述べよ
  37. 32でoutcomeを治療方法として多変数のロジスティック回帰分析を行え
  38. 37の傾向スコアを取り出して、元のデータベースに貼り付けよ
  39. 38の傾向スコアを治療法毎にヒストグラムにせよ(重なりがわかるように)
  40. 傾向スコアのみを交絡因子として調整し、退院時後遺症の有無に関する治療方法のオッズ比を求めよ
  41. 傾向スコアマッチングに関して1:1などの比率・復元・キャリパの3つの言葉を用いて簡潔に述べよ
  42. 傾向スコアマッチングを行い、推定リスク比を算出せよ
  43. SDを用いてマッチング前後の両群のバランスチェックを行え
  44. マッチング後のデータを用いて39のヒストグラムを再度作成し、39と比較せよ
  45. マッチング前後で一般化妥当性が大きく損なわれていないか論じよ
  46. 傾向スコア逆確率による重み付けに関して、外れ値という言葉も用いて簡潔に特徴を述べよ
  47. 傾向スコア逆確率による重み付けを行い、治療効果のオッズ比を算出せよ
  48. 傾向スコア逆確率による重み付け前後のバランスチェックを43と同様に行え
  49. 40・42・47の結果の違いに関して考察せよ
  50. 32・49の結果を合わせてtable.3を作成し、wordファイルにて完成させよ

  1. introductionやdiscussionは自分で作らなあかんけど。↩

  2. 個人情報保護法などもあるし。↩

lasso回帰をする時に二値変数と連続変数を同時に扱うにはどうすれば良いのか？

R

と思って調べた。 英語ではいくつかあるが、日本語ではなかなか見当たらなかったし、周りの人1に聞いても意外と誰も答えられなかった。

英語では例えば

https://stats.stackexchange.com/questions/69568/whether-to-rescale-indicator-binary-dummy-predictors-for-lasso

である。

簡単に書くとまあ結論は出ていないので、ダミー変数にすれば良いのかな、と現状では理解した。ちなみに筆者は確かにglmnetを使用している。

係数が問題になりそうとの事だが、機械学習やっている人の意見では、予測精度が良ければそんなに気にしないのでは？との事だった。

特徴量の選択を目的に使用するならそれで良いのかなー、でもダミー変数にして片方だけ落ちたらどうしようとか思った。ダミー変数にする時にどうせスプラインとか書くし、総合評価かなーとも思う。

最終的にモデルをどのように評価するかに依るのかもしれないが、現段階ではノーフリーランチ定理を以て言及を避けるという結論に行き着いている。

  1. データ扱っているとこの教員の先生とか研究所の人とか↩

SPSSで年度別データをまとめる方法

SPSS

筆者はSPSSはあまり使わないが、仕事上どうしても必要になった事がある。

SPSSで各年度のデータをひとまとめにする(unionの方法に相当)方法を記載しておく。

まずはデータセットを３つ用意する

データセットの説明

  • id ：個人特定番号
  • va1 ：個人の情報（その年の所属など：年度中に変更の可能性あり）
  • va2 ：何年のデータか
  • va3 ：その年に何点だったか
  • va4：その人の属性（性別など）

1つ目 test3_1

f:id:yiliaojingji:20190507111914p:plain
SPSS_001_1

2つ目 test3_2

f:id:yiliaojingji:20190507112143p:plain
SPSS_001_2

3つ目

f:id:yiliaojingji:20190507112330p:plain
SPSS_001_3

これらをまとめて１行ずつに

ファイルの結合を選び、

f:id:yiliaojingji:20190507112437p:plain
SPSS_001_4

ケースの追加を選ぶと、

f:id:yiliaojingji:20190507112511p:plain
SPSS_001_5

どの変数をどのように残すかを聞かれるので、追加変数を設定する。

f:id:yiliaojingji:20190507130940p:plain
SPSS_001_9

あとはMatchsequenceを付けて1、一人一行にすれば良い。(distinct)

f:id:yiliaojingji:20190507130217p:plain
SPSS_001_7

重複があると勝手に横に並ぶようになる。

f:id:yiliaojingji:20190507130632p:plain
SPSS_001_8

  1. 別の記事で説明します↩

対数変換(log変換)して有意差がなくなる事もある？？

R

stepwise法で(AICでなく)p値を使って変数選択をする際に対数変換するとしないのとで有意差が変わってしまうのは問題ではないかとカンファレンスで指摘した事がある。

指摘しておいて不安になったので実際やってみた。

height <- c(158,162,177,173,166,168)
hei_log <- log(height)
length <- c(300,320,330,340,320,336)
data_set <- data.frame(HEIGHT=height,LOG=hei_log,GOAL=length)
summary(lm(length~height,data=data_set))

## Call:
## lm(formula = length ~ height, data = data_set)

## Residuals:
##       1       2       3       4       5       6 
##  -8.710   4.595 -10.515   6.181  -2.101  10.551 

## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  44.2219   102.0472   0.433   0.6871  
## height        1.6740     0.6094   2.747   0.0515 .
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## Residual standard error: 9.506 on 4 degrees of freedom
## Multiple R-squared:  0.6535,    Adjusted R-squared:  0.5669 
## F-statistic: 7.545 on 1 and 4 DF,  p-value: 0.05154

summary(lm(length~hei_log,data=data_set))

## Call:
## lm(formula = length ~ hei_log, data = data_set)

## Residuals:
##       1       2       3       4       5       6 
##  -8.292   4.631 -10.436   6.034  -2.274  10.336 

## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -1124.8      512.1  -2.196   0.0930 .
## hei_log        283.1      100.0   2.830   0.0473 *
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## Residual standard error: 9.321 on 4 degrees of freedom
## Multiple R-squared:  0.6669,    Adjusted R-squared:  0.5836 
## F-statistic: 8.008 on 1 and 4 DF,  p-value: 0.04735

まとめると、

変数 logなし(height) logあり(hei_log)
p value 0.0515 0.0473 *

という事で、「ありうるはありうる」という結論になった。

なぜそんな衣ついた表現なのかと言うと、ありうる状況を作るまでに身長の値を何度も変える必要があったからである。なので、ぎりぎりの有意差(0.05あたり)の変数がある場合は注意する、という程度で良いだろう。

確かに「その値が信頼区間の間で常に正かどうか」などを判断する場合には意味合いもそんなに変わらないのだろうという結論に至った。

心臓マッサージ(胸骨圧迫)をするべきかどうかで困るのでは？ (心肺停止の判断基準)

医療

今日は平和かなーと思っていたが、そう思った日に限って突然心肺停止などが起こるものである1。

私は目の前での心肺停止を結構対応しているので、皮肉を込めて「名探偵コナンみたいやな」と言われた事もある。
院外(駅)で1回、院内では5回くらいfirstで急変対応(Dr.callでの対応を除く)をしている。(一般内科の医師としては結構多い方)初めはパニックになることもあったが、回数を重ねるにつれて「自分マニュアル」ができてきて、落ち着いて行動できるようになってくる。逆にそこまで決めておかないと、いつまで経っても心肺停止に対してなんとなく恐怖心が残ってしまう。

複数の経験を経て辿り着いた個人的なポイントとしては、

  • 呼吸確認の際は下顎呼吸は意味ないので、「寝ている時のリズミカルな呼吸と明らかに違う」かどうか
  • 脈の確認は首の上の方(顎を片方ずつ掴む感じ)を触って触れるか否か & そもそも血圧が測れるか否か
  • 意識障害の確認は胸をグリグリしながらどんな反応するかをみる(普段普通に喋る人が、痛み刺激にやっと反応する程度なら胸骨圧迫した方が良い。私が今まで胸骨圧迫をし始めた人はたいてい痛み刺激の確認に反応はしていたが、蘇生後に聞いても心臓マッサージされてた事を全員覚えてない)

院内であれば、

  • とりあえず「おいも」(O2の投与・IV(静脈路)の準備・monitor管理(血圧とSpO2のモニタリング))を指示して、何時に何をしたかの記録をつけながら人を呼ぶ(いずれも悪いことはない)
  • バイタルが落ち着いても、とりあえず心疾患が背後にあるか否かは判断が難しく専門外の人間には最後までrule outが困難なので、わからんかったらバイタルを落ち着けて循環器内科callかどこかのERに送る。

このくらいでパニックにはならなくなる。あとはボスミンの投与間隔などは自分で覚えやすいように覚えておけば良い

以前、論文発表で「心臓マッサージのやり方を学んだ人の中で、50%以下の人しか現場で実践できていない」という発表を聞いた。発表者である救急医は驚いていたが、そんなもんだろうと思う。
非医療者が一番困るのは「この人は本当に心臓が止まっているのかどうか」という判断をどのように下すかではないだろうか。これは研修医くらいでも結構難しい。
よく見るパターンは基本的にパニックになっていて、「多分何かの間違いだろう」と思い、とりあえずなんとなく傍観してしまう。ベテランの看護師さんや誰かが胸骨圧迫(心臓マッサージ)をしだすと、慌てて色々と手伝うことになる。(自分の研修医時代のことでもある)

これが済んで落ち着けば、採血・心電図・心エコー・頭部＋胸腹部CT(解離とか疑う)くらいは念仏のようにとっておいても文句は言われない(心肺停止したくらいなんだから、明らかに餅詰めたとか窒息でなければむしろ堂々とここまでしましたと報告したら良いと思う）

  1. 外勤中だった↩