備忘録 as vet.

日々のアイデア、疑問など備忘録的に書きます。Scienceが好きです。

<論文感想>濃厚赤血球輸血を実施した犬における臨床スコア、臨床検査の関連性

pRBCを輸血した犬における血液検査・臨床スコアの関連性

Kisielewicz, C., Self, I. & Bell, R. Assessment of Clinical and Laboratory Variables as a Guide to Packed Red Blood Cell Transfusion of Euvolemic Anemic Dogs. J. Vet. Intern. Med. 28, 576–582 (2014).

https://onlinelibrary.wiley.com/doi/10.1111/jvim.12280

Intro

  • 輸血は副作用やコスト、献血検体の確保の難しさがあるため、輸血適用条件は厳密に検討する必要がある

  • しかしながら、輸血適用条件に関する統一見解は存在せず、Hb濃度を指標に主観的に輸血適否を決定しているのが現状である

  • Hbの他に乳酸濃度の増加とDO2低下が輸血適否のマーカーとして研究されている

  • 乳酸濃度上昇はDO2低下の結果として組織低酸素を示唆する所見

  • 中心静脈の酸素飽和度はDO2の近似指標と考えられるが、検査自体が侵襲的で難しい
  • 末梢静脈酸素濃度(CvO2)はDO2の代替指標となる可能性があるが、検証されていない

本研究の目的

pRBCを輸血された犬におけるHb、Ht、CvO2、乳酸、臨床スコアの変化量を評価することで、輸血実施条件を確立する基盤的情報を提供すること


Methods

  • 単施設Prospective study
  • inclusion: Ht37%>の犬 (n=24)
  • exclusion: 酸素投与、呼吸器・心血管系疾患、呼吸困難症状有り、歩様異常をきたす神経筋/整形学的疾患、脱水、循環血液量減少、高ビリルビン血症、治療介入が必要な輸血副反応が疑われる患者

  • 患者は急性貧血、慢性貧血に区分

  • 輸血必要性は担当獣医師がHt、臨床検査に基づき独立して判断
  • ただし、ADCAS、乳酸値、CvO2はブラインドで取得され、輸血実施判断には使用されない

    →ADCASは担当獣医師とは別の獣医師が評価

  • 輸血による目標PCVは担当獣医師によって独立して決定されるが、施設の推奨目標はPCV20-30%としている

  • 輸血の前投薬は実施せず。

Outcome

  • anemic dog clinical assessment score (ADCAS)

    →可視粘膜色、脈圧、HR、RR、運動不耐性を0-3の四段階、15点満点でスコア化

  • Ht,Hb,CvO2, Lactate


Results

  • n=24 (FS 14, F 1, MC 7, M 2; DEA1.1+ 10, - 12, unkown 2)
  • 6y, median 21kg
  • ADCASは輸血前後で減少(改善)
  • Ht, Hb, CvO2, Lactateは輸血後に改善

Discussion

  • 臨床的に再輸血が必要と思われた症例はADCASでも依然として高いスコアであったことから、ADCASは輸血反応性のスコアリングシステムとして有用かもしれない

  • Hbは輸血反応性を示す指標としてやはり使えるかも;とくにHb <5.8g/dLが閾値になるかも?

  • CvO2はHbとともにDO2の代替指標と考えられるが、今回の結果ではHb以上の情報を示さなかった

  • lactateは輸血後低下したが、顕著な減少ではなく、再輸血の必要性を判断することはできなかった

    →採血のタイミングによって変わるのではないか?


感想

  • 今回のPopulatiuonは、貧血以外の基礎疾患をなるべく除外した中型ー大型犬が対象。
  • 今回の結果では臨床スコアが結構いい感じで、Hbと組み合わせて有用な輸血トリガーになる予感。ただし、考察でも述べられていたが、今回の循環血液量減少や黄疸などがない今回のPopulationだからこそ正確にスコアリングできて有用だった可能性はあるので、外挿したときの妥当性は検証が必要だと思われる。

  • lactateは輸血反応性の指標として有望だと思っていたが、今回の研究ではそこまででもなかった。これもやはり、組織低灌流になりにくい今回のpopulationを反映しているのかも。病態に応じて適切な指標が変わり得るということなのかもしれない。

  • 輸血後の再評価はどのタイミングで行われたか正確な明示がなかった(輸血後すぐと書いてあったが、、)ので、その部分の提示をしてほしかった。

  • 臨床スコアを(おそらく)独立した別の獣医師が評価していたり、Exclusionを厳密に明示していたりと、研究デザインを対処している点もよかった

