はじめに

階層データに対してマルチレベルモデリングなどの線形混合モデルをあてはめる際、よく使われるアプリケーションにRのlme4パッケージというのがあるのですが、これはRの関数にしては珍しく、フィッティングしたモデルの係数についてp値を出力しないという特徴があります。いままで僕はこの特性について、作者が反p値主義者だからとかそういう「思想」的な産物によるものだと思っていましたが、少し調べてみるとどうやらそういった理由によるものではないらしいということが分かりました。以下、lme4パッケージの作者Douglas Batesがこの件について述べた簡単な文章の拙訳を掲載します。原文はここです。

lmer, p-values and all that (2006)

lmer¹ パッケージによる線形混合モデルのあてはめ結果にp値が含まれていないことに対して驚愕するユーザーも少なくない。同様に単一のlmerモデルに対する分散分析の結果ついても、固定効果の項それぞれの二乗和とそれに対応する分子の自由度は提供するが、分母の自由度とp値については提供しない。

ユーザーは分母の自由度とそれに対応するp値が簡単に計算できると思っているので、これはlmerを作った人（私）に欠陥があったからだと言っている。

ここで私は再び、なぜ（lmerに）p値を引用しないのか、つまり、なぜ私がSASの結果を再現するという「明らかに正しい」アプローチをとらないのかについて説明しよう。もちろんRプロジェクトの目的はSASによって出力されるp値の結果を再現することだと考える人もいるけれども、私はその一派ではない。もし彼らが私を異端児で、地獄の炎で焼きつくされるべきだと感じるなら、私は来月開催される2006年Rユーザーの会の懇親会でエキサイティングな最後を迎えるのを心待ちにしようと思う。

よく知られているように、固定効果モデル行列²の係数についてのt検定量は分子の自由度が1のF検定量の平方根である。したがって、一般性を失うことなく、我々は分散分析の結果に現れるF検定量に集中することができる。大昔に「分散分析」または「実験計画」という授業を受けたことのある人なら、観測平方和と期待平方和に基づく分散の要素の推定方法と"error strata"³に基づく検定の方法を習ったはずだ。（もし覚えるように言われていなければあなたはラッキーかもしれない）。よってlmerモデル（SASによるものも含めて）によって作られるF検定量はerror strataに基づくものである。しかしこれが簡単な話ではない。

lmerによるパラメータ推定は最尤法または**制限付き最尤法(RMEL)**によって行われるが、それらは観測二乗和・期待二乗和に基づくものでもerror strataに基づくものでもない。MLやRMELによる推定はlmerが、何重にもネストした要因や部分的に重なりのあるグループを持つ偏りのあるデザインを扱うことを可能にしたという意味で、良い方法であることは間違いない。そしてこうした状況は大規模な観察研究でよく起こるケースである。

ここに私にとって興味深く、何年もの間取り組んできた線形混合モデルの式とパラメータ推定に関する性質がある。しかしここではあまり立ち入らずに、単にモデルとデータが与えられて、パラメータの推定値が求められたとしよう。どのようにF検定量が計算されるのだろうか？　平方和と分子の自由度は線形モデルと同様に計算される。線形混合モデルには線形モデルの"effect"要素と似たスロット⁴があり、モデル行列への属性の割り当てと同時に分子のF比を提供するスロットがある。分母は制限付き残差平方和をREML自由度で割ったもので、nを観測数、pを固定効果に対するモデル行列の列ランクとするとn-pに等しい。

ここで"the penalized residual sum of squares（制限付き残差平方和）"の"the"に注目して最後の文章を読んでみよう。すべてのF比は同じ分母を共有している。繰り返して言うが、すべてのF比は「同じ」分母を使っている。これが私が前提に問題があると考える理由である。ここでいう前提とはF統計量の分布は可変な分母の自由度を除いて既知の自由度を持つF分布に従い、検定ごとに異なる分母の自由度を割り当てる公式を発見することに酔ってどのようにp値を計算するかという問いに答えられるという前提である。分母は変わらない。分母は変わらない。なぜ分母の自由度は変わるのか？

混合モデルの固定効果の検定に関する研究の大部分はこれらの検定量が既知の自由度をもつF分布に従うと仮定しており、分母の自由度の近似をもとめることを目的としている。私は賛成しかねるが。

私が思うに実りあると考え、求めつづけている方法が一つだけである。混合効果モデルの罰則付きの二乗誤差は実際の罰則付き残差二乗和において残差の二乗和となる。そしてこの中に罰則付き最小二乗問題の自由度のように振る舞う量がある。私はそれを計算するコードを埋め込み、シミュレーションで振る舞いをみるつもりだ。

当面の間は、係数の一つ一つを評価するのであれば、私はMCMCサンプリングをおすすめする。それぞれの項を評価するのは難しいが、常にF比と分母の自由度の下限を得ることができる。

もしどうしてもSASによって得られる「明らかに正しい」分母の自由度の計算の実装に貢献したい人がいれば協力するのも吝かではない。しかしlmerに必要な罰則付きの二乗和アプローチと疎な行列の計算方法にはかなりの量の数式の移植を必要とすることを警告しておく。これらの式はn行n列行列の逆行列の計算する必要があり、nが100万を超える場合にはあなたも絶対に、絶対に計算したくないはずだ。

感想

読んでみると、真面目な話の間にちょいちょいきつめの皮肉が込められていてやっぱり癖のある人なんだなと思いました。技術的なことに関して言うと、結局混合モデルの係数をF検定の枠組みで検定するのは難しいから、どうしても検定をしたいならMCMCを使えってことなんでしょうか。この文章は2006年に出されたものなので、その後混合モデルのパラメータ推定に関する研究がどう進んだのかも気になります。この分野に詳しい方がいたら教えて下さい。

Linear Mixed Eeffect model via REML methodの略 ↩
model matrix の略。モデルを行列で表すということに対してまだ自分自身よくわかっていないので誤訳かもしれません。 ↩
すみませんここの訳分かりませんでした。直訳すると「誤差層」になりますが…。 ↩
これは専門用語か一般用語かの判別が付きませんでした。一般的に言うと「枠」とかでしょうか。 ↩