前回に引き続き、段階的なカテゴリーデータのモデリングに用いられる累積ロジット（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の累積確率
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も
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
+ 4*(logistic(coef(M2)[4]*x - coef(M2)[3])) # ランク3
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
- 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
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なのだが…じつにややこしい。
0 件のコメント:
コメントを投稿