Breidt, F.J., Opsomer, J.D. (2017) Model-assisted survey estimation with modern prediction techniques. Statistical Science, 32(2), 190-205.
のメモ、最終回。いよいよ面白い話になってくる。GREG推定でいうところの「作業モデル」を機械学習で組んじゃうのだ。
8. 統計的学習テクニックに基づくモデル・アシステッド手法
8.1 モデル・カリブレーション・アプローチ
線形モデルの外側には、観察データの非線形結合による予測手続きがたくさんあります。それらを用いると、GREG型推定量のカリブレーションとウェイティングの特性はますます複雑になります。
Wu & Sitter (2001)は、モデル・アシステッド推定量に非線形予測を組み込む単純で効率的な方法を提案しています。作業モデルを$$ y_k = \beta z_k + \epsilon_k $$ $$ \{\epsilon_k\} \ \mathrm{uncorrelated} \ (0, \sigma^2) $$ とし、任意の予測子 \(\hat{m}(\mathbf{x}_k)\)を共変量\(z_k\)として用いるのです。
[えええ？ 頭が混乱してきた。そもそも\(m_N(\mathbf{x}_k)\)というのは作業モデルの第1項であっで、GREG推定なら本来は\(m_N(\mathbf{x}_k) = \mathbf{x}_k^\top B_N\)で、でも標本しか観察できないから代わりに\(\hat{m}(\mathbf{x}_k) = \mathbf{x}^\top_k \hat{B} \)を使う、という話だよね？ いっぽうここでは、作業モデルの第1項を\(\beta z_k = \beta m_N(\mathbf{x}_k) \) とし、代わりに\(\hat{B} \hat{m}(\mathbf{x}_k) \)を使う、ということだろうか。たぶん\(z_k\)は\(\mathbf{x}_k\)の非線形関数でスカラーなのだろう。もはや\(m\)は\(m\)でない別の記号を使って書いたほうがよくない？]
[なにがなんだかわかんなくなりつつあるので、ここからしばらく逐語訳…]
\(\beta\)の推定は、母集団レベルでは$$ B_N = \sum_{k \in U} z_k y_k (\sum_{k_U} z^2_k)^{-1} $$ となり、標本レベルでは $$ \hat{B} = \frac{\sum_{k \in s} z_k y_k / \pi_k}{\sum_{k \in s} z^2_k / \pi_k} = \frac{\hat{m}(\mathbf{x}_k) y_k / \pi_k}{\sum_{k \in s} \hat{m}(\mathbf{x}_k)^2 / \pi_k} $$ となります。これを差分推定量にプラグインして $$ DIFF(y, \hat{B}\hat{m}(\mathbf{x}_k)) = \sum_{k \in U} \hat{B} \hat{m} (\mathbf{x}_k) + \sum_{k \in s} \frac{y_k – \hat{B} \hat{m}(\mathbf{x}_k)}{\pi_k} $$ となります。Wu & Sitter はこれを「モデル・カリブレーション推定量」と呼んでいます。
[虚心に数式をコーディングしていてだんだん腑に落ちてきた。もはや \(\hat{m}_N(\mathbf{x}_k)\)は母集団での\(y_k\)の予測というより、母集団で\(y_k\)と比例するはずのなにか、なのだ。それがなんなのかはきっと別の問題だ]
のちにWu(2003)は、抽出デザインにおけるある種のregularity条件の下で、このモデル・カリブレーション推定量が最適であることを示してます。ここで最適という意味は、カリブレーション推定量のクラスを通じて、漸近デザイン分散の超母集団モデルに基づく期待値が最小であるということです。
8.2 カーネル法と局所回帰
カーネル法では、作業モデルが局所的に単純であり(たとえば定数とか線形であり)、大域的にはスムーズであると仮定して、カーネル・ウェイティング関数によって決められた点の周りだけについて局所回帰関数を推定します。本章では局所多項回帰とそのモデル・アシステッド推定への拡張に焦点を当てます。
局所多項回帰
Breidt & Opsomer (2000)は、作業モデルの\(\mu(\cdot)\)をスカラー\(x\)の平滑化関数にして局所多項回帰で推定するというのを考えています。\(x_k\)における平滑化関数を\(q\)次多項式で局所的に推定するというアイデアです。母集団レベルでのあてはめはWLS回帰で行います。ウェイトは\(x_k\)を中心にしたカーネル関数で決めます。$$ \mathbf{y}^\top_U = [y_1, y_2, \ldots, y_N] $$ $$ \mathbf{X}_{Uk} = [1 \ x_j – x_k \ \cdots \ (x_j-x_k)^q]_{j \in U} $$ $$ \mathbf{W}_{Uk} = \mathrm{diag} \left( \frac{1}{h} K \left(\frac{x_j – x_k}{h} \right) \right)_{j \in U} $$ として $$ m_N(x_k) = (1, 0, \ldots, 0) \cdot (\mathbf{X}^\top_{Uk} \mathbf{W}_{Uk} \mathbf{X}_{Uk})^{-1} \mathbf{X}^\top_{Uk} \mathbf{W}_{Uk} \mathbf{y}_U $$ となります。
これのデザイン・ウェイテッド版を考えます。$$ \mathbf{y}^\top_s = [y_j]_{j \in s} $$ $$ \mathbf{X}_{sk} = [1 \ x_j – x_k \ \cdots \ (x_j-x_k)^q]_{j \in s} $$ $$ \mathbf{W}_{sk} = \mathrm{diag} \left( \frac{1}{\pi_j h} K \left(\frac{x_j – x_k}{h} \right) \right)_{j \in s} $$ として $$ \hat{m}_{LPR}(x_k) = (1, 0, \ldots, 0) \cdot (\mathbf{X}^\top_{sk} \mathbf{W}_{sk} \mathbf{X}_{sk})^{-1} \mathbf{X}^\top_{sk} \mathbf{W}_{sk} \mathbf{y}_s $$ これを差分推定量にプラグインした\(DIFF(y, \hat{m}_{LPR}(x_k))\)をLPR(局所多項回帰)推定量と呼ぶことにします。
その挙動はどうなるかというと[…中略…]、というわけで、差分推定量と同様の漸近分散とその推定量を持ちます。また、\(y\)と独立なウェイトによるウェイティングの形で書くこともできます。GREG推定量とのちがいは、すべての\(k \in U\)について\(x_k\)が必要になるという点です。
LPRの拡張や手法比較研究として以下があります…[略]
加法モデル
Breidt, Opsomer, Johnson, & Ranalli (2007)は、作業モデルの\(\mu(\mathbf{x}_k)\)がセミパラメトリック加法モデルであるモデル・アシステッド推定量を考えました。$$ \mu(\mathbf{x}_k) = \mu_1(x_{1k}) + \cdots + \mu_q(x_{qk}) + \mathbf{x}^\top_k \beta$$ 彼らはこのセミパラメトリック加法モデルを効率的に推定するデザイン・ウェイテッドのbackfittingアルゴリズムを提案しています。[なにいってんだかわからん]
[…以下、簡潔な研究レビュー。メモは省略する]
8.3 スプライン
カーネル法の欠点は、複数の共変量(特にカテゴリカル共変量と連続的共変量の組み合わせ)があったときの難しさ、そしてスパースなデータの場合の計算の大変さです。これらの困難さは、たくさんのスプライン(ないしなんらかの基底関数)を導入し、モデルの複雑をコントロールするための選択や正規化・罰則を導入することによってかなり軽減されます。
罰則つきスプライン
Breidt, Claeskens, & Opsomer(2005)は以下を提案しています。作業モデルはスカラー\(x\)の平滑化関数とするのですが、\(x\)の範囲を決めておいて、それを既知のノット\( \{\kappa_j\}_{j=1}^K \)で分割しておいて、線形混合モデル$$ y_k = [1, x_k, \ldots, x_k^q] \beta + [(x_k – \kappa_1)^q + (x_k – \kappa_2)^q + \cdots + (x_k + \kappa_K)^q] \mathbf{b} + \epsilon_k $$ $$ E\left( \left[ \begin{array}{c} \mathbf{b} \\ \mathbf{\epsilon}_U \end{array} \right] \right) = \mathbf{0} $$ $$ Var \left( \left[ \begin{array} \mathbf{b} \\ \mathbf{\epsilon}_U \end{array} \right] \right) = \sigma^2 \left[ \begin{array}{cc} \lambda^{-1} \mathbf{I}_k & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_n \end{array} \right] $$
この推定量の挙動も差分推定量に似ていて…[略]
[以下、いろいろ書いてあるけどめんどくさいのでパス]
回帰スプライン
Goga(2004, 2005)は回帰スプラインに基づくモデル・アシステッド推定量について研究しています。罰則のない回帰スプラインを使い、ドメインを\(K\)個のノットで分けてBスプライン基底を構築するのですが、\(K \rightarrow \infty\)として十分に密にします。[なんだか怖ろしい話だなあ。以下略]
8.4 ニューラル・ネットワークと関連手法
今度は、オリジナルの\(\mathbf{x}_k\)の線形結合として新しい共変量を導出し、作業モデルにおいては平均反応をそれらの共変量の非線形関数とするタイプの予測手法に注目しましょう。HasieらのESL本の11章をみてください[ニューラル・ネットワークの章だ]。
ニューラル・ネットワーク
Montanari & Ranalli (2005)は、作業モデルがレイヤーをスキップする結合を持つフィードフォワード・ニューラル・ネットワークであるモデル・アシステッド推定量を提案しています。$$ \mu(\mathbf{x}_k) = \mathbf{x}^\top_k \beta + \sum_{j=1}^M \alpha_j a(\gamma^\top_j \mathbf{x}_k)$$ \(a(\cdot)\)は活性化関数(多くの場合シグモイド)、\( \{ \gamma \}_{j=1}^M \) は未知パラメータです。で、モデル・カリブレーションを通じてウェイトをつくります。この推定量はデザイン一致で漸近不偏であり、分散推定量も一致で、シミュレーションでは局所多項回帰推定量を上回る、のだそうです。
射影追跡と単一インデクスモデル
射影追跡はニューラル・ネットワークの\( \{ \alpha_j a(\cdot) \}_{j=1}^M \)を未知の平滑化関数\( \{ \alpha_j (\cdot) \}_{j=1}^M \)に置き換えたモデルです。これを使ったモデル・アシステッド推定量の提案は見当たらないのですが、Wang(2009)は\(M=1\)の場合について検討しています。つまり$$ \mu(\mathbf{x}_k = \alpha(\gamma^\top \mathbf{x}_k) $$ というモデルです。これを単一インデクスモデルといいます。Wangは、$$ \mu_*(z) = E(\mu(\mathbf{x}_k | \gamma^\top \mathbf{x}_k = z] $$ という近似を行い、\(\mu_*(\cdot)\)を多項スプラインで推定しました。[???? なにいってんだかよくわからん]
この推定量もデザイン一致、漸近正規、分散推定量が一致です。
8.5 K最近隣法
KNNとは、要素\(k\)と共変量が近隣であるような\( \{y_l\}_{l \in L_k \cap s}\)(ただし\(|L_k| = K\))を平均して\(y_k\)を予測する手法です。Baffetta et al.(2009)はKNNを使ってモデル・アシステッド推定量を提案しています。[…略…]
8.6 木ベースの手法
木ベースの手法についてはESL本の9.2節を見てください。
回帰木
CARTやMARSのような標本の再帰分割によって回帰木を使うという手法を、全面的に使ったモデル・アシステッド推定は、我々は見たことがないです。でもその漸近理論の多くはすでにToth & Eltinge (2011)が与えています。彼らは、サーベイ・ウェイテッドな回帰木をごく一般的な回帰関数の推定量としてみたとき、漸近的にデザイン\(L^2\)一致性があることを示しています。こうしたサーベイ・ウェイテッド回帰木はモデル・アシステッド推定に容易に組み込めるでしょう。
内生的事後層別
すでに述べたように、事後層別推定量(PSD)はGREGの特殊ケースで、インジケータがカテゴリカル共変量であるものです。\(U = \cup_{h=1}^H U_h\)を母集団の分割とし、\(\mathbf{x}_k = [ \mathbf{1}_{k \in U_h} ]_{h=1}^H\)として、$$ DIFF(y, \mathbf{x}_\top \hat{\mathbf{B}}) = \sum_{k \in U} \mathbf{x}^\top_k \hat{\mathbf{B}} + \sum_{k \in s} \frac{y_k – \mathbf{x}^\top_k \hat{\mathbf{B}} }{\pi_k} $$ です。[なるほど… たしかにGREGだ…]
PSEにおいては、インジケータは標本においては既知で、母集団についてはその合計 \( \sum_{k \in U} \mathbf{1}_{k \in U_h} = | U_h | = N_h \)のみ既知です。PSEにおける作業モデルとは、事後層内の平均で、定数です。
森林の研究では、衛星画像の分類に基づく事後層別の構築が求められます。画像の教師付き学習のためには正解データが必要ですので、サーベイから正解データを得ます。サーベイのデータで分類アルゴリズムを訓練して、分類したイメージでサーベイデータを事後層別するのです。これを内生的事後層別といいます。[へええええ？ よくわからんが面白いなあ]
画像の分類を、あてはめたモデルからの予測の分類として表現すると、既知の分割点\( \{\tau_h\}_{h=0}^H \)について $$ \tilde{x}_k = [\mathbf{1}_{\tau_{h-1} \lt \hat{m}(\mathbf{x}_k) \leq \tau_h}]_{h=1}^H$$ を得て、PSE推定量にプラグインします。これを内生的事後層別推定量(EPSE)といいます。
この手法は、共変量の空間を矩形に分割しているという点で木ベースの手法です。[…以下略。全然知らん世界の話で面白い]
8.7 モデル選択と縮小手法
共変量ベクトル\(\mathbf{x}_k\)の次元数が大きすぎ、高い相関があったりすることもあります。Silva & Skinner (1997)は、復元なし単純無作為標本についての最良サブセットと前向きステップワイズ回帰で有限母集団の量を推定し、いくつかの係数をうまくゼロに設定しています。[どんな状況の話なんだろうか]
係数をゼロに移動させる縮小手法は、予測の正確性を非常に効率的に改善し、よってGREGウェイトの性質を改善します。例としてリッジ・カリブレーションがあります。
lassoはモデル選択と係数の縮小を同時にやる手法です。McConville(2011)はサーベイ・ウェイテッドlasso回帰係数をプラグインしたモデル・アシステッド推定を提案しています。デザイン一致性があり、推定量についての中心極限理論があり、分散推定量もデザイン一致です。lassoサーベイ回帰ウェイトもつくれます。
adaptive lassoを使うという提案もあります。共変量に高い相関があるとき、GREGよりMSEが小さくなるそうです。
8.8 アンサンブル学習と関連手法
世の中にはアンサンブル学習の手法がいっぱいあります。ここでの文脈でいえば、予測子\(\hat{\mu}_j(\cdot)\)をM個つくるわけです。それぞれの予測子は一部のデータに基づく、なんらかの予測手法による予測子です。で、$$ \hat{\mathbf{x}}_k) = \sum_{j=1}^M \omega_j \hat{\mu}_j(\mathbf{x}_k) $$ とします。\(\omega_j\)をどうやって決めるかがポイントになります。
モデル平均
Li & Opsomer (2006)は、回帰推定量の構築において最良のモデルを選ぶ代わりにモデル平均を使うという方法について検討しています。シミュレーションによればある程度効率がよいらしいのですが、まだ研究が進んでいません。おそらく、今後は補足情報の量が大きくなりますから、こういう手法が大事になるでしょう。[なるほど]
ランダム・フォレスト
Tipson, Opsomer, & Moisen (2013)は内生的事後層別のなかでさらにランダム・フォレストを使っています。直接にランダムフォレストを使ったモデル・アシステッド推定量は見たことがありません。
バギング
Wang, Opsomer, & Wang (2014 Survey Methodology)はデザインベースのセッティングにおけるブートストラップ累積(バギング)について検討しています。彼らによれば、バギングを使った推定量はデザイン一致で、効率もよいそうです。また、彼らはreplicateウェイト(分散推定のためにサーベイデータについていることが多いです)を使ってバギング推定量を作る方法を示しています。[へええ… よくわからんが面白そう]
直接にバギングを使ったモデル・アシステッド推定量は見たことがありませんが、今後有望な領域だと思います。
9. 考察
[…]
本論文で紹介した新しい手法が発展するのと同時に、モデル・アシステッド推定についての考え方の変化が生じているように思います。単純化してしまうと、伝統的な焦点は、推定量を解釈可能な母集団コントロール合計に一致させることにありました。事後層別でも、比推定でも、ほとんどのGREGアプローチでもそうです。いっぽう最近では、推定量の正確性を増すためにどういう予測モデルを作るかという点に焦点が当てられています。このような、統計的手法を予測の道具としてみるというアイデアは、Hastieらのいう「統計的学習」というアイデアと関連しています。この分野の研究者・実務家にとっては、本論文で紹介した手法の多くはおなじみのものだったでしょう。
モデル・アシステッド推定の文脈では、統計的学習の道具の向き不向きを評価する際にはなにより予測能力が大事なわけですが、他の考慮事項として以下の2点があります。
- デザイン・ウェイトを組み込めるか。推定量のデザイン一致性を維持するために大事です。[つまりあれね、\(\hat{m}(\mathbf{x}_k)\)を標本で推定する際に\(\pi_k\)を使えるかってことね]
- 反応変数の幅広い範囲についてよい予測ができるか。モデル・アシステッド推定の応用場面の多くにおいて、調査ウェイトを作って全ての変数に適用するという方法が使われます。ですから、特定の変数についてはすごく正確に推定できるけどほかの変数についてはだめ、というのでは困るわけです。この点は、調査変数、補足変数、手法の相互作用に依存しますから、どの文脈でどの手法が良いという一般的な言明は難しいです。とはいえ、これらの手法を実務に適用するには、いろんな調査変数と目標変数を通じた頑健性の評価が必要になります。ウェイトの安定性は多くの場合この文脈で検討されています。変動の大きなウェイトは、ウェイトとあまり相関しない変数や負の相関を持つ変数について不正確な推定につながるからです。[な・る・ほ・ど・ね… 納得…]
ところで、調査データは有限母集団の量の推定に使われるだけではなく、モデルのあてはめのためにもつかわれます。[そうそうそう！ この論文の主旨じゃないことはわかってんだけど、どっちかというとそっちのニーズが高いということもあるんですよ！ バリバリに有意抽出な標本でモデルを作って消費者についての理解を深めたい、とかさ]
この文脈では、パラメータ推定値はふつう、有限母集団の推定方程式の標本ウェイト版を解くことで得られます。たとえば、サーベイ・ウェイテッドなスコア・ベクトルが0になると設定して疑似最尤推定量(PMLE)を得て、それを解いて未知パラメータを推定したりします。
PMLEのような手法において、サーベイ・ウェイトは本論文で紹介したモデル・アシステッドなウェイトで置き換えてもいいですし、そうすることで、結果として得られるパラメータ推定量の特性にも影響があるでしょう。パラメータ推定におけるモデル・アシステッド・ウェイトの効果を研究するためには影響関数というのを使うことができます。この論文の掲載号のLumley & Scottをみてください。これはいまだ未開拓な領域です。
[へええええ！ たとえば抽出確率が不均等な調査データを使ってSEMを推定するとき、Mplusとかだと調査ウェイトが使えるけれど、素直に抽出確率の推定値の逆数を使うのではなく、たとえばなにかの補助変数を使ってつくったGREGウェイトを使ったっていいよ、ってことね。まじ？ そうなの？ そうなのか… ]
本論文を通じてお伝えしたかったのは、いろんな手法はあるけれど、推定量を構築する際の基本原理と、その統計的特性の研究というのは驚くほどシンプルなのだということです。[ちょっとわらっちゃった。先生、ハイブローすぎます。漸近挙動の解析はシンプルとはいいがたいっすよ]
4章で述べたレシピを使えば、実務家はいろんな共変量と予測手法を組み込むことができます。また3章・5章で示した漸近的枠組みは、推論のための道具の統計的評価と開発の一般的な経路を示しています。予測手法の母集団へのあてはめに基づく差分推定量が、デザイン・ベースの要素とモデル・ベースの要素のをつなぐリンクとなるのです。
——-
いやあ、面白かったー！
振り返ると、GREG推定量という考え方にはこれまで何度も触れたことがあった。このブログの記事を振り返ると
- 読了: Kalton & Flores-Cervantes (2003) 「ウェイトバック集計」のためのウェイト値のつくりかた (2015/04/01)
- 読了: Mercer, et al. (2017) 因果推論から眺めた非確率的標本調査 (2022/08/03)
- 読了：Rao (2020) 非確率標本に基づく統計的推測手法レビュー (2022/09/29)
あたりがそうだし、以前かなり真剣に読んだ Rao & Molina (2015) “Small Area Estimation”にも出てきた。でも、自分の問題として捉えられなかったんですよね。標本とは別になにかリッチな母集団データ、たとえばセンサスの個票データとかがあるときの話でしょ？ 公的調査とかで使うんじゃないの、と。
この論文のように、有限母集団についてのデザインベース推論だけど、母集団のモデルとみなしうるような予測モデルを補助的に使うのだ、その予測モデルとは事後層別だろうがlassoだろうがランダム・フォレストだろうがなんでもいいのだ… という風にはこれまで考えたことがなかった。広く視野をとれば、実は自分の仕事の中でも、著者らのいうモデル・アシステッド推論を使う余地があったかもしれない。ううむ。