<論文感想>sparse data biasに対する対処法

Sparse data biasの検出方法と対処方法に関するレビュー

Greenland, S., Mansournia, M. A. & Altman, D. G. Sparse data bias: a problem hiding in plain sight. BMJ 352, i1981 (2016).

https://www.bmj.com/content/352/bmj.i1981

前回読んだLancetのEditorial論文で引用されていたSparse data biasに関する解説論文

多変量回帰分析を行うときなど、多数の共変量を組み込むことで、アウトカムの分布がsparseになる状況が生じ得る

そのような状況では、最尤推定法(MLE)による回帰係数の推定値がNon-Null方向にBiasされる(これをsparse data biasと呼ぶ)

このBiasは: - low Event per variable - アウトカムと共変量の分布が著しく分離している

などによって生じるが、これはトータルのサンプルサイズが大きくても生じ得る

このような問題を対処しないと、仮に交絡因子(共変量)を調整したとしても、交絡によるバイアスの補正よりもspase data biasが上回る可能性がある

sparse data biasの検出方法と対処方法

<検出方法>

通常のMLEによる推定結果と、penalised estimationによる結果を比較したとき、結果に乖離がある場合、biasが示唆される

<対処法とその欠点>

  1. Firth bias adjustment

    sparseなデータがあった場合、アウトカム、曝露の分割表の各セルに0.5を追加するという方法 ある種のpenalised estimationであり、MLEに比べてBiasを軽減することはできるが、強い仮定を強いているのと、調整しても比率に関する推定値では大きなBiasが残る可能性がある →他のPenalisationが推奨される

  2. Stepwiseによる変数選択

    最終モデルの推定量の信頼区間が過剰に狭くなる;推定値がインフレする;重要な交絡因子が除去される可能性がある

  3. Exact logistic regression

    サンプルサイズや共変量数が多いとき、計算負荷がかかる;RR,RDの算出には使えない;過剰に保守的(信頼区間が非常に広くなる)

  4. 曝露モデルを使用する(PS matching, IPTW)

    アウトカムがレアだがExposureが比較的多い場合には使える 一方で、Exposureの発生頻度がアウトカムより少ない、データの分布が極端なときはモデル設定の正確性がアウトカムモデルより劣る可能性あり

  5. Penalisation

    data augmentationという方法でもっともらしい事前分布(信頼区間)を仮定し、実際のデータセットに拡張して組み込むことで、事後的な分布 (信頼区間)のBiasを補正して妥当な推定量を得る方法


頻度論的な立場でpenalaisationを解説しているが、同時にベイジアンの解釈も取り入れており、なんとなくわかる一方、理解は不十分。 そのうちベイズ統計学も本腰を入れて取り組みたいが、まだまだFrequentistの沼は深い。。

実用的には

We thus strongly recommend that basic data numbers within treatment or exposure and outcome categories be examined and presented, and that adjustment methods such as penalisation be applied whenever the numbers of events per covariate fall below four or five. The weighting (degree of penalisation) for each variable is best determined so that the implied prior interval encompasses the full range of reasonable possibilities for the effect of the variable.

この部分を特に意識していきたい。

<論文感想>心臓外科患者における輸血リスクの予測は本当に患者を救うのか?

心臓外科における輸血の予測モデルは臨床的に有用かどうか議論したEditorial論文

Bartoszko, J. & Karkouti, K. Can predicting transfusion in cardiac surgery help patients? Br. J. Anaesth. 119, 350–352 (2017).

