はじめに

TRIPODでは予測モデル作成において、欠測に対して、多重補完することをすすめている。 実施にMICEから予測モデル作成までは一般的な解析Codeで実施可能であるが、予測性能の評価に関しては、どのタイミングで統合すべきか、実際のCodeはどうなるのかわからなかった。 このため調べてみることにした。

統合のタイミング

Steyerberg 2019によれば正規分布の測定値は、インプットされたデータセット間でRubinのルールと容易に組み合わせることができる。例えば、モデルの切片を使ったCalibration-in-the-largeの推定値や、線形予測変数の回帰係数を使ったCalibration Slopeなどはそのまま平均をとって統合可能である。

AUCROC、または一致統計量（c）は、0と1の間の範囲をとる指標であり、ロジット変換することで複数の補完されたデータセットから、またはメタアナリシスの文脈で、推定値の統合が可能である。方法を記載した。

補完後の統合の方法

パラメーター Rubin’s Rule
回帰係数 そのまま統合
検定の推定値 Wald検定
患者毎の線形予測因子 そのまま統合
Calibration in the large (Calibration plot) そのまま統合
Calibration slope そのまま統合
AUC or c-statistic logit(c)

実際にSteif 2022をみると各データセットでAUCを算出しそれを統合している (Steif 2022)。
https://journals.lww.com/pccmjournal/Fulltext/2022/01000/Prediction_Model_Performance_With_Different.16.aspx

R codeがGit上で公開されている。https://github.com/part-cw/PCCM_2021_ImputationModel
Codeを確認するとAUCを変換せずそのまま統合していそうであった。

NRI、IDIは？

Kristinらの研究ではNRIは各補完インプットされたデータセットで別々に決定され、平均値として結合したと記載されている。（Rubinの法則）（Kristin JCE 2017）。

IDIに関してはDiscrimination slopesの予測差（イベント群と非イベント群の平均予測確率の差の差）なので単純に平均で良さそうだ。

補完の課題

補完による統合の課題は、欠損値が補完されたときの性能の楽観性の推定にあるようだ。

内的検証ではクロスバリデーションやブートストラップを行い楽観性補正を行うが、どのタイミングで補完するかはまだ議論されている。 一つの方法は、各補完されたデータセット内で補完する方法である。これは一般的で簡単なアプローチで、補完されたデータセット内のモデル性能の測定値を決定し、その後、モデル性能の全体的な測定値に対してこれらをプールする方法である。
もう一つの方法は、各ブートストラップサンプル内で補完する方法である。このようなブートストラップ後の複数のインピュテーションは、計算量が多くなる。特に、多くの欠損値を持つ小さなデータセットでは、かなりの数のインピュテーションが必要になり、差が生じるようだ (Steyeberg 2019 p150, Simone BMC Medical Research Methodology 2016)。

psfmi package

このパッケージは、Rubin’s Rules (RR) を用いて、多重補完されたデータセットに対して、線形、ロジスティック、Cox回帰モデルのプーリング、後方選択、前方選択を適用する関数を提供するものである。

Martijn W Heymans (2021). psfmi: Prediction Model Pooling, Selection and Performance Evaluation Across Multiply Imputed Datasets. R package version 1.0.0. https://mwheymans.github.io/psfmi/

予測モデルの検証は、クロスバリデーションまたは多重補完されたデータセット間のブートストラップによって実行でき、AUC値、R二乗、Hosmer and Lemeshow検定、尺度付きBrierスコア、キャリブレーションプロットなどのプールされたモデルのパフォーマンス測定が生成できる。また、複数の補完されたデータセットにまたがるロジスティック予測モデルを外部で検証する機能や、多重補完されたデータにおけるモデルを比較する機能が利用できる。

以下の論文で検証されている。 Eekhout (2017) doi:10.1186/s12874-017-0404-7.
Wiel (2009) doi:10.1093/biostatistics/kxp011.
Marshall (2009) doi:10.1186/1471-2288-9-57.

psfmiの使用

以降は、miceとpsfmiを使った（ロジスティック）回帰予測モデルの開発、内部検証、外部検証の方法を示していく。

Step1 Instal packages

インストール方法

#install.packages("psfmi")
library(psfmi)
library(mice)

Step 2 - Impute the missing data with mice

psfmiパッケージに含まれるlbp_origデータセット（159名、15変数のデータ）を使用する。

