重複測量分析：為什麼要在主體因素中嵌套實驗因素？

考慮一個純重複測量設計，（假設）3 個實驗對象內因子 A、B 和 C，以及（為簡單起見）每個因子 2 個水平。所以我們每個對像有 222 = 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|主題）的隨機截距？還有，不會

最後，我的直覺告訴我，這些隨機截距應該至少部分解釋通過將實驗因素的隨機斜率輸入模型中所捕獲的相同信息。那是對的嗎？

這是一個非常有趣的問題。我一直在思考這個和相關的事情。

對我來說，理解這一點的關鍵是要認識到：**分組因子的隨機截距並不總是足以捕捉數據中超過殘差變化的隨機變化。**正因為如此，我們有時會看到具有隨機截距的模型，用於固定因子和分組變量之間的交互（甚至有時隨機截距僅針對固定因子）。通常，我們建議一個因素可以是固定的或隨機的（截距），但不能兩者兼而有之 - 但是有一些重要的例外情況，您的示例就是其中之一。

理解這一點的另一個障礙是對於像我這樣來自觀察性社會和醫學科學的多級建模/混合模型背景的人來說，我們經常陷入重複測量和嵌套與交叉隨機效應的思考，而沒有理解事物是實驗分析略有不同。稍後再詳細介紹。

從評論來看，我們都發現了同樣的事情。在重複測量方差分析的背景下，如果您想獲得相同的結果，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

引用自：https://stats.stackexchange.com/questions/285021