基于lmer输出推导ISA变化下窝卵数的95%置信区间（单位/百分比）

我正在分析不透水地表面积占比（ISA）与繁殖窝卵数（CS）的关联，使用的混合效应模型结构如下：

fit=lmer(CS~ISA+Mass_Fem*LayDate+Year+(1|Site),data=DataGT_noNA) 

变量说明：

• CS：窝卵数（响应变量）
• ISA：不透水地表占比（%，核心解释变量）
• Mass_Fem：雌性体重（克）
• LayDate：产卵日期序数
• Year：5水平分类因子
• Site：8水平随机效应（分组变量）

我已经通过effects包提取了特定ISA值对应的窝卵数估计值：

effects_ISA <- effects::effect(term= "ISA", mod= fit, xlevels=list(ISA=c(0,10,20,30,40,50,60,70)))
x_ISA <- as.data.frame(effects_ISA)

现在需要完成两个计算：

1. 给定ISA百分比变化时，窝卵数的绝对减少量及95%置信区间（例如ISA增加30%对应窝卵数减少0.54，需补充95%CI）
2. 对应窝卵数的百分比减少量及95%置信区间

方法1：基于已提取的效应值计算 #

既然已经有了不同ISA对应的预测值（包含置信区间），直接对目标区间的数值做差即可：

计算绝对减少量及95%CI #

以ISA从0%增加到30%为例：

# 提取ISA=0和ISA=30的行
isa0 <- x_ISA[x_ISA$ISA == 0, ]
isa30 <- x_ISA[x_ISA$ISA == 30, ]

# 绝对减少量：用ISA=0时的预测值减去ISA=30时的预测值
abs_decrease <- isa0$fit - isa30$fit
# 95%置信区间：ISA=0的置信上下限分别减去ISA=30的对应值
abs_ci_lower <- isa0$lower - isa30$upper
abs_ci_upper <- isa0$upper - isa30$lower

# 输出结果
cat(sprintf("ISA增加30%时，窝卵数绝对减少量：%.2f，95%%CI：[%.2f, %.2f]\n", abs_decrease, abs_ci_lower, abs_ci_upper))

计算百分比减少量及95%CI #

百分比减少量基于绝对减少量除以ISA=0时的窝卵数预测值：

# 百分比减少量
pct_decrease <- (abs_decrease / isa0$fit) * 100
# 95%置信区间：用绝对减少量的CI分别除以ISA=0的fit值
pct_ci_lower <- (abs_ci_lower / isa0$fit) * 100
pct_ci_upper <- (abs_ci_upper / isa0$fit) * 100

# 输出结果
cat(sprintf("ISA增加30%时，窝卵数百分比减少量：%.1f%%，95%%CI：[%.1f%%, %.1f%%]\n", pct_decrease, pct_ci_lower, pct_ci_upper))

方法2：直接从模型系数计算（更高效） #

ISA是连续变量，模型中ISA的系数就是每增加1%ISA时窝卵数的变化量，直接用系数乘以变化幅度即可：

提取模型系数及计算绝对变化 #

# 提取固定效应系数表
coef_table <- summary(fit)$coefficients
isa_coef <- coef_table["ISA", "Estimate"]
isa_se <- coef_table["ISA", "Std. Error"]

# 计算ISA增加30%的绝对变化：系数×30
abs_change_30 <- isa_coef * 30
# 95%CI：(系数±1.96×标准误)×30
abs_ci_30_lower <- (isa_coef - 1.96*isa_se) * 30
abs_ci_30_upper <- (isa_coef + 1.96*isa_se) * 30

# 因为是减少量，取绝对值（注意如果系数为负，CI上下限会反转）
abs_decrease_30 <- abs(abs_change_30)
abs_ci_30_lower_abs <- abs(abs_ci_30_upper)
abs_ci_30_upper_abs <- abs(abs_ci_30_lower)

计算百分比减少量及95%CI #

先获取ISA=0时的窝卵数预测值（控制其他协变量在均值/参考水平）：

# 构建参考水平数据：其他协变量取均值或第一水平
ref_data <- data.frame(
  ISA = 0,
  Mass_Fem = mean(DataGT_noNA$Mass_Fem, na.rm=TRUE),
  LayDate = mean(DataGT_noNA$LayDate, na.rm=TRUE),
  Year = levels(DataGT_noNA$Year)[1],
  Site = levels(DataGT_noNA$Site)[1]
)

# 预测ISA=0时的窝卵数（忽略随机效应，只看固定效应）
cs_isa0 <- predict(fit, newdata=ref_data, re.form=NA)

# 百分比减少量
pct_decrease_30 <- (abs_decrease_30 / cs_isa0) * 100
# 95%CI
pct_ci_30_lower <- (abs_ci_30_lower_abs / cs_isa0) * 100
pct_ci_30_upper <- (abs_ci_30_upper_abs / cs_isa0) * 100

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