0

我很难回到 emmeans 包给出的“CL”值。

我的第一个目标是进行一种功率研究：我正在研究植物，我想看看我可以在下个季节改进多少实验设计，以便在我的基因型之间产生更显着的差异。我的基因型是重复的，这被用作阻止因素。我想玩重复次数（假设均值和残差 SE 保持不变）。最后，我想在预算问题和统计能力之间找到一个折衷方案。

我正在考虑查看置信区间 (CI)，记住如果两种基因型的 CI 至少包含一个共同值，那么这两种基因型没有显着差异。所以我想看看在玩重复次数时我能减少多少 CI。

我对使用 Tukey 方法进行多重比较的 CI 计算的猜测是：

μi ± (σ/(2*√(ni)))*q

μi 是第 i 个基因型的平均值

ni 第 i 个基因型的观察次数，

σ 残差标准误差和

q 以 t（基因型的数量）和 dfE（残差的自由度）作为 α = 0.05 的参数的学生化范围统计

这相当于：

μi ± (SEi*q)/2

SEi 是第 i 个基因型的标准误

这是一个包含 6 个基因型和 3 个重复的小示例数据集以及我运行的代码：

library(emmeans)

# import DF
df <- structure(list(Geno = structure(c(6L, 3L, 2L, 1L, 4L, 5L, 2L, 
                                     4L, 5L, 1L, 3L, 6L, 2L, 3L, 1L, 5L, 6L, 4L),
                                   .Label = c("A", "B", "C", "D", "E", "F"),
                                   class = "factor"), variable = c(2.279628571, 3.925157143, 3.089, 2.26, 2.503,
                                                                   2.495114286, 2.867166667, 3.069238095, 3.884285714, 3.409595238, 
                                                                   3.710714286, 1.763142857, 2.865285714, 4.219214286, 3.263452381, 
                                                                   3.359428571, 2.335285714, 2.443), 
                 Rep = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L), 
                                 .Label = c("1", "2", "3"), class = "factor")), .Names = c("Geno", "variable", "Rep"),
            row.names = c(NA, 18L), class = "data.frame")

# Number of observation
N <- nrow(df)

# Number of treatments
t <- nlevels(df$Geno)

# Define the model
model <- lm(data=df, variable ~ Geno + Rep)

# Compute means and confidence intervals
E_means <- as.data.frame(emmeans(model, pairwise~Geno)$emmean)

# Extract Standard Errors
se <- E_means[,"SE"]

# Studentized range statistic
q <- qtukey(0.95, nmeans=t, df=E_means[,"df"])

expected_CI <- se*q/2
emmeans_CI <- ((E_means[,"upper.CL"] - E_means[,"lower.CL"])/2)

print(expected_CI)
print(emmeans_CI)

请注意，大部分时间我都在处理不平衡的数据（所以我对每种基因型都有不同的 SE）。我想首先说明一个简单（=平衡）的情况。

到目前为止，对于我所做的每个测试，expected_CI 和 emmeans_CI 总是不同的（差别不大，但仍然不同），所以我想我的计算方式与 emmeans 包不同。所以这是我的问题：它是如何在 emmeans 包中完成的？

任何帮助将非常感激 ！

4

1 回答 1

3

您似乎将 EMM 与 EMM 的差异混为一谈。此外，您使用的公式仅适用于平衡单向设计。

以您的示例为例，如果您尝试：

emm <- emmeans(model, pairwise ~ Geno)

然后做：

confint(emm$emmeans, adjust = "tukey")
confint(emm$contrasts, adjust = "tukey")

你会注意到一些事情。首先，EMM 本身的 CI 不使用指定的 Tukey 调整（它用 Sidak 替换它），因为该调整仅适用于成对比较。其次，两个摘要具有不同数量的结果（除非恰好有 3 个处理）和不同的标准误。例如，有 4 个处理，有 4 个 EMM，但有 6 个成对比较；对于不平衡设计，可能有 6 个不同的 SE 与这些比较相关联。

如果您想了解成对比较的 Tukey 调整，您需要使用适合成对比较的统计数据 - 第二个摘要。Tukey 调整的 CI 使用估计的差异加或减sqrt(0.5*qtukey(.95, nmeans, df)) * SE，其中SEs 来自成对差异的结果。

于 2018-04-20T21:57:48.690 回答