3月25日に社会心理学会で「第2回春の方法論セミナー：GLMMが切り開く新たな統計の世界」が上智大学で開催されます。
第2回春の方法論セミナーのWebページ（社会心理学会Webサイト内）
セミナーでは，一般化線形混合モデル（GLMM）について，どういう方法で，社会心理学でどのように活用できるのかについて解説される予定です。
しかし，セミナーでは，具体的にどうやってGLMMを実行するか，どのソフトウェアが対応しているのか，などについては話をしません。
そこでこの記事では，SAS無償版でGLMMを実行する方法についてまとめておきます。
◆SAS無償版
昨年に，SASは個人利用を目的とした無償版を提供しています。基本的にはSASの学習用のものですが，SASのすべての機能が利用可能です。
もし手元にGLMMを使った方がいいデータをお持ちの人は，無償ですから，是非一度試してみてください。
無償版SASはWindows,Mac,Ubuntuなど，どんな環境でも利用できるのも利点です。
◆GLMMをSASで実行する
GLMMは，いろんなソフトウェアで実行可能です。Rでも，SPSSでも，Mplusでも可能です。しかし，SASのGLIMMIXプロシージャはGLMMを実行するうえで，他のソフトウェアに対していくつかアドバンテージがあります。
１．無償である（SPSSやMplusは有料）
２．推定精度のよい積分計算法（ガウスエルミート求積法）を選択できる（SPSSはあまりよい推定ができない手法しか選べない）
３．多様な分布を扱える（ほかのソフトでは選べない分布がある）
４．複数の変量効果を指定できる（Mplusではネストされたデータのみで，一つの変量効果のみ）
５．誤差の共分散構造を指定できる（RやMplusではそれができない）
などなど，SASは他のソフトウェアに対してかなり網羅的な推定についての指定が可能です。しかも，（今は）無償です。なのでSASを使えるようになっておけば，GLMMは全く怖いものではありません。
◆今回使うデータ
今回使うダミーのデータセットは，以下のようなデータです。
10人がそれぞれ4回刺激を提示されて，それに対して2値の反応をとる，というデータです。たとえば，ある英単語を覚えた後，4つの単語について正確に日本語訳を回答できたか，といったものです。説明変数は１～５の範囲を取る連続変量だとします。
IDは個人を識別する変数で，tは試行(trial)，yは目的変数，xは説明変数です。とりあえず3人分のデータだけを下に挙げておきます。
ID t y x
1 1 0 1
1 2 0 1
1 3 0 1
1 4 0 1
2 1 0 1
2 2 1 1
2 3 0 1
2 4 0 1
3 1 1 2
3 2 0 2
3 3 1 2
3 4 0 2
IDが縦に同じものが並んでいるのは，同じ人が4回反応していることを意味しています。SASのGLIMMIXプロシージャは，このように縦並びのデータを用います。
このデータは回答が２値なので，二項分布（あるいはベルヌーイ分布）に従うと考えられます。二項分布は平均が決まると自動的に散らばりも決まるタイプの分布（パラメータが１つしかない）なので，説明変数以外の個人差の変動を推定することができません。
もし説明変数で十分に個人差を説明できない場合，おそらく過分散が生じることが予想されます。
過分散が生じると，多くの場合，推定値が一致推定量にならないばかりか，推定精度が過大評価されます。よって，全然正しい推定ができなくなる，ということです。
なので，個人差を変量効果として推定するために，GLMMが必要なデータであると言えます。
◆SASの基本的な使いかた
SASはSPSSのようにグラフィカルなユーザーインターフェイスはなく，コードを書いてプログラムを実行するタイプのソフトウェアです。なので，データの読み込みなどもコードを書いて行う必要があります。
SASのファイル読み込みは，それだけで1つの記事がかけてしまうので，今回はデータをSASに貼り付けて実行する方法を書きます。
多くの方は，データの管理はtxtファイルcsvファイル，あるいはExcelファイルで行っていると思います。そこで，スペース区切り，カンマ区切り，タブ区切りのそれぞれについて書いておきます。
まず，分析データに名前を付けます。ここでは"test"としておきましょう。
data test;
SASでは文末にセミコロンが要ります。
次に，読み込みファイルを指定するのですが，ここではSASに貼り付けるので"cards"と書きます。
infile cards;
cardsは呪文だと思ってください。上の場合スペース区切りじゃないと読めないので，カンマ区切りにしたい場合は，
infile cards dlm = ',';
と書きます。dlm=は何でデータを区切っているかの指定です。タブ区切りの場合は，
infile cards dlm='09'x;
と書きます。'09'xはタブ区切りを意味します。
続いて，input文で読み込み変数名を書きます。もちろんこの変数名はデータの並びと同じである必要があります。
input ID t y x;
このとき，変数の中で数値ではなく文字データ（男・女，とか，正解・不正解とか）の場合は，変数名のあとに$をつけます。もしyが"正解" と"不正解"という文字データの場合，
input ID t y $ x;
と入力します。
続いて，cards;というコードの下にデータを入力します。データの最後にもセミコロンがいるので注意してください。
最後にRUN;をつけると，そこまでのコードを一度切って走らせろ，ということになります。
これらのコードを通して入力すると（カンマ区切りの場合），
DATA test;
INFILE CARDS dlm=',';
INPUT ID t y x;
CARDS;
1,1,0,1
1,2,0,1
1,3,0,1
1,4,0,1
2,1,0,1
2,2,1,1
2,3,0,1
2,4,0,1
3,1,1,2
3,2,0,2
3,3,1,2
3,4,0,2
4,1,1,2
4,2,1,2
4,3,0,2
4,4,0,2
5,1,0,3
5,2,0,3
5,3,0,3
5,4,0,3
6,1,1,3
6,2,1,3
6,3,1,3
6,4,1,3
7,1,1,4
7,2,1,4
7,3,1,4
7,4,1,4
8,1,0,4
8,2,1,4
8,3,1,4
8,4,0,4
9,1,1,5
9,2,1,5
9,3,1,5
9,4,1,5
10,1,1,5
10,2,1,5
10,3,1,5
10,4,1,5
;
RUN;
というようになります。ここまでをSASにコピペして走らせてみましょう。
エラーは英語で出ますが，セミコロンが抜けてないか，全角文字が入ってないか，スペルミスがないかを注意してみてください。
◆GLIMMIXプロシージャの使いかた
それでは，個人から反復測定しているので個人差の程度を変量効果に入れた，GLMMを構成しましょう。まず，PROC GLIMMIXとかき，その後に推定方法についての指定をします。
PROC GLIMMIX method = quad;
method = quadは，ガウス－エルミート求積法を意味します。もし積分計算をラプラス近似にする場合は，method = laplaceとします。それ以外にも，SPSSが採用している擬尤度を使う場合は，method = RSPLとします。
積分計算の方法はどれがいいかわかりにくいかと思いますが，次の基準で選ぶといいでしょう。
１．基本はquadで指定。
２．変量効果の種類が二つ以上ある場合はlaplaceで指定（２つ以上だとquadが使えない）。
３．残差共分散の構造を指定したい場合は，RSPLで指定。
ただ，RSPLは尤度を使って推定していないので，情報量規準が利用できない，一致推定量にならないなどの致命的な欠点があります。なので，基本はquadにしましょう。laplaceはあくまで近似なので精度は落ちますが，あまりにモデルが複雑じゃなければquadとほとんど同じ結果になるようです。
続いて，カテゴリカル変数を指定します。
CLASS ID;
今回は個人差をモデルにいれるので，IDをカテゴリカル変数としてい指定します。もし説明変数も名義データの場合はCLASSに指定します。今回は説明変数が連続量を仮定しているので指定していません。
次に，一番重要なモデルを記述します。まず，MODELと書いて，続いて式を書きます。
MODEL y = x/ dist = binomial link = logit solution ddfm = bw;
これでOKです。これで二項分布に従うyをxで予測することを意味します。solutionは回帰係数を出力するためのオプションです。ddfmは自由度の計算方法ですが，今回のようなネスト構造のデータの場合はbw（Between-Within）とします。データがネストされていない構造なら何も書かなくていいです。
※注：SatterthwateやKenword-Rogerの方法は，quadやlaplaceでは使用できません。RSPLなら可。
dist = binomialは，二項分布に従うことを意味します。二項分布を指定するとデフォルトでリンク関数はlogitになるので，別に書かなくてもいいです。もしリンクを自分でいじりたければ，link = のあとを変えればOKです。
GLIMMIXで指定可能な分布は，
| gaussian | 正規分布 |
| binomial | 二項分布 |
| poisson | ポアソン分布 |
| negbinomial | 負の二項分布 |
| beta | ベータ分布 |
| binary | ベルヌーイ分布 |
| lognormal | 対数正規分布 |
| multinomial | 多項分布 |
などがあります。他にも指数分布や逆正規分布などもあります。
リンク関数には，以下の様なものがあります。
| identity | 恒等 |
| logit | ロジット |
| cumlogit | 累積ロジット |
| log | 対数 |
| probit | プロビット |
| cumprobit | 累積プロビット |
identityは正規分布のデフォルトのリンク関数で，変換なし，ということです。logitは二項分布とセットでロジスティック回帰になります。cumlogitとmultinomialを組み合わせると順序ロジスティックになります。logはpoissonのデフォルトで，この組み合わせでポアソン回帰になります。
今回は二項分布なのでリンクはデフォルトでlogitになるわけです。
モデルが指定終わったら，次に変量効果を指定します。RANDOMと書いて，
RANDOM intercept / subject = ID;
と書きます。これで，IDで区別された目的変数の個人差を変量効果として推定します。interceptは切片の意味です。intercept以外に，説明変数を指定すると，回帰係数の変量効果も推定することができます。これは回帰係数の個人差を推定することを意味しています。ただ，多くの場合回帰係数の変量効果を推定すると，計算量が多くなるので，むやみに指定すると不安定になることがあるので注意しましょう。
なお，subjectで指定できるのはカテゴリカル変数のみです。なので，前もってCLASSで指定しておかないといけない点は注意しましょう。
最後に，RUN;と書くとプロシージャを走らせて出力する，という指定になります。
これらをまとめて書くと，
PROC GLIMMIX method = quad;
CLASS ID;
MODEL y = x/ dist = binomial link = logit solution ddfm = bw;
RANDOM intercept / subject = ID;
RUN;
この結果は，以下のようになります。
| 適合度統計量 | |
|---|---|
| -2 対数尤度 | 37.32 |
| AIC (小さいほどよい) | 43.32 |
| AICC (小さいほどよい) | 43.99 |
| BIC (小さいほどよい) | 44.23 |
| CAIC (小さいほどよい) | 47.23 |
| HQIC (小さいほどよい) | 42.33 |
| 条件付分布の適合度統計量 | |
|---|---|
| -2 log L(y | r. effects) | 26.26 |
| Pearson カイ 2 乗 | 20.99 |
| Pearson カイ 2 乗 / DF | 0.52 |
| 共分散パラメータの推定 | |||
|---|---|---|---|
| 共分散パラメータ | サブジェクト | 推定値 | 標準誤差 |
| Intercept | ID | 2.2730 | 2.9794 |
| 固定効果の解 | |||||
|---|---|---|---|---|---|
| 効果 | 推定値 | 標準誤差 | 自由度 | t 値 | Pr > |t| |
| Intercept | -3.7373 | 1.9359 | 8 | -1.93 | 0.0896 |
| x | 1.4864 | 0.6940 | 8 | 2.14 | 0.0646 |
| 固定効果の Type III 検定 | ||||
|---|---|---|---|---|
| 効果 | 分子の自由度 | 分母の自由度 | F 値 | Pr > F |
| x | 1 | 8 | 4.59 | 0.0646 |
結果は，かろうじて有意ではなく，xの効果があるとは主張しがたいです。
ではここで，変量効果を指定しないモデルを走らせてみましょう。RANDOMステートメントを*マークでコメントアウトして走らせてみます。
PROC GLIMMIX method = quad;
CLASS ID;
MODEL y = x/ dist = binomial link = logit solution ddfm = bw;
*RANDOM intercept / subject = ID;
RUN;
すると以下のようになります。
| 適合度統計量 | |
|---|---|
| -2 対数尤度 | 39.83 |
| AIC (小さいほどよい) | 43.83 |
| AICC (小さいほどよい) | 44.15 |
| BIC (小さいほどよい) | 47.21 |
| CAIC (小さいほどよい) | 49.21 |
| HQIC (小さいほどよい) | 45.05 |
| Pearson カイ 2 乗 | 35.88 |
| Pearson カイ 2 乗 / DF | 0.90 |
| パラメータ推定値 | |||||
|---|---|---|---|---|---|
| 効果 | 推定値 | 標準誤差 | 自由度 | t 値 | Pr > |t| |
| Intercept | -2.7071 | 0.9919 | 38 | -2.73 | 0.0096 |
| x | 1.0543 | 0.3335 | 38 | 3.16 | 0.0031 |
| 固定効果の Type III 検定 | ||||
|---|---|---|---|---|
| 効果 | 分子の自由度 | 分母の自由度 | F 値 | Pr > F |
| x | 1 | 38 | 9.99 | 0.0031 |
なんと，推定値は小さくなったのに，余裕で有意になってしまいました。このように，独立にサンプリングされていないデータでは，個人差を考慮に入れないと，ぜんぜん違う結果になってしまいます。当然，前者の結果が正しいといえるでしょう。
なお，今回は４回の試行をそれぞれ２値データとして分析しましたが，４回中何回正解したか，という二項分布のデータとしても分析することができます。この場合でも実は結果はおなじになります。ただ，尤度は変わってきます。
どちらにしても個人差を考慮しないといけない点は，同じです。以下のコードを走らせます。ここでは，yには先程のyの個人事の合計値を入れています。
DATA test;
INFILE CARDS dlm=',';
n = 4;
INPUT ID y x;
CARDS;
1,0,1
2,1,1
3,2,2
4,2,2
5,0,3
6,4,3
7,4,4
8,2,4
9,4,5
10,4,5
;
RUN;
PROC GLIMMIX method = quad;
CLASS ID;
MODEL y/n = x/ dist = binomial link = logit solution ddfm = bw;
RANDOM intercept / subject = ID;
RUN;
さっきとちょっと違うのは，n=4;というのがINFILEのあとにあります。これはnという新しい変数を定義して，それに４を入れろ，ということを意味しています。
そして，GLIMMIXのMODELステートメントに，y/n = xと書いています。これは，４回中のyの比率を目的変数としたモデルを推定しろ，ということです。
すると以下の様な結果が得られます。
| 適合度統計量 | |
|---|---|
| -2 対数尤度 | 23.80 |
| AIC (小さいほどよい) | 29.80 |
| AICC (小さいほどよい) | 33.80 |
| BIC (小さいほどよい) | 30.71 |
| CAIC (小さいほどよい) | 33.71 |
| HQIC (小さいほどよい) | 28.80 |
| 条件付分布の適合度統計量 | |
|---|---|
| -2 log L(y | r. effects) | 12.73 |
| Pearson カイ 2 乗 | 3.16 |
| Pearson カイ 2 乗 / DF | 0.32 |
| 共分散パラメータの推定 | |||
|---|---|---|---|
| 共分散パラメータ | サブジェクト | 推定値 | 標準誤差 |
| Intercept | ID | 2.2730 | 2.9794 |
| 固定効果の解 | |||||
|---|---|---|---|---|---|
| 効果 | 推定値 | 標準誤差 | 自由度 | t 値 | Pr > |t| |
| Intercept | -3.7373 | 1.9359 | 8 | -1.93 | 0.0896 |
| x | 1.4864 | 0.6940 | 8 | 2.14 | 0.0646 |
| 固定効果の Type III 検定 | ||||
|---|---|---|---|---|
| 効果 | 分子の自由度 | 分母の自由度 | F 値 | Pr > F |
| x | 1 | 8 | 4.59 | 0.0646 |
一番最初の結果と，推定結果は同じです。ただ，尤度が異なります。
今回は変量効果が一つだけのモデルを組みましたが，次の記事では変量効果が二つの場合の例を挙げます。