2014年6月10日火曜日

累積ロジットとGLM二項分布の比較・再&続

(うっかり、同様の検証記事を消してしまったので、ついでにアップデートします)
前回に引き続き、段階的なカテゴリーデータのモデリングに用いられる累積ロジット(cumulative logit)、
例えば、悪い、ふつう、よい、のようなデータを関連しそうな要因で解析するような場合に用いる。

しかし、悪い=0、ふつう=1、よい=2、のように数値化してしまえば、二項分布型のGLMではダメなのだろうか?たぶん、間隔のいびつなデータでは累積ロジットの方が適しているのだろうが。

テストデータを用いて、両者の推定と推定値の求め方を比較する。

# まず解析用データの設定
logistic <- function(xx) 1 / (1 + exp(-xx))
N <- 3 # 最大3、つまり 0, 1, 2, 3 の値を取りうる
X <- rep(c(1:10), each=10)
Y1 <- rbinom(100, N, prob=logistic(-5 + 0.8*X)) # logisticの中身を基にして、二項乱数を発生
Y2 <- factor(Y1, ordered=T) # 累積ロジット用に、ランク化した応答変数を用意
D <- data.frame(X, Y1, Y2)

# 解析の実行
require(ordinal)
require(VGAM)
M1 <- glm(cbind(Y1, N-Y1) ~ X, family=binomial, data=D) # 二項分布のGLM
M2 <- clm(Y2 ~ X, data=D) # 累積ロジット
M3 <- vglm(Y2 ~ X, family=cumulative(parallel=T), data=D)
 # 比較用にvglm版、parallel=Tで切片のみ複数になる、Fにすると回帰係数もランクごとに推定(切片のみランク毎の場合、比例オッズモデルともいうようだ)

### GLM
summary(M1)
# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  -5.9324     0.6644  -8.929   <2e-16 ***
# X             0.9241     0.1011   9.141   <2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# (Dispersion parameter for binomial family taken to be 1)
#
#     Null deviance: 287.723  on 99  degrees of freedom
# Residual deviance:  84.861  on 98  degrees of freedom
# AIC: 139.14

### 累積ロジット(clm)
summary(M2)
#  link  threshold nobs logLik AIC    niter max.grad cond.H 
#  logit flexible  100  -67.06 142.12 5(0)  2.20e-08 2.6e+03
#
# Coefficients:
#   Estimate Std. Error z value Pr(>|z|)    
# X   1.2295     0.1659    7.41 1.26e-13 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Threshold coefficients:
#     Estimate Std. Error z value
# 0|1    5.955      0.913   6.522
# 1|2    8.082      1.130   7.153
# 2|3    9.716      1.283   7.575

### 累積ロジット(vglm)# 回帰係数の正負が逆になる
summary(M3)
# Coefficients:
#              Estimate Std. Error z value
# (Intercept):1   5.9548    0.92423  6.4430
# (Intercept):2   8.0821    1.14814  7.0393
# (Intercept):3   9.7158    1.31086  7.4117
# X                  -1.2295    0.16832 -7.3046 
#
# Number of linear predictors:  3 
#
# Names of linear predictors: logit(P[Y<=1]), logit(P[Y<=2]), logit(P[Y<=3])
#
# Dispersion Parameter for cumulative family:   1
#
# Residual deviance: 134.1208 on 296 degrees of freedom # 過小分散気味、clmには無い情報!
#
# Log-likelihood: -67.06038 on 296 degrees of freedom


### 推定値を得るには少し手間がかかる
pre1 <- round(3*fitted(M1)) # GLM:試行回数*確率、を整数値に丸める
pre2 <- as.numeric(predict(M2,type="class")$fit) - 1 
# 累積ロジット(clm):predict関数で推定ランクを求め、
# ラベルに合うように 1 を引く(ランクは1,2,3,4、元の値は0,1,2,3なので)
pre3 <- apply(fitted(M3), 1, which.max) - 1
# 累積ロジット(vglm):何番目のランクの確率が最大かを求め、
# ラベルに合うように 1 を引く