data(lbp_orig)
head(lbp_orig)
#see help(lbp_orig)

このデータ欠損値を処理するために、5つの補完データセットを生成する。
欠損値の量に応じて、補完するデータセットの数を増やすこと（欠損割合、デフォルトで200など）。

imp <- mice(lbp_orig, m=5, maxit=10, seed = 750)

miceパッケージのcomplete関数を使用して、補完されたデータセットを抽出する。
action = “long” and include = FALSEとすることで、インプットされたデータセットを下に重ねて一つの長いデータセットとし、（欠損値を含む）元のデータセットはその長いデータセットに含まれないようにすることができる。

data_comp <- complete(imp, action = "long", include = FALSE)

Step 3 - Develop the model

5つの多重補完されたデータセットに対して後方選択（p-valueが0.05で、選択手法D1）を実行するために、psfmiパッケージのpsfmi_lr関数を使用する。

この関数では、後方選択時に非線形関係の場合に三次スプライン係数を含めることが可能で、予測変数は設定 keep.predictors で含めることによって後方選択時にモデルに強制的に入れることができる。設定方向をFWに変更すると、前方選択が行われる

プールされた係数、標準誤差、95 信頼区間、p 値を導き出すための基本的なプーリング手順は、ルービンの法則（RR）である。しかし，RRは，モデルが連続変数または2値変数を含む場合にのみ可能である．
モデルがカテゴリ変数（> 2 カテゴリ）または制限付き3次スプライン変数も含む場合，特別な手順が利用可能である．

これらのプーリング方法は，以下のとおりである．

D1: 全共分散行列のプーリング
D2: はカイ2乗値のプーリング D3とD4: 尤度比統計量（Meng と Rubinの方法）のプーリング
MPR: 中央値p値のプーリング（MPRルール）である．

スプライン回帰係数は、rmsパッケージのrcs関数（restrictedcubic splines）を用いて定義。knotsで定義した3ノット以上が必要。

pool_lr <- psfmi_lr(data=data_comp, nimp=5, impvar=".imp", 
              formula=Chronic ~ Gender + Age + JobControl + Tampascale + 
              Pain + Radiation + JobDemands + SocialSupport + Smoking + 
              factor(Satisfaction) + factor(Carrying) + rcs(Function, 3), 
              keep.predictors = "Radiation", 
              p.crit = 0.157, #P-value  selection  criterium. 
              method="D1", 
              direction = "BW")
## Removed at Step 1 is - Gender
## Removed at Step 2 is - JobDemands
## Removed at Step 3 is - Smoking
## Removed at Step 4 is - Age
## Removed at Step 5 is - Tampascale
## Removed at Step 6 is - rcs(Function,3)
## Removed at Step 7 is - JobControl
## Removed at Step 8 is - SocialSupport
## 
## Selection correctly terminated, 
## No more variables removed from the model

モデルの確認

pool_lr$RR_model_final
## $`Step 9`
##                    term   estimate std.error  statistic        df      p.value
## 1           (Intercept) -5.6290593 1.1754061 -4.7890336  39.95530 2.318419e-05
## 2                  Pain  0.9354091 0.1851492  5.0521901  67.68397 3.532584e-06
## 3             Radiation  0.5186631 0.5151459  1.0068276  91.48901 3.166731e-01
## 4 factor(Satisfaction)2 -0.4282950 0.6531511 -0.6557364  44.49919 5.153699e-01
## 5 factor(Satisfaction)3 -2.6431297 0.9429219 -2.8031270  26.83928 9.283867e-03
## 6     factor(Carrying)2  1.2554087 0.6254666  2.0071556  78.40085 4.817958e-02
## 7     factor(Carrying)3  1.6871723 0.6473849  2.6061349 115.47451 1.036350e-02
##            OR    lower.EXP  upper.EXP
## 1 0.003591953 0.0003338808  0.0386429
## 2 2.548255652 1.7610642412  3.6873197
## 3 1.679780430 0.6037849102  4.6732905
## 4 0.651619185 0.1747816872  2.4293596
## 5 0.071138280 0.0102714430  0.4926917
## 6 3.509272384 1.0103539396 12.1887907
## 7 5.404177760 1.4991278521 19.4814186
pool_lr$multiparm_final
## $`Step 9`
##                       p-values D1 F-statistic
## Pain                 1.275830e-06   25.524625
## Radiation            3.148421e-01    1.013702
## factor(Satisfaction) 2.426042e-02    4.267156
## factor(Carrying)     3.261509e-02    3.483376

