線形混合モデルを理解したい

この記事は

R Advent Calendar 2021の12日目の記事になるらしいです。
昨日の記事はしょこさんのLasso回帰スクラッチ実装です。
これはガチガチにすごいので正則化回帰完全理解を試みる勢は必読。
そして明日もしょこさんがロジスティック回帰をスクラッチ実装するらしいです。スクラッチ実装はいいぞ。
僕はやってませんが。

• 一般化線形混合モデルをRで実施する方法を書いています。
• ランダム効果の推定について2021å¹´12月11日現在で僕が追いきれている部分まで記述しています。
• 混合モデルの拡張(正則化/非線形モデル)についても今後の調査・リサーチのためにメモしています。

お久しぶりです

5月以来ブログの更新をしていませんでした。
今年1年を振り返ることはまた改めて記事にしようかなと思いますが、
主に以下のようなことがありました。

• 転職活動をしていたが疲れた。
• そんなことをしていたら1年育てた後輩が引っこ抜かれた
• 採用活動をしている
• R以外を触ることが増えてきた(最近はパワーポイントとPythonとかRustとかを触っている)。

今日書くこと

2019年のアドベントカレンダーでは正規線形モデルから一般化線形混合モデルまでの変遷について、
適当に書いていました。
今年はその中でも一般化線形混合モデルを細かく書こうと思います。

なぜ？

面白いからです。

環境

今回のスクリプトは以下の環境で実行しています。スクリプトは以下。

github.com

version

# platform       x86_64-pc-linux-gnu         
# arch           x86_64                      
# os             linux-gnu                   
# system         x86_64, linux-gnu           
# status                                     
# major          4                           
# minor          1.1                         
# year           2021                        
# month          08                          
# day            10                          
# svn rev        80725                       
# language       R                           
# version.string R version 4.1.1 (2021-08-10)
# nickname       Kick Things 

R言語の良いところはOSへの依存がかなり低いことだと思います。
最近はローカルでも仮想環境を立てて環境依存を回避する動きも非常に重要ではありますが、
データ分析からプログラミングを始め、Rが母国語でRが好きな人は、
まずローカルの環境で一通りコードが書けることは重要なステップだと思われます。 ちなみに、フォルダの管理については過去に記事を書いたので、
軽く参考にしてもらえるならいいかなと思います。

データ分析をちゃんと管理しよう with R 【フォルダ構成編】

データ

palmerpenguinsで読み込めるペンギンデータを用います。
測定年、ペンギンの種類、居住地など、混合モデルを適用するにも適した良いデータです。
データ分析初心者！という人にも、欠損値やカテゴリ変数といった基本的な概念がしっかり含まれていて、
とてもいいなあと思います。

library(palmerpengins)
summary(penguins)

# species          island    bill_length_mm  bill_depth_mm  
# Adelie   :152   Biscoe   :168   Min.   :32.10   Min.   :13.10  
# Chinstrap: 68   Dream    :124   1st Qu.:39.23   1st Qu.:15.60  
# Gentoo   :124   Torgersen: 52   Median :44.45   Median :17.30  
# Mean   :43.92   Mean   :17.15  
# 3rd Qu.:48.50   3rd Qu.:18.70  
# Max.   :59.60   Max.   :21.50  
# NA's   :2       NA's   :2      
#
# flipper_length_mm  body_mass_g       sex           year     
# Min.   :172.0     Min.   :2700   female:165   Min.   :2007  
# 1st Qu.:190.0     1st Qu.:3550   male  :168   1st Qu.:2007  
# Median :197.0     Median :4050   NA's  : 11   Median :2008  
#  Mean   :200.9     Mean   :4202                Mean   :2008  
#  3rd Qu.:213.0     3rd Qu.:4750                3rd Qu.:2009  
#  Max.   :231.0     Max.   :6300                Max.   :2009  
#  NA's   :2         NA's   :2       

実装

