如何获取线性混合效应模型中随机效应deployid的p值？

获取随机效应p值的操作方法 #

首先需要明确：lme()默认不输出随机效应的p值，是因为随机效应的方差参数检验属于边界检验（方差取值不能小于0），常规检验统计量的分布假设不成立，需要通过以下方法计算：

方法1：似然比检验（最常用） #

核心逻辑是拟合包含该随机效应的模型，和去掉该随机效应的嵌套模型，比较两个模型的拟合优度差异是否显著：

1. 首先将原模型和无随机效应的对照模型都用最大似然法（ML）拟合（REML法仅适用于固定效应相同的模型比较，检验随机效应差异需要用ML）：

# 含deployid随机效应的ML拟合模型
model_ml <- lme(divedepth ~ oarea, random = ~1|deployid, data = GDataTimes, method = "ML")
# 无随机效应的对照模型，用nlme包的gls拟合，保证模型框架一致
model_no_random <- gls(divedepth ~ oarea, data = GDataTimes, method = "ML")

2. 用anova函数比较两个嵌套模型，输出的p值就是随机效应deployid对应的显著性检验结果：

anova(model_ml, model_no_random)

方法2：参数自助法（Bootstrap，结果更稳健） #

如果对似然比检验的近似分布准确度有顾虑（尤其你的数据分组只有9组，属于分组数较少的情况），可以用参数自助法估计p值，核心逻辑是通过多次重采样模拟原假设（随机效应方差为0）下的统计量分布，计算实际观测统计量对应的分位数作为p值。

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