Step 4 - Determine model performance

ステップ 3 で選択したモデルの性能は、pool_performance 関数を用いて求めることができる。

perf <- pool_performance(data=data_comp, nimp=5, impvar=".imp", 
    formula = pool_lr$formula_final[[1]], 
    cal.plot=TRUE, plot.method="mean", groups_cal = 10)

perf
## $ROC_pooled
##                     95% Low C-statistic 95% Up
## C-statistic (logit)  0.8144      0.8978 0.9462
## 
## $coef_pooled
##           (Intercept)                  Pain             Radiation 
##            -5.6290593             0.9354091             0.5186631 
## factor(Satisfaction)2 factor(Satisfaction)3     factor(Carrying)2 
##            -0.4282950            -2.6431297             1.2554087 
##     factor(Carrying)3 
##             1.6871723 
## 
## $R2_pooled
## [1] 0.5842591
## 
## $Brier_Scaled_pooled
## [1] 0.4882759
## 
## $nimp
## [1] 5
## 
## $HLtest_pooled
##         D2          p        df1        df2 
##  0.7163405  0.6758435  8.0000000 42.8715173 
## 
## $model_type
## [1] "binomial"

Step 5 - Internally validate the model

モデルの内部検証を行うために、psfmi_関数を使用します。
この関数では、MIデータでモデルを内部検証するために5つの異なる手法を使用することができます（Vignettesを参照 https://mwheymans.github.io/psfmi/articles/）。

ここでは、クロスバリデーション、より具体的にはcv_MI_RRという手法を使用する。

この手法では、クロスバリデーションの手順に後方選択を組み込むことができる。

psfmi_lr関数によって選択された最後のモデルは、内部的に検証できる。したがって、フルモデルからクロスバリデーション中にバックワード選択を適用したい場合、まず、変数選択なしでpsfmi_lr関数を適用しなければならない。予測モデルがオーバーフィットする主な理由の1つは変数選択であることが分かっているので、ここではこれを適用する。

pool_val <- psfmi_lr(data=data_comp, 
                     formula = Chronic ~ Gender + Age + JobControl + Tampascale + 
                       Pain + Radiation + JobDemands + SocialSupport + Smoking +
                       factor(Satisfaction) + 
                       factor(Carrying) + rcs(Function, 3), p.crit = 1,
                     direction="BW",nimp=5, impvar=".imp", method="D1")

set.seed(200)

res_cv <- psfmi_validate(pool_val, val_method = "cv_MI_RR", data_orig = lbp_orig,
                         folds = 5,p.crit=0.05, BW=TRUE, nimp_mice = 10, 
                         miceImp = miceImp, printFlag = FALSE)
## 
## fold 1
## Removed at Step 1 is - Age
## Removed at Step 2 is - Smoking
## Removed at Step 3 is - Gender
## Removed at Step 4 is - Tampascale
## Removed at Step 5 is - JobDemands
## Removed at Step 6 is - rcs(Function,3)
## Removed at Step 7 is - SocialSupport
## Removed at Step 8 is - JobControl
## Removed at Step 9 is - Radiation
## Removed at Step 10 is - factor(Carrying)
## 
## Selection correctly terminated, 
## No more variables removed from the model
## 
## fold 2
## Removed at Step 1 is - Gender
## Removed at Step 2 is - rcs(Function,3)
## Removed at Step 3 is - Radiation
## Removed at Step 4 is - Age
## Removed at Step 5 is - JobControl
## Removed at Step 6 is - Smoking
## Removed at Step 7 is - SocialSupport
## Removed at Step 8 is - JobDemands
## Removed at Step 9 is - Tampascale
## Removed at Step 10 is - factor(Carrying)
## Removed at Step 11 is - factor(Satisfaction)
## 
## Selection correctly terminated, 
## No more variables removed from the model
## 
## fold 3
## Removed at Step 1 is - Gender
## Removed at Step 2 is - JobDemands
## Removed at Step 3 is - Tampascale
## Removed at Step 4 is - Smoking
## Removed at Step 5 is - SocialSupport
## Removed at Step 6 is - JobControl
## Removed at Step 7 is - Age
## Removed at Step 8 is - rcs(Function,3)
## Removed at Step 9 is - Radiation
## 
## Selection correctly terminated, 
## No more variables removed from the model
## 
## fold 4
## Removed at Step 1 is - Smoking
## Removed at Step 2 is - Gender
## Removed at Step 3 is - JobControl
## Removed at Step 4 is - Radiation
## Removed at Step 5 is - JobDemands
## Removed at Step 6 is - Age
## Removed at Step 7 is - factor(Carrying)
## Removed at Step 8 is - SocialSupport
## Removed at Step 9 is - Tampascale
## 
## Selection correctly terminated, 
## No more variables removed from the model
## 
## fold 5
## Removed at Step 1 is - JobDemands
## Removed at Step 2 is - Gender
## Removed at Step 3 is - Smoking
## Removed at Step 4 is - Radiation
## Removed at Step 5 is - Tampascale
## Removed at Step 6 is - rcs(Function,3)
## Removed at Step 7 is - Age
## Removed at Step 8 is - SocialSupport
## Removed at Step 9 is - JobControl
## 
## Selection correctly terminated, 
## No more variables removed from the model

