回帰クリギング

応用統計学および地球統計学において、回帰クリギング（RK）は、従属変数を補助変数（デジタル標高モデリング、リモートセンシング/画像、主題図から得られるパラメータなど）に回帰させ、回帰残差を補間（クリギング）する空間予測手法です。これは、補助予測変数を直接クリギング重みの解法に用いる、ユニバーサルクリギングや外部ドリフト付きクリギングなどと呼ばれる補間法と数学的に等価です。 ^{[ 1 ]}

回帰クリギングは、空間データに対する最良線形不偏予測器（BLUP）の実装であり、空間変動の普遍モデルを仮定した最良線形補間器です。Matheron（1969）は、ある場所における目的変数の値は、決定論的成分と確率的成分の和としてモデル化できると提案しました。^{[ 2 ]}

${\displaystyle Z(\mathbf {s} )=m(\mathbf {s} )+\varepsilon '(\mathbf {s} )+\varepsilon ''}$

彼はこれを空間変動の普遍モデルと呼んだ。空間変動の決定論的要素と確率的要素はそれぞれ個別にモデル化できる。この2つのアプローチを組み合わせると、次の式が得られる。

${\displaystyle {\hat {z}}(\mathbf {s} _{0})={\hat {m}}(\mathbf {s} _{0})+{\hat {e}}(\mathbf {s} _{0})=\sum \limits _{k=0}^{p}{{\hat {\beta }}_{k}\cdot q_{k}(\mathbf {s} _{0})}+\sum \limits _{i=1}^{n}\lambda _{i}\cdot e(\mathbf {s} _{i})}$

どこ${\displaystyle {\hat {m}}(\mathbf {s} _{0})}$適合された決定論的部分であり、${\displaystyle {\hat {e}}(\mathbf {s} _{0})}$は補間された残差であり、${\displaystyle {\hat {\beta }}_{k}}$推定された決定論的モデル係数（${\displaystyle {\hat {\beta }}_{0}}$は推定切片です。${\displaystyle \lambda _{i}}$クリギング重みは残差の空間依存構造によって決定され、${\displaystyle e(\mathbf {s} _{i})}$位置における残差は${\displaystyle {\mathbf {s} }_{i}}$回帰係数${\displaystyle {\hat {\beta }}_{k}}$は、例えば通常の最小二乗法（OLS）や、最適には一般化最小二乗法（GLS）などの何らかの適合方法によってサンプルから推定できます。 ^{[ 3 ]}

${\displaystyle \mathbf {\hat {\beta }} _{\mathtt {GLS}}=\left(\mathbf {q} ^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {q} \right)^{-\mathbf {1} }\cdot \mathbf {q} ^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {z} }$

どこ${\displaystyle \mathbf {\hat {\beta }} _{\mathtt {GLS}}}$は推定回帰係数のベクトルであり、${\displaystyle \mathbf {C} }$は残差の共分散行列であり、${\displaystyle {\mathbf {q} }}$はサンプリング地点における予測変数の行列であり、${\displaystyle \mathbf {z} }$は、目的変数の測定値のベクトルです。回帰係数のGLS推定は、実際には地理加重回帰の特殊なケースです。この場合、残差間の空間的自己相関を考慮するために、重みが客観的に決定されます。

変動の決定論的部分（回帰部分）が推定されたら、残差をクリギングで補間して推定トレンドに追加できます。残差の推定は反復プロセスです。まず、変動の決定論的部分をOLSを使用して推定し、次に残差の共分散関数を使用してGLS係数を取得します。次に、これらを使用して残差を再計算し、更新された共分散関数を計算し、これを繰り返します。これは多くの地質統計学者によって適切な手順として推奨されていますが、Kitanidis（1994）は、OLS残差から導出された共分散関数（つまり1回の反復）の使用は、数回の反復後に導出された関数と十分に異ならないため、多くの場合満足できるものであることを示しました。つまり、最終的な予測にはあまり影響しません。MinasnyとMcBratney（2007）は同様の結果を報告しており、より高度な統計的手法を使用するよりも、より高品質のデータを使用することの方が重要であるようです。^{[ 4 ]}

行列表記では、回帰クリギングは一般的に次のように記述されます。^{[ 5 ]}

${\displaystyle {\hat {z}}_{\mathtt {RK}}(\mathbf {s} _{0})=\mathbf {q} _{\mathbf {0} }^{\mathbf {T} }\cdot \mathbf {\hat {\beta }} _{\mathtt {GLS}}+\mathbf {\lambda } _{\mathbf {0} }^{\mathbf {T} }\cdot (\mathbf {z} -\mathbf {q} \cdot \mathbf {\hat {\beta }} _{\mathtt {GLS}})}$

どこ${\displaystyle {\hat {z}}({\mathbf {s} }_{0})}$位置における予測値${\displaystyle {\mathbf {s} }_{0}}$、${\displaystyle {\mathbf {q} }_{\mathbf {0} }}$はベクトルです${\displaystyle p+1}$予測因子と${\displaystyle \mathbf {\lambda } _{\mathbf {0} }}$はベクトルです${\displaystyle n}$残差を補間するために使用されるクリギング重み。RKモデルは、空間データの最良の線形予測器と考えられています。^{[ 5 ] [ 6 ]}地理空間と特徴空間の両方で新しい場所の位置（外挿）を反映する予測分散があります。

${\displaystyle {\hat {\sigma }}_{\mathtt {RK}}^{2}(\mathbf {s} _{0})=(C_{0}+C_{1})-\mathbf {c} _{\mathbf {0} }^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {c} _{\mathbf {0} }+\left(\mathbf {q} _{\mathbf {0} }-\mathbf {q} ^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {c} _{\mathbf {0} }\right)^{\mathbf {T} }\cdot \left(\mathbf {q} ^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {q} \right)^{\mathbf {-1} }\cdot \left(\mathbf {q} _{\mathbf {0} }-\mathbf {q} ^{\mathbf {T} }\cdot \mathbf {C} ^{-\mathbf {1} }\cdot \mathbf {c} _{\mathbf {0} }\right)}$

どこ${\displaystyle C_{0}+C_{1}}$敷居の変動と${\displaystyle {\mathbf {c} }_{0}}$は、未訪問地点における残差の共分散ベクトルです。

多くの（地理）統計学者は、空間データに対する最良の線形不偏予測モデルは1つしかない（例えば回帰クリギング）と考えており、通常のクリギング、環境相関、ポリゴンごとの値の平均化、逆距離補間などの他のすべての手法は、その特殊なケースと見なすことができる。残差に空間自己相関がない場合（純粋なナゲット効果）、回帰クリギングは純粋な多重線形回帰に収束する。なぜなら、共分散行列（${\displaystyle \mathbf {C} }$）は単位行列になります。同様に、目的変数が補助予測変数と相関がない場合、決定論的部分が（グローバル）平均値と等しくなるため、回帰クリギングモデルは通常のクリギングモデルに帰着します。したがって、純粋クリギングと純粋回帰は、回帰クリギングの特殊なケースとしてのみ考慮する必要があります（図を参照）。

地質統計学の文献では、本質的に同じ、あるいは少なくとも非常に類似した手法に対して、さまざまな用語が用いられています。これは利用者を混乱させ、マッピングプロジェクトに適した手法の選択を妨げます。実際、ユニバーサルクリギング、外部ドリフト付きクリギング、回帰クリギングは、基本的に同じ手法です。

Matheron (1969) は当初この手法をLe krigeage universelと名付けましたが、この手法はトレンドが座標の関数としてモデル化されるクリギングの一般化されたケースとして意図されていました。そのため、多くの著者は、座標のみが予測変数として使用されるケースにユニバーサル クリギング(UK) という用語を留保しています。変動の決定論的部分 (ドリフト) が座標ではなく補助変数の線形関数として外部で定義される場合、外部ドリフト付きクリギング(KED) という用語が好まれます (Hengl 2007、「回帰クリギングについて: 方程式からケーススタディまで」による)。UK または KED の場合、予測はクリギングと同様に行われますが、残差の共分散行列が補助予測変数で拡張される点が異なります。ただし、ドリフトと残差は別々に推定してから合計することもできます。この手順は Ahmed et al. (1987) と Odeh et al. によって提案されました。 (1995) は後にこれを回帰クリギングと名付け、一方 Goovaerts (1997) はトレンドモデルを使用したクリギングという用語を補間法のファミリーを指すために使用し、RK を変動する局所平均を持つ単純クリギングと呼んでいます。 Minasny と McBratney (2007) は、この手法を単に経験的最良線形不偏予測器、すなわちE-BLUPと呼んでいます。^{[ 7 ] [ 8 ] [ 9 ] [ 4 ]}

KEDの場合、新しい場所での予測は以下のように行われます。

${\displaystyle {\hat {z}}_{\mathtt {KED}}(\mathbf {s} _{0})=\sum \limits _{i=1}^{n}w_{i}^{\mathtt {KED}}(\mathbf {s} _{0})\cdot z(\mathbf {s} _{i})}$

のために

${\displaystyle \sum \limits _{i=1}^{n}w_{i}^{\mathtt {KED}}(\mathbf {s} _{0})\cdot q_{k}(\mathbf {s} _{i})=q_{k}(\mathbf {s} _{0})}$

のために${\displaystyle k=1,\ldots ,p}$または行列表記では：

${\displaystyle {\hat {z}}_{\mathtt {KED}}(\mathbf {s} _{0})=\mathbf {\delta } _{\mathbf {0} }^{\mathbf {T} }\cdot \mathbf {z} }$

どこ${\displaystyle z}$は目的変数です。${\displaystyle q_{k}}$は予測変数、つまり新しい場所での値です。${\displaystyle ({\mathbf {s} }_{0})}$、${\displaystyle {\mathbf {\delta } }_{\mathbf {0} }}$は KED 重みのベクトルです (${\displaystyle w_{i}^{\mathtt {KED}}}$)${\displaystyle p}$は予測変数の数であり、${\displaystyle \mathbf {z} }$はベクトルです${\displaystyle n}$主要地点での観測。KED重みは拡張行列を用いて解かれる。

${\displaystyle \mathbf {\lambda } _{\mathbf {0} }^{\mathtt {KED}}=\left\{w_{1}^{\mathtt {KED}}(\mathbf {s} _{0}),\ldots ,w_{n}^{\mathtt {KED}}(\mathbf {s} _{0}),\varphi _{0}(\mathbf {s} _{0}),\ldots ,\varphi _{p}(\mathbf {s} _{0})\right\}^{\mathbf {T} }=\mathbf {C} ^{{\mathtt {KED}}-1}\cdot \mathbf {c} _{\mathbf {0} }^{\mathtt {KED}}}$

どこ${\displaystyle {\mathbf {\lambda } }_{\mathbf {0} }^{\mathtt {KED}}}$は、解決された重みのベクトルです。${\displaystyle \varphi _{p}}$ラグランジュ乗数は、${\displaystyle {\mathbf {C} }^{\mathtt {KED}}}$は残差の拡張共分散行列であり、${\displaystyle {\mathbf {c} }_{\mathbf {0} }^{\mathtt {KED}}}$これは、新しい位置における共分散の拡張ベクトルです。

KEDの場合、残差の拡張共分散行列は次のようになります（Webster and Oliver、2007; p.  183）：^{[ 10 ]}

${\displaystyle \mathbf {C} ^{\mathtt {KED}}=\left[{\begin{array}{ccccccc}C(\mathbf {s} _{1},\mathbf {s} _{1})&\cdots &C(\mathbf {s} _{1},\mathbf {s} _{n})&1&q_{1}(\mathbf {s} _{1})&\cdots &q_{p}(\mathbf {s} _{1})\\\vdots &&\vdots &\vdots &\vdots &&\vdots \\C(\mathbf {s} _{n},\mathbf {s} _{1})&\cdots &C(\mathbf {s} _{n},\mathbf {s} _{n})&1&q_{1}(\mathbf {s} _{n})&\cdots &q_{p}(\mathbf {s} _{n})\\1&\cdots &1&0&0&\cdots &0\\q_{1}(\mathbf {s} _{1})&\cdots &q_{1}(\mathbf {s} _{n})&0&0&\cdots &0\\\vdots &&\vdots &\vdots &\vdots &&\vdots \\q_{p}(\mathbf {s} _{1})&\cdots &q_{p}(\mathbf {s} _{n})&0&0&\cdots &0\end{array}}\right]}$

そして${\displaystyle \mathbf {c} _{\mathbf {0} }^{\mathtt {KED}}}$このような：

${\displaystyle \mathbf {c} _{\mathbf {0} }^{\mathtt {KED}}=\left\{C(\mathbf {s} _{0},\mathbf {s} _{1}),\ldots ,C(\mathbf {s} _{0},\mathbf {s} _{n}),q_{0}(\mathbf {s} _{0}),q_{1}(\mathbf {s} _{0}),\ldots ,q_{p}(\mathbf {s} _{0})\right\}^{\mathbf {T} };q_{0}(\mathbf {s} _{0})=1}$

したがって、KEDは通常のクリギングと全く同じように見えますが、共分散行列/ベクトルに補助予測変数の値が追加されている点が異なります。

KEDは一見RKよりも計算が簡単そうに見えますが、 KEDのバリオグラムのパラメータも回帰残差から推定する必要があるため、別の回帰モデリング手順が必要です。残差間に空間相関がある可能性が高いため、この回帰はGLSであるべきです。多くの分析者は代わりにOLS残差を使用しますが、これはGLS残差とそれほど違いがないかもしれません。しかし、空間相関がある場合は最適ではなく、実際、クラスター化されたサンプルポイントやサンプル数が比較的小さい場合、かなり異なる可能性があります（${\displaystyle \ll 200}$)

