在R中用optim实现带random intercepts的logistic regression：SE异常排查
2026-7-19
用optim实现带随机截距的Logistic回归：问题排查与修正
问题背景
尝试用R的optim函数估计带随机截距的Logistic回归模型，点估计结果尚可，但标准误（SE）表现异常（出现NaN），且运行时产生大量警告，需要排查代码错误并修正。
1. 模拟数据生成
#install.packages('arm') library(arm) # Gelman and Hill R package # 生成模拟数据 set.seed(1234) n <- 1000 # 观测数 groups <- 1:5 # 随机截距分组 my.obs <- sample(groups, n, replace = TRUE) # 为每个观测分配组别 my.cov <- rnorm(n, mean = 0.30, sd = 0.05) # 协变量数据 mu.a <- 0.30 # 随机截距均值 sigma.a <- 0.10 # 随机截距标准差 B0 <- rnorm(length(groups), mu.a, sigma.a) # 每组的随机截距真实值 B1 <- 0.45 # 协变量斜率真实值 my.inter <- B0[my.obs] # 每个观测对应的截距 lin.pred <- B0[my.obs] + B1 * my.cov # 线性预测值 logit.p <- exp(lin.pred) / (1 + exp(lin.pred)) # 观测概率 y <- rbinom(n, size=1, prob=logit.p) # 模拟二分类结果
2. 对比：固定截距Logistic模型
# 固定截距逻辑回归分析 mylogit <- glm(y ~ my.cov, family = "binomial") summary(mylogit)
该模型协变量斜率估计值（0.2947）与真实值（0.45）偏差较大：
Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.2670 0.3962 0.674 0.50 my.cov 0.2947 1.2988 0.227 0.82
3. 对比：arm包的随机截距模型
my.obs <- as.factor(my.obs) my.model <- glmer (y ~ my.cov + (1 | my.obs), family=binomial(link="logit")) summary(my.model) # 提取估计系数 fixed.effects <- fixef(my.model) fixed.effects random.effects <- ranef(my.model) random.effects my.coef <- coef(my.model) # 固定截距+随机截距的组合值 my.coef #$my.obs # (Intercept) my.cov #1 0.1683482 0.2352318 #2 0.3502571 0.2352318 #3 0.4548216 0.2352318 #4 0.2937610 0.2352318 #5 0.1377921 0.2352318 my.est.intercepts <- (as.numeric(fixed.effects[1]) + as.numeric(unlist(random.effects$my.obs))) my.est.intercepts #[1] 0.1683482 0.3502571 0.4548216 0.2937610 0.1377921 # 真实截距值 #[1] 0.4128684 0.2005167 0.3369223 0.2706445 0.3036248
4. 自定义optim实现的问题代码
logit.negLL = function(betas, y, my.obs, my.cov) { b1 = betas[1] # 随机截距均值 b2 = betas[2] # 随机截距标准差 b3 = betas[3] # 协变量斜率 b.obs1 <- rnorm(1, b1, b2) b.obs2 <- rnorm(1, b1, b2) b.obs3 <- rnorm(1, b1, b2) b.obs4 <- rnorm(1, b1, b2) b.obs5 <- rnorm(1, b1, b2) llh <- rep(0,length(y)) llh.total <- 0 for(i in 1:length(y)) { if(my.obs[i] == 1) my.prob = exp(b.obs1 + b3 * my.cov[i]) / (1 + exp(b.obs1 + b3 * my.cov[i])) if(my.obs[i] == 2) my.prob = exp(b.obs2 + b3 * my.cov[i]) / (1 + exp(b.obs2 + b3 * my.cov[i])) if(my.obs[i] == 3) my.prob = exp(b.obs3 + b3 * my.cov[i]) / (1 + exp(b.obs3 + b3 * my.cov[i])) if(my.obs[i] == 4) my.prob = exp(b.obs4 + b3 * my.cov[i]) / (1 + exp(b.obs4 + b3 * my.cov[i])) if(my.obs[i] == 5) my.prob = exp(b.obs5 + b3 * my.cov[i]) / (1 + exp(b.obs5 + b3 * my.cov[i])) llh[i] <- llh[i] + log(my.prob*y[i] + (1-my.prob)*(1 - y[i])) llh.total <- llh.total + llh[i] } -1 * llh.total } my.output <- optim(c(0.3, 0.1, .45), logit.negLL, y = y, my.obs = my.obs, my.cov = my.cov, method="BFGS", hessian=TRUE) # 出现50+警告（用warnings()查看前50条） # 点估计结果 R.est <- my.output$par R.est #[1] 0.21501634 0.02136591 0.47118554 true.values <- c(0.30, 0.10, 0.45) # 计算协方差矩阵与标准误 covmat = solve(my.output$hessian) # 标准误结果异常 R.SE <- diag(covmat)^0.5 R.SE #[1] NaN 0.001433485 0.002519182
错误分析与修正方案
核心错误点
- 似然函数逻辑错误：当前代码在似然函数中随机生成组截距，违背了随机效应模型的似然计算逻辑。随机截距是需要积分的潜在变量，而非每次计算似然时随机采样的变量，正确做法是计算边际似然（对随机截距的分布积分）。
- 标准差参数无约束：
b2作为随机截距的标准差必须大于0，但代码未施加约束，优化过程中可能出现b2<=0的情况，导致rnorm报错（这是大量警告的来源），同时破坏Hessian矩阵计算，最终出现NaN标准误。 - 代码效率低下：逐个观测判断组别的循环逻辑冗余，易出错且效率低。
修正后的代码
方法1：精确边际似然（适用于组数较少的场景）
logit.negLL_fixed = function(betas, y, my.obs, my.cov) { mu_a = betas[1] # 随机截距均值 sigma_a = exp(betas[2]) # 用exp确保标准差始终为正 beta1 = betas[3] # 协变量斜率 groups = unique(my.obs) total_ll = 0 for (g in groups) { # 提取当前组的所有观测 idx = my.obs == g y_g = y[idx] x_g = my.cov[idx] # 定义当前组的条件似然（对组截距a_g的函数） cond_ll = function(a_g) { lin_pred = a_g + beta1 * x_g prob = plogis(lin_pred) # plogis是logit逆函数，等价于exp(x)/(1+exp(x)) sum(dbinom(y_g, size=1, prob=prob, log=TRUE)) + dnorm(a_g, mu_a, sigma_a, log=TRUE) } # 数值积分计算当前组的边际似然 marginal_ll_g = integrate(Vectorize(cond_ll), lower=-Inf, upper=Inf)$value total_ll = total_ll + marginal_ll_g } -total_ll # 返回负对数似然 } # 初始值注意：betas[2]是sigma_a的对数，初始值为log(0.1) my.output_fixed <- optim(c(0.3, log(0.1), 0.45), logit.negLL_fixed, y = y, my.obs = my.obs, my.cov = my.cov, method="BFGS", hessian=TRUE) # 转换参数：sigma_a是exp(betas[2]) R.est_fixed = c(my.output_fixed$par[1], exp(my.output_fixed$par[2]), my.output_fixed$par[3]) R.est_fixed # 用delta方法转换标准差的标准误 covmat_fixed = solve(my.output_fixed$hessian) se_mu = sqrt(covmat_fixed[1,1]) se_sigma = exp(my.output_fixed$par[2]) * sqrt(covmat_fixed[2,2]) se_beta1 = sqrt(covmat_fixed[3,3]) R.SE_fixed = c(se_mu, se_sigma, se_beta1) R.SE_fixed
方法2：蒙特卡洛近似边际似然（适用于组数较多的场景）
logit.negLL_mc = function(betas, y, my.obs, my.cov, n_samples=20) { mu_a = betas[1] sigma_a = exp(betas[2]) beta1 = betas[3] groups = unique(my.obs) total_ll = 0 for (g in groups) { idx = my.obs == g y_g = y[idx] x_g = my.cov[idx] # 多次采样组截距，近似积分 mc_ll = replicate(n_samples, { a_g = rnorm(1, mu_a, sigma_a) lin_pred = a_g + beta1 * x_g prob = plogis(lin_pred) sum(dbinom(y_g, size=1, prob=prob, log=TRUE)) + dnorm(a_g, mu_a, sigma_a, log=TRUE) }) marginal_ll_g = log(mean(exp(mc_ll))) # 重要性采样的对数均值 total_ll = total_ll + marginal_ll_g } -total_ll } my.output_mc <- optim(c(0.3, log(0.1), 0.45), logit.negLL_mc, y = y, my.obs = my.obs, my.cov = my.cov, method="BFGS", hessian=TRUE, n_samples=50) # 转换参数 R.est_mc = c(my.output_mc$par[1], exp(my.output_mc$par[2]), my.output_mc$par[3]) R.est_mc # 用delta方法计算标准误 covmat_mc = solve(my.output_mc$hessian) se_mu_mc = sqrt(covmat_mc[1,1]) se_sigma_mc = exp(my.output_mc$par[2]) * sqrt(covmat_mc[2,2]) se_beta1_mc = sqrt(covmat_mc[3,3]) R.SE_mc = c(se_mu_mc, se_sigma_mc, se_beta1_mc) R.SE_mc
修正效果
- 无警告产生：通过
exp(betas[2])确保标准差始终为正。 - 似然计算符合逻辑：点估计更接近真实值。
- 标准误计算正常：Hessian矩阵有效，不会出现NaN。
内容的提问来源于stack exchange，提问作者Mark Miller
相关产品推荐
相关产品推荐
