1何のための手法か

応答が割合であるとき、予測値が区間の外に出ることを構造的に防ぐ。

応答 \(y_s\) が「集落 \(s\) で失われた語彙の割合」のように \((0,1)\) に閉じているとする。素朴に正規分布を仮定した線形回帰 \(y_s = \mathbf{x}_s^\top \boldsymbol{\beta} + \varepsilon_s\) を当てると、線形予測子は実数全体を動くので、予測値が 0 を下回ったり 1 を超えたりする。割合として解釈できない値が出る時点で、モデルは応答の定義域と矛盾している。

回避策として持ち出されがちな二つの手が、どちらも適切でない。ひとつは閾値で二値化して（たとえば 0.5 以上を 1 とする）ロジスティック回帰に持ち込む方法だが、これは連続量として測れているものを捨てている。もうひとつはロジスティック回帰をそのまま当てる方法だが、二項分布は 0 か 1 かの試行を数える分布であって、連続な割合のための分布ではない。

ベータ回帰は、\((0,1)\) 上の連続分布であるベータ分布を応答分布に据え、その平均をリンク関数で線形予測子につなぐ。定義域の問題は、分布の選択そのものによって解消される（Ferrari & Cribari-Neto 2004）。

2起源 ―― どの分野の誰が作ったか

割合データに回帰を当てる枠組みが欠けていた、という不足から生まれた。

Silvia Ferrari と Francisco Cribari-Neto が 2004 年に Journal of Applied Statistics に発表した定式化である。ブラジルの統計学者による仕事で、動機は明快だった。割合・比率のデータは分野を問わず現れるのに、それを応答にとる回帰の標準的な枠組みが整っていなかった。上に挙げた二つの回避策の不備が、そのまま出発点になっている。

発表後の普及は速く、割合データが扱われる領域に広く浸透した。ただし注意すべきなのは、この論文が新しい分布や新しい推定原理を発明したわけではないという点である。ベータ分布は古くから知られており、一般化線形モデルの枠組みも 1972 年に整っていた。Ferrari と Cribari-Neto の貢献は、既にある部品を組み合わせ、回帰モデルとして使える形に仕上げたことにある。次節の系譜がそれを示す。

3数学的にどこから来ているか

指数型分布族 → 一般化線形モデル → ベータ回帰、という三段の接続である。

出発点は、ばらばらに見える分布が共通の形式を持つという発見にある。Fisher（1934）らが独立に定式化した指数型分布族は、密度を次の形に書けるものの総称である。

指数型分布族の一般形
\[ f(y \mid \theta) \;=\; h(y)\,\exp\bigl[\,\eta(\theta)\,T(y) - A(\theta)\,\bigr] \]

正規分布・二項分布・ポアソン分布、そしてベータ分布がこの族に属する。分布の違いは \(h,\eta,T,A\) の違いに還元される。

Nelder & Wedderburn（1972）は、この族のどの分布を応答に選んでも同じ手続きで回帰ができる枠組みを与えた。線形予測子 \(\mathbf{x}^\top\boldsymbol{\beta}\) と、応答分布の平均 \(\mu\) を、リンク関数 \(g\) で結ぶ。

一般化線形モデルの接続
\[ g(\mu_i) \;=\; \mathbf{x}_i^\top \boldsymbol{\beta} \]

正規分布に恒等リンクを選べば線形回帰、二項分布に logit を選べばロジスティック回帰、ポアソン分布に log を選べばポアソン回帰。ベータ分布に logit を選べばベータ回帰になる。

ベータ回帰は、この枠組みの空席を埋めた仕事である

上の対応表を見ると、ベータ分布の行だけが長く空いていたことが分かる。Ferrari & Cribari-Neto（2004）は、ベータ分布を平均と精度で再パラメータ化して一般化線形モデルの形に乗せ、その空席を埋めた。だからこの手法の系譜は、統計学の内部で閉じている。他分野から借りてきたものではない。

4定義と仮定

平均と精度で書き直したベータ分布に、logit リンクで予測子を接続する。