https://doi.org/10.1093/bja/aex216

  • 心臓手術では20-60%の患者が周術期に輸血を受けており、輸血は患者アウトカムに悪影響を及ぼすと考えられている。

  • 輸血リスクの指標は、臨床試験での高リスク患者の選択や、メタ解析などでの患者集団のリスク特性の理解に有用である。

  • 術前の正確なリスク層別化は、積極的な術前貧血管理や術中止血管理などの適切な対応につながる可能性がある。

  • 過去10年に多くの輸血リスクスコアが提案されているが、広く臨床現場で使われるには至っていない。

  • ACTA-PORTは、心臓麻酔学会ACTAの協力で大規模多施設データから開発されたという強みがあるが、施設間の輸血方針の違いなどの課題もある。

  • 既存のスコアを比較検討し、臨床での実用性を高める工夫や、スコア使用が実際に臨床行動や患者アウトカムの改善につながるかを検証する必要がある。

  • 麻酔科領域では優れた臨床リスクスコアが多数あるが、それをどのように有効活用するかが今後の課題である。

輸血必要性の予測モデルを構築したところで、結局、臨床現場を変えることはできるのだろうか?という疑問は確かにある。

輸血が必要な患者だから輸血をするわけだし、輸血が必要な患者を高い精度で予測できたところで、そもそも治療方針は大きく変わらないのではないだろうか?

臨床医が事前に予測できないような意外な患者を予測することができるのであれば、有用性は高いと思うが、本論文でも紹介されている予測モデルの予測変数は輸血の予測に寄与しそうな(つまり臨床知見に基づく)因子を集めてモデルを構築しているので、結局臨床医と予測モデルと同じ思考過程でハイリスク患者を同定するだろう。

とすると、臨床医から見て輸血のハイリスクと思われる患者には、予測モデルの結果に関わらず凝固系の精査をするだろうし、事前に出血リスクが低減できるような準備をするだろうし、治療方針は変わらないんじゃないだろうか。

メリットがあるとすれば、熟練度の低い臨床医において、予測モデルを補助として使うことで、ハイリスク患者の見逃しを防ぐことができるのかもしれない。

ただ、そのようなシチュエーションは、例えば状況次第で若手医師しかいないような救急外来などで活用はするだろうが、心臓外科のような予定手術の場合、多くのベテラン医師が多数介在する現場では、そのような予測モデルが活躍する隙はあまりないのではないだろうか。

使われない予測モデルが大量に生み出され、捨てられている現代において、AUCを1%でも高めようとモデル改良をする努力よりも、予測モデルが真に役立つclinical settingを考えるほうが全然重要なんだと思う。

<論文感想>医学研究における統計解析結果の報告方法に関する推奨~Lancetより

医学研究における統計解析結果の報告方法に関する推奨

Mansournia, M. A. & Nazemipour, M. Recommendations for accurate reporting in medical research statistics. Lancet 403, 611–612 (2024).

https://doi.org/10.1016/S0140-6736(24)00139-9

Lancetの編集グループ発、過去3年間に渡る1000本を超える投稿論文のレビュー経験より、統計解析において過不足無く提示すべき情報、適切な報告に関するガイドラインをたった2枚のPDFにまとめて簡潔に提示したCorrespondence。 首肯するしかないもっともなデータの提示法から、知らなかった概念の提示までシンプルに纏めてくれている。各項目において重要な引用文献を提示してくれているので、この論文をアンカーにしてさらに深掘りできたのがとてもよかった。

  1. データの分布形状を確認して、要約統計量をmean+SDか、median+IQRで提示せよ

    Supplementには各変数のヒストグラムやTableを提示せよ

  2. 解析に用いたあらゆる統計モデルの仮定を吟味して、視覚的に妥当性を評価せよ

  3. p値は二分化せず、正確なp値を報告せよ(例:not p<0.05, but p=0.032)

  4. 適切に算出された効果量の95%信頼区間と臨床的重要性を見て効果の程度を解釈せよ;"no effect"は安易に使うな

  5. 交絡因子は有意差検定で選択せず、背景知識に基づき検討し、DAGでその思考を描出せよ

  6. 結果に影響しうるほどデータの欠測がある場合、逆確率重み付けか多重代入法による欠測補完をせよ

  7. 比率推定においてsparse-data biasの存在の評価と適切な対処をせよ

    Greenland et al., BMJ, 2016がこの辺りのわかりやすい解説論文を出している

  8. アウトカム発生頻度が高い場合、ORではなくRR,RDを報告せよ

    modified Poisson regression, regression standardizationでRR,RDを推定できるとのこと。知らなかったので深掘りしてみる。

    Zou, American Journal of Epidemiology, 2004; Greenland, American Journal of Epidemiology, 2004

  9. 統計モデルにおける交互作用項の評価では、乗法的(Odds or Rate)なスケールだけでなく加法的(Risk)で評価せよ

    ここも知らない概念だった。交互作用項には乗法的なもの、加法的なものがあるらしい。

    Knol et al., International Journal of Epidemiology, 2012 を参考にして深掘りしてみる。

