Неравная дисперсия в R-линейной смешанной модели с апостериорными средними значениями

Question

Я подогнал линейную смешанную модель к сельскохозяйственным данным, учитывая неравные различия между группами (культурными сортами), пропустив weights = varIdent(...)к lme. lsmeansпоказывает стандартные ошибки, которые не совпадают со значительными различиями.

Пример моего кода ниже. Данные для воспроизведения вывода можно найти здесь.

      library(nlme)
library(multcomp)
library(emmeans)

model <- lme(variable ~ cultivar*year, 
             random = ~1|block, 
             weights = varIdent(form = ~1|cultivar),  
             method = "REML", 
             na.action = na.omit, 
             data = ag.data)

Leastsquare <- lsmeans(model,"cultivar")
cld(Leastsquare, Letters = letters)

       cultivar  lsmean     SE df lower.CL upper.CL .group
 Golden      6.92  3.841  1    -41.9     55.7  a    
 Campfield  10.33  4.330  1    -44.7     65.4  a    
 Tom        17.50  0.167  1     15.4     19.6  a    
 Harrison   25.67 12.649  1   -135.1    186.4  ab   
 Puget      30.58 20.502  1   -229.9    291.1  ab   
 HVC        37.08  5.331  1    -30.7    104.8   b   
 COL        38.08  0.433  1     32.6     43.6   b   
 Brown      62.67 20.207  1   -194.1    319.4  ab   

Как может быть, что сорт Brownсущественно не отличается от Golden? Это приемлемо? Кто-нибудь видел результаты, подобные этому?

Answer 1

Я превратил ваш пример в репрекс . Пожалуйста, найдите мои комментарии ниже.

      library(tidyverse)

ag.data <- tibble::tribble(
       ~year,   ~cultivar, ~block,   ~variable,
  "nineteen",       "HVC",     1L, 14.33333333,
  "nineteen",       "HVC",     2L, 23.33333333,
  "nineteen",     "Puget",     1L, 2.333333333,
  "nineteen",     "Puget",     2L, 3.333333333,
  "nineteen", "Campfield",     1L,          NA,
  "nineteen", "Campfield",     2L,           4,
  "nineteen",       "Tom",     1L,          10,
  "nineteen",       "Tom",     2L,          10,
  "nineteen",     "Brown",     1L,          NA,
  "nineteen",     "Brown",     2L, 56.66666667,
  "nineteen",       "COL",     1L,          NA,
  "nineteen",       "COL",     2L, 51.66666667,
  "nineteen",    "Golden",     1L,           5,
  "nineteen",    "Golden",     2L, 1.666666667,
  "nineteen",  "Harrison",     1L, 52.33333333,
  "nineteen",  "Harrison",     2L, 4.333333333,
    "twenty",       "HVC",     1L, 45.66666667,
    "twenty",       "HVC",     2L,          65,
    "twenty",     "Puget",     1L, 17.33333333,
    "twenty",     "Puget",     2L, 99.33333333,
    "twenty", "Campfield",     1L, 11.66666667,
    "twenty", "Campfield",     2L, 21.66666667,
    "twenty",       "Tom",     1L, 25.33333333,
    "twenty",       "Tom",     2L, 24.66666667,
    "twenty",     "Brown",     1L, 45.33333333,
    "twenty",     "Brown",     2L,          92,
    "twenty",       "COL",     1L,          24,
    "twenty",       "COL",     2L,          25,
    "twenty",    "Golden",     1L,           3,
    "twenty",    "Golden",     2L,          18,
    "twenty",  "Harrison",     1L,          31,
    "twenty",  "Harrison",     2L,          15
  )

library(nlme)
library(emmeans)
library(multcomp)
library(multcompView)

model <- lme(
  variable ~ cultivar * year,
  random = ~ 1 | block,
  weights = varIdent(form = ~ 1 | cultivar),
  method = "REML",
  na.action = na.omit,
  data = ag.data
)

anova(model)
#>               numDF denDF   F-value p-value
#> (Intercept)       1    12 16502.874  <.0001
#> cultivar          7    12   193.823  <.0001
#> year              1    12   952.713  <.0001
#> cultivar:year     7    12   296.145  <.0001

ag.data %>%
  filter(!is.na(variable)) %>%
  ggplot(aes(y = variable, x = year)) +
  facet_wrap(vars(cultivar)) +
  geom_point() +
  stat_summary(fun = mean,
               color = "red",
               geom = "line",
               aes(group = 1)) +
  theme_bw()

      emm <- emmeans(model, ~ cultivar) %>% 
  cld(Letters = letters) %>% 
  as_tibble() %>% 
  mutate(cultivar = fct_reorder(cultivar, emmean))
#> NOTE: Results may be misleading due to involvement in interactions

ggplot(emm, aes(
  y = emmean,
  ymin = lower.CL,
  ymax = upper.CL,
  x = cultivar,
  label = str_trim(.group)
)) +
  geom_point() +
  geom_errorbar(width = 0.1) +
  geom_text(
    position = position_nudge(x = 0.1),
    hjust = 0,
    color = "red"
  ) +
  theme_bw()

^{Создано 24 января 2022 г. пакетом reprex (v2.0.1)}

Отвечая на ваш вопрос, почему сорта с самым низким и самым высоким emmeans существенно не отличаются: я думаю, глядя на второй график, становится ясно, что это связано с совершенно разной точностью, с которой оценивались emmeans. Это, в свою очередь, частично связано с разнородными отклонениями ошибок, которые вы допускаете в модели, а также с отсутствующими/несбалансированными данными. Я бы сказал, что вы «не привыкли видеть, что экстремальные значения статистически не отличаются друг от друга, когда промежуточные значения различаются», потому что обычно со сбалансированными данными и / или без гетерогенной дисперсии ошибок вы бы этого не сделали. Попробуйте запустить код без аргумент (т.е. со стандартной однородной дисперсией ошибки) - вы не найдете «проблему», с которой вы столкнулись, с этими результатами.

Кроме того, обратите внимание, что у вас на самом деле, кажется, есть взаимодействия сорта-года - см. и первый сюжет. Таким образом, просмотр средств сорта может ввести в заблуждение, поскольку примечание ниже говорит функция. Вместо этого вы можете попробовать исследовать emmean-сравнения в год с помощью .