在 R 中配置 Generalized Logit Model

已知資料 data.csv 如下:

treat, r.a, r.b, r.c
   t1,  10,  20,   5
   t2,  11,  22,   5
   t3,  21,  40,  11
   t4,  21,   5,  10

其中 treat 為具有 4 個水準的固定因子; r.a, r.b 及 r.c 分別為 3 個類型的反應類別, 其下為頻率.

以下為 R code:

# 系統第一次運作請先安裝 package nnet
install.packages("nnet")

# 引入 package nnet
require(nnet)

# 讀取資料
dat <- read.csv("data.csv")

# 建立 generalized logit model
m <- multinom(
  cbind(r.a, r.b, r.c) ~ treat,          # 其中 r.a 寫在最前面即可作為比較基準
  data=dat,
  contrasts=list(treat=contr.treatment), # 其中採用 treatment contrast 方便解釋
  Hess=T                                 # 其中採用 Hessian Matrix 求標準誤
)

# 係數與標準誤
(s.m <- summary(m))

# odds ratio
exp(s.m$coefficients)

# chi-square
(chisq <- (s.m$coefficients / s.m$standard.error)^2)

# p-value
pchisq(chisq, df=1, lower.tail=F)

可得到以下結果:

> (s.m <- summary(m))
Call:
multinom(formula = cbind(r.a, r.b, r.c) ~ treat, data = dat, 
    contrasts = list(treat = contr.treatment), Hess = T)

Coefficients:
    (Intercept)        treat2      treat3      treat4
r.b   0.6931308  2.158365e-05 -0.04877868 <strong>-2.12818025</strong>
r.c  -0.6931750 -9.525513e-02  0.04656602 -0.04876199

Std. Errors:
    (Intercept)    treat2    treat3    treat4
r.b   0.3872968 0.5351291 0.4718241 0.6305649
r.c   0.5477240 0.7687048 0.6622145 0.6690451

Residual Deviance: 346.6049 
AIC: 362.6049 


> # odds ratio
> exp(s.m$coefficients)
    (Intercept)   treat2    treat3    treat4
r.b   1.9999671 1.000022 0.9523919 <strong>0.1190537</strong>
r.c   0.4999861 0.909141 1.0476672 0.9524078


> # chi-square
> (chisq <- (s.m$coefficients / s.m$standard.error)^2)
    (Intercept)       treat2      treat3       treat4
r.b    3.202894 1.626795e-09 0.010688081 <strong>11.390878935</strong>
r.c    1.601630 1.535527e-02 0.004944712  0.005311928


> # p-value
> pchisq(chisq, df=1, lower.tail=F)
    (Intercept)    treat2    treat3      treat4
r.b  0.07350811 0.9999678 0.9176589 <strong>0.000738056</strong>
r.c  0.20567234 0.9013815 0.9439400 0.941899243

可以發現, 唯一達顯著之係數為 r.b 對比 r.a 時在 treat4 對比 treat1 的情況, 其係數 $\beta = -2.12818025$, 勝算比 $= \exp(-2.12818025) = 0.1190537$, 檢定量 $\chi^2_1 = 11.390878935$, $p = 0.000738056$, 顯示 treat4 比 treat1 的 r.b / r.a 達顯著差異. 由實際資料 $\frac{5/21}{20/10} = 0.1190476$ 也和模型所預測的勝算比 0.1190537 相當接近.