<論文感想>低次元データにおける予測モデル構築のための罰則付きパラメータ推定方法の性能比較

低次元二値アウトカム(アウトカム発生数がパラメータ数より十分多くない状況)の予測モデル構築におけるパラメータ推定法に関する解説と手法の比較

Pavlou, M., Ambler, G., Seaman, S., Iorio, M. D. & Omar, R. Z. Review and evaluation of penalised regression methods for risk prediction in low-dimensional data with few events. Stat. Med. 35, 1159–1177 (2016).

https://onlinelibrary.wiley.com/doi/10.1002/sim.6782

予測モデル構築において、予測因子の変数選抜や予測因子の回帰係数の推定の手法はいくつかある。医学分野で一番良く見るのは最尤推定でBackwards stepwiseによる変数選択法だが、問題点も多い。 アウトカム発生数が少ないようなシチュエーションでは、安定した推定自体が困難なときがあるしoverfittingにより予測性能が低下することもある。 パラメータ数に罰則を設け、推定した回帰係数をnull方向にshrinkageすることでoverfittingを低減する手法が提案されている。 様々なシチュエーションにおいて各種手法の性能を比較したReview論文。

shrinkageについては、以下の教科書がわかりやすくて参考になった。*1

Regression Modeling Strategies With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis link.springer.com

Intro

  • 予測モデル構築において、低次元データセット(特にパラメータ数>アウトカム発生数の状況)ではoverfittingが問題になる
  • 医学領域では、パラメータ数とイベント数の比の経験的な数字として10EPV(10 Event Per Variable)を一つの目標として使用されてきたが、これは科学的根拠に基づくわけではない
  • パラメータ選択の方法として従来stepwiseなどがよく用いられているが、変数選抜の再現性の低さ、不安定性などの問題点が多く存在する
  • さらに、パラメータ選択をしてもなお、EPVが低くなってしまう状況もありえる
  • そのようなoverfittingに対処する方法として、shrinkageという統計学的手法がある
  • 代表的なものにRidge、Lassoがあり、高次元データにおいてはbackwards 、univariate selectionに比べて性能が良いことが示されている
  • Lassoなどを拡張したelastic net, adaptive lasso, SCAD などは高次元データにおいて従来のRidge, Lassoより優れた性能を示した
  • しかし、低次元データセットにおける各種手法の優位性を比較した研究は少ない
  • 本論文では低次元データセットにおいて、頻度論的手法(MLE, penalized MLE, Ridge, Lasso, adaptive Lasso, Elastic net, Smoothly clipped absolute deviation[SCAD])とベイジアン手法(Bayesian ridge, Bayesian lasso, Stochastic search variable selection[SSVS])のそれぞれで予測性能を評価する
dataset
  • 陰茎がんの5年生存率に関するデータを使用
  • アウトカム(5年以内死亡)発生数は25。
  • 用いる予測候補因子は9つ:
    • 連続変数ー細胞増殖に関するバイオマーカー(Ki67、Mcm2、Ki67-g95), 診断時の年齢、組織浸潤の深度
    • 二値変数ーリンパ節の転移有無、リンパ管浸潤の有無、腫瘍拡大の程度、DNA倍数性

