ダミーコーディングとANOVAコーディング

対数線形モデルで各要因の係数を求める際，一般的には基準セルからの偏差を取る方法（ダミーコーディング）が多く用いられているようです。
しかし，心理学系の論文では推定値の合計が0になるよう標準化する方法（ ANOVAコーディング や Effectコーディング などと呼ばれる）が多く用いられています。

しかし，このANOVAコーディングによる係数の算出は，一般の統計ソフトで行うのはなかなか厄介です。できなくはないのですが，多くのソフトウェアはすべての係数を表示してくれないのです。

ダミーコーディングでは，特定のセル（ほとんどの場合，１番最初または１番最後）のセルを基準として係数を算出しますので，基準となったセルについては係数が算出されません。そのため，結果の出力では，それら基準となったセルについては情報は出力されません。

しかし，ANOVAコーディングでは，特定のセルを基準としているわけではないのですべてのセルに対して係数が算出されます。ところが，SPSSやSASといったメジャーどころを含むほとんどの統計ソフトは，ANOVAコーディングでの係数の表示の際にも，ダミーコーディングの場合と同じように，最初あるいは最後のセルの表示を省略してしまうのです（SPSSの場合やSASの場合はこちらを参照）。

ANOVAコーディングでは係数の合計が0になるように決定されるので，表示されていない部分はそれ以外の係数の合計を0から引けばも止まるのですが，いちいち手計算をするのも面倒ですし，要因の数が増えれば凡ミスをする可能性も高くなります（実際，こちらのSASの場合の資料では，「ｕ.4. の標準誤差は〜」のところで0.0814を0.0314と間違えて計算しています）。

Rにおける対数線形モデル

Rでは，対数線形モデルを実施できる関数が複数種類あり，その中には，loglm()のようにデフォルトでANOVAコーディングを用いているものもあります。

ここで，SASの場合の数値データを用いて，いくつかの方法で対数線形モデルを用いた分析を実施してみます。使用するデータは次の通りです。

sample_data <- '
性別 年齢   分類      Freq
男   ~19   1.親族       92
男   ~19   2.強盗殺人  110
男   ~19   3.その他    144
男   20~29 1.親族      175
男   20~29 2.強盗殺人  191
男   20~29 3.その他    389
男   30~39 1.親族      212
男   30~39 2.強盗殺人   85
男   30~39 3.その他    217
男   40~   1.親族      263
男   40~   2.強盗殺人   34
男   40~   3.その他    188
女   ~19   1.親族       18
女   ~19   2.強盗殺人    2
女   ~19   3.その他      8
女   20~29 1.親族       68
女   20~29 2.強盗殺人    7
女   20~29 3.その他     22
女   30~39 1.親族       69
女   30~39 2.強盗殺人    4
女   30~39 3.その他     17
女   40~   1.親族       46
女   40~   2.強盗殺人    3
女   40~   3.その他      2'

df <- read.table(textConnection(sample_data),head=T)

`loglm()`による対数線形モデル分析

loglm()はMASSパッケージに含まれています。loglm()は，度数を集計済みのデータからでも度数集計前のデータからでもぶんせきをおこなうことができます。度数集計済みのデータの場合には，formulaの部分が度数 ~ 要因A + 要因Bのような形になり，度数集計前のデータの場合には，formulaは~ 要因A + 要因Bという形になります。

サンプルのデータはすでに度数集計済みですので，分析用のスクリプトは次のようになります。ここでは，各要因の主効果と交互作用のすべてを含む飽和モデルの結果を求めています。

require(MASS)

res <- loglm(Freq ~ 性別*年齢*分類, data=df)

結果の表示についてはsummary(res)としたいところですが，summary()関数ではまともな結果は表示されません。代わりに，coef(res)で係数の推定値を表示させることができます。なお，モデル全体の分析結果についてはresあるいはanova(res)で表示できますが，飽和モデルの場合には測定値はすべて説明され尽くしますので，尤度比もカイ２乗も0になります。

loglm()で得られた推定値をみてみましょう。coef(res)とすると，次のような出力が得られます（一部のみ抜粋）。性別，年齢，分類のそれぞれで係数の合計が0になり，ANOVAコーディングされていることがわかります。

> coef(res)
$`(Intercept)`
[1] 3.711677

$性別
       女        男 
-1.295557  1.295557 

$年齢
       ~19      20~29      30~39        40~ 
-0.4025032  0.5611728  0.2271145 -0.3857841 

$分類
     1.親族  2.強盗殺人    3.その他 
 0.76181577 -0.83096082  0.06914505

係数だけが必要ならこれでよいのですが，大抵の場合には標準誤差や$z$が必要になるでしょう。残念ながら，loglm()ではそれらの値を得ることができません。また，loglm()では各要因の主効果や交互作用についての情報も得られないので，モデル選択をしたい場合にも適しません。