ベータ分布は通常 \(\mathrm{Beta}(a,b)\) と2つの形状パラメータで書かれるが、回帰に使うには平均を直接扱えたほうがよい。\(\mu = a/(a+b)\)、\(\phi = a+b\) と置き直すと次の形になる。

平均 \(\mu\) と精度 \(\phi\) による再パラメータ化
\[ f(y \mid \mu,\phi) \;=\; \frac{\Gamma(\phi)}{\Gamma(\mu\phi)\,\Gamma\bigl((1-\mu)\phi\bigr)}\; y^{\mu\phi-1}\,(1-y)^{(1-\mu)\phi-1} \]

\(\mu \in (0,1)\) … 分布の平均。これを線形予測子につなぐ
\(\phi > 0\) … 精度。大きいほど平均のまわりに集中する
分散 … \(\mathrm{Var}[y] = \mu(1-\mu)/(1+\phi)\)。平均に依存して変わる。等分散は仮定していない

この \(\mu\) を logit で線形予測子に接続する。米原調査の主モデルは、集落 \(s\) の語彙消失率について次の形をとる。

この研究の主モデル
\[ \operatorname{logit}\bigl(\mathrm{loss}_s\bigr) \;=\; \beta_0 + \beta_1\,\mathrm{growth}_s + \beta_2\,\log N^{\text{ms}}_s + \beta_3\,\log N^{\text{gp}}_s + u_{\text{district}[s]} \]

\(\mathrm{growth}_s\) … 人口成長率（2010年人口 / 2005年人口）を標準化したもの
\(\log N^{\text{ms}}_s,\ \log N^{\text{gp}}_s\) … 中学生・高年層それぞれの調査人数の対数。標本努力の統制項
\(u_{\text{district}[s]}\) … 中学校区のランダム切片（5区）

分散が平均に依存することは、仮定ではなく帰結である

上の分散の式が示すとおり、ベータ分布では平均が 0 や 1 に近づくほど分散が小さくなる。割合データで実際に観察される「端に寄るほどばらつきが小さい」という性質が、分布の側から自動的に出てくる。正規分布を当てた場合にこれを表現するには不等分散の追加構造が要るが、ベータ回帰では分布を選んだ時点で織り込まれている。

5端点の扱い ―― 0 と 1 がちょうど出たとき

ベータ分布は開区間の分布である。端点は密度を持たない。

定義域は \((0,1)\) であって \([0,1]\) ではない。\(y=0\) や \(y=1\) がデータに含まれると尤度が定義できず、当てはめは失敗する。割合データでは端点は珍しくないので、実務では何らかの処理が要る。

標準的な対処は Smithson & Verkuilen（2006）が与えた線形圧縮で、標本サイズ \(n\) に対して次のように区間の内側へ押し込む。