手法サマリ

  • MLE

    Yをアウトカム、Xを予測変数のベクトル、βを回帰係数のベクトルとするとき、ロジスティク回帰モデルをlogit(E(Y|X)) = βTX

    MLEでは対数尤度比: l(β)が最大となる値βを求めに行く

  • 罰則付きMLE

    l(β) − λpen(β), pen=罰則項、λ=チューニングパラメータ

  • Ridge

    多重共線性に対応するために開発。相関の高い説明変数の回帰係数をNull方向に縮約する

  • Lasso

    回帰係数をNull方向に縮約するが、完全に0になることもあるため、ある種の変数選択とみなせる

  • Elastic net

    RidgeとLassoの組み合わせバージョン。罰則項はridgeとlassoを足し合わせたもの

  • Adaptive lasso

    lassoの各罰則項にweightをかけたもの。weightはridge regressionにおける回帰係数の逆数にすることが多く、結果として回帰係数が大きい(予測寄与度が大きい)因子のshrinkageは小さく、予測寄与度が小さい因子のshrinkageは大きくなる

以下詳細省略(理解できていません) * Smoothly clipped absolute deviation(SCAD) * Bayesian approach * Bayesian ridge * Bayesian lasso * Stochastic search variable selection(SSVS)

results

1)実際の低次元データセットを用いた予測性能比較
  1. calibration slope 予測リスク(確率)を横軸、実際のリスクを縦軸に取り直線を引いたとき、完璧な予測であればslope=1となるが、new dataに対してoverfittingするとslopeは1から遠ざかる。 MLE、BEはともにshrinkageを行わない推定方法であり、slopeは明らかに1を下回っている。他のshrinkageの方法は全体的にslopeを保持しており、overfittingに対応できていると言える。 Ridgeはslope=1.2に近く、ややunderestimateな感じ。

  2. c-statistics (discrimination) いわゆるROC曲線のAUC。真のHigh risk患者を予測モデルでHigh riskと分類できる確率のこと。理想的にはc=1となり、c=0.5だと1/2の確率で分類を間違えることと同じ。 MLEに比べてやや改善しているが、手法によってはMLEとほぼ同性能のものもある(Adaptive Lasso, SCADなど)

  3. RPMSE 患者ごとの予測確率と実測確率の差の二乗和の平均の平方根。つまり予測確率と実測確率にどれくらい差があるかを見ているので、0に近いほど予測性能が高い。 これもMLEに比べていずれのshrinkage手法も良好な性能を示している。

Prediction ability comparison

2)推定リスクの被覆確率比較

MLEに比べて、ridge, lasso, Bayesian Lasso and SSVSのいずれも被覆確率は向上している。 実測確率における被覆確率では、RidgeはHigh riskになるほど被覆確率が低下しているが、Lassoは比較的保持(90%)されている。 Bayesian Lasso and SSVSはHigh riskではLassoと同じような動きだが、Low riskでは顕著に被覆確率が低下していることが特徴的。

coverage probability

3)回帰係数の推定量のバイアス比較

各回帰係数の真の推定量に対して、手法ごとの推定量がどれくらいバイアスを有しているのかを評価。 罰則付き回帰手法は、各回帰係数の推定量をnull方向にshrinkageすることでOverfittingを抑えていたが、言い換えれば回帰係数の推定量はバイアスを受けることと同じ。

真の回帰係数からの変動%で表しているが、そもそも真の回帰係数をどうやって求めたのか論文読んでも分からず。 LassoやAdaptive Lassoでは回帰係数が小さい(予測寄与度合いの低い)因子ではバイアスが大きくshrinkageの程度が大きく、予測寄与度合いの高い因子ではバイアスが小さくなるのがわかる。 一方でRidge, Erastic netなどは比較的バイアスの程度が小さく、予測因子間でもそこまで大きな変動がない。

Bias estimation

この結果のところでなるほど、と思ったのは以下の文:

While unbiased estimation of coefficients is important when the aim is to investigate associations, bias is considered to be a less important issue for risk prediction studies where the predictive performance of the model is of main interest.

いわゆるリスク因子解析のように各共変量自体の回帰係数の推定量に注目できるのは、あくまで不偏推定量のときのみであり、今回の手法のように予測妥当性を担保するために推定量にあえてバイアスをかけるモデルでは、各予測因子の回帰係数を解釈することは重視しない(そもそも解釈が困難だしミスリーディングになりそう)。

