1何のための手法か

観測が集団に入れ子になっているとき、独立を仮定した回帰では何が壊れるのか。

米原〜関ケ原の調査では、各集落の語彙消失率が応答であり、集落は中学校の学区に属している。同じ学区の集落は通学圏を共有し、同じ学校で同じ世代がことばを交わす。したがって、学区内の集落が互いに独立な観測だという前提は、そのままでは成り立たない。

独立を仮定した回帰をそのまま当てると何が起きるか。観点2の枠組みでは、推定のばらつきは情報量 \(I(\theta)\) の逆数として決まり、信頼区間はそこから出る。観測が重複した情報を含んでいるのに独立として数えると、実際より多くの情報があることになり、\(I(\theta)\) を過大に見積もる。結果として標準誤差が小さく出て、信頼区間が狭くなり、検定は有意を出しやすくなる。「手続きを繰り返したときの当たり率」という頻度論の約束が、そこで破れる。

対処は2通りある。ひとつは学区ごとに別々の切片を固定効果として推定する方法、もうひとつは学区ごとの切片のばらつきを分散パラメータ1つで表す方法である。線形混合モデルは後者にあたる。

2ランダム効果とは何をモデルに入れているのか

集団ごとに別々の値を推定するのではなく、値の散らばりの幅だけを推定する。

固定効果で書く場合
\[ y_i \;=\; \beta_0 \;+\; \sum_{d=2}^{D} \gamma_d\,\mathbf{1}\!\left[d(i)=d\right] \;+\; \mathbf{x}_i^{\top}\boldsymbol{\beta} \;+\; \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}\!\left(0,\sigma_\varepsilon^2\right) \]

\(\gamma_d\) … 学区 \(d\) の切片。基準の1区を除いて \(D-1\) 個のパラメータが増える
\(d(i)\) … 集落 \(i\) が属する学区。この研究では \(D=5\)

ランダム効果で書く場合
\[ y_i \;=\; \beta_0 \;+\; u_{d(i)} \;+\; \mathbf{x}_i^{\top}\boldsymbol{\beta} \;+\; \varepsilon_i, \qquad u_d \sim \mathcal{N}\!\left(0,\sigma_u^2\right),\quad \varepsilon_i \sim \mathcal{N}\!\left(0,\sigma_\varepsilon^2\right) \]

\(u_d\) … 学区 \(d\) の切片のずれ。これ自体はパラメータではなく、正規分布から引かれた確率変数と見る
\(\sigma_u^2\) … 推定するのはこの1つだけ。学区が5つでも50でも増えない

差は「何を推定するか」にある。固定効果版は学区の数だけパラメータが増え、各学区は自分の観測だけで決まる。ランダム効果版が推定するのは散らばりの幅 \(\sigma_u^2\) ひとつで、個々の \(u_d\) はその幅のもとで全体平均のほうへ引き寄せられる。観測の少ない学区ほど強く引き寄せられるこの性質を縮約（部分プーリング）という。

独立をどこで使っているか ―― 周辺共分散
\[ \mathrm{Cov}\!\left(y_i,\,y_j\right) \;=\; \begin{cases} \sigma_u^2 & d(i)=d(j),\ i \neq j \\[4pt] 0 & d(i) \neq d(j) \end{cases} \qquad\Longrightarrow\qquad \mathrm{ICC} \;=\; \frac{\sigma_u^2}{\sigma_u^2+\sigma_\varepsilon^2} \]

誤差 \(\varepsilon_i\) は互いに独立だが、同じ学区の観測は \(u_d\) を共有するため、周辺的には相関する。その大きさが級内相関 ICC である。独立の仮定は、誤差の層に押し込められただけで消えていない。

「独立をどの行で使ったか」を指せることが要点である

ふつうの回帰が置く独立の仮定は、観測どうしが無関係だという主張である。混合モデルはそれを弱め、学区を条件づけたうえでの独立に置き換える。仮定が消えたのではなく、どこに置かれているかが変わった。統計学Ⅱで「独立を式のどの行で使ったか」を毎回指し示しているのは、この種の読み替えを自分でできるようにするためである。学区をまたぐ相関（たとえば隣接する学区どうしが似る）は、この式では依然として 0 と置かれている。そちらを扱うには別の構造が要る（第5節）。

