三路重复测量方差分析的 lme4::lmer 等效项是多少？

机器算法验证 r 方差分析 混合模式 重复测量 lme4-nlme
2022-03-29 01:13:28

我的问题基于此响应，该响应显示了哪个lme4::lmer模型对应于双向重复测量方差分析：

require(lme4)
set.seed(1234)
d <- data.frame(
    y = rnorm(96),
    subject = factor(rep(1:12, 4)),
    a = factor(rep(1:2, each=24)),
    b = factor(rep(rep(1:2, each=12))),
    c = factor(rep(rep(1:2, each=48))))

# standard two-way repeated measures ANOVA:
summary(aov(y~a*b+Error(subject/(a*b)), d[d$c == "1",]))

# corresponding lmer call:
anova(lmer(y ~ a*b+(1|subject) + (1|a:subject) + (1|b:subject), d[d$c == "1",]))

我现在的问题是如何将其扩展到三向方差分析的情况：

summary(aov(y~a*b*c+Error(subject/(a*b*c)), d))
## [...]
## Error: subject:a:b:c
##           Df Sum Sq Mean Sq F value Pr(>F)
## a:b:c      1  0.101  0.1014   0.115  0.741
## Residuals 11  9.705  0.8822

自然扩展及其版本与 ANOVA 结果不匹配：

anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject), d))
## [...]
## a:b:c  1 0.1014  0.1014  0.1500

anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject) + 
               (1|a:b:subject) + (1|a:c:subject) + (1|b:c:subject), d))
## [...]
## a:b:c  1 0.1014  0.1014  0.1539

请注意，之前已经提出了一个非常相似的问题。但是，它缺少示例数据（在此处提供）。

2个回答

你的问题的直接答案是你写的最后一个模型，

anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject) + 
           (1|a:b:subject) + (1|a:c:subject) + (1|b:c:subject), d))

我相信“原则上”是正确的，尽管这是一个奇怪的参数化，在实际实践中似乎并不总是很好。

至于为什么你从这个模型得到的输出与aov()输出不一致，我认为有两个原因。

  1. 您的简单模拟数据集是病态的，因为最佳拟合模型是暗示负方差分量的模型，混合模型lmer()（和大多数其他混合模型程序）不允许使用。
  2. 即使使用非病理数据集，如上所述，您设置模型的方式在实践中似乎并不总是有效，尽管我必须承认我并不真正理解为什么。在我看来，这通常也很奇怪，但那是另一回事了。

让我首先在您的初始双向 ANOVA 示例中演示我更喜欢的参数化。假设您的数据集d已加载。您的模型（请注意，我从虚拟代码更改为对比代码）是：

options(contrasts=c("contr.sum","contr.poly"))
mod1 <- lmer(y ~ a*b+(1|subject) + (1|a:subject) + (1|b:subject),
         data = d[d$c == "1",])
anova(mod1)
# Analysis of Variance Table
#     Df  Sum Sq Mean Sq F value
# a    1 2.20496 2.20496  3.9592
# b    1 0.13979 0.13979  0.2510
# a:b  1 1.23501 1.23501  2.2176

它在这里工作得很好，因为它与aov()输出相匹配。我喜欢的模型涉及两个变化：手动对因子进行对比编码，以便我们不使用 R 因子对象（我建议在 100% 的情况下这样做），并以不同的方式指定随机效应：

d <- within(d, {
  A <- 2*as.numeric(paste(a)) - 3
  B <- 2*as.numeric(paste(b)) - 3
  C <- 2*as.numeric(paste(c)) - 3
})
mod2 <- lmer(y ~ A*B + (1|subject)+(0+A|subject)+(0+B|subject),
             data = d[d$c == "1",])
anova(mod2)
# Analysis of Variance Table
# Df  Sum Sq Mean Sq F value
# A    1 2.20496 2.20496  3.9592
# B    1 0.13979 0.13979  0.2510
# A:B  1 1.23501 1.23501  2.2176

logLik(mod1)
# 'log Lik.' -63.53034 (df=8)
logLik(mod2)
# 'log Lik.' -63.53034 (df=8)

这两种方法在简单的 2-way 问题中是完全等价的。现在我们将转向一个三向问题。我之前提到过您提供的示例数据集是病态的。因此，在处理您的示例数据集之前，我想做的是首先从实际方差分量模型（即，非零方差分量内置到真实模型中）生成数据集。首先，我将展示我的首选参数化如何比您建议的更好。然后我将演示另一种估计方差分量的方法，它不会强加它们必须是非负的。然后我们将能够看到原始示例数据集的问题。

除了我们将有 50 个主题之外，新数据集的结构将相同：