結果は以下にまとまる。　　 統合された結果として出てくる。　　 　　

res_cv
## $stats
##                  Train      Test
## AUC          0.8916129 0.8410027
## Brier scaled 0.4737013 0.3051174
## Rsq          0.5631277 0.4927999
## 
## $slope
##  Intercept      Slope 
## -0.0943762  0.8374978

Step 6 - Shrink pooled coefficients and determine new intercept

新しいデータでモデルがオーバーフィットしないように、前のステップで推定されたスロープ値0.849148を縮小係数として使用することができる。
これは、プールされた係数に収縮係数を掛け、さらに収縮した係数と一致する新しい切片の値を決定することで行う。

このためにpool_intadj関数を使用し、調整された係数と新しい切片の値を算出する。

pool_select <- psfmi_lr(data=data_comp, nimp=5, impvar=".imp", 
                        formula = pool_lr$formula_final[[1]],  
                        p.crit = 1, method="D1")

res <- pool_intadj(pool_select, shrinkage_factor = 0.8449148)

res$int_adj
## [1] -4.795145
res$coef_shrink_pooled
##                  Pain             Radiation factor(Satisfaction)2 
##             0.7903410             0.4382261            -0.3618728 
## factor(Satisfaction)3     factor(Carrying)2     factor(Carrying)3 
##            -2.2332194             1.0607134             1.4255169

Step 7 - Externally validate the model

欠損データにより5回インプットされた外部データセットでモデルを検証する。 lbpmi_extvalがすでに作成済みの検証用データセットとして提供されている.

head(lbpmi_extval)

関数はmivalext_lr: External Validation of logistic prediction models in multiply imputed datasets

Developした式にcal.plotで係数を指定して使用する.

res_extval <- mivalext_lr(data.val = lbpmi_extval, nimp = 5, impvar = "Impnr",
                          formula = pool_lr$formula_final[[1]], 
                          lp.orig = c(res$int_adj, res$coef_shrink_pooled), 
                          #Numeric vector of the original coefficient values
                          cal.plot = TRUE, plot.method = "mean")
## 
## Pooled performance measures over m = 5 imputed external validation datasets
##             correctly estimated

res_extvalにはROCの統合結果も格納されているのでこれを示す

res_extval$ROC
##                     95% Low C-statistic 95% Up
## C-statistic (logit)  0.7382      0.8482 0.9172
res_extval$R2
## [1] 0.45587
res_extval$HLtest
##          D2           p         df1         df2 
##   0.7147875   0.6781956   8.0000000 149.5561959
res_extval$LP_pooled_ext
## intercept     slope 
##  -0.01812   0.88741

外部データセットでは、AUCは0.8482、R2乗は0.45587である。
Hosmer and Lemeshow検定のp値は0.67820で、これは適合が満足のいくものであることを意味している。
これは、多重補完されたデータセットで平均化されたキャリブレーション曲線が最適（破線）線に近いというキャリブレーションプロットでも確認できる。

傾き（回帰係数）は0.88741の値を持ち、1の値からわずかにずれています。これは開発データと外部データセットで係数の値が異なることを意味する。