3観点2との対応 ―― 最尤で何を最大にしているか

ランダム効果は積分で消し、残った尤度を最大にする。

\(u_d\) はパラメータではなく確率変数なので、尤度に残すことができない。そこで積分して消す。残るのは固定効果 \(\boldsymbol{\beta}\) と分散 \(\sigma_u^2,\sigma_\varepsilon^2\) についての関数であり、これを最大にする。

周辺尤度 ―― ランダム効果を積分消去する
\[ L\!\left(\boldsymbol{\beta},\sigma_u^2,\sigma_\varepsilon^2\right) \;=\; \int p\!\left(\mathbf{y}\mid \mathbf{u}\right)\, p\!\left(\mathbf{u}\right)\, d\mathbf{u} \qquad\Longrightarrow\qquad \mathbf{y} \sim \mathcal{N}\!\left(X\boldsymbol{\beta},\; ZGZ^{\top}+R\right) \]

\(Z\) … どの観測がどの学区に属するかを示す計画行列
\(G=\sigma_u^2 I\)、\(R=\sigma_\varepsilon^2 I\) … ランダム効果と誤差の分散
線形かつ正規のときに限り、この積分は解析的に解けて多変量正規分布になる。ベータ回帰のような非正規の応答では解けず、数値積分か MCMC が要る。

ここから先は観点2の定石どおりである。対数尤度を微分してゼロと置き、数値的に解く。推定のばらつきは情報量から出て、そこから標準誤差と信頼区間と検定が組み立てられる。パラメータは固定した未知の値であり、確率変動するのは標本のほうだと考える。

REML ―― 分散推定の偏りを直す代わりに、比較できる範囲が狭まる

最尤で分散成分を推定すると、固定効果を同時に推定した分だけ下向きに偏る。制限最尤（REML）は固定効果の影響を除いた尤度を使ってこの偏りを補正する。ただし REML の尤度は固定効果の計画行列に依存するため、固定効果の構造が違うモデルどうしでは、REML の尤度も AIC も比較できない（Pinheiro & Bates 2000）。説明変数の組を選ぶために AIC を使うなら、最尤で当てはめ直す必要がある。

4同じ構造をベイズで推定すると、何が変わるのか

構造は同じ。違うのは \(\sigma_u\) に点を与えるか、分布を与えるかである。

この研究では、同じ米原〜関ケ原のデータが後に階層ベイズで解き直されている。ベイズ推定のページが扱っているのがそれで、モデル式の (1|district) は本ページの \(u_{d(i)}\) と同じものである。学区をランダム効果に置くという構造自体は変わっていない。変わったのは推定の枠組みと、応答に置く分布である。

　線形混合モデル（本ページ）階層ベイズ（ベイズ推定）
階層構造学区をランダム効果学区をランダム効果（同じ）
応答の分布正規分布ベータ分布・logit リンク
\(\sigma_u\) の扱い点推定を1つ返す事後分布を返す
係数の不確実性信頼区間（手続きの当たり率）確信区間（パラメータの確率）
出力できる量p 値\(P(\beta>0)\)
選択の道具AIC・BICLOO 近似交差検証

同じデータ・同じ階層構造に対する2つの枠組みの対比。応答の分布が変わっている点についてはベータ回帰のページを参照。

クラスタ数が少ないとき、この違いが効いてくる

学区は5つ、集落は25である。\(\sigma_u^2\) は「5個の値の散らばり」から推定されることになり、情報がきわめて少ない。最尤はこの状況で \(\hat{\sigma}_u = 0\) という境界解を返すことがある。分散は負になれないので、尤度の最大が母数空間の端に来てしまう場合である。推定値が 0 だからといって学区差が存在しないわけではなく、「このデータからは幅を推定できなかった」という意味にすぎないが、出力される数値は同じ 0 である。

ベイズでは \(\sigma_u\) に事前分布を置くため、事後分布が 0 に一点集中することはない。区間として返るので、「0 に接しているが幅がある」ことが読める。実際、同じ研究の空間分散についてICAR のページが報告している sdcar は 0.04、95%確信区間 [0.00, 0.16] で、下限が 0 に接しながら区間として表示されている。最尤ならここは 0 という1つの数字になり、接しているのか幅があるのかは区別できない。