set.seed(9852903)
d2 <- expand.grid(A=c(-1,1), B=c(-1,1), C=c(-1,1), sub=seq(50))
d2 <- merge(d2, data.frame(sub=seq(50), int=rnorm(50), Ab=rnorm(50),
  Bb=rnorm(50), Cb=rnorm(50), ABb=rnorm(50), ACb=rnorm(50), BCb=rnorm(50)))
d2 <- within(d2, {
  y <- int + (1+Ab)*A + (1+Bb)*B + (1+Cb)*C + (1+ABb)*A*B +
    (1+ACb)*A*C + (1+BCb)*B*C + A*B*C + rnorm(50*2^3)
  a <- factor(A)
  b <- factor(B)
  c <- factor(C)
})

我们要匹配的 F 比率是：

aovMod1 <- aov(y ~ a*b*c + Error(factor(sub)/(a*b*c)), data = d2)
tab <- lapply(summary(aovMod1), function(x) x[[1]][1,2:4])
do.call(rbind, tab)
#                          Sum Sq Mean Sq F value
# Error: factor(sub)       439.48    8.97        
# Error: factor(sub):a     429.64  429.64  32.975
# Error: factor(sub):b     329.48  329.48  27.653
# Error: factor(sub):c     165.44  165.44  17.924
# Error: factor(sub):a:b   491.33  491.33  49.694
# Error: factor(sub):a:c   305.46  305.46  41.703
# Error: factor(sub):b:c   466.09  466.09  40.655
# Error: factor(sub):a:b:c 392.76  392.76 448.101

这是我们的两个模型：

mod3 <- lmer(y ~ a*b*c + (1|sub)+(1|a:sub)+(1|b:sub)+(1|c:sub)+
  (1|a:b:sub)+(1|a:c:sub)+(1|b:c:sub), data = d2)
anova(mod3)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# a      1  32.73   32.73  34.278
# b      1  21.68   21.68  22.704
# c      1  12.53   12.53  13.128
# a:b    1  60.93   60.93  63.814
# a:c    1  50.38   50.38  52.762
# b:c    1  57.30   57.30  60.009
# a:b:c  1 392.76  392.76 411.365

mod4 <- lmer(y ~ A*B*C + (1|sub)+(0+A|sub)+(0+B|sub)+(0+C|sub)+
  (0+A:B|sub)+(0+A:C|sub)+(0+B:C|sub), data = d2)
anova(mod4)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# A      1  28.90   28.90  32.975
# B      1  24.24   24.24  27.653
# C      1  15.71   15.71  17.924
# A:B    1  43.56   43.56  49.694
# A:C    1  36.55   36.55  41.703
# B:C    1  35.63   35.63  40.655
# A:B:C  1 392.76  392.76 448.101

logLik(mod3)
# 'log Lik.' -984.4531 (df=16)
logLik(mod4)
# 'log Lik.' -973.4428 (df=16)

正如我们所看到的，只有第二种方法与 的输出相匹配aov()，尽管第一种方法至少在大致范围内。第二种方法也实现了更高的对数似然。我不确定为什么这两种方法会给出不同的结果，因为我再次认为它们“原则上”是等效的，但也许是出于一些数值/计算的原因。或者也许我弄错了，即使在原则上它们也不等同。

现在我将展示另一种基于传统 ANOVA 思想估计方差分量的方法。基本上，我们将为您的设计采用预期的均方方程，代入均方的观测值，并求解方差分量。为了获得预期的均方，我们将使用我几年前编写的 R 函数，称为EMS()，在此处记录。下面我假设该函数已经加载。

# prepare coefficient matrix
r <- 1 # number of replicates
s <- 50 # number of subjects
a <- 2 # number of levels of A
b <- 2 # number of levels of B
c <- 2 # number of levels of C
CT <- EMS(r ~ a*b*c*s, random="s")
expr <- strsplit(CT[CT != ""], split="")
expr <- unlist(lapply(expr, paste, collapse="*"))
expr <- sapply(expr, function(x) eval(parse(text=x)))
CT[CT != ""] <- expr
CT[CT == ""] <- 0
mode(CT) <- "numeric"
# residual variance and A*B*C*S variance are confounded in
# this design, so remove the A*B*C*S variance component
CT <- CT[-15,-2]
CT
#        VarianceComponent
# Effect  e b:c:s a:c:s a:b:s a:b:c c:s b:s a:s b:c a:c a:b s   c   b   a
#   a     1     0     0     0     0   0   0   4   0   0   0 0   0   0 200
#   b     1     0     0     0     0   0   4   0   0   0   0 0   0 200   0
#   c     1     0     0     0     0   4   0   0   0   0   0 0 200   0   0
#   s     1     0     0     0     0   0   0   0   0   0   0 8   0   0   0
#   a:b   1     0     0     2     0   0   0   0   0   0 100 0   0   0   0
#   a:c   1     0     2     0     0   0   0   0   0 100   0 0   0   0   0
#   b:c   1     2     0     0     0   0   0   0 100   0   0 0   0   0   0
#   a:s   1     0     0     0     0   0   0   4   0   0   0 0   0   0   0
#   b:s   1     0     0     0     0   0   4   0   0   0   0 0   0   0   0
#   c:s   1     0     0     0     0   4   0   0   0   0   0 0   0   0   0
#   a:b:c 1     0     0     0    50   0   0   0   0   0   0 0   0   0   0
#   a:b:s 1     0     0     2     0   0   0   0   0   0   0 0   0   0   0
#   a:c:s 1     0     2     0     0   0   0   0   0   0   0 0   0   0   0
#   b:c:s 1     2     0     0     0   0   0   0   0   0   0 0   0   0   0
#   e     1     0     0     0     0   0   0   0   0   0   0 0   0   0   0