予測研究をする上で、結局多変量解析するなら、リスク因子ごとの寄与度(OR)も算出されるから、それぞれのORについて臨床的な解釈も同時に行って、リスク因子解析&予測性能評価の両方やったら一石二鳥じゃん!って思っていたが、そういうわけにはいかなそう。興味のある結果(研究目的)にターゲットを絞って論文構成を考えないといけない理由は、論理の整然化だけでなく、解析上の問題もあるということか。

4)低次元データセットに様々なダミーデータを加えたシミュレーションによる予測性能比較
  1. noise variable 予測に全く寄与しない(回帰係数=0)のnoise variableを加えた状況でシミュレーション
  2. correlated predictors 変数同士が強く相関する多重共線性を持った変数を加えた状況でシミュレーション

結論

ベストな手法があるわけではない。Pros/Consを理解して、状況に応じてベターな手法を選びましょう。という結論。 とはいっても、以下のサマリーを見る限り、Elastic NetかSSVSを選んでおけば良さそうな印象。頻度論的にアプローチするなら、Elastic Net、Bayesian的にいくならSSVSという感じだろうか。

Summary

*1:ちなみに探せば無料でPDFダウンロードできる

Rで予測モデルの性能評価(Cross ValidationとBootstrapの比較)

予測モデルの性能を評価するときは、データの過学習(Overfitting)を考慮する必要がある。

モデル構築で使用したデータに過剰に最適化され、新規データの予測性能が低下する問題が生じるため。

過学習は予測性能を過剰に見積もるバイアスの原因なので、さまざまな方法で過学習の程度を評価し、バイアスを補正した形で予測性能を見る必要がある。

その方法として代表的にはCross Validation(CV)とBootstrapがある。

CV

CVでは、まずデータをモデル構築用のデータと、性能検証用のデータに分割する。

モデル構築用データを元に作成した予測モデルを使い、モデル構築には使用しなかった性能検証用のデータを代入して予測性能を評価する方法。

当然、一回だとデータの分割の仕方によって大きく影響を受けるので、分割領域を次々変えて多数のモデル構築→性能検証を繰り返し、最後に得られた予測指標の平均値を算出する。

欠点は、データを分割しているので、アウトカム発生数が少ないときには予測因子の投入数が制限されることなど。

以下の模式図が分かりやすい。

Iwagami, M. & Matsui, H. Introduction to Clinical Prediction Models. Ann. Clin. Epidemiology 4, 72–80 (2022).

Bootstrap

Bootstrap samplingは、あるオリジナル集団から、重複を許してランダムにサンプリングを行い、オリジナルデータと同数の新しいリサンプリング集団を作出すること。

これを何度も繰り返すことで、オリジナル集団と同数の集団を多数作成できる。

このリサンプリング集団を用いて、予測モデルを作成する。このときリサンプリング回数(B)の数だけ予測モデルが作成されることになる。

構築した予測モデルで、元のリサンプリング集団における予測指標(AUCなど)を算出すると、B個分のAUCが生成される(AUC1_i; i:1~B)。このAUCは当然、同じデータを使ってモデル構築とモデル性能評価も行っているので、過学習されていると考えることができる。

どの程度過学習されたのかを見るために、各予測モデルを使ってオリジナル集団における予測指標を算出する。ここでも予測モデルの数だけAUCが得られる(AUC2_i; i:1~B)ことになる。AUC2は、少なくともAUC1よりは過学習の影響を受けていないはずだから、AUC1とAUC2の差が過学習の程度(optimism)であり、推定値のバイアス分になる。

optimism=AUC1-AUC2もB個分存在するので、平均化したものを算出する。

最後に、オリジナル集団を使って作成した予測モデルを、オリジナル集団に当てはめたときのAUCを算出し、先程のoptimismの平均値を引いたものが、bias corrected AUC estimationになる。

CVではモデル構築に使うデータと性能評価のデータを分けることで過学習に対処していたが、bootstrapではオリジナル集団のサンプルサイズを保ちつつ、過学習の程度を定量化し、過学習によるバイアスを補正することで対処している点が異なるアプローチということになる。