2つの解析の数値を混同してはいけない

本ページ第7節の表に載せた AIC や係数は、正規分布を仮定した最尤推定の出力である。ベイズ推定のページが報告する \(\beta = +0.231\)、\(P(\beta>0)=0.993\) は、ベータ分布を仮定した事後分布から得た別の量である。応答の分布もリンク関数も推定原理も違うため、係数の大きさを直接比べることはできない。同じ研究の同じデータでありながら、2つの数値は別のモデルの出力である。

5空間的自己相関を頻度論で扱う ―― 空間誤差モデル

学区で括れない「隣と似る」を、誤差の相関としてパラメータ化する。

第2節の周辺共分散は、学区をまたぐ相関を 0 と置いていた。しかし地理上に並んだ集落では、学区が違っても隣り合えば似る。この構造を扱うのが空間モデルである。揖斐川上流域の18地点を対象とした解析では、空間誤差モデル（spatial error model）が用いられている。

空間誤差モデル ―― 誤差そのものが空間的に相関する
\[ \mathbf{y} \;=\; X\boldsymbol{\beta} + \mathbf{u}, \qquad \mathbf{u} \;=\; \lambda W \mathbf{u} + \boldsymbol{\varepsilon} \qquad\Longrightarrow\qquad \mathbf{u} \;=\; \left(I-\lambda W\right)^{-1}\boldsymbol{\varepsilon} \]

\(W\) … 近隣行列。どの地点とどの地点が隣かを定める
\(\lambda\) … 空間自己相関の強さ。\(\lambda=0\) なら通常の回帰に戻る
ランダム効果と違い、空間構造は誤差の相関として書かれる。推定は最尤で行う。

ICAR との対比 ―― 同じ問題を、別の層で書いている

ICAR は各地点に空間ランダム効果 \(\phi_i\) を置き、その \(\phi\) に「隣どうしは近い値を取る」という事前分布を与える。空間誤差モデルは \(\phi\) を明示せず、誤差の共分散構造をパラメータ \(\lambda\) ひとつで表す。前者はベイズの枠組みで空間構造を事前分布として書き、後者は頻度論の枠組みで相関構造として書く。扱っている問題は同じで、どの層に構造を置くかが違う。

近隣の定義は、ここでも設計物である

この解析の近隣リストは、直線距離ではなく集落間の歩行コスト距離から作られている。山間の渓谷に位置する集落では、直線で近くても尾根を越える必要があれば実際の往来は遠い。歩行コスト関数（Aiken 1977、Langmuir 1984）を使うのはその補正である。ただし、どの関数を使うか、どこで切って隣とするかは分析者が決める。ICAR のページで「何を隣とするかはデータから決まらない」と述べたのと同じ論点が、頻度論側でもそのまま当てはまる。

6AICによるモデル選択 ―― 観点5への橋

当てはまりから複雑さを差し引く。情報理論の量が、頻度論の推定の上に乗る。

赤池情報量規準
\[ \mathrm{AIC} \;=\; -2\,\ell\!\left(\hat{\theta}\right) + 2k \]

\(\ell(\hat{\theta})\) … 最大対数尤度。第3節で最大化したもの
\(k\) … 推定したパラメータの数
第1項は当てはまりの良さ、第2項は複雑さの罰則である（Akaike 1974）。

AIC は観点5（情報理論）の量である。真の分布をモデルで表したときの無駄な符号長、すなわちカルバック・ライブラー情報量を、手元のデータから推定する形になっている。最尤で当てはめ、情報理論の基準で選ぶ。この2つの観点は対立せず、推定は頻度論、選択は情報理論という分担で組み合わされている。同じ観点5にある正規化圧縮距離が「圧縮できた量で近さを測る」のに対し、AIC は「符号長でモデルを比べる」という位置にある。

米原〜関ケ原では、全部分集合の総当たりで選んでいる

説明変数は経度・緯度・面積・標高・土地利用（建物用地・森林・田）・人口・人口密度・人口増加率・高齢化率・15歳未満率と多い。GISA2016 の解析では MuMIn::dredge によって全部分集合を総当たりし、AIC の低い上位を取り出したうえで、\(\Delta\mathrm{AIC} < 4\) のモデルを平均している。総合・関西系・東京系・美濃系それぞれの gain / loss について、同じ手続きが実施された。この探索の危うさについては第8節で述べる。

