2022-02-24

繰り返し計算による回帰係数の決定

ロジスティック回帰の係数は最尤法により決まるわけだが，繰り返し計算によって収束し値が求まる．

2021-07-28

オーバーフィットとはなんぞや

回帰モデルに投入する変数が多いとオーバーフィットする，という．ではオーバーフィット（過剰適合）とは何を指すのか？

よく例に出されるのはn個の点とそれを通るn-1次式である．全ての点を正確に通るが補外はおろか補間もできそうにない． f:id:sintheta:20210728160046p:plain

n <- 10
x <- 1:n
y <- rnorm(n)
df <- tibble(y, x)
model <- lm(y ~  poly(x,n-1), df)
newx <- seq(1,n,length.out = 10*n)
predy <- predict(model, newdata = tibble(x=newx))
plot(newx, predy, type = "l")
points(x, y, col="red")

これと説明変数が多いときの過剰適合は同じだろうか？

以下の論文は他の何かを検索していたときに，Q&Aサイトで紹介されていた．引用数はPubMedだと396．この人はHarrellの弟子なのでしょうか．

Babyak MA. What you see may not be what you get: a brief, nontechnical introduction to overfitting in regression-type models. Psychosom Med. 2004 May-Jun;66(3):411-21. doi: 10.1097/01.psy.0000127692.23278.a9. PMID: 15184705.

アブストラクトから．

オーバーフィッティングの概念は，利用可能なデータから多くのことを求めるという観点から提示されます。あるデータセットに一定の数の観測値がある場合、許容できる不確実性の度合いで導き出すことのできるモデルの複雑さには上限があります。複雑さは、データ分析のどの段階においても、同じデータセットに対して費やされる自由度の数（相互作用や非線形項などの複雑な項を含む予測変数の数）の関数として生じます。

イントロダクションから．こういうときの optimisticは楽観的という訳でいいのだろうか．

それは、回帰型モデルにおいて、サンプルの特異な特性を利用する問題（オーバーフィットと呼ばれる）です。オーバーフィッティングは、過度に楽観的なモデル結果をもたらします。オーバーフィッティングされたモデルに現れる「所見」は、実際には母集団に存在しないため、再現されない。

本文4ページから．

利用可能なサンプルサイズに対して自由度が高すぎる回帰モデルでは，繰り返しのサンプルで重みが大きく変動する傾向があります．サンプル間の変動が大きいと、あるサンプルで回帰重みの一部が非常に大きくなる可能性が高くなり、過剰に楽観的な適合になってしまいます。

そうですね．

少し違った見方をすると，統計学者は，10個の予測変数と20個のオブザベーションを持つ回帰を推定することは，ある意味で，それぞれがN = 2の標本サイズを持つ10個の別々の1-予測変数の回帰を推定することと同等であるとよく指摘する．回帰係数の不安定さの問題は、どの回帰重みに注意を払うかを先験的に選択するのではなく、それらのより大きな回帰重みに注目してしまう（さらに悪いことに、最終的なモデルのために、より大きな、またはより有意な重みを選んでしまう）我々の自然な傾向によってさらに悪化します。

なるほど．自由度が高いほどサンプルの変動によって大きな値をもつ回帰係数が生じてしまう．そして回帰係数の大きさ・重要さにだけ注目するので，サンプルの変動の影響で大きな回帰係数をもってしまった変数を選んでしまう，ということか．

2021-07-19

ロスマンと多重比較

ロスマンは多重比較を考慮しないそうだ．

Rothman KJ. No adjustments are needed for multiple comparisons. Epidemiology. 1990 Jan;1(1):43-6. PMID: 2081237.

引用数は4000以上．並の疫学者はひれ伏すであろう．ロスマンが言っているから多重比較は考慮しません，と言ってみたい．

アブストラクトでは以下のように言っている．

多重比較の調整を日常的に行うことを推奨する理論的根拠は、「偶然」が観察された現象の第一次説明として機能するという「普遍的な帰無仮説」です。この仮説は、自然は規則的な法則に従っており、観察によって研究することができるという、経験的研究の大前提を覆すものである。多重比較の調整を行わない方針の方が、評価対象のデータが乱数ではなく、実際の自然界の観察結果である場合には、解釈の誤りが少なくて済むので望ましい。

