GLMM（一般化線形混合モデル）をSASで実行する方法２

 

昨日アップした，GLMMをSASで実行する方法１に続いて，変量効果を２種類実行する方法を書きます。この記事の前に，１のほうを先に見てください。

 

◆個人差以外の変動を変量効果として推定する

先ほどの分析では，個人差の変動を変量効果として推定しましたが，ここで提示する刺激がすべての参加者で同じ場合，刺激による回答の変動も目的変数の分散を説明する可能性があります。

例えば，すべての回答者に同じ英単語４つ提示する場合，単語ごとの覚えやすさが異なりますから，その影響をモデルに組み込む必要が出てきます。もちろん，その単語の覚えやすさにほとんど違いがない場合もあるでしょう。それらの違いはAICなどの情報量規準を参考にすれば，その変動をモデルに入れたほうが予測誤差が小さくなるかどうかの判断ができます。

 

◆SASのGLIMMIXプロシージャで２つ以上の変量効果を推定する

個人差以外に，提示した刺激差という違う変動を変量効果に組み入れる場合，変量効果の種類が２つとなります。

この場合，前の記事にも書いたように，ガウス-エルミート求積法（method = quad）が使えません。なので，変量効果が２種類以上ある場合はラプラス近似による積分計算を用います。ただ，同じ個人変動について切片と回帰係数という二つの変量効果を推定する場合は，1種類と判断されます。要は，個人差と刺激差というように，発生源が異なるものが二つ以上ある場合はラプラス近似を使わないといけない，ということです。

刺激の種類はtという変数で識別されるので，前もってCLASSにtも追加しておきます。

するとコードは，

PROC GLIMMIX method = laplace;
CLASS ID t;
MODEL y = x/ dist = binomial link = logit solution ddfm = bw;
RANDOM intercept / subject = ID;
RANDOM intercept / subject = t;
RUN;

というように，RANDOMステートメントを二つ書きます。

これを走らせると，以下のような結果になります。

適合度統計量
-2 対数尤度	37.44
AIC (小さいほどよい)	45.44
AICC (小さいほどよい)	46.59
BIC (小さいほどよい)	37.44
CAIC (小さいほどよい)	41.44
HQIC (小さいほどよい)	37.44

条件付分布の適合度統計量
-2 log L(y | r. effects)	25.70
Pearson カイ 2 乗	20.29
Pearson カイ 2 乗 / DF	0.51

共分散パラメータの推定
共分散パラメータ	サブジェクト	推定値	標準誤差
Intercept	ID	2.2501	3.2710
Intercept	t	0.06917	0.7470

固定効果の解
効果	推定値	標準誤差	自由度	t 値	Pr > |t|
Intercept	-3.8041	2.0912	3	-1.82	0.1665
x	1.5100	0.7605	27	1.99	0.0573

固定効果の Type III 検定
効果	分子の自由度	分母の自由度	F 値	Pr > F
x	1	27	3.94	0.0573

xの推定値が前の記事のものとは少し変わり，やや大きくなっています。しかし，ｔの変量効果の分散成分は標準誤差に比べてかなり小さいことがわかります。

加えて，AICをみると45.44と，前の43.32よりも大きくなっています。このことからも，ｔの変動はモデルに入れなくてもいい，という判断ができそうです。

 

◆ネストされたデータ構造ではないけど，GLMMを使う必要がある場合

GLMMはサンプリングの非独立性があるようなデータにだけ使うものだ，と思う方もいるかもしれませんが，そうではありません。

例えばポアソン分布や二項分布は正規分布と異なり，推定するパラメータが一つだけです。正規分布は平均と分散の二つを推定しますが，ポアソンや二項分布は平均だけなのです。つまり，これらの分布は平均が決まると，分散が自動的に決まってしまうような分布なのです。

しかし，心理学にかぎらず，データが全てその分布に綺麗に当てはまるわけではありません。緑本にあるように，どうしても既存のきれいな分布からはみ出る，説明変数では説明しきれない個体差や個人差がデータには含まれます。

このような，分散パラメータが平均によって決まるタイプの分布をモデルに当てはめるとき，説明しきれない残差変動が生じることがあります。これを過分散（Over-dispersion）と呼びます。過分散が生じると，推定の一致性が損なわれたり，標準誤差が小さく推定されたり，正しい推定ができなくなります。

例えば次のような負の二項分布に従うデータyを考えます。負の二項分布は，ポアソン分布に比べて分散が大きい分布なので，ほぼ過分散が生じます。また，その乱数と弱い関連がある別の変数xを作りました。

ID    y    x
1    1    4
2    5    3
3    4    6
4    6    5
5    7    5
6    1    5
7    7    3
8    3    4
9    1    4
10    9    7
11    0    3
12    5    4
13    1    2
14    0    4
15    0    5