`glm()`による対数線形モデル分析

ダミーコーディングによる推定

その名前が示す通り，対数線形モデルも線形モデルの一種ですから，glm()で実施することもできます。glm()で対数線形モデルによる分析を実施するには，family = poisson(link="log")あるいはfamily = "poisson"と指定します。

res <- glm( Freq ~ 性別*年齢*分類, family='poisson', data=df)

基本的な結果はsummary()で得ることができます（一部を抜粋）。

> summary(res)

Call:
glm(formula = Freq ~ 性別 * 年齢 * 分類, family = "poisson", 
    data = df)

Deviance Residuals: 
 [1]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

Coefficients:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                      2.89037    0.23570  12.263  < 2e-16 ***
性別男                           1.63142    0.25773   6.330 2.45e-10 ***
年齢20~29                        1.32914    0.26507   5.014 5.32e-07 ***
年齢30~39                        1.34373    0.26467   5.077 3.83e-07 ***
年齢40~                          0.93827    0.27802   3.375 0.000739 ***
分類2.強盗殺人                  -2.19722    0.74536  -2.948 0.003200 ** 
分類3.その他                    -0.81093    0.42492  -1.908 0.056335 .

さて，この結果では「性別男」の係数は表示されていますが，「性別女」がありません。また，年齢では「年齢~19」が，分類では「分類1.親族」がありません。さらに，それぞれの係数は，先ほどのloglm()のものとは異なっています。これは，glm()による対数線形モデルの推定値がダミーコーディングによって算出されているためです。つまり，性別，年齢，分類についての各水準の係数を算出する際に，ここでは表示されていない「性別女」，「年齢~19」，「分類1.親族」のそれぞれを基準としているのです。基準となったセルの係数は0ですので，結果では省略されます。

ANOVAコーディングによる推定

さて，ANOVAコーディングでの推定値を得るにはどうすればよいでしょうか。ANOVAコーディングによる推定値は，次のようにしてcontrastを指定することで求められます。

# ANOVAコーディング
res <- glm( Freq~性別*年齢*分類, family='poisson', data=df,
				contrasts=list(
					性別=contr.sum, 
					年齢=contr.sum, 
					分類=contr.sum))

結果は次のようになります。

> summary(res)

[ 略 ]

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)        3.71168    0.06127  60.583  < 2e-16 ***
性別1             -1.29556    0.06127 -21.147  < 2e-16 ***
年齢1             -0.40250    0.11653  -3.454 0.000552 ***
年齢2              0.56117    0.08223   6.825 8.80e-12 ***
年齢3              0.22712    0.09240   2.458 0.013969 *  
分類1              0.76182    0.06663  11.434  < 2e-16 ***
分類2             -0.83096    0.10212  -8.137 4.05e-16 ***

[ 略 ]

この結果で表示されている係数の推定値は，ANOVAコーディングを用いているloglm()のものと同じであることが確認できます。

また，carパッケージのAnova()を使用すれば，それぞれの主効果，交互作用についての検定も可能です。

> require(car)
> Anova(res)
Analysis of Deviance Table (Type II tests)

Response: Freq
               LR Chisq Df Pr(>Chisq)    
性別            1616.40  1  < 2.2e-16 ***
年齢             198.59  3  < 2.2e-16 ***
分類             263.06  2  < 2.2e-16 ***
性別:年齢         21.85  3  7.025e-05 ***
性別:分類        166.26  2  < 2.2e-16 ***
年齢:分類        195.51  6  < 2.2e-16 ***
性別:年齢:分類    10.40  6     0.1086    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

しかしsummary()による結果には，ダミーコーディングの場合と同様に，各要因の最後の水準の結果が表示されていません。SASやSPSSと同様に，Rのglm()もANOVAコーディングで計算はしてくれるのですが，結果をきちんと表示してくれないのです。

ただし，glmオブジェクトからは残りの推定値の算出に必要な値を得ることができますので，これらを用いて係数や標準誤差を求めることは可能です。残りの係数を算出するにはcoef()で得られる係数行列を，残りの標準誤差を算出するにはvcov()で得られる分散共分散行列を使用します。

たとえば，summary()の結果では年齢は年齢1，年齢2，年齢3の３つしかありません（loglm()の係数と比較してみると，表示されていないのは40~の係数です）が，残りの係数（40〜の係数）は次の手順で求めることができます。


# 係数行列を取得（summary()の結果に対してcoef()を使用する）
> res.coef <- coef(summary(res)) 

# 係数を確認（年齢の係数は3行目〜5行目）
> res.coef
                      Estimate Std. Error      z value     Pr(>|z|)