Bootstrapの方法では、AUC2はAUC1に比べて過学習のバイアスは少ないが過学習がゼロになっているわけではないので、得られたoptimismだけでバイアスの補正は十分なのか?という部分が気になる。つまり、bootstrapではオリジナル集団から重複を許したリサンプリングなので、モデル構築で使用したデータと性能評価で使用したデータの構成要素は同じであり、完全に別個のデータではない。AUC2に残っている過学習バイアスを取り除く方法は無いのだろうか。宿題。

CVをRで実装してみる

#k-fold CV 
k <-5 
AUC_CV <-data.frame(matrix(ncol=2,nrow=k)) 
form<-as.formula(
  DAY30 ~ rcs(AGE, 4) + DIA + HYP + HRT + ANT + PMI + 
    FAM + TTR + rcs(WEI, 4) + AGE + ANT*AGE + DIA*HYP + HYP*FAM + HRT*ANT)

#Randomly shuffle the data 
data2r<-data2[sample(nrow(data2)),] 
 
#Create k equally size folds 
folds <- cut(seq(1,nrow(data2r)),breaks=k,labels=FALSE) 
 
#Perform k fold cross validation 
for(i in 1:k){  

  #Segment data2r by fold using the which() function  
testIndexes <- which(folds==i,arr.ind=TRUE) 
testData <- data2r[testIndexes, ] 
trainData <- data2r[-testIndexes, ] 
fit1<-glm(formula = form, family = binomial(link = "logit"), data = trainData)
testData$fitted<-predict(fit1,testData,type="response") 
ROC <- roc(testData$DAY30, testData$fitted) 
AUC_CV[i,2] <- ROC$auc 
AUC_CV[i,1] <- i 

#print(paste(i,"th cross-validation iteration is completed.", sep="")) 
}  
names(AUC_CV)[2] <- "AUC_cv" 
names(AUC_CV)[1] <- "k" 
print(summary(AUC_CV$AUC_cv)) 

CVによるAUCは平均0.8269

Bootstrap samplingをRで実装

B<-100 #number of resampling 
N<-dim(data2)[1] #data2のdimentionを返す→2973行x33列→そのうち一番目(2973)を取得しNとする
AUC.i1<-AUC.i2<-numeric(B) 

for(i in 1:B){  

  #bootstrap random sampling 
  bs.i<-sample(1:N,N,replace=TRUE) #重複を許して1~2973までの数を2973回サンプリングする
  
  #bootstrap samplingに基づき、data2から集団を再構成 →data2.iとする
  data2.i<-data2[bs.i,] 
  
  #bootstrapしたデータセットごとで回帰モデル構築→各bootstrap dataからリスク予測確率の算出→data2.1のprob_1列に格納
  gm1.i<-glm(formula = form, family = binomial(link = "logit"), data = data2.i)
  data2.i$prob_1<-predict(gm1.i,data2.i,type="response") 
  
  #AUCの算出;Bootstrap回数(B=100)分のAUCを算出して結果を格納
  AUC.i1[i]<-roc(DAY30~prob_1,data=data2.i)$auc 
  
  #bootstrapで作成した各回帰モデルを使って、オリジナルデータ(data2)からリスク予測確率を算出→data2のprob_2列に格納;AUC算出
  data2$prob_2<-predict(gm1.i,data2,type="response") 
  AUC.i2[i]<-roc(DAY30~prob_2,data=data2)$auc #ここで得られるAUCはbootstrap glmを用いてoriginal data2に外挿したときのAUC
  
  #print(paste(i,"th bootstrap iteration is completed.",sep="")) 
}  

#optimism(過学習によるバイアス)
opt1<-AUC.i1-AUC.i2 
summary(opt1) 

hist(opt1) 
lam1<-mean(opt1) #estimate of the optimism 

#original dataを使って、originalのstepwise modelからorignialのAUCを算出→AUC1
gm1<-glm(formula = form, family = binomial(link = "logit"), data = data2)
data2$prob1<-predict(gm1,data2,type="response") 
AUC1<-roc(DAY30~prob1,data=data2)$auc #original dataにおけるstepwise model(original model)のAUC

AUC.c1<-AUC1-lam1 #bias corrected AUC estimate
AUC.c1

bias corrected AUC estimateは0.8292でありcross validationのAUC(0.8269)とほぼ一致!

