R线性混合效应模型：求解单个固定效应方差及占比问题

我是一名正在学习R语言的用户，已使用lme4包构建线性混合效应模型，尝试求解单个固定效应的方差。参考Nakagawa & Schielzeth 2013的方法，我通过固定效应设计矩阵乘以固定效应估计向量，再计算拟合值的方差来估计单个固定效应的方差，但结果有时出现负值，不确定是否正确；同时想了解如何将这些值转换为方差解释百分比，附上示例代码及结果，恳请解答。

示例数据 #

data
     Ho          SP        SL          DP         DL            species
1 0.432 -0.62240975 0.9615917 -0.20535833 -0.8038894 Accipiter gentilis
2 0.355 -0.85728592 1.0821313 -0.07917467 -0.8055025 Accipiter gentilis
3 0.370 -0.91221808 0.7616323 -0.17403906 -0.8041293 Accipiter gentilis
4 0.447 -1.15202246 0.9326085 -0.07640844 -0.8031998 Accipiter gentilis
5 0.452  0.12952583 0.7566257  0.35393076 -0.8058789 Ales Alces
6 0.406 -0.06462976 0.6679433  0.86746704 -0.7952550 Ales Alces

示例代码及结果 #

lme <- lmer(data$Ho ~ data$SP + data$DP + data$SL + data$DL + (1|data$spp)) 
vector <- fixef(lme)
vector <- vector[-1]
datanew <- data.frame(data$SP, data$DP, data$SL, data$DL)
matrix <- datanew * vector
variance <- data.frame(var(matrix)[1], var(matrix)[2], var(matrix)[3], var(matrix)[4])

运行结果：

var.matrix..1. var.matrix..2. var.matrix..3. var.matrix..4.
1    0.002225362   1.712515e-05   -0.001885247   -0.001308941

1. 负值出现的原因及错误根源 #

你当前的计算逻辑存在核心问题：直接将自变量数据框与固定效应系数逐元素相乘，得到的并非单个固定效应的边际贡献拟合值，这种操作忽略了固定效应间的相关性，再加上样本量过小（示例仅6个样本），就会出现不合理的负值。

Nakagawa & Schielzeth 2013的方法要求：单独提取每个固定效应对应的设计矩阵列，乘以其系数得到该效应的拟合值后，再计算方差——不是将所有变量与系数相乘后直接取列方差。

2. 正确计算单个固定效应方差的方法 #

修正后的代码如下，核心是逐个处理每个固定效应：

library(lme4)

# 修正模型公式（建议用规范写法，避免$data$的冗余引用）
lme <- lmer(Ho ~ SP + DP + SL + DL + (1|species), data = data)

# 提取固定效应系数
fixefs <- fixef(lme)

# 初始化结果容器
single_effect_var <- numeric(length(fixefs)-1)
names(single_effect_var) <- names(fixefs)[-1]

# 逐个计算单个固定效应的方差
for (var_name in names(single_effect_var)) {
  # 提取对应变量的列（连续变量直接取原始数据；分类变量需用model.matrix提取设计矩阵列）
  var_col <- data[[var_name]]
  # 计算该效应的拟合值：系数 × 变量值
  effect_fit <- fixefs[var_name] * var_col
  # 计算方差
  single_effect_var[var_name] <- var(effect_fit)
}

# 查看结果
single_effect_var

关键提示： #

• 如果固定效应是分类变量，不能直接用原始数据列，必须通过model.matrix(lme)提取对应的设计矩阵列后再计算拟合值。
• 小样本下方差估计仍可能不稳定，但不会出现负值；若仍有负值，需检查变量是否存在极端共线性。

3. 转换为方差解释百分比 #

要计算单个固定效应对总方差的解释比例，需先明确模型的总方差构成（固定效应方差+随机效应方差+残差方差），代码如下：

# 提取模型方差分量
var_comp <- VarCorr(lme)
# 随机效应方差（物种水平截距方差）
rand_var <- as.numeric(var_comp$species[1])
# 残差方差
resid_var <- attr(var_comp, "sc")^2
# 所有固定效应的总方差
total_fixed_var <- sum(single_effect_var)
# 模型总方差
total_var <- total_fixed_var + rand_var + resid_var

# 计算单个固定效应的方差解释百分比
var_explained_pct <- (single_effect_var / total_var) * 100
var_explained_pct

说明： #

• Nakagawa & Schielzeth的方法中，还区分了边际R²（固定效应解释的方差占总方差比例）和条件R²（固定+随机效应解释的方差占总方差比例），单个固定效应的解释比例是边际R²的细分。
• 如果自变量未标准化，建议先对变量做标准化处理再建模，避免尺度差异影响方差解释比例的可比性。

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