EPV≧10 (c) - Logistic回帰に関するメモ

続いて，2019年の論文である． Google scholarでは引用数150． Crossrefでは83．Web of Scienceでは74．

van Smeden M, Moons KG, de Groot JA, et al. Sample size for binary logistic prediction models: Beyond events per variable criteria. Statistical Methods in Medical Research. 2019;28(8):2455-2474. doi:10.1177/0962280218784726

ありがたいことにフリーである．アブストラクトには以下のように書いてある．

EPVは予測性能の指標とは強い関係がなく、（二値）予測モデル開発研究の適切な基準とはならないことがわかった。我々は、予測変数の数、全サンプルサイズ、イベントの割合を考慮することで、サンプル外の予測性能を近似できることを示す。予測モデルの新しいサンプルサイズ基準の開発は、これら3つのパラメータに基づいて行われるべきであることを提案し、サンプルサイズ決定を改善するための提案を行う。

結論には以下のように書いてある．

サンプルサイズは、rMPSEやMAPEなどの意味のあるサンプル外予測性能の尺度に基づいて決定されるべきです。

ここで，

$\begin{aligned} {\rm MSPE} = \frac{1}{N}\sum_i (\pi_i - \hat{\pi}_i) ^2,\ {\rm MAPE} = \frac{1}{N}\sum_i |\pi_i - \hat{\pi}_i| \end{aligned}$

である． 5.3節には以下のように書いてある．

臨床予測モデルは新規患者の確率を推定するために使用されるので、rMSPEとMAPEは予測モデルを開発する際に直接的な関連性があります。

具体的にどうやってMSPEとMAPEを計算するかというと，

サンプル外のrMPSEおよびMAPEは、本論文で行ったように、シミュレーションによって近似することができます。シミュレーションコードは、GitHubで公開しています。また、rMPSEとMAPEは、我々のメタモデルの結果で近似することもできます（表4）。

とある．シミュレーションコードを読んでみたが，これからMSPEやMAPEを近似計算するのは無理ではなかろうか．真のモデルが分からないし，予測因子の分布も分からない．何種類も試せば近似がよくなるのだろうか．そうは思えないのだが．

メタモデルの結果で近似するのが現実的だろう．近似計算の具体例は論文の続きに書いてあり，後半で述べる． MSPEやMAPEを予測変数の数・サンプルサイズ・イベント割合で予測するわけである．シミュレーションではAUCや予測変数の相関係数などもハイパーパラメータとして変化させているが，これらはMSPEやMAPEの予測には効果がないという結果であった．

データが許容できる予測変数の数の求め方は以下のようになろう．

MSPEやMAPEの上限を定める．
サンプルサイズ・イベント割合はデータで定まるので，メタモデルからMSPEやMAPEの値が上限以下となる予測変数の数を求める．

なおいつものことだが，メタモデルの結果は検証していないし，補外するなと書いてある．

予測モデルで許容可能なMSPEやMAPEは合理的に決められるだろうか？小さければ小さいほどよいのだろうが.. また決められたとして近似するメタモデルが信じられるだろうか？コードを読んだ感想としては信じがたい．全体としての感想は，この論文も参考程度にしかできない，というものである．

Rのコードがgithubで手に入るのはよい．

github.com

実行してみたがエラーで止まった．以下のようにsim_LRMs.Rを変更したら動くようになった．修正が正しいかどうかは知らない．

  # scope <- attr(terms(working), "term.labels")
  scope <- terms(working)
      # working <- update(working, formula = newform, pl = FALSE,data=object$data)
      working <- update(working, formula = newform, pl = FALSE, data = object$model)
      working <- update(working, formula = newform, data = object$model)
    # working <- update(working, formula = newform,data=object$data)

シミュレーション用のデータを作成するのだが，AUCも事前に設定する．特定のAUCの値になるようにデータを生成できるとは知らなかったので，データ生成のコードを眺めてみたが，設定用のリストを見てもそれらしきパラメータは見当たらない．他の条件とパラメータから残り物を調べると，dgm.parがそれらしい．そこで，データを生成してみたが，AUCは設定値とはならなかった.. サンプルサイズが大きくなれば設定したAUCの値となるのかと思い試してみると，うまくいった！

condno <- 225#1,113
condID <- names(simlists)[condno]
simlist.validation <- simlists[[condID]]$validation
SIMx <- do.call(what=simlist.validation$datagen,args=simlist.validation$args)

design.mat <- as.matrix(cbind(1,SIMx))
dgm.par <- simlist.validation$dgm.par
p <- exp(design.mat%*%dgm.par)/(1+exp(design.mat%*%dgm.par))
p <- p %>% 
  as.vector()
y <- rbinom(length(p),1,p) %>% 
  as.vector()
pROC::roc(y ~ p)

$\begin{aligned} X \sim {\rm MN}(\mu, \Sigma)\\ p = \frac{1}{1+\exp(- ^t X \beta)}\\ Y \sim {\rm Ber}(p) \end{aligned}$

という状況下では，サンプルサイズを増加させるとAUCは一定値に収束するように思える． $\mu, \Sigma$ は事前に設定するので， $\beta$ によってAUCを決めるようだ． $\beta$ はグリッドサーチで求めるのかもしれない．サンプルサイズが小さなときはAUCは設定値と異なるのだが，それはよしと考えるのだろう，多分．