CVとBootstrapでぜんぜん違う過学習に対するアプローチを取っているけども、ほぼ同じような推定値になっているのはすごい。

実際に手を動かしてみるとCVやBootstrapの概念を深く理解できるようになったと思う。やはり論文や教科書を読むだけじゃなくて自分でやってみるのが一番勉強になる。

Rでサンプルサイズと検出力、SDの関係をシミュレーションしてみた

大学院の授業でRを用いたサンプルサイズ(必要症例数)の算出方法のシミュレーションを学んだので、拡張して遊んでみた。

サンプルサイズの求め方はいくつかあるが、一例として

α(有意水準)、β(検出力)、SD(標準偏差)、興味のある効果量の差、を決めることで算出することができる。

通常はα=5%,β=80-90%として、SDは先行研究や手元の自施設データから決めていく事が多い。興味のある効果量の差は臨床的に重要と思われる最小の効果量を定めることで決定していく。

実際にはシミュレーションベースで必要症例数を決定していくことも多いらしい。

ということで症例数とSDを色々変えてみたときのβの変化の様子をシミュレーションで表現してみた。

 

<シチュエーション>

等分散の正規分布のデータを2つ用意。

2群の差を10で固定し、有意水準α=5%としたときの2群の平均値差が0か否か、つまり

H0: 平均値差=0

H1: 平均値差≠0

という仮説検定を考える。この場合はt検定が適用になる。

 

差を10で固定した2つの正規分布に従うデータをランダムに生成するので、当然2群間に差はあるのだが、同じ試行を何百回と繰り返していくと、確率的なゆらぎによって有意差を検出できない状況も生じてくる。

仮に100回検定をした場合に、80回は差を検出できたとき、検出力は80%といえる。

サンプルサイズとデータが従う正規分布標準偏差(データのばらつき)を変えたときに、検出力βがどう変わるのか、をシミュレーションしていく。

 

下準備。

正規分布に従うデータを2個生成し、差が10となるようにする。

上記2群のデータから1000回t検定を繰り返すシミュレーションを実施し、有意水準5%で有意差を検出できた回数/全試行回数からβ(power)を算出するような自作関数を作成。

sim<-
  function(x,y=10){ #x:症例数,y:SD, default=10

set.seed(1)
count<-0
n<-1000

for(i in 1:n){
  a1<-rnorm(x,mean=120,sd=y)
  a2<-rnorm(x,mean=130,sd=y)
  p<-t.test(a1,a2,var.equal=T)$p.value
  p01<-ifelse(p<0.05,1,0)
  count<-count+p01
}

power<-100*(count/n)
return(power)
  
  }

 

次に、症例数、SD(データのばらつき)を変えたときに検出力がどのように変化するのか可視化。今回は症例数を2~50までの範囲で設定。

#症例数と検出力、SDの関連性をグラフ化

res_sd<-function(y){ #y:SD
  
result<-
  matrix(nrow=49,ncol=3,
       dimnames=list(NULL,c("n","power(%)","SD"))
              )

for(x in 2:50){
result[x-1,1]<-x
result[x-1,2]<-sim(x,y)
result[x-1,3]<-y
}

return(result)

}



#SDごとに結果を結合
results<-
  res_sd(1) %>% 
  rbind(
    res_sd(5),res_sd(10),res_sd(20),res_sd(50),res_sd(100)
  )


#結果の図示
results %>% 
  as.data.frame() %>% 
  mutate(SD=as.factor(SD)) %>% 
  ggplot(aes(x=n,y=`power(%)`,color=SD))+
  geom_point()+
  geom_line()

横軸にサンプルサイズ(n)、縦軸に検出力(power)としてSDごとに色分けしてグラフを図示。

  • サンプルサイズが増えるにつれ検出力も増大
  • SDが小さいほうが検出力が上がりやすい

ということがわかる。

概ね予想通りの見た目であったが、想像よりもSDの影響が大きいんだなという印象。

SDが2群の差よりも大きい(SD: 20~)と、検出力を保つためには一気にサンプルサイズがたくさん必要になるっぽい。

簡単なシミュレーションだったが、R codingの良い練習になった。