順序がある
カテゴリカル変数の
解析メモ
よくある状況だけどめんどくさい問題
これまで連続変数を説明変数とする解析が多くなってきていたが、今一度ここで、カテゴリカル変数の解析に思いをはせてみよう。特に一元配置分散分析の時を思い出してもらうとよいが、例えば、3水準あったときにこのカテゴリカル変数には順序の制約がなかった。薬剤無し、薬剤X、薬剤Yみたいな、3水準の場合、通常、これらの間に、効果の強さという意味で、薬剤X > 薬剤Y > 薬剤無しという順序を考える意味はないだろう。なぜなら、3水準は全く別の処理で、同じ軸上で考える意味はないはずだからだ。
次のような状況はどうだろう。ある植物の多さが、照度に与える影響を調べている。けれども、植物の個体数を完全に数え上げることは難しい。そこで、「全くいない」、「少しいる」、「多い」の3段階でざっくり評価することにした。おおざっぱだが、でも確実にこの3水準には順序が存在する。もちろん、がんばって数えなさいよ、というツッコミはあるだろうけど、たくさんのデータをとるためにはやむを得ないこともあるし、まあまあよくある状況である。こんなとき、逆にこれらの3水準を別の処理として扱うのは不自然だ。こんなとき、どうやって対処するべきか。本項では、R内のlm関数が順序制約のあるカテゴリカル変数をデフォルトでどう扱っているかの理解を深めつつ、いくつかの解釈を広げる術を紹介しよう。
順序制約のあるカテゴリカル変数
まず、R内で順序制約のあるカテゴリカル変数がどう表されているかを確認してみよう。以下のように、自分で乱数を生成した。x0は普通のカテゴリカル変数、x1には順序制約を付けた。
------------------------------------------------------
library(plotn)
lv <- c("A", "B", "C", "D")
r <- 30
x <- rep(lv, each = r)
x0 <- factor(x, ordered = FALSE)
x1 <- factor(x, ordered = TRUE)
y <- rep(c(2, 3, 5, 6), each = r) + rnorm(length(x), 0, 0.5)
d <- data.frame(x0, x1, y)
d$x0
## [1] A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A B B B B B B B
## [38] B B B B B B B B B B B B B B B B B B B B B B B C C C C C C C C C C C C C C
## [75] C C C C C C C C C C C C C C C C D D D D D D D D D D D D D D D D D D D D D
## [112] D D D D D D D D D
## Levels: A B C D
d$x1
## [1] A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A B B B B B B B
## [38] B B B B B B B B B B B B B B B B B B B B B B B C C C C C C C C C C C C C C
## [75] C C C C C C C C C C C C C C C C D D D D D D D D D D D D D D D D D D D D D
## [112] D D D D D D D D D
## Levels: A < B < C < D
head(d)
## x0 x1 y
## 1 A A 1.3964671
## 2 A A 2.1387146
## 3 A A 2.5422206
## 4 A A 0.8271511
## 5 A A 2.2145623
## 6 A A 2.2530279
boxplotn(y ~ x0, data = d) #図1の描画
------------------------------------------------------
図1 データの様子
x1の中身を確認してみると、順序としてA < B < C < Dという制約がついている。とりあえず、天下り的に順序制約のないx0と制約のあるx1をそれぞれ説明変数としたときの線形回帰を行ってみよう。
------------------------------------------------------
fit1 <- lm(d$y ~ d$x0)
summary(fit1)
##
## Call:
## lm(formula = d$y ~ d$x0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.02464 -0.29404 -0.08036 0.22515 1.35613
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.85179 0.08351 22.174 < 2e-16 ***
## d$x0B 0.87240 0.11810 7.387 2.5e-11 ***
## d$x0C 3.21751 0.11810 27.244 < 2e-16 ***
## d$x0D 4.20398 0.11810 35.597 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4574 on 116 degrees of freedom
## Multiple R-squared: 0.9348, Adjusted R-squared: 0.9331
## F-statistic: 554 on 3 and 116 DF, p-value: < 2.2e-16
fit2 <- lm(d$y ~ d$x1)
summary(fit2)
##
## Call:
## lm(formula = d$y ~ d$x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.02464 -0.29404 -0.08036 0.22515 1.35613
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.92526 0.04175 94.007 < 2e-16 ***
## d$x1.L 3.34450 0.08351 40.049 < 2e-16 ***
## d$x1.Q 0.05703 0.08351 0.683 0.496
## d$x1.C -0.63311 0.08351 -7.581 9.21e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4574 on 116 degrees of freedom
## Multiple R-squared: 0.9348, Adjusted R-squared: 0.9331
## F-statistic: 554 on 3 and 116 DF, p-value: < 2.2e-16
------------------------------------------------------
fit1の中身は通常の線形回帰の結果である。モデル式は、dB、dC、dDをダミー変数とし、推定したいパラメータをβA、βB-A、βC-A、βD-Aとすれば下記のとおりである。βA、βB-A、βC-A、βD-Aはそれぞれ、「カテゴリAの平均の推定値」「カテゴリBとAの平均の推定値の差」「カテゴリCとAの平均の推定値の差」「カテゴリDとAの平均の推定値の差」になる。
それと比べると順序制約のあるx1で線形回帰した結果は、それとは異なった結果となっている。これは何を表しているのだろうか。その詳細に立ち入る前に、対比contrastという考え方を紹介する必要がある。
対比contrastは「カテゴリ間の比較を表す方法」
まず、対比の定義は下記のようにあらわす。
ここで非常にわかりにくいのがciである。これは一度具体的に当てはめて考えたほうがわかりやすい。例えば、カテゴリA~Dの平均をmA~mDとあらわし、それに対応する「係数」（推定するパラメータと区別して、こう表記する）をcA~cDとあらわすとしよう。(cA, cB, cC, cD) = (1, 0, 0, 0)のとき、対比はmAとなり、カテゴリAの平均そのもの、つまりβAである。では、(cA, cB, cC, cD) = (-1, 1, 0, 0)のときには、対比はmB - mAとなり、「カテゴリBとカテゴリAの平均値の差」、つまりβB-Aである。以上を考慮すれば、この「係数」の行列と各カテゴリの平均を使えば、推定するパラメータは以下のようにあらわせる。
このときの「係数」の行列の逆行列は掃き出し法を使って簡単に求めることができて、パラメータに逆行列を掛ければ、各カテゴリの平均値を計算することができる。この逆行列は丁度、ダミー変数を行列としてまとめたものと同じになっていることがわかるだろう。
こういうスタンダードな方法以外にも、例えば、(cA, cB, cC, cD) = (1/2, 1/2, -1, 0)とすれば、「（カテゴリAとBの平均）とカテゴリCの平均の差」ということも表現可能である。特に、「係数」の和＝Σci= 0の時に限って、対比と呼ぶことが多い（粕谷先生のページ）。
もう少し、この対比という値について深掘りしよう。今、「係数」と各カテゴリの平均の共分散（をカテゴリの数だけかけたもの）を考えてみると、その値は対比と一致することがわかる。
さらに上記の式をciの分散で割ることを考える。この時、もし、「係数」の平方和＝Σci^2 = 1となるように標準化していたとすると、対比は、cを説明変数、mを被説明変数としたときの線形回帰の傾きと一致することがわかる。
上記の議論に立ち戻って、対比は平均値の差を表していたことからも、この結果は妥当なものである。別の見方を考えれば、対比はciとmiの共分散を表す指標でもあるわけなので、「平均値の挙動が、与えられた「係数」とどれくらい相関しているかの指標」、でもあるわけである。この考えが、次の順序制約のあるカテゴリカル変数の解釈に必要になる。
順序制約カテゴリカル変数の解析は多項式対比
では、順序制約があるときのカテゴリカル変数の解析結果は、どんなことをやっているのだろうか。それは、多項式対比polynomial contrastと呼ばれる方法である。「係数」を定めることで、対比の形を決めることができたので、ciの決め方から紹介するのが良いだろう。
ざっくりいうと、Σci= 0およびΣci^2 = 1の制約を守りつつ、1次関数、2次関数、3次関数……に従うciを定めている。Rではそのような「係数」をcontr.poly関数で生成できる。これを使って図示してみよう。この関数は、引数にカテゴリの数をとる。今回は、とりあえず4を代入する。
------------------------------------------------------
contr.poly(4)
## .L .Q .C
## [1,] -0.6708204 0.5 -0.2236068
## [2,] -0.2236068 -0.5 0.6708204
## [3,] 0.2236068 -0.5 -0.6708204
## [4,] 0.6708204 0.5 0.2236068
plotn(contr.poly(4)[,1]) #図2の描画
plotn(contr.poly(4)[,2]) #図3の描画
plotn(contr.poly(4)[,3]) #図4の描画
------------------------------------------------------
図2 1次の「係数」
図3 2次の「係数」
図4 3次の「係数」
生成された「係数」の列に関して、L、Q、Cはそれぞれ1次～3次を表している（linear, quadric, cubic）。それぞれの「係数」が1次～3次式に従っているのだが、ちょっとわかりにくい。とりあえず、20を代入して確認してみよう。
------------------------------------------------------
plotn(contr.poly(20)[,1]) #図5の描画
plotn(contr.poly(20)[,2]) #図6の描画
plotn(contr.poly(20)[,3]) #図7の描画
plotn(contr.poly(20)[,4]) #図8の描画
------------------------------------------------------
図5 1次の「係数」
図6 2次の「係数」
図7 3次の「係数」
図8 4次の「係数」
たしかに、何やらよくはわからないが、1～3次関数に従った値が生成されているようだ。さらに、その「係数」の和は0で、平方和は1となる。
------------------------------------------------------
sum(contr.poly(4)[,1])#ほぼ0
## [1] 0
sum(contr.poly(4)[,2])#ほぼ0
## [1] 0
sum(contr.poly(4)[,3])#ほぼ0
## [1] -2.775558e-17
sum(contr.poly(4)[,1]^2)
## [1] 1
sum(contr.poly(4)[,2]^2)
## [1] 1
sum(contr.poly(4)[,3]^2)
## [1] 1
------------------------------------------------------
この「係数」と、各カテゴリの平均を使って、パラメータは下記のようにあらわされている。
実際、下記のように計算可能である。
------------------------------------------------------
fit2 <- lm(d$y ~ d$x1)
summary(fit2)
##
## Call:
## lm(formula = d$y ~ d$x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.02464 -0.29404 -0.08036 0.22515 1.35613
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.92526 0.04175 94.007 < 2e-16 ***
## d$x1.L 3.34450 0.08351 40.049 < 2e-16 ***
## d$x1.Q 0.05703 0.08351 0.683 0.496
## d$x1.C -0.63311 0.08351 -7.581 9.21e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4574 on 116 degrees of freedom
## Multiple R-squared: 0.9348, Adjusted R-squared: 0.9331
## F-statistic: 554 on 3 and 116 DF, p-value: < 2.2e-16
sum(tapply(d$y, d$x0, mean) * 1/4 * c(1,1,1,1))
## [1] 3.925261
sum(tapply(d$y, d$x0, mean) * contr.poly(4)[,1])
## [1] 3.344496
sum(tapply(d$y, d$x0, mean) * contr.poly(4)[,2])
## [1] 0.05703022
sum(tapply(d$y, d$x0, mean) * contr.poly(4)[,3])
## [1] -0.6331095
------------------------------------------------------
対比とは「係数」とカテゴリの平均の相関の程度を表すことを思い出せば、1～3次の対比は、それぞれ1～3次関数との相関の程度を表していると言える。つまり、n個の順序制約付きのカテゴリカル変数を説明変数に持つ線形回帰で推定されるパラメータは、各カテゴリの平均と1～n-1次関数との相関の程度を推定している、ということである。
この解析は、各カテゴリが等間隔に並んでいるとき、順序が大きくなっていくときのトレンドを検出している、と理解できる。例えば、1次対比の値が正であれば、順序が大きくなっていくときに平均値が上昇する傾向を検出できる。負であれば逆に減少する傾向を検出できるだろう。2次対比が正であれば、真ん中のカテゴリが小さく、小さいor大きいカテゴリは大きな値を示す傾向、負であればその逆の傾向を検出できる。3次対比が正であれば、順序の小さい方から小→大→小→大という風に変動していく傾向を……というように、カテゴリ間の平均の変化を多項式で表現しているわけである。
ちなみに「係数」の行列の逆行列を求めれば、パラメータから各カテゴリの平均値を逆算できる。下記のように計算するとわかるが、「係数」の行列の1行目を1におきかえ、その行列を転置したものになっている。
------------------------------------------------------
x <- rbind(1/4 * c(1,1,1,1), t(contr.poly(4)))
solve(x)
## .L .Q .C
## [1,] 1 -0.6708204 0.5 -0.2236068
## [2,] 1 -0.2236068 -0.5 0.6708204
## [3,] 1 0.2236068 -0.5 -0.6708204
## [4,] 1 0.6708204 0.5 0.2236068
------------------------------------------------------
しかしながら、多項式対比は各カテゴリが等間隔に並んでいることを仮定している。線形のトレンドを検出できる点で便利であるが、単純に隣接するカテゴリ間だけで大小関係を知りたい場合もあるだろう。そこで、次のような方法を使えることを紹介しよう。
累積型ダミー変数による解析
Dobson and Barnett (2018)による"An Introduction to Generalized Linear Models Fourth Edition"のp.43には下記のようなダミー変数を使った解析が紹介されている。
上段の行列はダミー変数をまとめた行列、下段の行列は対比の「係数」をまとめた行列である。ダミー変数をまとめた行列から掃き出し法を使えば、その逆行列、つまり「係数」の行列を作成することができる。上記のモデル構造はA→B→C→Dという順序がある中で、その隣との差を検出するものになっていることがわかるだろう。
このモデル構造を指定しながら、解析をするには以下のように行う。オブジェクトmmのなかは長いので各自確認してみるとよいが、x0がカテゴリBの値のときは、B列に1、カテゴリCの時はBとC列に1、というようにデータが格納されている。なお、カテゴリAを指定してないが、これはlm関数で自動で切片項が挿入され、そこで推定されるため問題ない。
------------------------------------------------------
mm <- matrix(data = 0, nrow = nrow(d), ncol = length(unique(d$x0))-1)
mm[d$x0 == "B", 1] <- 1
mm[d$x0 == "C", 1:2] <- 1
mm[d$x0 == "D", 1:3] <- 1
colnames(mm) <- c("B", "C", "D")
fit3 <- lm(d$y ~ mm)
summary(fit3)
##
## Call:
## lm(formula = d$y ~ mm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.02464 -0.29404 -0.08036 0.22515 1.35613
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.85179 0.08351 22.174 < 2e-16 ***
## mmB 0.87240 0.11810 7.387 2.50e-11 ***
## mmC 2.34511 0.11810 19.857 < 2e-16 ***
## mmD 0.98646 0.11810 8.353 1.64e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4574 on 116 degrees of freedom
## Multiple R-squared: 0.9348, Adjusted R-squared: 0.9331
## F-statistic: 554 on 3 and 116 DF, p-value: < 2.2e-16
------------------------------------------------------
この結果のパラメータの推定値は、上記で紹介したとおり、今注目しているカテゴリの平均と直前のカテゴリの平均の差になっている。
------------------------------------------------------
sum(tapply(d$y, d$x0, mean) * c(1,0,0,0))
## [1] 1.851788
sum(tapply(d$y, d$x0, mean) * c(-1,1,0,0))
## [1] 0.8724033
sum(tapply(d$y, d$x0, mean) * c(0,-1,1,0))
## [1] 2.345109
sum(tapply(d$y, d$x0, mean) * c(0,0,-1,1))
## [1] 0.9864638
------------------------------------------------------
この解析では、カテゴリ間の距離について等しいことを仮定していないので、平均値が一貫して増加しているかなどを確認することが可能である。多項式対比のように、線形トレンドという1つの数値で要約はされないため解釈は難しい面もあるが、覚えておくと便利な解析の一つである。