cbind(X, Y1, pre1, pre2, pre3) 
# 推定は完全ではないが、GLMと累積ロジットの推定結果は近いものになった。
     X Y1 pre1 pre2 pre3
1    1  0    0    0    0
2    1  0    0    0    0
3    1  0    0    0    0
4    1  0    0    0    0
5    1  0    0    0    0
(中略)
51   6  1    1    1    1
52   6  1    1    1    1
53   6  2    1    1    1
54   6  2    1    1    1
55   6  2    1    1    1
56   6  1    1    1    1
57   6  2    1    1    1
58   6  1    1    1    1
59   6  1    1    1    1
60   6  0    1    1    1
61   7  3    2    2    2
62   7  1    2    2    2
63   7  1    2    2    2
64   7  3    2    2    2
65   7  3    2    2    2
66   7  1    2    2    2
67   7  3    2    2    2
68   7  2    2    2    2
69   7  1    2    2    2
70   7  3    2    2    2
71   8  3    2    3    3
72   8  2    2    3    3
73   8  3    2    3    3
74   8  3    2    3    3
75   8  2    2    3    3
76   8  1    2    3    3
77   8  1    2    3    3
78   8  3    2    3    3
79   8  3    2    3    3
80   8  3    2    3    3
81   9  2    3    3    3
82   9  1    3    3    3
83   9  3    3    3    3
(後略)


### では、推定確率の曲線を描くにはどうしたらよいか?
# これがけっこう面倒くさい。まずはvglmを用いてチェックする。

# vglmをfitted()すると、ランク毎の確率が出てくる
head(fitted(M3)) # head()で頭出しする
                       0                     1                      2                  3
1   0.991209875 0.007734525 0.0008493615 0.0002062383
2   0.991209875 0.007734525 0.0008493615 0.0002062383
3   0.991209875 0.007734525 0.0008493615 0.0002062383
4   0.991209875 0.007734525 0.0008493615 0.0002062383
5   0.991209875 0.007734525 0.0008493615 0.0002062383
6   0.991209875 0.007734525 0.0008493615 0.0002062383

# clmでは、
M2@fitted.values で同様にランク毎の確率が得られる。



# まず、この出来合いの推定値を図示してみる
plot(X, fitted(M3)[,1], xlim=c(0,10), ylim=c(0,1), col="lightblue", ylab="probability") # 0
points(X, fitted(M3)[,2], col="turquoise") # 1
points(X, fitted(M3)[,3], col="royalblue") # 2
points(X, fitted(M3)[,4], col="darkblue") # 3

# モデルの構造を考えると推定確率曲線は次の計算でいいはず
#(coef(M3)[4]は回帰係数だが、正負が逆なので - を付ける)
# 累積ロジットの名前通り、累積で表されている
curve(1 - logistic(-coef(M3)[4]*x - coef(M3)[1]), 
  add=T, lwd=2, col="lightblue") # ランク0: 1 - ランク1の累積確率
curve(logistic(-coef(M3)[4]*x - coef(M3)[1]) - logistic(-coef(M3)[4]*x - coef(M3)[2]), 
add=T, lwd=2, col="turquoise") # ランク1: ランク1の累積確率 - ランク2の累積確率
curve(logistic(-coef(M3)[4]*x - coef(M3)[2]) - logistic(-coef(M3)[4]*x - coef(M3)[3]),
add=T, lwd=2, col="royalblue") # ランク2: ランク2の累積確率 - ランク3の確率
curve(logistic(-coef(M3)[4]*x - coef(M3)[3]), 
  add=T, lwd=2, col="darkblue") # ランク3: ランク3の確率(もはや累積でない)

curve(logistic(coef(M1)[1] + coef(M1)[2]*x), add=T, col="tomato", lwd=2) # 比較用にGLMも


# 色が薄い〜濃いにかけて、それぞれ、0, 1, 2, 3 になる確率、赤はGLMの確率