端点を内側へ寄せる変換（この研究では使っていない）
\[ y' \;=\; \frac{y\,(n-1) + 1/2}{n} \]

\(n\) が大きいほど移動量は小さい。ただしこれはデータを動かす操作であり、端点が多いデータでは結果が変換に依存しうる。端点そのものに意味がある場合は、0 と 1 に確率質量を置く zero-one inflated の定式化を検討すべき場面である。

この研究では変換を必要としなかった

米原調査の消失率・獲得率は元データで 0〜100 のパーセント表記であり、100 で割って \([0,1]\) に直している。その結果として得られた 25 自治会の値は端点を1つも含まなかったため、圧縮変換を適用していない。これは手法の一般的性質ではなく、このデータがたまたま満たした条件である。端点を含むデータに同じモデルを当てる場合には、上の処理か別の定式化が必要になる。

6起源は最尤、実装はベイズ

この区別は、結果の読み方を変える。

Ferrari & Cribari-Neto（2004）の定式化は最尤推定である。一般化線形モデルの枠組み自体、反復重み付き最小二乗で最尤推定量を得る手続きとして整えられた（Nelder & Wedderburn 1972）。この系譜において、パラメータは固定された未知の値であり、標本のほうが確率変動する。

一方この研究の確定版モデルは、brms を通じた Stan によるベイズ推定で当てはめている。パラメータに事前分布を置き、事後分布を MCMC で得る。報告される区間は信用区間であって信頼区間ではなく、\(P(\beta > 0)\) はパラメータについての確率として読める。

観点2に置いているのは、定式化の系譜による

このページを観点2（頻度論）に置いているのは、ベータ回帰という手法が指数型分布族と一般化線形モデルという最尤の系譜から生まれたものだからである。推定の枠組みそのものはベイズである。同じ尤度に事前分布を掛けて事後分布を得ているのであって、標本のばらつきで判断しているのではない。

この区別は名目上のものではない。後掲の表で base の係数に添えた 0.993 は、事後分布において \(\beta>0\) である確率であり、p 値ではない。頻度論の枠組みで同じ数字を出すことはできない。逆に、頻度論のベースラインとして併走させた betareg の結果には p 値が付いており、こちらはパラメータについての確率を与えない。7つの観点を式で読むで述べているとおり、尤度はベイズでも使われる。最尤は推定原理のひとつであって、頻度論と同義ではない。

7この研究での使い方

人口成長率が語彙消失率を押し上げるか、を標本努力と空間構造を統制して問う。

  1. 応答の構成　集落 \(s\) の項目 \(j\) について、高年層（60歳以上）の使用割合を \(g_{sj}\)、中学生の使用割合を \(m_{sj}\) とし、\(\mathrm{loss}_s = (1/J)\sum_j \max(g_{sj}-m_{sj},\,0)\) で消失率を定める。獲得率 \(\mathrm{gain}_s\) は差の向きを逆にしたもの。
  2. 頻度論のベースライン　betareg で logit リンクのベータ回帰を当て、標準化した人口規模・人口成長率・面積を予測子にする。小標本（n=25）で先に素描を取っておく段階。
  3. 階層ベイズへの移行　brms ＋ cmdstanr/Stan で、中学校区のランダム切片を入れた階層ベータ回帰を当てる。標本努力の統制項として中学生・高年層それぞれの調査人数の対数を入れる。
  4. 事前分布　傾きに \(\mathrm{Normal}(0,1)\)、切片に \(\mathrm{Normal}(0,2)\)、精度 \(\phi\) に \(\mathrm{Gamma}(0.1,0.1)\) の弱情報事前。
  5. サンプリング　4 鎖・4,000 反復（うち 2,000 をウォームアップ）、adapt_delta = 0.97、max_treedepth = 12。採用モデルは \(\hat{R}=1.00\)、bulk-ESS > 3,000、発散遷移 1 以下。
  6. 空間の確認　3 km 以内を隣接とする二値の隣接行列を作り、ICAR 空間ランダム効果を入れたモデルと入れないモデルを LOO 近似交差検証で比較する。
モデル人口成長率の係数95% 信用区間\(P(\beta>0)\)
base（標本努力を統制）+0.231[+0.059, +0.397]0.993
+ covariates（面積・標高・建物用地）+0.229[+0.032, +0.422]0.987
+ spatial（ICAR、3 km 隣接）+0.228[+0.045, +0.402]0.992

階層ベイズベータ回帰（n=25、中学校区ランダム切片）。係数は logit スケール、予測子は標準化済み。数値はリポジトリの loss_growth_coefs.csv と作業ログに記録された実行結果の引用。

頻度論ベースラインとの対照

betareg による当てはめでは、人口成長率の標準化係数が +0.2522（標準誤差 0.0582、\(z=4.334\)、\(p=1.5\times10^{-5}\)）。同じモデルに含めた人口規模は \(-0.1274\)（\(p=0.101\)）、面積は \(+0.0295\)（\(p=0.690\)）で、いずれも 5% 水準に達しない。獲得率を応答にすると人口成長率の係数は \(-0.139\)（\(p=0.267\)）となり、符号が反転したうえで有意でなくなる。頑健に出ているのは消失率の側だけである。

ベイズの階層モデルで得た \(+0.231\) と、頻度論の \(+0.2522\) は近い。小標本で事前分布が弱情報であれば両者が接近するのは想定どおりで、事前分布が結果を作っているわけではないことの確認になる。

空間ランダム効果は採用していない

LOO 比較では、ICAR を入れたモデルは入れないモデルよりわずかに悪い（elpd 差 \(-0.7 \pm 0.3\)）。加えて n=25 に対して 25 個の空間ランダム効果を置くのは過剰であり、サンプリングにも treedepth の問題が出た。採用モデルは非空間の階層ベータである。空間項は「必要かどうかをデータで確かめる」ために置いたものであって、補正として入れたものではない。

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

分布の選択が解決するのは定義域の問題だけである。

応答の定義域が正しくなっても、因果は出ない

ベータ回帰が保証するのは、予測値が \((0,1)\) に収まることと、分散が平均に依存する形で扱われることだけである。人口成長率と消失率の関連が \(+0.23\) で安定していることは、人口成長が消失を引き起こしている証拠ではない。観察データであり、処置の割り当てがない。ここで統制した標本努力・面積・標高・建物用地・空間構造は、あくまで測定できた交絡候補にすぎない。

小標本であることは、モデルの選択では解消しない

消失率が計算できたのは 25 自治会である。中学校区は 5 つしかなく、ランダム切片の分散を推定する材料は乏しい。弱情報事前は不確実性を適正に表現するための装置であって、標本を増やす代わりにはならない。実際、1 集落を除く leave-one-out の感度分析では、高年層の回答者が 1 名しかいない集落を除いた場合に係数が \(+0.12\)（\(p=0.25\)）まで落ちる。この関連は、最も情報の乏しい観測を残したときに最も強く出る。

応答の作り方そのものが結果を左右する

\(\mathrm{loss}_s\) は世代間の使用割合の差を項目について平均したものだが、差を \(\max(\cdot,0)\) で切り落としている。消失と獲得を別々の量として扱うための定義であり、この切り落としは負の側の情報を捨てる操作でもある。また各集落の割合は少数の回答者から推定されており、その不確実性は集計の時点で消えている。この点は、項目×集落×世代の使用者数を直接モデル化するカウントレベルのベータ二項モデルで別途確認されており、そちらでは世代×成長の交互作用が \(-0.17\)、95% 信用区間 \([-0.50,\,0.18]\) と 0 をまたぐ。集計を前提にしたときにだけ関連が明確に出ることは、読み手に伝えるべき性質である。

9出典とデータ

滋賀県米原市から岐阜県関ケ原町にかけての5中学校区を対象とした、3世代方言語彙調査の二次解析である。調査票は 40 の質問概念で 106 の方言語彙変異を引き出しており、関西・近江系 45、東京・標準系 31、美濃・尾張系 26 に若者語・新方言 4 項目を加えて事前分類されている。原調査は 2015 年にはびろネットが実施し、本研究はその匿名化データを解析したもの。解析は全工程がスクリプト化されている。

手法の文献

  • Ferrari, S. & Cribari-Neto, F. (2004) Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7), 799–815. doi:10.1080/0266476042000214501
  • Nelder, J. A. & Wedderburn, R. W. M. (1972) Generalized linear models. Journal of the Royal Statistical Society, Series A, 135(3), 370–384. doi:10.2307/2344614
  • Smithson, M. & Verkuilen, J. (2006) A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables. Psychological Methods, 11(1), 54–71. doi:10.1037/1082-989X.11.1.54
  • Bürkner, P.-C. (2017) brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1–28. doi:10.18637/jss.v080.i01

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

ベータ回帰は、指数型分布族という共通形式と、一般化線形モデルというリンク関数の枠組みから導かれる。この系譜は最尤推定を中心に据えた観点2（頻度論）のものだが、7つの観点を式で読むで述べているとおり、尤度はベイズでも使われる。この研究の実装は事前分布を置いた事後推定であり、同じ尤度が観点3（ベイズ）の側から読まれている。ひとつの式が複数の観点で読めることの実例にあたる。

このページは、小野原が実際の研究で使っている手法を、同じ分野の読み手に向けて整理したものです。数値はすべて解析スクリプトの出力と作業記録から引用し、推定値や概算は記載していません。
手法の位置づけの全体像は 7つの観点を式で読む、どの授業回がどの観点に当たるかは 手法と数学の地図 にあります。