多分「理論とかどうでもいいのでコード知りたい！」という人もいるんだろうなあと思います。
以降に僕なりに理論を噛み砕いているのでわからなくなったらまた読んでください。
Rで線形混合モデルを実装できるパッケージはいくつかあります。細かい点で違いはありますが、
この記事では、僕が一番好きな lme4 パッケージでデモンストレーションをします。
大変すばらしいことに、Rで回帰分析を知っている人であれば問題なく実行ができると思います。
知る必要のあるコードは「ランダム効果」を推定したい変数に関する部分です。
具体的には、ランダム効果を推定したい変数を(変数名 | グループ変数名)の形で書きます。
切片を加えたい場合は1を含みます。
今回はペンギンのくちばしの長さ(bill_length_mm)をペンギンの翼の長さ(flipper_length_mm)で回帰することを考えます。
ランダム効果については、ペンギンの種類（species）によって切片と翼の長さが異なるだろうという仮説に基づき(1 + flipper_length_mm | species)という形で記述します。

# define data
# to simple, we omit NA rows.
usedata <- penguins %>% 
  na.omit()

# glmm model run (for comparison)
glmm_model <- lme4::lmer(bill_length_mm ~ flipper_length_mm + (1 + flipper_length_mm|species),
                         data = usedata)

モデリングの部分はほんの2行、1関数で済みます。とても簡単。
結果を見てみます。

summary(glmm_model)
# Linear mixed model fit by REML ['lmerMod']
# Formula: bill_length_mm ~ flipper_length_mm + (1 + flipper_length_mm |      species)
# Data: usedata
# 
# REML criterion at convergence: 1596.2
# 
# Scaled residuals: 
#   Min      1Q  Median      3Q     Max 
# -2.6469 -0.6595  0.0297  0.6924  4.9957 
# 
# Random effects:
#   Groups   Name              Variance  Std.Dev. Corr 
# species  (Intercept)       2.6567656 1.62996       
# flipper_length_mm 0.0009521 0.03086  -1.00
# Residual                   6.7064017 2.58967       
# Number of obs: 333, groups:  species, 3
# 
# Fixed effects:
#   Estimate Std. Error t value
# (Intercept)        1.25267    4.34480   0.288
# flipper_length_mm  0.21789    0.02767   7.873
# 
# Correlation of Fixed Effects:
#   (Intr)
# flppr_lngt_ -0.886
# optimizer (nloptwrap) convergence code: 0 (OK)
# unable to evaluate scaled gradient
# Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

今回はいくつかwarningsがでてしまいましたが、この部分については一旦脇に置きます*1。
Fixed effectsにある部分は通常の回帰係数の解釈が可能です。Random effectsの部分はVariance(分散)の推定値のみが算出されています。係数をはじき出すにはranef()関数を使うことで可能です。

ranef(glmm_model)
# $species
#            (Intercept) flipper_length_mm
# Adelie      1.4836890      -0.028092178
# Chinstrap  -1.7997485       0.034099870
# Gentoo      0.3160595      -0.006007692

切片と翼の長さに関しては、固定効果とランダム効果の合計で「ペンギンの種類別の係数」が算出できます。
モデリング自体は非常に簡単です。ただモデリングそのものの過程は非常に難しいのです。
例えば今回は「ペンギンの種類によって、各部位の大きさ/長さは異なるだろう」という仮説をおいて、どうやらそれはある程度の妥当性を持って評価できそうだ、ということが分かっていますが、
現実問題でこういった仮説を妥当性を持つ形で表現できない場合、一般化線形混合モデルは「複雑過ぎる」モデルになってしまうでしょう。これは良くない。
ですのでちゃんと仮説を立てて分析しましょう。とは言いつつ使わないと使えるようにはならないので、
一つのやり口としては「ぼくのかんがえた一番複雑なモデル」を定義して、
そこから「もっとシンプルに解けないか」と下っていくというスタイルも(邪道ですが)いいんじゃないでしょうか。
R Advent Calendar要素はこれで終わりなので以降は蛇足です。

一般化線形混合モデルとは

一般化線形混合モデルは広く言えば線形回帰分析の1ファミリーです。
ただ、通常の回帰分析よりも要請する仮定が緩いことが特徴です。