# ちゃんと関数による推定プロットと、自前で作った推定曲線が一致した。計算の仕方はこれでよさそうだ。

# ちなみにclmの場合はこう計算する(回帰係数の前の - が不要)
curve(1 - logistic(coef(M2)[4]*x - coef(M2)[1]), add=T, lwd=2, col="lightblue") # 0
curve(logistic(coef(M2)[4]*x - coef(M2)[1]) - logistic(coef(M2)[4]*x - coef(M2)[2]), 
add=T, lwd=2, col="turquoise") # 1
curve(logistic(coef(M2)[4]*x - coef(M2)[2]) - logistic(coef(M2)[4]*x - coef(M2)[3]), 
add=T, lwd=2, col="royalblue") # 2
curve(logistic(coef(M2)[4]*x - coef(M2)[3]), add=T, lwd=2, col="darkblue") # 3


### 次に、推定値の曲線を図示してみる(上のは推定確率でした)
# 比較対照用にGLMの曲線と比較する
# 元データ
plot(Y1 ~ X, data=D, xlim=c(0,10), ylim=c(0,3))
 # GLM
curve(3*logistic(coef(M1)[1] + coef(M1)[2]*x), add=T, col="red", lwd=2)

# 累積ロジット(全部足し合わせるような累積構造になる)
curve(
   1*(1 - logistic(coef(M2)[4]*x - coef(M2)[1])) # ランク0
+ 2*(logistic(coef(M2)[4]*x - coef(M2)[1]) - logistic(coef(M2)[4]*x - coef(M2)[2])) # ランク1
+ 3*(logistic(coef(M2)[4]*x - coef(M2)[2]) - logistic(coef(M2)[4]*x - coef(M2)[3])) # ランク2
+ 4*(logistic(coef(M2)[4]*x - coef(M2)[3])) # ランク3
 - 1, # 0, 1, 2, 3なので、1,2,3,4から1を引く
add=T, col="blue", lwd=2)
# 赤:GLM、青:累積ロジット。この単純なケースではほとんど推定値は変わらない。


#### ついでにベイズ版でも計算
# 下準備、データを累積に変換。1以上、2以上、3以上にまとめる
rank <- matrix(0, nrow=100, ncol=3) # 3はランクの数-1、100はいわゆるN数
for (r in 1:3) rank[Y1 >= r, r] <- 1 
# GLM用の用意していた数値データY1を利用、ランク1~3以上の場合に各列に1を入れる

model <- function() {
  for (j in 1:3) { # ランクの数-1
  for (i in 1:100) { # いわゆるN数
  rank[i,j] ~ dbern(p[i,j]) # それぞれベルヌーイ分布
 logit(p[i,j]) <- alpha[j] + beta*X[i] # 切片だけランク毎になってる
 }
  alpha[j] ~ dnorm(0, 1E-6) # 切片をランク毎に推定
  }
  beta ~ dnorm(0, 1E-6) # 回帰係数は共通
}

data <- list(X=X, rank=rank)
parm <- list(beta=0, alpha=rnorm(3))
source("WBUGS.R") # 自作ラッパー関数
out <- wbugs(data, parm, model)

# 3 chains, each with 11000 iterations (first 1000 discarded), n.thin = 10
# n.sims = 3000 iterations saved
#                 mean    sd    2.5%     25%     50%     75%   97.5%  Rhat n.eff
# beta         1.328 0.166   1.022   1.217   1.318   1.434   1.680 1.001  2200
# alpha[1]  -6.377 0.889  -8.217  -6.960  -6.327  -5.770  -4.703 1.002  1100
# alpha[2]  -8.645 1.142 -11.040  -9.370  -8.563  -7.851  -6.528 1.002  1800
# alpha[3] -10.518 1.336 -13.331 -11.373 -10.420  -9.585  -8.139 1.001  3000
# deviance 154.787 2.934 151.200 152.600 154.050 156.200 162.000 1.002  1500
# DIC info (using the rule, pD = Dbar-Dhat)
# pD = 4.0 and DIC = 158.8