これをSASでポアソン回帰をしてみましょう。ただのポアソン回帰ならGLMプロシージャでもいいのですが，あとでGLMMをするのでGLIMMIXプロシージャを使います。

DATA test2;
INFILE CARDS  dlm=',';
INPUT ID y x;
CARDS;
1,1,4
2,5,3
3,4,6
4,6,5
5,7,5
6,1,5
7,7,3
8,3,4
9,1,4
10,9,7
11,0,3
12,5,4
13,1,2
14,0,4
15,0,5
;
RUN;

PROC GLIMMIX method = quad;
CLASS ID;
MODEL y = x/ dist = poisson link = log solution;
RUN;

このコードを走らせると，

適合度統計量
-2 対数尤度	75.62
AIC (小さいほどよい)	79.62
AICC (小さいほどよい)	80.62
BIC (小さいほどよい)	81.04
CAIC (小さいほどよい)	83.04
HQIC (小さいほどよい)	79.61
Pearson カイ 2 乗	33.90
Pearson カイ 2 乗 / DF	2.26

パラメータ推定値
効果	推定値	標準誤差	自由度	t 値	Pr > |t|
Intercept	-0.00465	0.5329	13	-0.01	0.9932
x	0.2698	0.1093	13	2.47	0.0283

固定効果の Type III 検定
効果	分子の自由度	分母の自由度	F 値	Pr > F
x	1	13	6.09	0.0283

xの係数は0.2698で，p値は<.05で有意でした。

しかし，これはもともとポアソン分布に従わないデータですから，本来はポアソン回帰をするのは誤りです。ではどうすればいいかといえば，三つの選択肢があります。

１．負の二項回帰分析を使って過分散を推定する
２．GLMMを使って，個体差の分散を正規分布で推定する
３．ロバストな標準誤差を使って，過分散で生じた標準誤差のバイアスを補正する（次のセクション）

今回は負の二項分布に従うデータを発生させているわけですから，負の二項分布回帰を使うのが一番いいでしょう。

PROC GLIMMIX method = quad;
CLASS ID;
MODEL y = x/ dist = negbinomial link = log solution;
RUN;

上のコードで負の二項分布回帰を実行すると，AICは79.62から74.39と小さくなり，xの推定値は0.2343でp値は0.225と有意ではありませんでした。また負の二項分布の分散成分は0.5889と推定されました。このように，ポアソン回帰だと過分散が生じていた，ということがわかります。

次に，GLMMを使ってポアソン回帰で説明できない個人差を別に正規分布で推定する，ということをやってみましょう。コードは以下です。なお，上ではありませんでしたがddfm=bwが付いている点に注意です。これがないとxの自由度が0として計算されてしまいます。

PROC GLIMMIX method = quad;
CLASS ID;
MODEL y = x/ dist = poisson link = log solution　ddfm=bw;
RANDOM intercept / subject = ID;
RUN;

今回のようにネストされていないデータであっても，変量効果の推定ができます。この結果，AICは75.40.とポアソンの79.62よりは小さくなりました。xの効果も非有意となりました。そして，推定された個人差の分散は0.5368でした。もちろん負の二項分布回帰に比べれば，確かにAICは悪いですが，この方法はポアソンでなくても二項分布など，他の分布に対しても使えるのが利点です。

 

◆ロバスト標準誤差による補正

最後に，ロバスト標準誤差を用いることで，過分散によるバイアスを補正する方法を紹介します。

次のように，EMPIRICALというオプションを加える事で，分布の仮定からの逸脱による標準誤差のバイアスを調整することができます。

PROC GLIMMIX method = quad EMPIRICAL;
CLASS ID;
MODEL y = x/ dist = poisson link = log solution ;
RUN;

これを走らせると，xの推定値は0.2698と同じですが，標準誤差が0.1093から0.1473と大きくなりました。これである程度は分布からの逸脱を補正することができます。

ただ，これはあくまで補正なので，データの生成メカニズムそのものを表現しているわけではありません。今回の例では負の二項分布から発生するデータなので負の二項分布回帰を使うべきですが，データが特定の分布と個体差という二つの変動が組み合わさったデータが手に入ることがほとんどです。そういう場合は，GLMMを使ってデータの発生をちゃんとモデリングするほうがいいでしょう。

あと，当然GLMMでもロバスト標準誤差を利用できます。次のように書けば，過分散を正規分布で推定しながら，それでもポアソン分布から逸脱した部分を考慮に入れて標準誤差を補正します。

PROC GLIMMIX method = quad EMPIRICAL;
CLASS ID;
MODEL y = x/ dist = poisson link = log solution ddfm = bw;
RANDOM intercept / subject = ID;
RUN;

このように，人事を尽くしてロバスト標準誤差，ではないですが，可能な限りモデリングをした上で，最後に補正する，というのがいいかなーというのが僕の考えです。

今回はSASに限った話をしましたが，次はRでGLMMの記事を書こうと思います。