三向嵌套随机效应ANOVA分析结果疑问求助

以下是模拟数据集构建与分析代码：

set.seed(100)
ps=factor(paste("countries_",rep(letters[1:3],each=200),sep=""))
vs=factor(paste("cities_",rep(rep(letters[1:4],each=50),3),sep=""))
qs=factor(paste("neighborhoods_",rep(rep(rep(letters[1:5],each=10),4),3),sep=""))
x=vector(mode="numeric",length=600)
for (i in 1:600) x[i]=rnorm(1)*2+as.integer((i-1)/10)
dataset=data.frame(countries=ps,cities=vs,neighborhoods=qs,values=x)
attach(dataset)

三个因子均为随机效应，城市嵌套于国家，社区嵌套于城市，使用lme4构建模型：

library(lme4)
mod=lmer(x~(1| countries)+(1| countries: cities)+(1| countries: cities: neighborhoods))

用lmerTest的rand()函数检验随机效应：

library(lmerTest)
rand(mod)

得到结果：

ANOVA-like table for random-effects: Single term deletions

Model:
1 | countries) + (1 | countries:cities) + (1 | countries:cities:neighborhoods)
                                     npar  logLik    AIC     LRT Df Pr(>Chisq)    
<none>                                  5 -1370.4 2750.8  
(1 | countries)                         4 -1378.0 2764.0  15.151  1  9.924e-05 ***
(1 | countries:cities)                  4 -1422.4 2852.8 103.999  1  < 2.2e-16 ***
(1 | countries:cities:neighborhoods)    4 -1440.5 2888.9 140.085  1  < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

存在两个疑问：

• 预期国家效应最强、城市次之、社区最弱，但LRT结果显示趋势相反；
• 所有卡方检验自由度均为1，无法理解。

1. LRT统计量趋势与预期相反的原因 #

你的模拟数据生成逻辑和预期效应强度完全相反。观察x的生成代码：

for (i in 1:600) x[i]=rnorm(1)*2+as.integer((i-1)/10)

as.integer((i-1)/10)的取值规律：

• 每10个样本（同一社区）对应同一个整数值，社区内部无组间差异；
• 每50个样本（同一城市）包含5个社区，对应5个不同的整数值，城市内的差异全部来自社区层级；
• 每200个样本（同一国家）包含4个城市，对应20个不同的整数值，国家内的差异全部来自城市层级。

这意味着社区层级的组间差异最大，城市次之，国家最小。LRT统计量反映的是去掉该随机效应后模型拟合度下降的幅度：去掉社区效应时模型拟合下降最多（LRT=140），去掉城市效应次之（LRT=104），去掉国家效应下降最少（LRT=15），这完全匹配你的数据生成逻辑，和你预期的效应强度相反是因为模拟数据的设置问题，而非分析错误。

如果需要生成国家效应最强的模拟数据，应让国家层级的组间差异最大，示例代码如下：

# 重新生成x：国家间差异>城市间>社区间
set.seed(100)
country_eff = rep(rnorm(3, 0, 10), each=200)  # 国家效应方差最大
city_eff = rep(rep(rnorm(4, 0, 5), each=50), 3)  # 城市效应方差次之
neigh_eff = rep(rep(rep(rnorm(5, 0, 2), each=10),4),3)  # 社区效应方差最小
x = rnorm(600, 0, 2) + country_eff + city_eff + neigh_eff

2. 卡方检验自由度为1的原因 #

rand()函数执行的是单一项删除的似然比检验：每次从全模型中移除一个随机效应项，对比移除后模型与全模型的拟合差异。每个随机截距项对应一个方差参数，移除一个随机效应等价于固定该方差为0，本质是检验“该随机效应的方差是否大于0”——这是单个参数的检验，因此卡方分布的自由度为1（等于全模型与简化模型的参数数量之差：全模型有5个参数，即固定截距+3个随机方差+残差方差；移除一个随机效应后模型有4个参数，5-4=1）。

内容的提问来源于stack exchange，提问作者Eric Wajnberg