なんというか、たまにSNSとかwebの記事なんかをみると、大企業のデータサイエンティストなる華やかな人々がビジネスへの貢献について華やかに語っておられて、彼我のちがいにちょっと目眩がすることがある。ああいう人たちってふだんなに食ってんだろうか。ステーキとかかな。なんか知らん横文字の料理とかかな。すくなくとも私みたいに冷やご飯にのりたま振って流しのまえで立ち食いしたりはしないんだろうな。知らんけど。

　まあとにかく、きっと皆さん私の知らないことをたくさん知っているので、たとえば目的変数がなにかの割合であるようなデータを渡されて回帰分析する羽目になったときも(突然に卑近な話になる)、きっとなにか私の知らない先端的な手法を使うのだろうなあと思う(いやいや、アシスタントに丸投げするんでしょうね)。いっぽう私はそのたびにこうジクジクと悩むわけです。毎度毎度binomial-logitでGLMしてていいの？ たまにはなんかこう気の利いた誤差分布とかないわけ？ 元の観察数がわかんなかったらロジット変換してOLSでいいの？ なんかもっとパンクな手法はないわけ？ とかなんとか。あーあ、残念な人生だ。

Cribari-Neto, F., Zeleis, A. (2010) Beta Regression in R. Journal of Statistical Software, 34(2), 1–24.

　仕事の都合でざっと目を通した奴。実際に読んだのは上記文献ではなく、その改訂版らしき R のbetaregパッケージのvignetteである。ちょっと都合があって、betaregを実戦投入しようかと思ったことがあったので。

　本文のメモの前におさらいしておくと、(第1種)ベータ分布の密度関数は $$ f(y; p,q) = \frac{1}{B(p, q)} y^{p-1} (1-y)^{q-1} $$ である。ベータ関数\(B(p, q)\)というのが出てくるけど、要は積分 $$ B(p,q) = \int^1_0 y^{p-1} (1-y)^{q-1} dy$$ で割りましたという話である。なおガンマ関数で書き直すならば $$ B(p,q) = \frac{\Gamma(p)\Gamma(q)}{\Gamma(p+q)}$$ である。ガンマ関数とは… なんだかしらんが階乗みたいなもんであろう。

　いわく。
　ベータ回帰とはFerrari & Cribari-Neto (2004)が提案したもので、\((0,1)\)に落ちる目的変数の回帰モデルである。\((a,b)\)に落ちる変数の場合は\( (y-a)/(b-a) \)と変換してから回帰すれば良い。\([0,1]\)に落ちる変数の場合は \( (y(n-1)+0.5)/n \) と変換してから回帰すれば良い。[←へー]
　ベータ分布の密度関数は$$ f(y; p,q) = \frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)} y^{p-1} (1-y)^{q-1} $$これをそのまま使ってもいいのだが、Ferrariたちは \(\mu = p/(p+q), \phi = p+q\)と置いて$$ f(y; \mu,\phi) = \frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1-\mu)\phi)} y^{\mu\phi-1} (1-y)^{(1-\mu)\phi-1} $$としている。こうすると平均は\(\mu\), 分散が\(\mu(1-\mu)/(1+\phi)\)となり、\(\phi\)は精度だということになる。以下ではこの分布を\(\mathcal{B}(\mu, \phi)\)と書く。
　ベータ回帰とはこういうモデルだ。無作為標本\(y_1, \ldots, y_n\)について \(y_i \sim \mathcal{B}(\mu_i, \phi)\)と仮定する。で、\((0,1)\)を実数にマップするリンク関数を\(g(\cdot)\)として$$ g(\mu_i) = x_i^T \beta $$とする (簡略のため \(x_{i1} = 1\)として切片項も含める)。
　リンク関数を噛ますのは\(\mu_i\)が\((0,1)\)でしか動かないからだが、ユーザに柔軟性を与えるためでもある。ロジットリンク、プロビットリンク、clog-logリンク, log-logリンク, コーシーリンク \(g(\mu) = tan\{\pi (\mu -0.5)\}\)などが使える。
　ベータ回帰モデルでは\(Var(y_i)\)が\(\mu\)の関数になる、つまり自然に分散不均一になるというのがポイント。
　最尤推定できる。対数尤度関数は…[メモは省略するけど、素直に密度関数の対数をとったのが個体尤度になる]

　variable dispersionベータ回帰という提案もある。精度パラメータも定数ではなくして$$ g_1(\mu_i) = x_i^T \beta$$ $$ g_2(\phi_i) = z_i^t \gamma$$ とモデル化する。さらに、この2本の式を非線形にするという提案もある(betaregパッケージは未対応だけど)。

　ところで、このモデルでは残差 \(y_i – \hat{\mu}_i \)はもともと不均一なので、それを $$ \hat{Var}(y_i) = \frac{ \hat{\mu_i}(1-\hat{\mu_i}) }{ 1+\hat{\phi}_i }$$ の平方根で割ったのを使う。これをピアソン残差とか標準化通常残差という。ほかにも標準化加重残差2という提案があって…[面倒くさいので略]

　このようにベータ回帰は、率や割合やジニ係数のような\((0,1)\)の連続変数のモデル化手法なのだが、試行数に占める成功割合のモデル化に使っても良い(試行数が十分に大きければ)。この場合は二項変数のGLMモデルのような使い方になるが、もっと柔軟である。dispersionを固定すればquasi-binomialモデルに相当するモデルになるし(しかしquasi-binomialモデルと違って完全にパラメトリック)、dispersionが変動するモデル化も容易である。[なるほど、成功数/試行数型のデータであっても、GLMより過分散を扱いやすいというメリットがあるってことね]

というわけで、betareg::betareg()のご紹介。主な引数は次のとおり。

• formula. y ~ x1 + x2というふうに書くと\(\phi\)は定数になる。y ~ x1 + x2 | z1 + z2という風に書くと\(g_2(\phi) = z_i^T \gamma\)となる。なお、デフォルトでは\(g_2(\cdot) = \log(\cdot)\)だから、y ~ x1 + x2 | 1と書くと\(\phi\)じゃなくて\(\log(\phi)\)を推定することになる。
• data, subset, na.action, weight, offset [書いてないけどlm()とかと同じであろう]
• link. \(g_1(\cdot)\)を指定する。デフォルトは”logit”.　ほかにいろいろ使える。
• link.phi. \(g_2(\cdot)\)を指定する。”identity”, “log”, “sqrt”とかが使える。
• control = betareg.control(...). optim()に渡される。

　…ここからは事例を5つ紹介。メモは省略するけど、最後の例は、時系列回帰で係数が時変しないかどうかしらべる strucchange::gefp() という関数の引数としてbetaregを与えるという話であった。へええ。
　最後にちらっと触れられていたけど、ベータARMAモデルという提案もあるのだそうだ。へえええ。

　というわけで、ざーっと読んだだけだけど、ベータ回帰モデルの主旨がわかったのでよしととしよう。誤差項が正規分布でも二項分布でもなくてベータ分布だと考えるってことね。だからほっといても分散不均一になり、かつ分散のモデル化も柔軟になるわけね。でも期待値に対して適切なリンク関数を考えないといけないという点はGLMと同じで、つまり別に仕事が楽になるわけじゃないのね。ちょっとがっくり。