2021-07-18

『完全無欠の賭け』

『感染の法則』が面白かったので同じ著者の本を借りた．やはり軽快な文章で面白い．この本の訳では一人称は私．脚注番号があるのに本の中に脚注がない．巻末に「本書の原注は，下記のURLよりPDFファイルをダウンロードしてご覧ください．」とあった．今どきはこんな分割方法もあるのですね．薄い本の方が好まれるということだろうか．

2021-07-15

EPV≧10 (d)

次は，

Riley RD, Snell KI, Ensor J, et al. Minimum sample size for developing a multivariable prediction model: PART II - binary and time-to-event outcomes [published correction appears in Stat Med. 2019 Dec 30;38(30):5672]. Stat Med. 2019;38(7):1276-1296. doi:10.1002/sim.7992

である．

2021-07-15

EPV≧10 (c)

続いて，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は多くなる．これを指して低いイベント割合では識別と較正が貧弱になると述べているようだ．

2021-07-14

EPV≧10 (b)

ロジスティック回帰においてEPV≧10の起源は以下の論文のようである．まだ25年前だから新しいと感じる．ネット上にPDFが転がっているが真っ当なものかどうかはあやしい．

Peduzzi P, Concato J, Kemper E, Holford TR, Feinstein AR. A simulation study of the number of events per variable in logistic regression analysis. J Clin Epidemiol. 1996 Dec;49(12):1373-9. doi: 10.1016/s0895-4356(96)00236-3. PMID: 8970487.

現時点ではGoogle scholarで引用6330，ELSEVIERのCiting articlesが4124と凄まじい．（しかし差が2000件以上あるのはどういうわけだろうか？）使用データはネットから入手できなさそうである．

この研究はシミュレーション研究である． EPVとモデル性能の関係を調べるために，ある実データからEPVを変えたデータを作ってモデルを作成しその性能を評価する．

サンプルサイズ673，イベント数252の実データを元データに使う．予測変数は7つで，EPV=36である．元データから作成したモデルを参照モデルとする．

元データからEPV=2,5,10,15,20,25となるデータを作成する．サンプルサイズはオリジナルと同じ673とする．例えばEPV=2の場合は，2*7=14の生存サンプルを参照モデルから計算した生存確率をもとに抽出する．このときサンプルが実際は死亡であっても生存として扱うようだ．また同じサンプルを繰り返し選んでもよい．死亡の659=673-14サンプルの選び方も同様である．各EPVの値に対してそれぞれ500セットのデータを作成し，それぞれのデータに対してモデルを作成する．ただし収束しなかったモデルは性能評価に利用しない．このようにしてデータを作成すると，データの真のモデルの回帰係数が参照モデルと同じ値となる．（論文中に証明されている．）このことを使ってEPVの値に対するモデルの性能を評価する．

読んだ感想としては，一つの特定のデータを使っているだけので，これから何かを主張するには根拠が弱い気がする．またモデルが収束しなかった場合にそのモデルの指標を使わないのもまずい気がする．著者もlimitationsで，

単一の「実在」データ
調査した変数はすべて離散的で、
有病率や転帰との関連性（すべて正）も中程度の範囲しかなかった。
連続尺度や変数間の相互作用は含まれていなかった。
その他（略）

と述べている．アブストラクトの結論では

EPV値が10以上の場合、大きな問題は発生しませんでした。しかし，EPV値が10未満の場合は（略）。ロジスティックモデルの有効性には，他の要因（イベントの総数やサンプルサイズなど）も影響すると思われますが，今回の結果は，EPVが低いと大きな問題が生じる可能性があることを示しています。

と慎重に述べているので，EPV<10はまずそうですよ，というのが結論となりそうだ． EPV≧10なら大丈夫とは言っていない．

EPV≧10は経験則として利用しやすいのが広まった理由だろうか．この論文が発表されるまでモデルに投入する変数の数に制限がなかったのなら大きな進歩であったろう．とはいえ根拠は薄弱である．よく使われるようになった現在，意味は失われ形式だけが残るのである．

Logistic回帰に関するメモ

Logistic回帰についてまとめます．