# BUGS版では、切片が正負が逆になった
# きちんとコードを書いた結果がこれなので、符号が逆なのはvglmやclmなのだが…じつにややこしい。

2014年6月9日月曜日

累積ロジットの汎用Rパッケージ {ordinal}

累積ロジット(cumulative logit model)を使う際に、今ひとつ使い勝手のいいRパッケージがないのが気になっていた。
(cf. 累積ロジット:よい、ふつう、わるい、のような段階的な現象についての推定に用いる。等間隔に近ければ二項分布のGLMで構わないだろうけれど、そうでない場合にはこれが適しているようだ)

例えば、vglmではランダム効果が入れられない。

mixcatでは、ランダム効果が入れられるが、対数尤度までは出力できても、AICなど情報量規準は出してくれない(手計算すれば良い話ではあるが…)。

最近、ordinalというパッケージを見つけた。1つのパッケージでGLM版、GLMM版の両方が含まれているし、使用法もglm()やglmer()と同じようにしてくれていて使いやすい。

ちなみに、GLM版にstepAICは使えたが、dredgeはダメだった。
GLMM版では、モデル選択関数はdrop1が使用できた。

require(ordinal)
example(ordinal) # パッケージで用意されている例を出力

# まずはGLM版
ordinl> ## A simple cumulative link model:
ordinl> fm1 <- clm(rating ~ contact + temp, data=wine)

ordinl> summary(fm1)
formula: rating ~ contact + temp
data:    wine

 link  threshold nobs logLik AIC    niter max.grad cond.H
 logit flexible  72   -86.49 184.98 6(0)  4.01e-12 2.7e+01

Coefficients:
           Estimate Std. Error z value Pr(>|z|)    # 回帰係数部分(この例ではカテゴリカルだが)
contactyes   1.5278     0.4766   3.205  0.00135 **
tempwarm     2.5031     0.5287   4.735 2.19e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:           # 各ランクの切片
    Estimate Std. Error z value
1|2  -1.3444     0.5171  -2.600
2|3   1.2508     0.4379   2.857
3|4   3.4669     0.5978   5.800
4|5   5.0064     0.7309   6.850


# 比較対照用にvglmでの結果
require(VGAM)
vm1 <- vglm(rating ~ contact + temp, family=cumulative(parallel=T), data=wine)

Coefficients:
Estimate Std. Error z value
(Intercept):1 -1.3444 0.50850 -2.6438
(Intercept):2 1.2508 0.43908 2.8487
(Intercept):3 3.4669 0.59711 5.8061
(Intercept):4 5.0064 0.72906 6.8669
contactyes -1.5278 0.47362 -3.2258 # 回帰係数部分の正負がclmと逆なことに注意
tempwarm -2.5031 0.53199 -4.7052


# ordinalに戻って、こちらはGLMM版
ordinl> ## A simple cumulative link mixed model:
ordinl> fmm1 <- clmm(rating ~ contact + temp + (1|judge), data=wine) # glmer()同様の構造

ordinl> summary(fmm1)
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: rating ~ contact + temp + (1 | judge)
data:    wine

 link  threshold nobs logLik AIC    niter    max.grad cond.H # ちゃんとAICも算出する
 logit flexible  72   -81.57 177.13 331(996) 1.05e-05 2.8e+01

Random effects:                                     # ランダム効果
 Groups Name        Variance Std.Dev.
 judge  (Intercept) 1.279    1.131
Number of groups:  judge 9

Coefficients:
           Estimate Std. Error z value Pr(>|z|)  
contactyes   1.8349     0.5125   3.580 0.000344 ***
tempwarm     3.0630     0.5954   5.145 2.68e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Threshold coefficients:
    Estimate Std. Error z value
1|2  -1.6237     0.6824  -2.379
2|3   1.5134     0.6038   2.507
3|4   4.2285     0.8090   5.227
4|5   6.0888     0.9725   6.261