假设我正在使用multcomp
软件包对R(Bretz et al. (2011) Multiple Comparisons with R. Chapman & Hall/CRC)进行剂量反应分析。 (模拟的)数据如下:
dat <- data.frame(Group = rep(c(0, 0.3, 0.7, 1.2, 1.8, 2.5), each = 5),
Response = c(rnorm(5, 20, 1.2),
rnorm(5, 19.5, 1.2),
rnorm(5, 19, 1.2),
rnorm(5, 15, 1.2),
rnorm(5, 12, 1.2),
rnorm(5, 11, 1.2)
)
)
第一步,我想确定反应是否存在下降趋势(即平均反应水平是否随剂量增加而下降)。这可以通过一项特别强大的测试-Williams测试(Williams (1971) A Test for Differences between Treatment Means When Several Dose Levels are Compared with a Zero Dose Control. Biometric 27:103-117)来完成。这是R代码:
dat$Group = factor(dat$Group)
M <- lm(Response ~ Group, data = dat)
trend = glht(M, linfct = mcp(Group = "Williams"), alternative = "less")
summary(trend)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Williams Contrasts
Fit: lm(formula = Response ~ Group, data = dat)
Linear Hypotheses:
Estimate Std. Error t value Pr(<t)
C 1 >= 0 -7.8117 0.7587 -10.296 < 1e-08
C 2 >= 0 -7.3490 0.6571 -11.184 < 1e-08
C 3 >= 0 -6.5282 0.6195 -10.538 < 1e-08
C 4 >= 0 -5.1096 0.5998 -8.519 < 1e-08
C 5 >= 0 -4.3293 0.5877 -7.366 3.01e-08
(Adjusted p values reported -- single-step method)
测试对比中最低的p值
但是,下一步,我想确定所谓的“最小有效剂量(MED)”,即对我的反应变量有重大影响的最低剂量水平。威廉姆斯(Williams)在他的原始出版物中建议依次应用许多t样测试(即比较要控制的第二高剂量,然后比较要控制的第三高剂量,依此类推),并在第一个无关紧要的结果处停止该程序。先前的重要结果将对应于MED。不幸的是,尽管与t统计量相似,但Williams提出的检验统计量在无剂量反应的无效假设下并未遵循标准的t分布。在他的原始论文中,作者确实列出了其测试统计量的一些关键值。
但是,我想知道这样的顺序Williams测试是否有R实现。可以使用
multcomp
包以某种方式完成(例如,通过以某种方式指定对比度)吗?我花了很多时间在网上寻找答案,但不得不放弃。任何帮助将不胜感激。 最佳答案
您可以使用Dunnett对比来找到MED。
也许是“加强威廉姆斯”的对比:
# step-up Williams contrast matrix
n <- tapply(dat$Group, dat$Group, length)
k <- length(n)
CM <- c()
for (i in 1:(k - 1)) {
help <- c(-1, n[2:(i + 1)] / sum(n[2:(i + 1)]), rep(0 , k - i - 1))
CM <- rbind(CM, help)
}
rownames(CM) <- paste("C", 1:nrow(CM))
CM
# supply to glht()
summary(glht(M, linfct = mcp(Group = CM), alternative = "less"))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: User-defined Contrasts
Fit: lm(formula = Response ~ Group, data = dat)
Linear Hypotheses:
Estimate Std. Error t value Pr(<t)
C 1 >= 0 0.1535 0.7630 0.201 0.7214
C 2 >= 0 -0.2032 0.6608 -0.308 0.5259
C 3 >= 0 -1.5409 0.6230 -2.473 0.0219 *
C 4 >= 0 -3.2164 0.6032 -5.332 <0.001 ***
C 5 >= 0 -4.2203 0.5910 -7.141 <0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)