性能評価は三つの観点から行う．

識別能：△AUC（AUCの損失）
較正：CS（較正傾きの中央値），CIL（較正切片の中央値）
予測誤差：Brierスコア，MSPE（予測誤差2乗平均），MAPE（予測誤差絶対値の平均）

$\begin{aligned} {\rm CIL} = \frac{1}{N}\sum_i (Y_i - \hat{\pi}_i),\\ {\rm Brier} = \frac{1}{N}\sum_i (Y_i - \hat{\pi}_i) ^2,\ \end{aligned}$

5000回実行するのは時間がかかった．並列処理を行った．

library(parallel)
detectCores()
detectCores(logical = FALSE)

fctnames <- c(
                'simlists',
                'BackwardFR',
                'bias',
                'single.run',
                'random.sampling',
                'gen.MVNX',
                'gen.binY',
                'firth.lrm.plusbw',
                'tryCatch.W.E',
                'heuristicshrink.lrm',
                'penal.cv.lrm2',
                'stats.firth.lrm',
                'lp.firth.lrm',
                'predict.lrm',
                'c.stat2',
                'brier',
                'MSPE',
                'MAPE',
                'cond.error',
                'calibration',
                'coef.firth.lrm',
                'pseudo.Rsqrs',
                'LL',
                'stats.glmnet.lrm',
                'lp.penal.lrm',
                'coef.penal.lrm',
                'shrinkage.factor'
              )


cl <- makeCluster(7)
clusterEvalQ(cl, {
  require(MASS) # MVN function
  require(glmnet) # Ridge and Lasso LRMs
  require(logistf) # ML and Firth's LRMs
})
clusterExport(cl, c(fctnames, "condno", "condID", "simlist", "simlist.validation", "seeds"))
system.time(
  parLapply(cl, 1:5000, function(iter) {
    condno <- 1
    condID <- names(simlists)[condno]
    simlist <- simlists[[condID]]$development
    simlist.validation <- simlists[[condID]]$validation
    seeds <- sample(1:1e6,3)   # note: in our simulation the seeds were pre-defined
    single.run(simlist,simlist.validation,seeds,condno,condID,iter)
    }
  )
)
stopCluster(cl)

4.1節のボックスプロットを描いてみた． f:id:sintheta:20210727130153p:plain

図2のEPV=5のMAPEと同じはずだがかなり違う．図2だとMAPEの最大値が0.2以下のようだが，この図だと0.5以上となっている．またステップワイズを実行したほうがMAPEが小さくなっている．私が何か間違えているのだろうか．

Rコードは以下の通り．

library(tidyverse)
reslists <- list()
for (iter in 1:5000) {
  filename <- paste0("./OUT/main00001/main00001_", sprintf("%05.0f",iter), ".rds")
  res <- readRDS(file=filename)
  reslists <- append(reslists, list(res))
}

reclist <- function(x, FUN, index) {
  if (FUN(x)) (x[index]) 
  else if (is.list(x)) do.call("c", lapply(x, reclist, FUN, index))
}

plot_box <- function(reslists, index) {
  resv <- reslists %>% 
    reclist(function (x) {is.list(x) & (index %in% names(x))}, index) %>% 
    unlist()
  res <- tibble(name=names(resv), value=resv) %>% 
    filter(str_detect(name, "validate")) %>% 
    mutate(name = str_replace_all(name, "logistf|validate|hs|glmnet", "")) %>% 
    mutate(name = str_replace_all(name, index, ""))
  res %>% 
    ggplot(aes(x=name,y=value)) +
    geom_boxplot() +
    labs(y=index) +
    coord_flip() +
    theme_classic()
}
index <- "MAPE"
plot_box(reslists, index)

5.3節のrMPSEとMAPEの計算に戸惑ったので書いておこう．

N <- 400
p <- 8
EF <- 1/4
lin <- -0.59 -1.06*log(N) + 0.36*log(EF) + 0.94*log(p)
sqrt(exp(lin))
lin <- -0.93 -0.88*log(N) + 0.50*log(EF) + 0.49*log(p)
sqrt(exp(lin))
lin <- -0.48 -0.53*log(N) + 0.31*log(EF) + 0.48*log(p)
(exp(lin))
lin <- -0.61 -0.45*log(N) + 0.36*log(EF) + 0.26*log(p)
(exp(lin))

それから

$\begin{aligned} {\rm EPV} = \frac{N\times ({\rm event\ fraction})}{P}=400\cdot \frac{1}{4}\cdot\frac{1}{8} = 12.5 \end{aligned}$

なのだが本文に ${\rm EPV} = 20$ とあるのは間違いだろうか．

また図4の見方が一見して分からなかった．サンプルサイズに目盛りがないのは結果が相対的になるからだろうか． ${\rm EPV} = 10$ ならば，

$\begin{aligned} N= \frac{{\rm EPV}\times P}{{\rm event\ fraction}} = \frac{10 P}{{\rm event\ fraction}} \end{aligned}$

である．つまり双曲線である．他の曲線もラベルの値が一定となるようなEvents fractionに対する $N$ の値を図示しているようだ． MSPE, MAPE, Brierスコアはイベント割合が1/2から1/16になると必要なサンプルサイズは少なくなるが， ΔAUCとCSは多くなる．これを指して低いイベント割合では識別と較正が貧弱になると述べているようだ．