logit尺度二项式预测转换回响应尺度的时机与方法

问题背景 #

处理巢成功二项式数据（成功=1/失败=0，判定标准为至少孵化1只幼鸟），采用带随机截距（bird ID）的glmer模型，设置family=binomial(link="logit")，核心需求是将logit尺度的预测值转换回响应尺度（成功概率）。

使用predict(scale="response")得到的结果明显偏低，因此手动模拟两种转换逻辑：

• 先转换再平均（对应示例中example()函数的fit1）：对每个随机截距样本，先计算logit逆转换（sigmoid函数）得到个体成功概率，再取所有样本的均值
• 先平均再转换（对应示例中example2()函数的fit2）：先对logit尺度的预测值（含随机截距）取平均，再做logit逆转换，结果与predict(scale="response")默认输出一致

两种方法结果差异显著，需要明确规则、正确方法及原理。

数学原理：Jensen不等式的影响 #

差异的核心源于logit逆函数（sigmoid函数）的非线性，结合Jensen不等式：

• sigmoid函数 ( \sigma(x) = \frac{e^x}{1+ex} ) 是凸函数（二阶导数恒正）
• 根据Jensen不等式，对于凸函数有：( E[\sigma(X)] > \sigma(E[X]) )
  其中 ( X = \text{固定效应} + u )，( u \sim N(0, \sigma_u^2) ) 是随机截距，( E[\cdot] ) 表示期望（即模拟中的样本均值）

这直接解释了两种方法的差异：

• fit1计算的是 ( E[\sigma(X)] )：对每个随机截距对应的个体概率取平均，是总体平均成功概率（边际均值）
• fit2计算的是 ( \sigma(E[X]) )：先对logit尺度的预测值取平均（此时随机截距的期望为0，结果等于固定效应的线性组合），再转换，是随机截距为0的"典型个体"的成功概率（条件均值，当随机截距取期望时）

正确方法选择 #

根据研究目标选择对应方法：

1. 目标：估计总体平均成功概率

   • 选择fit1的逻辑：先对每个随机截距样本计算sigmoid，再取平均
   • 在R中无需手动模拟，可直接用emmeans包计算边际均值：

     library(emmeans)
     # 假设模型名为model
     emmeans(model, ~ age, type = "response")
     

   • 或用predict函数指定忽略随机效应的边际预测：

     predict(model, newdata = your_newdata, type = "response", re.form = NA)
     

2. 目标：估计"典型个体"（随机截距为0）的成功概率

   • 选择fit2的逻辑，对应predict(model, type = "response", re.form = ~0)的结果，这是固定效应对应的基准个体概率

3. 为什么默认predict结果偏低？

   • 默认predict.glmer输出的是条件均值：如果是已有数据，会使用拟合得到的随机截距估计值；如果是新数据，默认re.form=~0（即随机截距取0），因此结果等于sigma(固定效应)，也就是"典型个体"的概率，而非总体平均。

示例代码验证 #

优化后的模拟代码清晰展示了两种逻辑的差异：

# 先转换再平均：模拟总体平均成功概率（边际均值）
example <- function(int, slope, slope2, slope3, var_random_intercept, maxage, minage, replicates){
  age <- seq(minage, maxage, 1)
  y <- array(data = NA, dim = c(replicates, length(age))) 
  for(i in 1:replicates) {
    u <- rnorm(1, 0, sqrt(var_random_intercept)) # 单次抽样随机截距，避免重复抽样
    linear_pred <- int + u + slope*age + slope2*(age^2) + slope3*(age^3)
    y[i,] <- exp(linear_pred)/(1+exp(linear_pred))
  }
  return(colMeans(y))
}

# 先平均再转换：模拟典型个体的成功概率（条件均值）
example2 <- function(int, slope, slope2, slope3, var_random_intercept, maxage, minage, replicates){
  age <- seq(minage, maxage, 1)
  y <- array(data = NA, dim = c(replicates, length(age))) 
  for(i in 1:replicates) {
    y[i,] <- int + rnorm(1,0, sqrt(var_random_intercept)) + slope*age + slope2*(age^2) + slope3*(age^3)
  }
  return(exp(colMeans(y))/(1+exp(colMeans(y))))
}

# 生成预测结果
example_pred <- data.frame(
  "age" = c(2:38), 
  "fit1" = example(int=-2.267e+00, slope=6.150e-02, slope2=3.296e-03, slope3=-1.487e-04, 
                   var_random_intercept=0.3068, maxage=38, minage=2, replicates=1000000),
  "fit2" = example2(int=-2.267e+00, slope=6.150e-02, slope2=3.296e-03, slope3=-1.487e-04, 
                    var_random_intercept=0.3068, maxage=38, minage=2, replicates=1000000)
)

# 可视化差异
library(ggplot2)
ggplot() +
  geom_line(data=example_pred, aes(x=age, y=fit1, color="总体平均概率")) +
  geom_line(data=example_pred, aes(x=age, y=fit2, color="典型个体概率")) +
  labs(x="年龄", y="巢成功概率", color="预测类型")

注：优化了原example()函数，避免同一随机截距重复抽样（原代码中分子分母各抽一次随机截距，属于错误，应单次抽样后复用）。

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