R语言从零实现线性混合模型极大似然函数遇optim参数报错求助

错误原因 #

optim的核心要求是：目标函数的第一个参数必须是待优化的参数向量par。你定义的mylmax函数把y,x,group放在参数列表最前面，完全没有接收optim传入的par参数，因此触发"unused argument (par)"错误。

同时你的代码还有几个关键问题需要修正：

• 待优化参数（beta、sigma2、tau2）被硬编码为固定值，没有通过par传入，完全没起到优化作用
• 方差协方差矩阵v的计算错误（随机效应部分的结构不符合线性混合模型要求）
• 函数没有返回对数似然值（缺少return(l)）
• 直接优化方差值可能出现负数，建议用对数转换约束参数为正

修正后的完整代码 #

library(lme4) # 加载sleepstudy数据集
data(sleepstudy)

x <- sleepstudy$Days
group <- sleepstudy$id
y <- sleepstudy$Reaction

# 定义对数似然函数，第一个参数为待优化的par
mylmax <- function(par, y, x, group) {
  n <- length(y)
  # 从par中拆分参数：前两个是beta，后两个是对数转换的sigma2和tau2（保证方差为正）
  beta <- par[1:2]
  log_sigma2 <- par[3]
  log_tau2 <- par[4]
  sigma2 <- exp(log_sigma2)
  tau2 <- exp(log_tau2)
  
  # 构建固定效应设计矩阵X（含截距项）和随机效应设计矩阵Z（随机截距模型）
  X <- cbind(1, x)
  Z <- model.matrix(~factor(group) - 1) # 分组指示矩阵
  
  # 计算方差协方差矩阵V
  V <- sigma2 * diag(n) + tau2 * Z %*% t(Z)
  # 计算V的逆矩阵
  inv_V <- solve(V)
  
  # 计算对数似然（极大似然）
  reml <- FALSE
  l <- -0.5 * (log(det(V)) + t(y - X %*% beta) %*% inv_V %*% (y - X %*% beta))
  
  # 如果是REML，调整似然值
  if (reml) {
    l <- l - 0.5 * log(det(t(X) %*% inv_V %*% X))
  }
  
  # 返回似然值（optim默认最小化，通过fnscale=-1转为最大化）
  return(as.numeric(l))
}

# 初始参数：beta初始值设为0，log_sigma2和log_tau2设为log(1000)
init_par <- c(0, 0, log(1000), log(1000))

# 运行optim进行极大似然估计
opt_result <- optim(par = init_par, 
                    fn = mylmax, 
                    y = y, x = x, group = group, 
                    method = "L-BFGS-B", 
                    hessian = TRUE,
                    control = list(fnscale = -1))

# 查看优化结果
opt_result

关键调整说明 #

1. 参数顺序修正：把par作为mylmax的第一个参数，确保optim能正常传入待优化参数
2. 参数拆分与约束：从par中拆分出beta、对数转换的sigma2和tau2，用exp()还原为正数方差，避免优化时出现无效的负值
3. 协方差矩阵修正：正确构建随机效应的设计矩阵Z，计算符合随机截距线性混合模型的协方差矩阵V
4. 返回值修正：添加return(as.numeric(l))，确保函数返回有效的对数似然值
5. 初始参数设置：给待优化参数设置合理的初始值，保证优化过程收敛

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