通常の回帰モデル

例えば日本在住の人たちからいい感じに標本を抽出して、
所得を年齢や性別、職業など(以後長いので「個人属性」といいます)の説明変数(${X}$で表記)で回帰するシンプルな回帰モデルを考えます。 $$Income_i = \sum_j^{J}\beta_j X_{i,j} + \epsilon_i$$ カッコよく行列化して $${Income} = {X}{\beta} + {\epsilon}$$ を考えます。
統計モデルにはいろいろと前提を設けます。上記のモデルは全国で評価することを前提にしています。つまり、日本全体における 所得と個人属性との間にある「全国平均」の関係性をモデル化しています。この結果も無意味ではありませんが、
例えば就業構造基本調査のデータなどを見ていくと「居住する都道府県別に個人属性の回帰係数は異なるのではないか？」
という仮説がでてくるかもしれません。
ただ、上記のモデルでは「都道府県別に回帰係数が変わる」ようなモデルを作ることは難しいです。どのように組み直せば良いでしょうか。

線形混合モデル

都道府県別に回帰係数が異なるということは、例えば年齢という変数1つに対して、
47個の回帰係数を得るようなモデルを考えられそうです。
そして、その47個の回帰係数が標本数だけ得られれば良いので、例えば以下のような形の書き換えができるかもしれません。

$${Income} = {X}{\beta} + {X}{b} + {\epsilon}$$

ここで新しく入ってきた${b}$は行は人の数、列が都道府県の数だけあるベクトルです。
これによって都道府県別に職業の影響を評価するというモデルが、(一般化)線形混合モデルの基本的な構造です。
都道府県別でなくても、他の例として、店舗の時系列データを使った売上予測モデルに線形混合モデルを使う場合は店舗別に変動しうる説明変数を${b}$にかかる変数行列(計画行列)${Z}$として入力する、すなわち $${Income} = {X}{\beta} + {Z}{b} + {\epsilon}$$ の形式で表現することが多いです。

固定効果とランダム効果

回帰係数について、${\beta}$を「固定効果(fixed effect)」、${b}$を「変量効果」あるいは「ランダム効果(random effect)」と呼びます。
なぜこれらのように呼ぶのかについてはWikipedia先生でも簡便な記載がありますが、$\beta$は定数である一方、${b}$は確率変数として扱うことが理由です。

の推定

さて、${b}$は回帰係数の行列である、つまりパラメータであるので、
データの構造から推定される必要があります。
一方でこれが「ランダム効果」と言われるように、こちらは確率変数としての扱いをする必要があります。
この扱いによってパラメータ推定が確率的に可能になります。
理論と実装はギャップがあるものなので、この辺はBates, et. al(2015)をベースに整理していきます。

定式化

以下のような線形混合モデルを考えます。
$$Y = {X}{\beta} + {Z}{b} + {\epsilon}$$ 目的変数$Y$が特定の分布(ここでは正規分布ということにします)に従うことを仮定します。

$$ Y \sim N(\mu, \sigma ^2)$$

さて、上記の仮定のもとで$\mu$を推定すると、通常の線形回帰のパラメータ推定と同じ定式化になります。
今回は混合モデルで、${b}$を推定したいので、$b$で条件付けられた$Y$の分布として考えようと思います。すなわち

$$(Y|B = b) \sim N(X \beta + Z b + \epsilon, \sigma ^2 W^{-1})$$

