推 cywhale: very helpful to me, thanks for sharing! 09/01 21:50
※ 編輯: andrew43 (122.117.37.172), 09/01/2016 23:15:46
→ andrew43: 做了一點小改動，但不影響主要內容。 09/01 23:16
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() 對組別的名字有點敏感，
在上述的例子中已有提醒，請注意。
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 122.117.37.172
※ 文章網址: https://www.ptt.cc/bbs/R_Language/M.1472728842.A.ED1.html