Как вычисляются доверительные интервалы значений при корректировке по методу Тьюки?

Question

У меня проблемы с возвратом к значениям "CL", заданным пакетом emmeans.

Моей первой целью было провести своего рода силовое исследование: я работаю над растениями и хочу посмотреть, насколько я смогу улучшить свой экспериментальный дизайн на следующий сезон, чтобы иметь более существенные различия между моими генотипами. Мои генотипы повторяются, и это используется в качестве блокирующего фактора. Я хочу поиграть с количеством повторений (при условии, что средства и остаточный SE останутся прежними). В конце я хочу найти компромисс между бюджетными вопросами и статистической властью.

Я подумывал взглянуть на доверительные интервалы (КИ), имея в виду, что если КИ двух генотипов охватывают хотя бы одно общее значение, эти два генотипа существенно не отличаются. Поэтому я хочу посмотреть, насколько я могу уменьшить 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 для каждого генотипа). Я хочу сначала прояснить для простого (= сбалансированного) случая.

Ожидаемые значения_CI и emmeans_CI всегда различны для каждого теста, который я до сих пор проводил (не сильно отличается, но все же отличается), поэтому я полагаю, что я не делаю вычисления таким же образом, как пакет emmeans. Итак, вот мой вопрос: как это делается в пакете emmeans?

Любая помощь приветствуется!

Answer 1

Похоже, вы путаете EMM с различиями EMM. Кроме того, используемые вами формулы применимы только к сбалансированным односторонним конструкциям.

На вашем примере, если вы попробуете:

emm <- emmeans(model, pairwise ~ Geno)

а затем сделать:

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

Вы заметите несколько вещей. Во-первых, CI для EMM сами по себе не используют указанную настройку Tukey (она заменяет ее на Sidak), потому что эта настройка подходит только для парных сравнений. Во-вторых, две сводки имеют разное количество результатов (кроме случаев, когда имеется ровно 3 обработки) и разные стандартные ошибки. Например, с 4 обработками, есть 4 EMM, но 6 парных сравнений; и с несбалансированным дизайном, может быть 6 разных SE, связанных с этими сравнениями.

Если вы хотите понять настройку Тьюки для парных сравнений, вам нужно использовать статистику, подходящую для парных сравнений - это второе резюме. Скорректированные на Тьюки КИ используют расчетные различия плюс или минус sqrt(0.5*qtukey(.95, nmeans, df)) * SE, где SEs приходят из результатов для парных различий.