(Intercept)        3.711677011 0.06126573  60.58325301 0.000000e+00
性別1             -1.295557355 0.06126573 -21.14652725 2.970190e-99
年齢1             -0.402503224 0.11653036  -3.45406321 5.522079e-04
年齢2              0.561172832 0.08222508   6.82483787 8.802467e-12
年齢3              0.227114505 0.09239540   2.45807163 1.396853e-02
分類1              0.761815767 0.06662498  11.43438722 2.815172e-30
分類2             -0.830960816 0.10212235  -8.13691454 4.054787e-16

[ 略 ]

# 年齢40~の係数は，(-1) * (他の年齢水準[3行目~5行目]の総和) で求められる
> 年齢40.coef <- (-1) * sum(res.coef[3:5, 1])  
> 年齢40.coef
[1] -0.3857841

# 分散共分散行列を取得（glmオブジェクトにvcov()を使用）
> res.vcov <- vcov(res)

# 標準誤差は，(他の年齢水準の分散と共分散の総和) の正の平方根
> 年齢40.stderr <- sqrt(sum(res.vcov[3:5,3:5]))
> 年齢40.stderr
[1] 0.1271404

# zは係数/標準誤差
> 年齢40.z <- 年齢40.coef / 年齢40.stderr
> 年齢40.z
[1] -3.034316

# Pr(>|z|)は両側確率
>  pnorm(abs(年齢40.z),lower.tail=F)*2
[1] 0.002410814

これならば，結果の数値を手入力して計算し直すよりは凡ミスの可能性は低くなるでしょう。しかしなかなか面倒です。主効果でこれですから，交互作用になったら大変です。たとえば，glm()の結果では「年齢40~」と「分類」の要因の交互作用についての情報は表示されません。この部分の係数そのものは先ほどの手順と同様にして単純に求められますが，標準誤差を計算するためには，すでにある「年齢1」，「年齢2」，「年齢3」と「分類1」,「分類2」の部分の分散共分散行列を合成して「年齢4」と「分類1」,「分類2」の分散共分散行列を作成し，そしてその結果から「年齢4」と「分類3」の標準誤差を求めるといった手順が必要になります。以下に，実際の作業手順を示しておきます。ここでは，年齢*分類の交互作用について，まず年齢4と分類1~2の分散と共分散を求めて完全な分散共分散行列を作成し，それらを用いて年齢4と分類3の分散共分散行列を作成しています。標準誤差は作成した分散共分散行列の１番右下にある分散の正の平方根です。

#### 年齢:分類の係数の算出
age.type.vcov1<- res.vcov[13:15,13:15] # 年齢:分類1 * 年齢:分類1
age.type.vcov1<- cbind(age.type.vcov1, apply(age.type.vcov1,1,function(x) (-1)*sum(x))) # 共分散行列の最後の列は，(-1) * 他の列の合計
age.type.vcov1<- rbind(age.type.vcov1, apply(age.type.vcov1,2,function(x) (-1)*sum(x))) # 共分散行列の最後の行は，(-1) * 他の行の合計

age.type.vcov2<-res.vcov[16:18,16:18] # 年齢:分類2 * 年齢:分類2
age.type.vcov2<- cbind(age.type.vcov2, apply(age.type.vcov2,1,function(x) (-1)*sum(x))) # 共分散行列の最後の列は，(-1) * 他の列の合計
age.type.vcov2<- rbind(age.type.vcov2, apply(age.type.vcov2,2,function(x) (-1)*sum(x))) # 共分散行列の最後の行は，(-1) * 他の行の合計

age.type.vcov3<-res.vcov[13:15,16:18] # 年齢:分類1 * 年齢:分類2
age.type.vcov3<- cbind(age.type.vcov3, apply(age.type.vcov3,1,function(x) (-1)*sum(x))) # 共分散行列の最後の列は，(-1) * 他の列の合計
age.type.vcov3<- rbind(age.type.vcov3, apply(age.type.vcov3,2,function(x) (-1)*sum(x))) # 共分散行列の最後の行は，(-1) * 他の行の合計

## 全部を合計
## ※ res.vcov[16:18,13:15]の結果はres.vcov[13:15,16:18]と同じなので，ここではres.vcov[13:15,16:18]の結果を２倍して合計している

age.type.vcov4<- age.type.vcov1 + age.type.vcov2 + (age.type.vcov3 *2)

age.type.stderr <- sqrt(age.type.vcov4[4,4]) # 標準誤差は，末尾の分散の正の平方根

こうなると，いくら結果の数値を手入力しなくて済むにせよ，やはり作業途中で凡ミスをやらかしてしまう可能性が高くなります。そこで，次のエントリではこれらの作業を自動化できるようにしてみます。