揖斐川上流域では、空間項の要否を AIC で判断している

空間誤差モデルの AIC は \(-52.203\)、同じ説明変数で空間項を持たない線形モデルの AIC は \(-50.466\) と報告されている。差は約 1.7 で、空間項を入れたほうが小さい。ただし AIC は相対指標であり、この差の大きさが何を意味するかは別途の判断になる。

7この研究での使い方

2つの地域、2つの応答、2つの空間の扱い。

(A) 米原〜関ケ原 ―― 語彙の消失率・獲得率と環境・社会的要因

  1. 応答の構成　滋賀県米原市と岐阜県関ケ原町の集落について、40項目106語の調査から、60歳以上世代（142人）から中学生世代（544人）へかけての語の得失率の平均値を集落ごとに算出する。関西系・美濃尾張系・東京系の3タイプそれぞれについて、消失率と獲得率を分けて作る。
  2. 説明変数　土地利用（建物用地・森林・田の割合）は国土交通省「国土数値情報 土地利用細分メッシュ」（平成18年度）から、人口・人口密度・人口増加率・高齢化率・15歳未満率・面積は総務省統計局 e-Stat の平成22年国勢調査（小地域）から取る。
  3. 標準化　説明変数を scale で標準化し、係数の大きさを比較できるようにする。
  4. モデル　中学校の学区をランダム効果とする線形混合モデルを当てる。
  5. 選択　全部分集合を総当たりし、AIC の最も低いモデルを各応答について選ぶ。
応答（AIC 最小モデル）AICBIClogLikdeviancedf.resid
関西系　消失率70.0576.53−30.0360.0522
東京系　消失率63.3074.98−22.6545.3018
美濃系　消失率75.2784.35−30.6461.2720
関西系　獲得率76.0381.22−34.0268.0323
東京系　獲得率66.1275.19−28.0652.1220
美濃系　獲得率69.8280.19−26.9153.8219

GISA2016 発表要旨の表1より、各応答について AIC が最小となったモデルの適合指標。係数については同表に11変数分が掲載されているが、ここでは本文が明示的に述べている知見のみを下に記す。

要旨本文が述べている知見は次の3つである。どのタイプでも、人口が多い集落ほど消失率・獲得率がともに低い。消失率に特有の要素は人口増加率であり、獲得率に特有の要素は面積である。すなわち人口が増えている場所では語が失われやすく、集落の面積が広ければ新しい語を得やすい。さらに東京系・美濃尾張系の項目では、経度が大きい側（関ケ原側）で消失率・獲得率がともに高い。

要旨は27集落と書くが、モデルに入ったのは25集落と読める

上の表の df.resid は、標本数から固定効果の数（切片を含む）を引いた値として読める。関西系獲得率は説明変数1つ＋切片で df.resid 23、東京系消失率は説明変数6つ＋切片で 18 である。いずれも \(n=25\) と置いたときに一致する。要旨は対象地域を27集落と記しているが、消失率・獲得率が得られた25集落がモデルに投入されたと考えるのが整合的である。作業記録も n=25 と記しており、後のベイズ再解析も25自治会を用いている。

(B) 揖斐川上流域 ―― 系統推定の合成変数と集落の条件

  1. 対象　岐阜県西北部、揖斐川上流域の18地点。徳山村8地点、藤橋村2地点、坂内村3地点、久瀬村5地点である。
  2. 言語データ　山田（1978, 1980）が報告した47項目の基礎語彙を用いる。
  3. 系統推定　基礎語彙について系統推定を行う。枝の信頼度はブートストラップ法1000回で評価し、徳山村とそれ以外が 68.3%、徳山村・坂内村とそれ以外が 95.2% の信頼度で分岐すると報告されている。
  4. 合成変数の作成　系統推定で算出された距離に多次元尺度構成法を適用し、集落間距離の分散を最大にする第1軸（MDS1）を目的変数とする。
  5. 説明変数　戸数、人口（1875年〜1920年）、田の面積、畑の面積、石高、標高を用いる。いずれも揖斐郡志と角川日本地名大辞典に依る。
  6. 空間誤差モデル　歩行コスト距離から近隣リストを作り、空間誤差モデルを最尤推定で当てはめる。
