请求协助将剂量反应实验似然函数转R代码以求解参数r

核心思路 #

剂量反应实验中，每个剂量组的感染数Y_i服从二项分布Binomial(N_i, p_i)，其中p_i是该剂量下的感染概率。假设采用常用的指数剂量反应模型：p_i = 1 - exp(-r*D_i)，则似然函数为所有组的二项概率乘积，我们通过最大化似然（或最小化负对数似然）来估计参数r。

数据集加载 #

先加载你提供的数据集：

Di <- c(4.2,8,10.6,11.9,12.2,12.6,13,13.2,28.3,28.7,68.2,69,71.1,83,91.7,95.5,106.5,136.5,171.2,186.9,309.3,1557.1)
Ni <- c(1,1,1,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
Yi <- c(0,0,1,1,2,1,0,1,0,1,1,0,1,1,1,1,1,1,1,1,1,1)
df <- data.frame(Di, Ni, Yi)

MLE求解代码 #

步骤1：定义负对数似然函数 #

neg_log_likelihood <- function(r, data) {
  D <- data$Di
  N <- data$Ni
  Y <- data$Yi
  
  # 按模型计算感染概率
  p <- 1 - exp(-r * D)
  
  # 边界处理避免log(0)报错
  p <- pmax(pmin(p, 0.9999), 0.0001)
  # 计算对数似然并取负
  log_lik <- sum(dbinom(Y, size = N, prob = p, log = TRUE))
  
  return(-log_lik)
}

步骤2：调用optim()求解最优参数 #

# 设置参数初始值（r必须为正，设为0.01）
initial_r <- 0.01

# 用L-BFGS-B方法限制r>0，确保参数符合实际意义
result <- optim(par = initial_r, 
                fn = neg_log_likelihood, 
                data = df,
                method = "L-BFGS-B",
                lower = 0)

# 输出求解结果
print(result)

步骤3：提取估计结果 #

estimated_r <- result$par
cat("估计的参数r值：", estimated_r, "\n")

min_neg_log_lik <- result$value
cat("最小负对数似然值：", min_neg_log_lik, "\n")

自定义模型说明 #

如果你的似然方程对应的感染概率模型不是指数形式，只需修改neg_log_likelihood()函数中p的计算表达式即可，替换为你给定的方程对应的感染概率公式。

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