重复测量分析：为什么要在主体因素中嵌套实验因素？

机器算法验证 r 混合模式 lme4-nlme 重复测量 随机效应模型
2022-03-23 07:45:10

考虑一个纯重复测量设计，（假设）3 个实验对象内因子 A、B 和 C，以及（为简单起见）每个因子 2 个水平。所以我们每个对象有 2*2*2 = 8 次测量。

现在我想用线性混合效应模型测试固定效应。我已经阅读了几个来源（例如 Andy Field 的书“使用 R 发现统计”，以及这个网站：http ://www.jason-french.com/tutorials/repeatedmeasures.html ），对于 lme，应该使用以下内容句法：

model <- lme(dv ~ A*B*C, random = ~1|id/A/B/C)

但是，我不明白您为什么要将主题内的因素“嵌套”在模型的随机部分中，而不仅仅是使用 (1|id)。这有什么意义，它有什么作用？

从概念上讲，我不明白为什么要将实验固定因素嵌套在随机主题因素中。到目前为止，我理解嵌套的方式，您只会用它来解释某些较低因子水平仅存在于某些较高因子水平内的事实 - 例如城市中学校内班级内的学生等。这如何适用于重复测量具有完全交叉的学科内因素的设计？

从数学上讲，我理解这一点的方式是，这样的模型将首先估计每个受试者的随机截距，捕捉受试者之间因变量平均值的随机差异。因此，在（假设）20 个受试者的情况下，我们得到 20 个不同的随机截距。然后，显然，该模型估计每个受试者组合和因素 A 水平的随机截距（导致 40 个随机截距），然后对于受试者、因素 A 和因素 B 的每个组合（80 个随机截距），一直到最具体的水平，我们得到与测量值一样多的估计随机截距（160）。这有什么意义，为什么我们不仅要估计每个主题（1|主题）的随机截距？还有，不会

最后，我的直觉告诉我，这些随机截距应该至少部分解释通过将实验因素的随机斜率输入模型中所捕获的相同信息。那是对的吗？

1个回答

这是一个非常有趣的问题。我一直在思考这个和相关的事情。

对我来说，理解这一点的关键是要认识到：分组因子的随机截距并不总是足以捕捉数据中超过残差变化的随机变化。正因为如此，我们有时会看到具有随机截距的模型，用于固定因子和分组变量之间的交互（甚至有时随机截距仅针对固定因子）。通常，我们建议一个因素可以是固定的或随机的（截距），但不能两者兼而有之 - 但是有一些重要的例外情况，您的示例就是其中之一。

理解这一点的另一个障碍是对于像我这样来自观察性社会和医学科学的多级建模/混合模型背景的人来说，我们经常陷入重复测量和嵌套与交叉随机效应的思考，而没有理解事物是实验分析略有不同。稍后再详细介绍。

从评论来看，我们都发现了同样的事情。在重复测量方差分析的背景下，如果您想获得相同的结果，lmer那么您适合：

 y ~ A + B + (1|id) + (1|id:A) + (1|id:B)

我在C不失一般性的情况下丢弃了因子。

有些人指定的原因1|id/A/B是他们正在使用nmle:lme而不是lme4:lmer. 我不确定为什么需要这样做，lme()但我相当确定要复制重复测量方差分析 - 每个组合id和因素都有变化 - 然后你将上面的模型拟合到lmer(). 请注意，这(1|id/A/B)似乎相似，但它是错误的，因为它也适合(1|id:A:B)与剩余方差无法区分的情况（如您的评论中所述）。

重要的是要注意（因此值得重复），我们只适合这种类型的模型，我们有理由相信每个组合id和因素都有变化。通常对于混合模型，我们不会这样做。我们需要了解实验设计。一种常见的实验类型是所谓的裂区设计，其中也使用了块。这种类型的实验设计采用不同“水平”的随机化 - 或者更确切地说是不同的因素组合，这就是为什么对此类实验的分析通常包括乍一看似乎很奇怪的随机截距项。然而，随机结构是实验设计的一个属性，如果不了解这一点，几乎不可能选择正确的结构。

所以，关于你的实际问题，实验有重复因子设计，我们可以使用我们的朋友，模拟，进一步调查。

我们将为模型模拟数据：

 y ~ A + B + (1|id)

和

 y ~ A + B + (1|id) + (1|id:A) + (1|id:B)

看看当我们使用这两个模型来分析两个数据集时会发生什么。

set.seed(15)
n <- 100 # number of subjects
K <- 4 # number of measurements per subject

# set up covariates
df <- data.frame(id = rep(seq_len(n), each = K),
                 A = rep(c(0,0,1,1), times = n),
                 B = rep(c(0,1), times = n * 2)
)