項目推定値標準誤差z 値p 値
切片0.108437710.039873662.71950.0065374
戸数0.003521940.001069483.29310.0009908
人口−0.000563700.00017643−3.19510.0013979
田−0.001341130.00067133−1.99770.0457472
畑−0.001156080.00037578−3.07650.0020944
石高0.000293810.000180601.62690.1037645
標高−0.000431520.00013314−3.24110.0011909

空間誤差モデルの偏回帰係数。発表資料に掲載された errorsarlm の出力からの引用。観測数18、推定パラメータ数9。\(\lambda = -1.2505\)、LR 検定統計量 3.7368（p = 0.053227）、Wald 統計量 24.445（p = 7.648e−07）、AIC −52.203（線形モデルでは −50.466）。

資料はこの結果を、MDS1 が増加すると戸数と石高が増加し、田・畑・標高が減少する関係として読んでいる。石高については p 値が 0.104 で 0.05 をやや上回るが、有意傾向とみなせるレベルだと述べている。

\(\lambda\) が負である

推定された \(\lambda\) は \(-1.2505\) と負である。空間モデルで通常想定されるのは「隣どうしは似る」という正の空間自己相関であり、負の \(\lambda\) はその逆、隣どうしが異なる方向の構造を意味する。山間の渓谷に沿って並ぶ集落について、この符号が何を反映しているのかは資料には記されていない。あわせて LR 検定の p 値は 0.053 で、空間項そのものの必要性は明確には支持されていない。観測数18に対してパラメータ9という規模も、この検定の解釈を慎重にさせる条件である。

8限界 ―― この手法で言えないこと

「ランダム効果の分散がゼロ」は、ふつうの検定にかけられない

\(H_0:\sigma_u^2=0\) を尤度比検定で調べたくなるが、この帰無仮説は母数空間の境界にある。分散は負になれないため、帰無仮説のもとでパラメータが内点にあるという通常の正則条件が満たされない。このとき尤度比統計量は自由度1のカイ二乗分布に従わず、Self & Liang (1987) が一般の枠組みを与え、Stram & Lee (1994) が縦断混合モデルについて \(\chi^2_0\) と \(\chi^2_1\) の 50 対 50 の混合になることを示している。標準的なカイ二乗分布を当てると保守的な側に外れる、すなわち有意になりにくい方向に誤る。

AIC は相対指標であり、同じデータ・同じ応答でなければ比べられない

AIC の絶対値そのものには意味がない。上の表の 70.05 や \(-52.203\) は、それ単体では何も語らない。比較できるのは、同一のデータに対して同一の応答を当てはめたモデルどうしの差だけである。欠測の扱いが違って行数が変わっていれば、その時点で比較は成立しない。統計学Ⅰ 第11回でも同じ注意を扱っている。

加えて、第3節で述べたとおりREML で当てはめたモデルの AIC は、固定効果が違うモデルのあいだで比較できない。説明変数の組を選ぶ用途に使うなら最尤での当てはめが必要になる。今回参照した2つの資料には、どちらで推定したかの記載がない。

全部分集合の総当たりは、標本が小さいほど危うい

説明変数が11あれば部分集合は2000通りを超える。そのすべてに AIC を計算して最小を選ぶ手続きは、\(n=25\) という標本では過学習と偽陽性の危険が大きい。選ばれたモデルは「このデータに最もよく当てはまった組み合わせ」であって、母集団で最も良いモデルとは限らない。加えて、選択という操作そのものを経た後の係数や p 値は、選択を経ていない場合と同じようには読めない。

この点は外部からの批判ではなく、研究自身の作業記録が「n=25 で全探索は過学習・偽陽性リスク大」と記しているものである。同じ記録は、割合データに正規分布を仮定していたことも分布の誤指定として挙げ、ベータ回帰への移行を課題としている。

応答が割合であることを、正規分布は扱えていない

消失率・獲得率は \([0,1]\) に閉じた割合である。正規分布を仮定した線形混合モデルは、予測値が区間の外に出ることを構造的に防げない。標本の範囲内で当てはまりが良くても、定式化としては誤っている。この問題に対処するために応答分布をベータ分布へ置き換えたのがベータ回帰であり、同じ研究の後続の解析がその形を採っている。本ページが扱っているのは、その前段にあたる解析である。