# get mean squares
(MSmod <- summary(aov(y ~ a*b*c*factor(sub), data=d2)))
#                   Df Sum Sq Mean Sq
# a                  1  429.6   429.6
# b                  1  329.5   329.5
# c                  1  165.4   165.4
# factor(sub)       49  439.5     9.0
# a:b                1  491.3   491.3
# a:c                1  305.5   305.5
# b:c                1  466.1   466.1
# a:factor(sub)     49  638.4    13.0
# b:factor(sub)     49  583.8    11.9
# c:factor(sub)     49  452.2     9.2
# a:b:c              1  392.8   392.8
# a:b:factor(sub)   49  484.5     9.9
# a:c:factor(sub)   49  358.9     7.3
# b:c:factor(sub)   49  561.8    11.5
# a:b:c:factor(sub) 49   42.9     0.9
MS <- MSmod[[1]][,"Mean Sq"]

# solve
ans <- solve(CT, MS)
cbind(rev(ans[c(grep("e",names(ans)),grep("s",names(ans)))])/
        c(1,2,2,2,4,4,4,1))
# s     1.0115549
# a:s   1.5191114
# b:s   1.3797937
# c:s   1.0441351
# a:b:s 1.1263331
# a:c:s 0.8060402
# b:c:s 1.3235126
# e     0.8765093
summary(mod4)
# Random effects:
#  Groups   Name        Variance Std.Dev.
#  sub      (Intercept) 1.0116   1.0058  
#  sub.1    A           1.5191   1.2325  
#  sub.2    B           1.3798   1.1746  
#  sub.3    C           1.0441   1.0218  
#  sub.4    A:B         1.1263   1.0613  
#  sub.5    A:C         0.8060   0.8978  
#  sub.6    B:C         1.3235   1.1504  
#  Residual             0.8765   0.9362  
# Number of obs: 400, groups:  sub, 50

好的，现在我们将回到原来的例子。我们试图匹配的 F 比率是：

aovMod2 <- aov(y~a*b*c+Error(subject/(a*b*c)), data = d)
tab <- lapply(summary(aovMod2), function(x) x[[1]][1,2:4])
do.call(rbind, tab)
#                       Sum Sq Mean Sq F value
# Error: subject       13.4747  1.2250        
# Error: subject:a      1.4085  1.4085  1.2218
# Error: subject:b      3.1180  3.1180  5.5487
# Error: subject:c      6.3809  6.3809  5.2430
# Error: subject:a:b    1.5706  1.5706  2.6638
# Error: subject:a:c    1.0907  1.0907  1.5687
# Error: subject:b:c    1.4128  1.4128  2.3504
# Error: subject:a:b:c  0.1014  0.1014  0.1149

这是我们的两个模型：

mod5 <- lmer(y ~ a*b*c + (1|subject)+(1|a:subject)+(1|b:subject)+
  (1|c:subject)+(1|a:b:subject)+(1|a:c:subject)+(1|b:c:subject),
  data = d)
anova(mod5)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# a      1 0.8830  0.8830  1.3405
# b      1 3.1180  3.1180  4.7334
# c      1 3.8062  3.8062  5.7781
# a:b    1 1.5706  1.5706  2.3844
# a:c    1 0.9620  0.9620  1.4604
# b:c    1 1.4128  1.4128  2.1447
# a:b:c  1 0.1014  0.1014  0.1539

mod6 <- lmer(y ~ A*B*C + (1|subject)+(0+A|subject)+(0+B|subject)+
  (0+C|subject)+(0+A:B|subject)+(0+A:C|subject)+
  (0+B:C|subject), data = d)
anova(mod6)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# a      1 0.8830  0.8830  1.3405
# b      1 3.1180  3.1180  4.7334
# c      1 3.8062  3.8062  5.7781
# a:b    1 1.5706  1.5706  2.3844
# a:c    1 0.9620  0.9620  1.4604
# b:c    1 1.4128  1.4128  2.1447
# a:b:c  1 0.1014  0.1014  0.1539

