重回帰のときに必要なデータの数はパラメータの数の10倍です

MikuHatsune2016-09-28

読んだ。
J Clin Epidemiol. 1995 Dec;48(12):1503-10.
パラメータ数が $m$ のデータを標本数 $n$ 個観測して重回帰なりなんらかの回帰分析を行うのだが、「標本数はいくつあったらいいんですか?」と聞かれることが多々ある。そういうとき、パラメータの10倍(10 $m$ )あったらいいんじゃないっすかね〜(適当と答えることが多いのだが、パラメータ数に対する標本数の割合(EPV, event per variables)が10 ならいいんじゃないっすかね(適当と言っている論文。
　
論文中では673人の患者に対して7つのパラメータが測定されていて、252件の死亡についてCox 回帰を行うが、データセットをすべて使えば252/7 = 36 EPV のところを、2, 5, 10, 15, 20 とEPV の割合を変えてデータセットをリサンプリングしてCox 回帰を繰り返し、パラメータについて推定される係数の分布、つまりEPV が不十分だと推定値がめちゃくちゃになるのではないかと思って解析している。
　
結論からいうとEPV=10 くらいがいいんじゃないかと言っている。これをR でやってみよう。
データは重症患者に右心カテ(RHC) をするのがいいのか悪いのかというのを傾向スコアでやったデータセットrhc というのをパクってくる。パラメータは連続量と男女の01 に限定し、死亡の2値を予測するロジスティック回帰を想定すると、23パラメータ、5735症例になる。死亡が3722、生存が2013なので、EPV は87 くらいまでとれるが、シミュレーション上は50くらいまでとする。つまり、生存と死亡がそれぞれ1130ずつのデータセットを作る(死亡生存の比率は残さず、簡単に1:1 としておく)。
相関や多重共線性をみるためのcoplot は割愛するが、相関は総じて0.3 以下で、パラメータ間の強い相関はなかった。
　
Fig.1 に相当する、EPV を変えてロジスティック回帰(論文ではCox回帰)したときの、推定されたパラメータの係数の値の分布。中央の破線はすべてのデータを使った時の推定値である。EPV が小さい、例えば3-5くらいのときは、中央値はすべてのデータを使ったときの推定値とほぼ一致するが、ばらつきが正負両方向に大きい。
EPV を大きくしていけば、ばらつきはちいさくなる。
今回はage のパラメータだけプロットしたが、ほとんどこの傾向になる。一部、二峰性になるものもあったが。
論文ではCHF (慢性心不全)についてのみプロットしてある。

　
Fig.2 に相当する、すべてのデータセットを使って推定した真の値(推定という時点で真ではないが)から、EPV を変化させた時の推定値が相対的にどれくらい差がでてしまうかのプロット。論文ではEPV が10を超えるあたりから、相対差が±10% にはいるので、これくらいならprudent じゃないか、と言っている。
prudent と言っているのは、本当なら死ぬほど症例集めたいけど、現実的には難しいのでこれくらい、という程度である。しかしこれでも多くの研究では難しいことがほとんどであるが…
図では収束性がなかった数個のパラメータを除いてプロットしている。

　
というのが論文ではシミュレーションを用いて説明されている。
しかし、モデルの作成に用いた7つのパラメータの選び方が、死亡にp<0.1 で関連のあったものを選んでいるらしい。このパラメータセットを使ったモデルをゴールデンスタンダードとしているため、実際の答えは神のみぞ知る、という感じ。EPV をいじるのに、サンプル数しかいじっていないので、パラメータ数をいじった場合はどうなるか。
　
EPV については、こちらでも数行だけ言及されている。

数学いらずの医科統計学第2版

作者: 津崎晃一
出版社/メーカー: メディカルサイエンスインターナショナル
発売日: 2011/03/24
メディア: 単行本
購入: 1人クリック: 37回
この商品を含むブログ (2件) を見る

何も考えないなら10 $m$ , 変数選択をするならば20 $m$ , とかいろいろ書いてあるが、データはあればあるだけいいのは当然である。かといって、10 $m$ を満たせないときの解析は無駄かというと、そうでもない。というのは、モデルがそもそも正しいかという問題があるし、測定されているデータがものすごい精度で、ばらつきがほぼない、死亡と生存の2群をほぼ正確にわけている、とかなら少なくてもいいと思う。
本気で正統派な回帰分析をするなら、あるパラメータがどの程度死亡に寄与しそうか、効果量と検出力を決めてやれば必要なサンプル数は出てくるはずだが、回帰分析を「意味のありそうな、有意になるパラメータを探す」という探索的な解析と、「観測できるパラメータたちを使って、今後の患者の予後予測をする」というモデル構築的な解析とではたぶん本当のところは違うような気がするのだが、いろいろ考えると難しい。
　
当然ながら、EPV = 10 じゃなくてもいいんじゃないか、という話はある。
Am J Epidemiol. 2007 Mar 15;165(6):710-8.

library(abind)

# ネットが繋がっていれば
rhc <- read.csv("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.csv")
rhc0 <- rhc[, c(9, 22:27, 29:44, 46)] # NA とか変なカテゴリパラメタなどを削除したとして
EPV <- c(2:10, 12, 15, 20, 25, 30, 50)

iter <- 1000 # シミュレーション回数
g <- glm(death ~ ., data=rhc0, family=binomial(link="logit")) # すべての症例を使ったフルモデル
res <- fitted.values(g) > 0.5

res <- mapply(function(j){
  n <- (ncol(rhc0)-1) * j # EPV に対応した症例数
  c0 <- mapply(function(i){
    idx <- mapply(function(z) sample(seq(z), size=n), table(rhc0$death), SIMPLIFY=FALSE)
    rhc1 <- do.call(rbind.data.frame, mapply(function(z1, z2) z1[z2,], split(rhc0, rhc0$death), idx, SIMPLIFY=FALSE))
    g0 <- glm(death ~ ., data=rhc1, family=binomial(link="logit"))$coefficients
    return(g0)
  }, seq(iter))
  return(c0)
}, EPV, SIMPLIFY=FALSE)

res <- do.call(abind, c(res, along=3))


cols <- rainbow(length(EPV))
par(mfrow=c(2, 2))
for(i in 1:nrow(res)){
  xl <- range(res[i, , -1])
  d0 <- apply(res[i,,], 2, density)
  yl <- c(0, max(mapply(function(z) z$y, d0)[,-1]))
  plot(0, type="n", xlim=xl, ylim=yl,
       xlab="coefficients", ylab="Density", main=paste("Parameter:", names(g$coefficients)[i]))
  for(j in 1:length(EPV)){
    lines(d0[[j]]$x, d0[[j]]$y, lwd=2, col=cols[j])
  }
  abline(v=g$coefficients[i], lty=3, col="blue")
  legend("topright", legend=paste0("EPV = ", EPV), col=cols, lty=1, lwd=3)
}


# relative bias
hoge <- mapply(function(i) apply((res[i,,]-g$coefficients[i])/g$coefficients[i], 2, median), 1:nrow(res))
matplot(EPV, hoge[,c(1:3,5:7,9,12:14,19,24)], type="o", ylim=c(-1, 1), pch=16, lty=1, lwd=3, xlab="Event per variable", ylab="Relative bias")