考虑随机效应不确定性的边际均值

Question

当我们对实验单元进行重复测量时，通常不能将这些单元视为“独立”的，需要以我们对标准误差进行有效估计的方式进行建模。

当我比较通过使用混合模型（将单位视为随机效应）计算处理的边际均值而获得的区间时，在另一种情况下，首先对单位进行平均，然后在平均响应，我得到完全相同的不确定区间。

我们如何将单位测量的不确定性纳入我们认为治疗效果的不确定性中？

为了真正传播所有的不确定性，我们不应该看看处理的样子，对一个单元的“所有可能的测量值”进行平均吗？

``` r

    library(dplyr) 
    #> 
    #> Attaching package: 'dplyr'
    #> The following objects are masked from 'package:stats':
    #> 
    #>     filter, lag
    #> The following objects are masked from 'package:base':
    #> 
    #>     intersect, setdiff, setequal, union
    library(emmeans) 
    library(lme4) 
    #> Loading required package: Matrix
    library(ggplot2)
    
    tmp <- structure(list(treatment = c("A", "A", "A", "A", "A", "A", "A", 
    "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B", "B", "B", 
    "B", "B", "B", "B"), response = c(151.27333548, 162.3933313, 
    159.2199999, 159.16666725, 210.82, 204.18666667, 196.97333333, 
    194.54666667, 154.18666667, 194.99333333, 193.48, 191.71333333, 
    124.1, 109.32666667, 105.32, 102.22, 110.83333333, 114.66666667, 
    110.54, 107.82, 105.62000069, 79.79999821, 77.58666557, 75.78666928
    ), experimental_unit = c("A-1", "A-1", "A-1", "A-1", "A-2", "A-2", 
    "A-2", "A-2", "A-3", "A-3", "A-3", "A-3", "B-1", "B-1", "B-1", 
    "B-1", "B-2", "B-2", "B-2", "B-2", "B-3", "B-3", "B-3", "B-3"
    )), row.names = c(NA, -24L), class = c("tbl_df", "tbl", "data.frame"
    ))
    
    
    ### Option 1 - Treat the experimental unit as a random effect since there are 
    ### 4 repeat observations for the same unit  
    
    lme4::lmer(response ~ treatment + (1 | experimental_unit), data = tmp) %>% 
      emmeans::emmeans(., ~ treatment) %>% 
      as.data.frame() 
    #>   treatment   emmean       SE df  lower.CL upper.CL
    #> 1         A 181.0794 10.83359  4 151.00058 211.1583
    #> 2         B 101.9683 10.83359  4  71.88947 132.0472
      #ggplot(.,aes(treatment, emmean)) + 
      #geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL))  
    
    
    
    
    
    ### Option 2 - instead of treating the unit as random effect, we average over the 
    ### 4 repeat observations, and run a simple linear model  
    
    tmp %>%
      group_by(experimental_unit) %>%
      summarise(mean_response = mean(response)) %>%
      mutate(treatment = c(rep("A", 3), rep("B", 3))) %>%
      lm(mean_response ~ treatment, data = .) %>%
      emmeans::emmeans(., ~ treatment) %>%
      as.data.frame() 
    #>   treatment   emmean       SE df  lower.CL upper.CL
    #> 1         A 181.0794 10.83359  4 151.00058 211.1583
    #> 2         B 101.9683 10.83359  4  71.88947 132.0472
      #ggplot(., aes(treatment, emmean)) +
      #geom_pointrange(aes(ymin = lower.CL, ymax = upper.CL))  
    
    
    
    
    ### Whether we include a random effect for the unit, or average over it and THEN model it, we find no difference in the 
    ### marginal means for the treatments 
    
    ### How do we incoporate the variation of the repeat measurments to the marginal means of the treatments? 
    ### Do we then ignore the variation in the 'subsamples' and simply average over them PRIOR to modeling? 
    
    
    <sup>Created on 2021-07-31 by the [reprex package](https://reprex.tidyverse.org) (v2.0.0)</sup> 

Answer 1

emmeans()确实考虑了随机效应的误差。这是我删除复杂的管道序列时得到的结果：

> mmod = lme4::lmer(response ~ treatment + (1 | experimental_unit), data = tmp)
> emmeans(mmod, "treatment")
 treatment emmean   SE df lower.CL upper.CL
 A            181 10.8  4    151.0      211
 B            102 10.8  4     71.9      132

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

如图所示。如果我拟合一个将实验单位作为固定效应的固定效应模型，我会得到：

> fmod = lm(response ~ treatment + experimental_unit, data = tmp)
> emmeans(fmod, "treatment")
NOTE: A nesting structure was detected in the fitted model:
    experimental_unit %in% treatment
 treatment emmean   SE df lower.CL upper.CL
 A            181 3.25 18    174.2      188
 B            102 3.25 18     95.1      109

Results are averaged over the levels of: experimental_unit 
Confidence level used: 0.95 

后一个结果的 SE 相当低，这是因为 experimental_unit 中的随机变化被建模为固定变化。

显然，您所做的管道说明了随机效应的变化，并包括 EMM 中的变化。我认为这是因为您为每个实验单元分别做了一些事情，并以某种方式将这些结果结合起来。我对 7 步长的管道序列不太满意，我不明白为什么这只会导致一组方法。

我建议最后不要使用as.data.frame()。这消除了有助于理解你所拥有的东西的注释。如果您这样做是为了获得更高的数字精度，我会声称这些数字是您不需要的，它只是夸大了您有权要求的精度。

一些后续cmets的说明

随后，我确信我们在 OP 第二部分的管道操作中看到的确实包括计算每个 EU 的平均值，然后对其进行分析。

让我们在正式模型的上下文中来看看。我们有（抱歉，MathJax 在 stackoverflow 上不起作用，但我还是会把标记留在那里）

$$ Y_{ijk} = \mu + \tau_i + U_{ij} + E_{ijk} $$

其中 $Y_{ijk}$ 是第 i 个处理中的第 k 个响应测量值和第 i 个处理中的第 j 个 EU，rhs 项分别代表总体平均值、（固定）处理效果、（随机）EU 效果，以及（随机）误差效应。我们假设随机效应都是相互独立的。通过平衡设计，EMM 只是边缘手段：

$$ \bar Y_{i..} = \mu + \tau_i + \bar U_{i.} + \bar E_{i..} $$

在哪里有一个“。”下标意味着我们对该下标进行平均。如果每个处理有 n 个 EU，每个 EU 有 m 个测量值，我们就得到了

$$ Var(\bar Y_{i..} = \sigma^2_U / n + \sigma^2_E / mn $$

现在，如果我们提前汇总欧盟的数据，我们将从

$$ \bar Y_{ij.} = \mu + U_{ij} + \bar E_{ij.} $$

但是，如果我们然后通过对 j 求平均来计算边际均值，我们得到的结果与我们之前用 $\bar Y_{i..}$ 所做的完全相同，并且方差与已经显示的完全一样。这就是为什么我们是否先聚合并不重要。