KEDの限界の一つは、共変量が空間的に滑らかに変化しない場合に拡張行列が不安定になることである。RKの利点は、トレンド推定と残差の空間予測を明示的に分離することで、KEDで使用できる単純な線形手法ではなく、任意の複雑な形式の回帰を使用できる点にある。さらに、2つの補間成分を個別に解釈することも可能となる。回帰を重視する理由は、変動の決定論的部分（回帰）の適合が、確率的部分（残差）の適合よりも最終的なマップの品質向上にしばしば有利となるためでもある。

回帰クリギングは、例えばR統計計算環境において、gstatパッケージやgeoRパッケージなどを用いて自動化できます。典型的な入力/出力は以下のとおりです。

入力:

• 補間セット（ポイントマップ）—${\displaystyle z(\mathbf {s} _{i})}$${\displaystyle i=1,\ldots ,n}$主要な拠点において。
• 最小値と最大値の期待値および測定精度（${\displaystyle \Delta z}$);
• 連続予測変数（ラスターマップ）—${\displaystyle q(\mathbf {s} )}$; 未訪問の新たな場所で
• 離散予測因子（ポリゴンマップ）
• 検証セット（ポイントマップ）—${\displaystyle z*(\mathbf {s} _{j})}$${\displaystyle j=1,\ldots ,l}$（任意）
• ラグ間隔と制限距離（バリオグラムに適合させるために必要）

