1. 套件名称：
multcomp_1.4-6
multcompView_0.1-7
2. 套件主要用途：
multcomp 的主要工作是：
在建立模型后，
利用调整某因子的编码，
达到特定的对比组合或所有两两对比（即事后比较）。
其中可选用不同的 p-value 校正方式。
multcompView的主要工作是：
所有两两对比完成后，
利用所截取出的 p-value 或是否显著的 boolean 值自动进行分群，
分群的结果是给予各组 abcde……之类的分群编码。
3. 套件主要函数列表：
a. multcomp::glht()
在模型中套用特定的对比组合或所有两两对比
b. multcompView::multcompLetters()
按 p-value 或是否显著 boolean 进行分群
4. 分享内容：
## 建立一个 one-way Poisson regression 模型
## TRT 为自变量，共四组，平均分别为 1、1.01、2、4
dt <- data.frame(
FREQ = c(rpois(10, 1), rpois(10, 1.01), rpois(10, 2), rpois(10, 4)),
TRT = gl(4, 10, labels = c("GroupA", "GroupB", "GroupC", "GroupD"))
)
poisson.model <- glm(FREQ ~ TRT, data = dt, family = poisson)
######## CASE 1
## 利用 multcomp::glht() 进行自己设计的
## (GroupB - GroupC = 0)
## (GroupD - GroupA = 1)
## (GroupC = 0)
## 三组对比与检验，并未考虑 p-value 校正
## 如果想采用更复杂的对比，例如 A-(B+C)/2 = 0，请见 help(glht)
## 如果模型中有其它自变量，会一并纳入考虑，并改变指定因子的对比
poisson.model.p1 <- glht(
poisson.model,
linfct = mcp(TRT = c(
"GroupB - GroupC = 0",
"GroupD - GroupA = 1",
"GroupC = 0"
))
)
summary(poisson.model.p1)
######## CASE 2
## 利用 multcomp::glht() 达到所有两两事后比较
## 注意！经验告诉我此时 TRT 的组别名称不要有 space 或怪怪的字符，
## 纯英数最好，不然会卡关
poisson.model.pairwise <-
glht(poisson.model, linfct = mcp(TRT = "Tukey"))
## 配合 Bonferroni p-value 校正
## 校正选项可另见 help(p.adjust)
poisson.model.pairwise.summary <-
summary(
poisson.model.pairwise, test = adjusted("bonferroni")
)
print(poisson.model.pairwise.summary)
## 把已校正的 p-value 截取出来
p.vals <- poisson.model.pairwise.summary$test$pvalues
## 注意！p.vals 的名字不要有 space 或怪怪的字符，纯英数最好，不然会卡关
## 组和组之间只有一个 "-" 来隔开就好
names(p.vals)
names(p.vals) <- gsub(" ", "", names(p.vals))
names(p.vals)
## 利用 multcompView() 以 p < 0.05 进行分群
multcompLetters(p.vals)
## 利用 multcompView() 以 p < 0.01 进行分群
## 此时输入的是 boolean
multcompLetters(p.vals < 0.01)
## 同上，但借由改变 threshold
multcompLetters(p.vals, threshold = 0.01)
5. 备注
a. multcomp 和 multcompView 都是可以独立进行的，不相依。
只不过为了演示，我同时在此一并介绍。
b. glht() 和 multcompView() 对组别的名字有点敏感，
在上述的例子中已有提醒，请注意。