在R中使用optim直接估计Beta回归的alpha与beta参数遇阻求助

我在R中尝试用optim估计Beta回归的参数，当采用betareg包的mu-phi参数化方式时能成功估计，但无法直接用optim估计年度alpha和beta参数，怀疑是编码问题而非统计问题。

1. 模拟数据生成 #

先生成包含2年、每年1000个生存比例样本的数据集，将生存比例建模为年份的函数：

set.seed(1234)

# 创建数据
year     <-  1:2
years    <- length(year)  # 年份数
n.trials <-  1000         # 每年样本量

# 定义线性预测器参数
B0       <-  0.7
B1       <- -0.05

# 计算mu和phi
my.means <- exp(B0 + B1 * year) / (1 + exp(B0 + B1 * year))
my.means

my.phi <- 100

# 将mu和phi转换为年度alpha和beta
ab.fun <- function(mu, phi)
{
  a <- mu * phi
  b <- phi - mu * phi
  ab <- data.frame(a, b)
  return(ab)
}

alphabeta <- ab.fun(my.means, my.phi)
alphabeta
#          a        b
# 1 65.70105 34.29895
# 2 64.56563 35.43437

my.alpha <- alphabeta$a
my.beta  <- alphabeta$b

# 生成生存比例的随机样本
my.survival <- rep(NA, (years*n.trials))
my.year     <- rep(NA, (years*n.trials))

iter <- 1

for(i in 1:years) {
     for(j in 1:n.trials) {
          my.survival[iter] <- rbeta(1, my.alpha[i], my.beta[i])
          my.year[iter] <- i
          iter <- iter + 1
     }
}

head(my.survival)
head(my.year)

2. 用betareg包估计参数 #

使用betareg包得到的估计结果如下：

library(betareg)

gy <- betareg(my.survival ~ my.year)
summary(gy)
# Call:
# betareg(formula = my.survival ~ my.year)
# 
# Standardized weighted residuals 2:
#     Min      1Q  Median      3Q     Max 
# -3.4235 -0.6520 -0.0099  0.6561  3.5155 
# 
# Coefficients (mean model with logit link):
#              Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  0.721712   0.014738   48.97  < 2e-16 ***
# my.year     -0.067093   0.009293   -7.22 5.19e-13 ***
# 
# Phi coefficients (precision model with identity link):
#       Estimate Std. Error z value Pr(>|z|)    
# (phi)   100.69       3.17   31.77   <2e-16 ***
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
# 
# Type of estimator: ML (maximum likelihood)
# Log-likelihood:  3268 on 3 Df
# Pseudo R-squared: 0.02544
# Number of iterations: 8 (BFGS) + 2 (Fisher scoring) 
# Warning message:
# In deparse(x$call, width.cutoff = floor(getOption("width") * 0.85)) :
#   invalid 'cutoff' value for 'deparse', using default

3. 用optim估计mu-phi参数化的模型 #

用optim实现的mu-phi参数化模型，得到的估计结果和betareg几乎一致：

# 用optim估计参数
beta.reg = function(betas, my.survival, my.year){

     phi   = betas[1]
     B0    = betas[2]
     B1    = betas[3]

     n <- length(my.survival)
     llh <- rep(0, n)

     for(i in 1:n){

          mu <- exp(B0 + B1 * my.year[i]) / (1 + exp(B0 + B1 * my.year[i]))

          alpha <- mu * phi
          beta  <- phi - mu * phi

          llh[i] = llh[i] + dbeta(my.survival[i], alpha, beta, log = TRUE)
     }

     -sum(llh)
}

alpha.init <- 2
beta.init  <- 2
phi.init   <- 25
B0.init    <- 0
B1.init    <- 0

my.model <- optim(c(phi.init, B0.init, B1.init), 
                 beta.reg, my.survival = my.survival, my.year = my.year, method = "BFGS", hessian=TRUE)


my.model$par
#[1] 100.71561736   0.72166040  -0.06705041

my.model$value
#[1] -3268.381

covmat <- solve(my.model$hessian)

mySE <- diag(covmat)^0.5
mySE
#[1] 3.170301989 0.014736516 0.009291484

4. 尝试直接估计alpha和beta参数的问题 #

我尝试直接用optim估计年度alpha和beta参数，但代码运行失败，得到的结果不合理，还出现了奇异矩阵错误。以下是尝试的代码：