出力:

• 予測結果と相対予測誤差のマップ。
• 最適な予測因子サブセットと相関の有意性（調整済みR二乗値）
• バリオグラムモデルパラメータ（例：${\displaystyle C_{0}}$、${\displaystyle C_{1}}$、${\displaystyle R}$）
• GLSドリフトモデル係数
• 検証ポイントにおける予測精度：平均予測誤差（MPE）および二乗平均平方根予測誤差（RMSPE）。

回帰クリギングは、気象学、気候学、土壌マッピング、地質マッピング、種分布モデリングなど、さまざまな応用分野で使用されています。回帰クリギングを通常のクリギングなどと比較した場合の唯一の要件は、1つ以上の共変量レイヤーが存在し、それらが対象となる特徴と有意に相関していることです。回帰クリギングの一般的な応用例をいくつか挙げます。

• 地統計学的マッピング：回帰クリギングにより、ハイブリッド地統計学的手法を使用して、例えば土壌特性の空間分布をモデル化することができます。
• マップのダウンスケール：回帰クリギングは、既存のさまざまなグリッドマップをダウンスケールするためのフレームワークとして使用できます。この場合、共変量レイヤーは、元のポイントデータよりも高い解像度（サンプリング強度に対応）で利用可能である必要があります。^{[ 11 ]}
• 誤差伝播：回帰クリギングモデルを用いて生成されたシミュレーションマップは、シナリオテストや伝播する不確実性の推定に使用できます。

