已知資料 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 相當接近.