これはBates et. al(2015)の式(2)に相当する変換です。ランダム効果${b}$によって条件付けられた我々は$Y$が正規分布に従うとうれしいというモデルを考えます。その時の平均$\mu$が回帰式で求められます。${W}$は対角行列で、分散の重みを意味します。今回はあんまりでてこないので、「分散を制御する重みだなあ」という感じで見てください。
今回の関心は${b}$を推定することにあるのですが、これを1つ1つ丁寧に推定するのは大変に骨が折れます。都道府県別のパラメータを推定する場合、変数1つにつき47個必要で、それが標本数分あるのですから、1つ1つ推定しようとすると容易に計画行列がランク落ちしてしまいます。
そのためには推定するパラメータを少なくできると良いです。例えば${b}$が何らかの分布に従っていることにできると、その分布パラメータを決定するだけで良くなるので、うれしいです。
例えば多変量正規分布に従ってくれるととても扱いやすくなります。とりあえず${b}$を確率変数と捉えて、以下のように仮定するといい気がします。
$$B\sim{N}({0}, {\Sigma})$$ ここで、$\Sigma$は分散共分散行列です。これは目的変数の分布ではなく、ランダム効果パラメータの分布に対する仮定なので、目的変数の分布に関しては一般化線形モデルの範疇で(つまり指数型分布族に含まれる分布であれば)拡張が可能です。実際lme4では、glmer関数でロジスティック回帰とポアソン回帰が可能です。
さて、これで${\Sigma}$を推定することができればランダム効果${b}$が確率的に生成できるらしい、というところまで落ち着きました。
ここで素敵なことは、${b}$の成分1つ1つを推定しなくても、多変量正規分布の分散成分が推定できればよいという話にすり替えられたことにあります。
ここから推定に入るのですが、コンピュータがこれを計算できるようにするには、この${\Sigma}$を分解してわちゃわちゃする必要があるようなのですが、
僕がここから理解できていないので理解できたら追記します。締切優先です*2。

効果の有意性について

隠れた秘密ですが、lme4では回帰係数の有意性を算出できません。というか算出「しない」というポリシーで動いています。理由はシンプルに「計算が大変である」ということです。
有意性について考えます。
一方でlmerTestというライブラリは、この自由度の計算を近似計算で実現しています。
ただ、こちらは固定効果に対する有意性だけで、混合効果の有意性は評価できなかった、はずです。

統計的学習への拡張

ここからは統計的学習の文脈での現時点での僕の関心の話です。
詳細は完全に理解してから記事を書こうと思いますが、モチベーションとしては 混合モデルの正則化や、決定木をはじめとする非線形モデルへの拡張を目論みたくなるところにあります。
一般化線形モデルの範疇であれば、正則化を実施することが可能であることはTibshiraniらによって知られています。
Rのアドベントカレンダーなので、Rの話をしますが、例えばLasso/Ridge/ElasticNetについては、glmnetパッケージで容易に実装することができます。
一方で、混合モデルのパラメータは固定効果と変量効果に分かれており、glmnetでは難しい気がしてきます。
とりわけ変量効果については、変量効果として必要な変数がなんなのかをヒューリスティックに決定したくない場合がままあるかなと思います。つまりは変数選択の需要が相対的に高いというわけです。こうした正則化つき一般化線形混合モデルについてもglmmLassoパッケージを利用することで可能です。
決定木系で混合モデルを行う場合、最近はGPBoostというモデルが提案されています。

github.com

LightGBMをベースに混合効果部分を追加しているので、LightGBMのパラメータを経由して正則化項を追加できます。
もちろんRでも利用でき、gpboostとして実装されています。

参考文献

Baayen,R.H., Davidson, D.J., and Bates, D.M. (2008) Mixed-effects modeling with crossed random effects for subjects and items, Journal of Memory and Language
Bates,Douglas., Mächler, Martin., Bolker,Ben., and Walker, Steve.,Fitting Linear Mixed-Effects Models Using lme4
Groll, A. and G. Tutz (2014).Variable selection for generalized linear mixed models by L 1-penalized estimation
Sigrist, Fabio., 2020, Gaussian Process Boosting
Tibshirani, Robert., 1996, Regression Shrinkage and Selection via Lasso

lme4でp値が出力されない件と、lmerTestのp値 について

*1:モデルのパラメータ推定が収束していないので、無視しちゃいけないっちゃいけないですが

*2:ちなみにちょっと話すと$\Sigma$をYの分布パラメータと$\theta$及び$\Lambda$を用いて $\sigma ^2 \Lambda(\theta) \Lambda(\theta) ^T$と分解することで
推定をしやすくしているようですが、なぜこれでいいのかはわかりません。