回帰クリギングに基づくアルゴリズムは、可能な共変量の数が日々増加しているため、地球統計学においてますます重要な役割を果たしています。^{[ 1 ]}例えば、DEM は現在、多くのソースから入手可能です。地形の詳細かつ正確な画像は、SPOTやASTERなどのリモートセンシングシステムから注文できます。SPOT5 は高解像度ステレオ (HRS) スキャナを提供しており、これを使用して最大 5 m の解像度で DEM を作成できます。^{[ 12 ]}航空レーザースキャナを使用すると、より細かい標高の差も取得できます。データのコストは、技術の進歩に伴い無料または価格が下がっています。NASA は、2000 年にシャトル レーダー地形ミッションで世界の地形の大部分を記録しました。 ^{[ 13 ]} 2004 年の夏以降、これらのデータは、ほぼ全世界について約 90 m の解像度 (北米大陸については約 30 m の解像度) で (たとえばUSGS ftpを介して)入手可能になっています。同様に、MODISマルチスペクトル画像は250mの解像度で無料でダウンロードできます。Landsat画像の膨大な無料リポジトリも、Wayback MachineのGlobal Land Cover Facility (GLCF)アーカイブ（2012年6月22日）からダウンロードできます。

• 第2章、回帰クリギング、トミスラフ・ヘングル（2009）『地質統計マッピングの実践ガイド』、291ページ、ISBN 978-90-9024981-0。
• Hengl T.、Heuvelink GBM、Rossiter DG (2007)「回帰クリギングについて：方程式からケーススタディまで」Computers & Geosciences . 33 (10): 1301–1315 . Bibcode : 2007CG.....33.1301H . doi : 10.1016/j.cageo.2007.05.001 .

• Gstatパッケージ（KEDを実装）
• GeoRパッケージ（KEDを実装）