logLik(mod5)
# 'log Lik.' -135.0351 (df=16)
logLik(mod6)
# 'log Lik.' -134.9191 (df=16)

在这种情况下，这两个模型产生的结果基本相同，尽管第二种方法的对数似然略高。两种方法都不匹配aov()。但是，让我们看看当我们像上面那样求解方差分量时得到什么，使用不将方差分量约束为非负的 ANOVA 程序（但它只能用于没有连续预测变量且没有缺失数据；经典的方差分析假设）。

# prepare coefficient matrix
r <- 1 # number of replicates
s <- 12 # number of subjects
a <- 2 # number of levels of A
b <- 2 # number of levels of B
c <- 2 # number of levels of C
CT <- EMS(r ~ a*b*c*s, random="s")
expr <- strsplit(CT[CT != ""], split="")
expr <- unlist(lapply(expr, paste, collapse="*"))
expr <- sapply(expr, function(x) eval(parse(text=x)))
CT[CT != ""] <- expr
CT[CT == ""] <- 0
mode(CT) <- "numeric"
# residual variance and A*B*C*S variance are confounded in
# this design, so remove the A*B*C*S variance component
CT <- CT[-15,-2]

# get mean squares
MSmod <- summary(aov(y ~ a*b*c*subject, data=d))
MS <- MSmod[[1]][,"Mean Sq"]

# solve
ans <- solve(CT, MS)
cbind(rev(ans[c(grep("e",names(ans)),grep("s",names(ans)))])/
        c(1,2,2,2,4,4,4,1))
# s      0.04284033
# a:s    0.03381648
# b:s   -0.04004005
# c:s    0.04184887
# a:b:s -0.03657940
# a:c:s -0.02337501
# b:c:s -0.03514457
# e      0.88224787
summary(mod6)
# Random effects:
#  Groups    Name        Variance  Std.Dev. 
#  subject   (Intercept) 7.078e-02 2.660e-01
#  subject.1 A           6.176e-02 2.485e-01
#  subject.2 B           0.000e+00 0.000e+00
#  subject.3 C           6.979e-02 2.642e-01
#  subject.4 A:B         1.549e-16 1.245e-08
#  subject.5 A:C         4.566e-03 6.757e-02
#  subject.6 B:C         0.000e+00 0.000e+00
#  Residual              6.587e-01 8.116e-01
# Number of obs: 96, groups:  subject, 12

现在我们可以看到原始示例的病态之处。最佳拟合模型是暗示几个随机方差分量为负的模型。但是lmer()（和大多数其他混合模型程序）将方差分量的估计限制为非负。这通常被认为是一个合理的约束，因为方差当然永远不会真正为负。然而，这种约束的结果是混合模型无法准确表示具有负类内相关性的数据集，即来自同一集群的观察结果较少的数据集（而不是更多）平均相似于从数据集中随机抽取的观察值，因此集群内方差大大超过集群间方差。这样的数据集是在现实世界中偶尔会遇到（或意外模拟！）的完全合理的数据集，但它们不能用方差分量模型来合理地描述，因为它们暗示负方差分量。但是，如果软件允许，它们可以由此类模型“不明智地”描述。aov()允许它。lmer()才不是。

a, b,是c固定效应还是随机效应？如果它们是固定的，你的语法将只是

summary(aov(y~a*b*c+Error(subject), d))
Error: subject
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 11  13.47   1.225               

Error: Within
          Df Sum Sq Mean Sq F value  Pr(>F)   
a          1   1.41   1.408   1.730 0.19235   
b          1   3.12   3.118   3.829 0.05399 . 
c          1   6.38   6.381   7.836 0.00647 **
a:b        1   1.57   1.571   1.929 0.16889   
a:c        1   1.09   1.091   1.339 0.25072   
b:c        1   1.41   1.413   1.735 0.19168   
a:b:c      1   0.10   0.101   0.124 0.72518   
Residuals 77  62.70   0.814  

library(lmerTest)
anova(lmer(y ~ a*b*c+(1|subject), data=d))
Analysis of Variance Table of type 3  with  Satterthwaite 
approximation for degrees of freedom
      Sum Sq Mean Sq NumDF  DenDF F.value   Pr(>F)   
a     1.4085  1.4085     1 76.991  1.7297 0.192349   
b     3.1180  3.1180     1 76.991  3.8291 0.053995 . 
c     6.3809  6.3809     1 76.991  7.8363 0.006469 **
a:b   1.5706  1.5706     1 76.991  1.9289 0.168888   
a:c   1.0907  1.0907     1 76.991  1.3394 0.250716   
b:c   1.4128  1.4128     1 76.991  1.7350 0.191680   
a:b:c 0.1014  0.1014     1 76.991  0.1245 0.725183
其它你可能感兴趣的问题
上一篇Scikit 二项式偏差损失函数 下一篇自举是否适合这种连续数据？