# 
df$y <- df$A + 2 * df$B + 3 * df$id + rnorm(n * K)

m0 <- lmer(y ~ A + B + (1|id) , data = df)
m1 <- lmer(y ~ A + B + (1|id) + (1|id:A) + (1|id:B), data = df)
summary(m0)

m0是这些数据的“正确”模型，因为我们只是创建y了固定效应id（我们用随机截距捕获）和单位方差。这有点滥用，但它很方便并且可以满足我们的要求：

 Groups   Name        Variance Std.Dev.
 id       (Intercept) 842.1869 29.0205 
 Residual               0.9946  0.9973 
Number of obs: 400, groups:  id, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept) 50.47508    2.90333   17.39
A            1.01277    0.09973   10.15
B            2.06675    0.09973   20.72

正如我们所看到的，我们恢复了残差的单位方差和固定效应的良好估计。然而：

> summary(m1)

Random effects:
 Groups   Name        Variance  Std.Dev.
 id:B     (Intercept) 0.000e+00  0.0000 
 id:A     (Intercept) 8.724e-03  0.0934 
 id       (Intercept) 8.422e+02 29.0204 
 Residual             9.888e-01  0.9944 
Number of obs: 400, groups:  id:B, 200; id:A, 200; id, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept) 50.47508    2.90334   17.39
A            1.01277    0.10031   10.10
B            2.06675    0.09944   20.78

这是一个奇异的拟合 - 对术语方差的估计为零id:B并且接近于零id:A- 我们碰巧知道在这里是正确的，因为我们没有模拟这些“水平”的任何方差。我们还发现：

> anova(m0, m1)

m0: y ~ A + B + (1 | id)
m1: y ~ A + B + (1 | id) + (1 | id:A) + (1 | id:B)
   npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
m0    5 1952.8 1972.7 -971.39   1942.8                     
m1    7 1956.8 1984.7 -971.39   1942.8 0.0052  2     0.9974

意味着我们非常喜欢（正确的）模型m0

因此，现在我们也在这些“级别”上模拟具有变化的数据。由于m1是我们要模拟的模型，我们将其设计矩阵用于随机效应：

# design matrix for the random effects
Z <- as.matrix(getME(m1, "Z"))

# design matrix for the fixed effects
X <-  model.matrix(~ A + B, data = df)

betas <- c(10, 2, 3) # fixed effects coefficients
D1 <- 1 # SD of random intercepts for id
D2 <- 2 # SD of random intercepts for id:A
D3 <- 3 # SD of random intercepts for id:B

# we simulate random effects
b <- c(rnorm(n*2, sd = D3), rnorm(n*2, sd = D2), rnorm(n, sd = D1))
# the order here is goverened by the order that lme4 creates the Z matrix

# linear predictor
lp <- X %*% betas + Z %*% b

# add residual variance of 1
df$y <- lp + rnorm(n * K)

m2 <- lmer(y ~ A + B + (1|id) + (1|id:A) + (1|id:B), data = df)
m3 <- lmer(y ~ A + B + (1|id), data = df)
summary(m2)

'm2' 是这里的核心模型，我们得到：

Random effects:
 Groups   Name        Variance Std.Dev.
 id:B     (Intercept) 6.9061   2.6279  
 id:A     (Intercept) 4.4766   2.1158  
 id       (Intercept) 2.9117   1.7064  
 Residual             0.8704   0.9329  
Number of obs: 400, groups:  id:B, 200; id:A, 200; id, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept)  10.3870     0.3866  26.867
A             1.8123     0.3134   5.782
B             3.0242     0.3832   7.892

截距的 SDid有点高，但除此之外，我们对随机和固定效应有很好的估计。另一方面：

> summary(m3)

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 6.712    2.591   
 Residual             8.433    2.904   
Number of obs: 400, groups:  id, 100

Fixed effects:
            Estimate Std. Error t value
(Intercept)  10.3870     0.3611  28.767
A             1.8123     0.2904   6.241
B             3.0242     0.2904  10.414

虽然固定效应点估计还可以，但它们的标准误差更大。当然，随机结构是完全错误的。最后：

> anova(m2, m3)

Models:
m3: y ~ A + B + (1 | id)
m2: y ~ A + B + (1 | id) + (1 | id:A) + (1 | id:B)
   npar    AIC    BIC   logLik deviance Chisq Df Pr(>Chisq)    
m3    5 2138.1 2158.1 -1064.07   2128.1                        
m2    7 1985.7 2013.7  -985.87   1971.7 156.4  2  < 2.2e-16 ***

表明我们非常喜欢m2

其它你可能感兴趣的问题
上一篇Kolmogorov-Smirnov 检验：p 值和 ks 检验统计量随着样本量的增加而减少 下一篇如何选择 GAM 参数的类型