# 用optim估计参数
beta.reg = function(betas, my.survival, my.year){

     B0           = betas[1]
     B1           = betas[2]
     alpha1       = betas[3]
     alpha2       = betas[4]
     beta1        = betas[5]
     beta2        = betas[6]
     phi          = betas[7]

     n <- length(my.survival)
     llh <- rep(0, n)

     # 硬编码年份值为1:2简化mu计算
     mu1  <- exp(B0 + B1 * 1)  / (1 + exp(B0 + B1 * 1))
     mu2  <- exp(B0 + B1 * 2)  / (1 + exp(B0 + B1 * 2))

     alpha1   = mu1 * phi
     beta1    = phi - mu1 * phi
     alpha2   = mu2 * phi
     beta2    = phi - mu2 * phi

     iter <- 1

     for(i in 1:years){
          for(i in 1:n.trials){
               if(i ==  1) llh[iter] = llh[iter] + dbeta(my.survival[iter], alpha1, beta1, log = TRUE)
               if(i ==  2) llh[iter] = llh[iter] + dbeta(my.survival[iter], alpha2, beta2, log = TRUE)
               iter <- iter + 1
          }
     }
     -sum(llh)
}

B0.init    <- 0
B1.init    <- 0
alpha.init <- rep(2, years)
beta.init  <- rep(2, years)
phi.iter   <- 100

my.model <- optim(c(B0.init, B1.init, alpha.init, beta.init, phi.iter), 
                 beta.reg, my.survival = my.survival, my.year = my.year, method = "BFGS", hessian=TRUE)


my.model$par
#[1]   1.0013198  -0.2096154   2.0000000   2.0000000   2.0000000   2.0000000 338.7167214

my.model$value
#[1] -9.090913

covmat <- solve(my.model$hessian)
# Error in solve.default(my.model$hessian) : 
#  Lapack routine dgesv: system is exactly singular: U[3,3] = 0

mySE <- diag(covmat)^0.5
mySE

问题原因分析 #

这段代码有两个核心问题：

1. 参数冗余：同时传入了B0/B1、alpha1/alpha2/beta1/beta2和phi，但这些参数之间存在严格的函数关系（alpha=mu*phi、beta=phi-mu*phi，而mu由B0/B1决定），导致参数空间完全共线性，optim无法优化这种冗余模型，最终出现奇异矩阵错误。
2. 循环变量冲突：内层循环和外层循环都用了i作为变量，导致逻辑错误，无法正确匹配年份对应的样本。

5. 正确的直接估计alpha和beta的方案 #

如果要直接估计年度的alpha和beta参数，不需要保留B0/B1和phi，直接将alpha1、alpha2、beta1、beta2作为待估参数即可，因为模型本质是按年份分组的Beta分布，每个年份对应一组alpha和beta。

修改后的代码如下：

# 直接估计alpha和beta的对数似然函数
beta_reg_alpha_beta <- function(params, y, year) {
  # params顺序：alpha1, alpha2, beta1, beta2
  alpha1 <- params[1]
  alpha2 <- params[2]
  beta1 <- params[3]
  beta2 <- params[4]
  
  # 按年份分配参数
  alpha <- ifelse(year == 1, alpha1, alpha2)
  beta <- ifelse(year == 1, beta1, beta2)
  
  # 计算对数似然并取负（因为optim求最小值）
  -sum(dbeta(y, shape1 = alpha, shape2 = beta, log = TRUE))
}

# 初始化参数（用模拟的真实值作为初始值，加快收敛）
init_params <- c(alphabeta$a[1], alphabeta$a[2], alphabeta$b[1], alphabeta$b[2])

# 用optim优化
model_alpha_beta <- optim(init_params, 
                          beta_reg_alpha_beta, 
                          y = my.survival, 
                          year = my.year,
                          method = "BFGS",
                          hessian = TRUE)

# 查看估计结果
model_alpha_beta$par
# 应该接近模拟的真实值：65.70, 64.57, 34.30, 35.43

# 计算标准误
cov_mat <- solve(model_alpha_beta$hessian)
se <- sqrt(diag(cov_mat))
se

扩展：带线性趋势的alpha-beta模型 #

如果需要将alpha/beta和年份关联成线性模型（比如alpha = a0 + a1*year），可以重新构建参数化方式，确保参数无冗余：

# 带线性趋势的alpha-beta模型
beta_reg_trend <- function(params, y, year) {
  # params顺序：a0, a1, b0, b1（alpha = a0 + a1*year，beta = b0 + b1*year）
  alpha <- params[1] + params[2] * year
  beta <- params[3] + params[4] * year
  
  # 确保alpha和beta为正（Beta分布参数必须大于0）
  if(any(alpha <= 0) || any(beta <= 0)) return(Inf)
  
  -sum(dbeta(y, shape1 = alpha, shape2 = beta, log = TRUE))
}

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