2

我正在尝试复制用于单向混合效应方差分析的 Tukey 测试的教科书示例（来自 Statistics，William L. Hays p 581-583），但我使用 lme 和 glht 获得的 p 值没有意义

该研究重复测量了四种治疗方法和 10 名受试者

数据

subject=c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 
    6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10)

treatment=c("a1", "a2", "a3", "a4", "a1", "a2", "a3", "a4", "a1", "a2", "a3", 
    "a4", "a1", "a2", "a3", "a4", "a1", "a2", "a3", "a4", "a1", "a2", 
    "a3", "a4", "a1", "a2", "a3", "a4", "a1", "a2", "a3", "a4", "a1", 
    "a2", "a3", "a4", "a1", "a2", "a3", "a4")

response=c(11, 9, 5, 17, 14, 12, 10, 18, 15, 7, 6, 21, 17, 10, 13, 22, 15, 7, 6, 
    15, 14, 8, 13, 22, 9, 6, 7, 15, 17, 11, 10, 19, 10, 13, 14, 23, 12, 
    8, 11, 20)

dataFrame=data.frame(subject, treatment, response)

该模型

library(nlme)
model = lme(response~ treatment,random=~1|subject,data=dataFrame)
anova(model)

            numDF denDF  F-value p-value
(Intercept)     1    27 375.9198  <.0001
treatment       3    27  43.4507  <.0001

这个 F 值与 Hay 的 ( ) 足够接近F = 43.41，我很确定我的模型很好。

图基测试

library(multcomp)
glht.out =glht(model, mcp(treatment="Tukey"))
summary(glht.out)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lme.formula(fixed = response ~ treatment, data = dataFrame, random = ~1 | 
    subject)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)    
a2 - a1 == 0   -4.300      1.006  -4.276   <0.001 ***
a3 - a1 == 0   -3.900      1.006  -3.879   <0.001 ***
a4 - a1 == 0    5.800      1.006   5.768   <0.001 ***
a3 - a2 == 0    0.400      1.006   0.398    0.979    
a4 - a2 == 0   10.100      1.006  10.044   <0.001 ***
a4 - a3 == 0    9.700      1.006   9.647   <0.001 ***

这与书本不一致。Hay's 只报告了一些比较，并给出了 HSD 和均值，而不是 p，但对于 a3-a1 的对比，发现与 an 没有显着差异，HSD = 4.02并且mean = 3.9根据我的计算有p.value = 1-ptukey(3.9*sqrt(10/8.27),4,9)=0.05719563.

R 输出也没有意义，因为应该控制多重比较的 Tukey 检验 p 值远小于配对 t 检验的 p 值（p=0.0142，使用t.test(c(11, 14, 15, 17, 15, 14, 9, 17, 10, 12),c(5, 10, 6, 13, 6, 13, 7, 10, 14, 11), paired=TRUE).

知道我做错了什么以及如何在 R 中正确执行 Tukey 测试吗？

4

1 回答 1

0

我调查了你的问题，这些是我的发现。不幸的是，我没有海伊的书，所以我只关注你发布的信息。

1）关于海伊的书报告的HSD = 4.02。该lme函数使用RELM最大化受限对数似然的拟合方法，而该LM方法最大化对数似然，?lme以获取更多信息。与 ML 方法拟合的模型显示相同的 F 值。

model <- lme(response~treatment, random=~1|subject, data=dataFrame, method="ML")
anova.lme(model)
            numDF denDF  F-value p-value
(Intercept)     1    27 375.9227  <.0001
treatment       3    27  43.4506  <.0001

使用 glht 进行多重比较时，z value更接近书 4.088 报告的结果，但对比仍然发现差异显着。

glht.out <- glht(model, mcp(treatment="Tukey"), alternative="two.sided")
summary(glht.out)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lme.formula(fixed = response ~ treatment, data = dataFrame, random = ~1 | 
    subject, method = "ML")

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)    
a2 - a1 == 0  -4.3000     0.9539  -4.508  < 1e-04 ***
a3 - a1 == 0  -3.9000     0.9539  -4.088  0.00027 ***
a4 - a1 == 0   5.8000     0.9539   6.080  < 1e-04 ***
a3 - a2 == 0   0.4000     0.9539   0.419  0.97520    
a4 - a2 == 0  10.1000     0.9539  10.588  < 1e-04 ***
a4 - a3 == 0   9.7000     0.9539  10.168  < 1e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

2）关于p值的计算。前面的输出显示 p 值是使用single-step method. 该?summary.glht描述解释了在其中实现adjusted并且p.adjust可以使用的任何方法。以下是如何更改它的示例：

summary(glht.out, test = adjusted(type = "single-step"))
summary(glht.out, test = adjusted(type = "none"))

不幸的是，我无法从 glht 输出中重现 p 值。所以我检查了这个TukeyHSD功能。它需要aov模型输出。

aovmodel <- aov(model)
summary(aovmodel)
            Df Sum Sq Mean Sq F value   Pr(>F)    
treatment    3  659.0  219.67   26.95 2.57e-09 ***
Residuals   36  293.4    8.15                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

运行 TukeyHSD 时，我们得到以下输出

TukeyHSD(aovmodel, "treatment")
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = model)

$treatment
      diff       lwr        upr     p adj
a2-a1 -4.3 -7.738482 -0.8615177 0.0093920
a3-a1 -3.9 -7.338482 -0.4615177 0.0210165
a4-a1  5.8  2.361518  9.2384823 0.0003378
a3-a2  0.4 -3.038482  3.8384823 0.9891641
a4-a2 10.1  6.661518 13.5384823 0.0000000
a4-a3  9.7  6.261518 13.1384823 0.0000000

至少对我而言，该代码比来自 multcomp 的代码更容易理解。p 值可以如下重现：

MSE <- 8.15
center <- -3.9/sqrt(MSE/10)  # -4.32002
ptukey(abs(center), 4, 36, lower.tail=FALSE)
[1] 0.02101645

请注意，HSD 现在为 -4.32002，自由度为 Nk = 40 - 4 = 36。p 值为 0.021，略高于配对 t 检验的值。这是重现输出的要点。

希望这可以帮助

于 2018-02-21T10:51:20.783 回答