近隣の定義と、クラスタ数の少なさ

空間誤差モデルの近隣行列も、混合モデルのクラスタ（学区）も、分析者が決める枠組みである。歩行コスト関数の選択、距離の切り方、どの単位を集団と見なすか。いずれもデータからは決まらない。同じ論点はICAR のページで、同じ研究群において隣接の定義を変えると空間自己相関の有無について異なる指標が出た事実として記録されている。

クラスタ数についても同様である。学区は5つ、地点は18ないし25である。この規模で分散成分を推定するとき、推定値は不安定になり、検定の検出力も低い。「学区差が検出されなかった」ことと「学区差がない」ことは、この標本では区別しにくい。

9出典とデータ

2つの解析を扱った。(A) は滋賀県米原市から岐阜県関ケ原町にかけての集落で行われた3世代方言語彙調査に基づく。原調査は2015年度に米原市「まいばら協働事業提案制度」のもとではびろネット（滋賀県米原市柏原）が実施し、調査総数は809人である。(B) は岐阜県揖斐川上流域18地点を対象とし、山田（1978, 1980）が報告した基礎語彙を用いている。

手法の文献

  • Akaike, H. (1974) A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6), 716–723. doi:10.1109/TAC.1974.1100705
  • Self, S. G. & Liang, K.-Y. (1987) Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82(398), 605–610. doi:10.1080/01621459.1987.10478472
  • Stram, D. O. & Lee, J. W. (1994) Variance components testing in the longitudinal mixed effects model. Biometrics, 50(4), 1171–1177.
  • Pinheiro, J. C. & Bates, D. M. (2000) Mixed-Effects Models in S and S-PLUS. Springer.

研究の文献

  • Bromham, L., Hua, X., Fitzpatrick, T. G. & Greenhill, S. J. (2015) Rate of language evolution is affected by population size. Proceedings of the National Academy of Sciences, 112(7), 2097–2102. doi:10.1073/pnas.1419704112
  • Aoki, K., Lehmann, L. & Feldman, M. W. (2011) Rates of cultural change and patterns of cultural accumulation in stochastic models of social transmission. Theoretical Population Biology, 79(4), 192–202.
  • Trudgill, P. (2011) Sociolinguistic Typology: Social Determinants of Linguistic Complexity. Oxford University Press.
  • 山田達也 (1978) 揖斐川上流域における基礎語彙について. 名古屋市立大学教養部紀要 人文社会研究, 22, 33–66.
  • 山田達也 (1980) 福井県今庄、池田方言及び滋賀県木之本方言の基礎語彙について ―― 徳山村方言に関連して. 名古屋市立大学教養部紀要 人文社会研究, 24, 23–34.
  • Aiken, R. (1977) Wilderness areas in Scotland. Unpublished Ph.D. thesis, University of Aberdeen.　／　Langmuir, E. (1984) Mountaincraft and leadership. The Scottish Sports Council/MLTB, Cordee.　（歩行コスト関数の出典として (B) が引くもの。(B) の本文は「Aitken」、参考文献欄は「Aiken」と表記が揺れている）

この手法がどの観点に属するか

線形混合モデルは、パラメータを固定された未知の値と見て、ランダム効果を積分消去した周辺尤度を最大にする。返すのは点推定と信頼区間であり、確率がかかっているのは区間のほうである。7つの観点を式で読むの観点2（頻度論）に属する。一方、モデルを選ぶ道具として用いる AIC は観点5（情報理論）の量であり、同じ観点には正規化圧縮距離がある。推定は頻度論、選択は情報理論という組み合わせで動いている点が、この手法の位置づけである。同じデータを観点3（ベイズ）で解き直した場合についてはベイズ推定とICAR 空間ランダム効果を参照。

このページは、小野原が実際の研究で使っている手法を、同じ分野の読み手に向けて整理したものです。数値はすべて発表資料に印刷された値および解析の作業記録から引用し、推定値や概算は記載していません。確認できなかった書誌情報は記載していません。
手法の位置づけの全体像は 7つの観点を式で読む、どの授業回がどの観点に当たるかは 手法と数学の地図 にあります。