前の二つの記事では、20人のCD4データを人単位で選び直し、平均変化の標準誤差と信頼区間を求めました。ところが実際の解析では、一人分の行を選べば済む場面ばかりではありません。回帰分析なら行を選ぶのか残差を選ぶのか、多施設研究なら施設を選ぶのか施設内で選ぶのか、時系列なら1時点なのか連続する期間なのかという指定が加わります。
同じ「ブートストラップを行いました」という一文でも、この指定が違えば、近似している標本変動が別のものになります。名前ではなく、何を独立な抽出単位と考えたかが方法の中身です。
ここでは、回帰での行の抽出と残差の抽出、不均一分散を残すワイルド法、観測値ではなくモデルから生成するパラメトリック法、施設や個人をまとめて持ち込むクラスター抽出、そして連続する時点を残すブロック抽出を順に扱います。信頼区間の作り方そのものは、五つの作り方で端点はどこまで動くかで扱いました。
回帰分析で何を再抽出するか
回帰分析では、対象者の行を選び直す方法と、説明変数を固定してアウトカムを作り直す方法で、近似している標本変動が異なります。同じ「回帰係数のブートストラップ」でも、どちらを使ったかを方法に記す必要があります。
この区別について、Freedmanは、回帰に対するブートストラップで、仮定したデータ生成の仕組みを再抽出が反映する必要性を論じています。分布を指定しない方法という印象だけで、行の抽出と残差の抽出を入れ替えると、近似している標本変動が変わります。
対象者の行を選ぶ方法
pairs bootstrap、case bootstrapでは、アウトカムと説明変数を持つ一人分の行を、そのまま復元抽出します。たとえば、年齢、治療群、アウトカムを別々に選ぶのではなく、その組を同じ対象者から取ります。
そのため、対象者ごとの共変量とアウトカムの組を独立な標本として扱える観測研究では、この方式が候補になります。そして、各再抽出標本で回帰モデルを当てはめ直し、目的の係数や予測値を取り出します。
ここで、元標本で当てはめた回帰係数を固定し、抽出された行に対する予測だけを繰り返しても、係数の推定の揺れは反映されません。どこまでを推定し直すかが、再抽出単位と同じくらい重要になります。
残差を選ぶ方法
説明変数を固定し、回帰モデルの平均構造が正しく、誤差が独立で同じ分布を持つと仮定する場合には、残差ブートストラップが候補になります。
説明変数ベクトルを$\mathbf x_i$、推定係数を$\widehat{\boldsymbol\beta}$、中心化などを行った残差から抽出する値を$e_i^*$とすると、作り直すアウトカムは次の形です。
$$Y_i^*=\mathbf x_i^\top\widehat{\boldsymbol\beta}+e_i^* \tag{1}$$
式(1)は、元モデルが予測する平均に、残差から選んだばらつきを加えています。$\top$は転置で、説明変数の値は固定したまま、アウトカム側の変動を作っています。
ただし、説明変数によって残差分散が異なる場合に、全員の残差を一つに集めて配り直すと、その違いを壊してしまいます。また、生の残差はモデル当てはめによる収縮を含むため、中心化や自由度、レバレッジへの対応が必要になる場合があります。
なお、回帰の平均構造や残差の状態は、回帰分析の仮定や不均一分散性と合わせて確認すると、どの揺れを再現したいかを説明しやすくなります。
不均一分散とワイルド法
ワイルド・ブートストラップでは、各観測の残差にランダムな重みを掛け、残差の大きさと説明変数の関係を残しながらアウトカムを作る方法を使います。
調整した残差を$\widetilde e_i$、平均0、分散1の重みを$V_i$とすると、基本的な形は次のようになります。
$$Y_i^*= \mathbf x_i^\top\widehat{\boldsymbol\beta} +\widetilde e_iV_i \tag{2}$$
式(2)は、残差を別の対象者へ配り直す代わりに、各人の残差の大きさを利用して、符号や大きさをランダムに動かしています。重みの分布や残差調整は、使う方法に合わせて指定します。
このような残差に基づく再抽出は、不均一分散への対応を含む回帰の再標本化研究につながっています。ただし、平均0、分散1という条件は重みに課される最低限の要請で、それだけで任意のモデルに使えるかが決まるわけではありません。
そのうえ、施設内に相関があるなら、対象者ごとに独立の重みを与える方式とは別に、施設単位の重みを使う方法などが必要です。不均一分散と施設内依存は異なる問題なので、ワイルド法という名称だけでは、その両方を扱ったことにはなりません。
不均一分散のあるデータで確かめる
説明変数が大きいほど残差のばらつきも大きくなるデータを乱数で生成し、行を選ぶ方法と残差を選ぶ方法で傾きの分布を比べます。200例について、投与量が1から5の範囲に一様に分布し、残差の標準偏差が投与量の4乗に比例する設定です。
set.seed(42)
dose <- runif(200, min = 1, max = 5)
response <- 2 + 0.5 * dose + rnorm(200, sd = 0.02 * dose^4)
fitted_model <- lm(response ~ dose)
adjusted <- residuals(fitted_model) / sqrt(1 - hatvalues(fitted_model))
slope_case <- replicate(2000, {
rows <- sample.int(200, 200, replace = TRUE)
coef(lm(response[rows] ~ dose[rows]))[2]
})
slope_residual <- replicate(2000, {
drawn <- fitted(fitted_model) + sample(adjusted, 200, replace = TRUE)
coef(lm(drawn ~ dose))[2]
})
round(c(case = sd(slope_case), residual = sd(slope_residual)), 4)#> case residual
#> 0.3533 0.2928行を選んだ場合の標準誤差は0.3533、残差を選んだ場合は0.2928でした。後者は、通常のlm関数が返す標準誤差0.2926とほぼ同じ値です。残差を全員分まとめて配り直した時点で、残差の大きさと投与量の関係が消え、等分散を仮定した式と同じ答えに戻っています。
ここで読み違えやすいのは、狭い方が精度の高い推定だという読み方です。この設定では、残差の分散が投与量によって変わることを壊した結果として幅が縮んでいるので、狭さは情報の増加を表していません。
パラメトリック法との違い
ここまでは観測値そのものを選び直してきましたが、パラメトリック・ブートストラップでは、推定した分布や回帰モデルからデータを生成します。たとえば、推定した正規分布から新しい値を生成したり、ロジスティック回帰の推定確率から二値アウトカムを生成したりします。
ここで生成する値は、実際に追加収集した測定値ではありません。仮定したモデルの下で研究を繰り返した場合の推定の揺れを調べるためのデータです。また、各生成データへモデルを当てはめ直すことで、元の推定手順の変動を近似します。
観測範囲を超える値
ノンパラメトリックな行の復元抽出では、個々の測定値は元データにある値から選ばれます。一方、パラメトリック法では、モデルが許す範囲で元標本の最小値や最大値を超える値も生成できます。
この違いから読み取れるのは、情報の多さではなく、仮定の置き方の差です。観測していない範囲の値を生成できるのは、その範囲について分布の形を仮定しているためです。そして、裾の形を誤ると、整った再抽出分布を得られても、区間の被覆率が合わなくなる可能性があります。
そして、ノンパラメトリックという言葉も、「仮定がない」とは読みません。独立な抽出単位、母集団との対応、推定量の近似可能性などの条件は残ります。また、標本から観測できなかった裾の情報は、単純な復元抽出だけでは追加できません。
推定するモデルと検定用のモデル
パラメトリック・ブートストラップで信頼区間を作る場合と、帰無仮説の下で検定統計量を作る場合では、生成に使うモデルの制約が異なることがあります。
たとえば、治療効果0を検定するために制約付きモデルからデータを生成する方法は、推定した治療効果を残したままデータを生成する方法とは目的が異なります。何を再現する分布が必要なのかを決めてから、生成モデルを指定します。
多施設データと反復測定
1施設に100人ずつ、10施設から計1,000人を集めても、同じ施設の対象者に共通する変動があれば、1,000人全員を独立として抽出し直す方法では、その構造を再現できません。また、同じ人の10回の測定を10人として選び直す操作でも、個人内の関連を壊してしまいます。
そこで、クラスター・ブートストラップでは、独立な抽出単位と考える施設や個人を選び、その中の観測をまとめて持ち込みます。FieldとWelshは、階層データに対する複数の再抽出方式について、成立条件が仮定するモデルに依存することを論じています。
施設を選ぶ場合と施設内で選ぶ場合
施設を母集団からの標本と考え、選んだ施設の全対象者を観測している設定なら、施設ごとに抽出する方式が候補になります。施設が10なら各回に10施設分を復元抽出し、同じ施設が複数回現れることもあります。
ただし、施設ごとの人数が異なると、選ぶ施設の数は10でも、その標本に含まれる対象者の総数は元の1,000人から変わる場合があります。対象者数が変わったからという理由で、後から行を切り詰める操作はしません。
反対に、研究対象の施設を固定し、その施設内の対象者抽出だけを考えるなら、施設内での再抽出が候補になります。また、施設と施設内の対象者の両方が標本化される研究では、二段階の抽出を反映する方法を考えます。
もっとも、どの階層も自動的に抽出し直せばよい、という選び方にはしません。施設を固定した推論なのか、別の施設へ一般化したい推論なのかで、再現すべき研究が異なるためです。また、施設の平均を等しく重み付けする対象と、対象者全体の平均を知りたい対象も区別します。
重複施設の識別と層別指定
同じ施設が2回抽出された時、二つのコピーを一つの施設IDのまま混合モデルへ渡すと、意図した別々のクラスターとして扱われない場合があります。クラスターのコピーごとに新しい識別子を付け、個人のコピーについても対応するIDを管理します。
また、boot関数のstrataへ施設IDを指定するだけでは、施設そのものを抽出する方式にはなりません。通常は各層の中で行を抽出する指定なので、施設内の再抽出とクラスター抽出を混同しないようにします。
なお、再抽出単位と解析モデルの関係は、frailtyモデルと多施設データとも共通する問題です。ただし、混合モデルやfrailtyモデルを使ったという情報だけで、再抽出方式が一つに決まるわけではありません。
少数クラスターの情報量
1,000人がいても、施設間の変動を10施設から学習する状況では、その変動に関する情報量は対象者数だけで決まりません。10施設を10,000回選び直しても、観測した施設の種類は10のままです。
そして、30施設や50施設なら一律に信頼できるという境界もありません。施設サイズの不均衡、施設レベルの説明変数、治療群への施設の偏りなどによって有限標本での性能が変わるため、少数クラスターでは小標本調整やモデルに基づく推論も含めて比較します。
また、施設を一つずつ除いた時に主要な推定値が大きく変わるなら、その事実を反復回数の多さで覆い隠さず、どの施設の情報に結果が依存しているかを確認します。
時系列とブロック抽出
日ごとの入院件数や、同じ装置から連続して得た測定では、隣り合う時点が関連する場合があります。観測時点を独立に選び直すと、その関連を消してしまい、平均や回帰係数の標準誤差が元の系列に対応しなくなることがあります。
そこで、ブロック・ブートストラップでは、連続する複数時点をひとまとまりとして抽出し、ブロック内の関連を残します。Künschの1989年の論文は、定常系列に対するこのような再標本化の理論を扱っています。
ただし、ブロックを長くすると一つのブロックに含む依存を残しやすい反面、元データから作れる実質的に異なるブロックは少なくなります。短いブロックでは長く続く依存を捨てる可能性があるため、ブロック長は計算の都合だけで決めません。
実装としては、移動ブロックや長さをランダムにする方式などがあり、bootパッケージには時系列用のtsboot関数もあります。
また、流行、測定法の変更、制度変更などで系列の分布が時期によって変わる場合には、単純な定常性の仮定が合わない可能性があります。その場合も、ブロック抽出という名称だけで非定常性まで処理したことにはなりません。
単位の指定を方法に書く
ここまで並べた方式は、どれも計算量の大小で選ぶものではありません。行、残差、施設、ブロックのどれを独立に得られた単位と考えるかは、研究の集め方から決まります。
私なら、コードを書く前に、この研究をもう一度やり直すとしたら何が入れ替わるのかを先に決めます。入れ替わるものが施設なら施設を選び、対象者なら対象者を選ぶという対応が付けば、関数の引数は後から決まります。
次の記事では、こうした単位を正しく指定しても近似が成り立たない場面を扱います。分布の上限を標本最大値で推定する場合や、イベントが観測されなかった場合には、回数を増やしても埋まらない制約が残ります。
参考文献
Freedman DA. Bootstrapping regression models. The Annals of Statistics. 1981;9(6):1218-1228.
Wu CFJ. Jackknife, bootstrap and other resampling methods in regression analysis. The Annals of Statistics. 1986;14(4):1261-1295.
Field CA, Welsh AH. Bootstrapping clustered data. Journal of the Royal Statistical Society: Series B. 2007;69(3):369-390.
Künsch HR. The jackknife and the bootstrap for general stationary observations. The Annals of Statistics. 1989;17(3):1217-1241.
Loy A, Korobova J. Bootstrapping clustered data in R using lmeresampler. The R Journal. 2022;14(4):103-120.
