R中用optim()求解二分类Logit模型MLE遇非数值二元操作符错误的求助

• 直接触发报错的诱因：你手动将因子类型的sector经as.matrix()转换后，矩阵存储的是分类因子的类型标记，并非可参与数值运算的0/1哑变量，当数值型的param[3]和因子类型的x2执行乘法时，R判定这是不兼容的非数值运算，直接抛出错误。
• 未触发但会导致结果错误的隐藏问题：
  • 即使强制把x2转成数值类型，手动编码的分类变量和R GLM内部的参照组哑变量编码规则不匹配，最终算出的参数会和GLM的参照结果存在系统性偏差
  • 你提取的因变量disease是factor类型，直接转矩阵后存储的是因子水平编码（1、2）而非似然计算需要的0/1取值，会导致负对数似然的计算完全偏离正确值

核心是让自定义似然函数用到的变量编码、取值和GLM内部的计算逻辑完全对齐，避免手动拆转换变量带来的类型错误和编码偏差：

• 直接用model.matrix()从已拟合的oreduced GLM对象中提取设计矩阵，这个矩阵已经自动包含了截距项、age的连续值列、sector对应的0/1参照组哑变量列，和GLM内部运算用的矩阵完全一致，不需要手动逐列提取变量
• 把因变量转成标准的0/1数值向量，匹配二项分布似然的计算要求
• 用矩阵乘法构造线性预测项，搭配R内置的plogis()函数计算概率，比手写逻辑斯蒂公式的数值稳定性更高

d <- read.delim("http://dnett.github.io/S510/Disease.txt")
d$disease = factor(d$disease)
d$ses = factor(d$ses)
d$sector = factor(d$sector)

# GLM参照拟合
oreduced <- glm(disease ~ age + sector, family = binomial(link = logit), data = d)
cat("GLM参照估计值：\n")
print(coef(oreduced))

# 提取对齐GLM逻辑的设计矩阵
X <- model.matrix(oreduced)
# 提取0/1编码的因变量
y <- as.numeric(d$disease) - 1

# 自定义负对数似然函数
nlldbin <- function(param) {
  eta <- X %*% param
  p <- plogis(eta)
  -sum(y * log(p) + (1 - y) * log(1 - p), na.rm = TRUE)
}

# 优化求解MLE
MLE_estimates <- optim(
  par = c(0.1, 0.1, 0.1),
  fn = nlldbin,
  hessian = TRUE
)
cat("\noptim求解的MLE估计值：\n")
print(MLE_estimates$par)

• 修复后不会再出现类型不匹配的报错
• 两组估计值的差异会落在optim的默认收敛精度（1e-4）范围内，完全可以作为对照验证
• 搭配hessian矩阵计算的标准误，也可以和GLM输出的系数标准误做交叉验证

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