R语言mhurdle包IHS选项报错，寻求IHS双hurdle模型实现方案

我之前也碰到过mhurdle包选择dist="IHS"时弹出object 'lin' not found的报错，这大概率是包内部的小bug——处理IHS分布的代码分支里，开发者遗漏了对内部变量lin的定义，而其他分布（如box-cox、log-normal）的分支都正确做了这个定义。下面给你两个可行的解决思路：

方案1：临时修复mhurdle包的IHS分支代码 #

你可以用trace()函数临时修改mhurdle.fit的内部逻辑，补上缺失的lin定义：

# 先加载mhurdle包
library(mhurdle)

# 自动定位IHS分支并插入变量定义
trace(mhurdle.fit, 
      at = which(grepl('dist == "IHS"', body(mhurdle.fit)))[1] + 2, 
      expr = quote(lin <- TRUE),
      print = FALSE)

修改完成后再运行你的IHS双hurdle模型，应该就能避开这个报错了。注意这个修改是临时的，重启R会话后需要重新执行。

方案2：手动实现IHS双hurdle模型（更稳定可控） #

如果不想依赖包的临时修复，你可以用maxLik包手动编写似然函数来实现模型，灵活性更高也更稳定。以下是完整示例：

步骤1：定义IHS变换及逆变换 #

# IHS变换函数（针对正数值）
ihs <- function(y) {
  log(y + sqrt(y^2 + 1))
}

# IHS逆变换（用于后续预测）
ihs_inv <- function(z) {
  (exp(z) - exp(-z)) / 2
}

步骤2：编写双hurdle模型的对数似然函数 #

假设你的模型设定是：

• 第一阶段（参与决策）：用probit模型估计y>0的概率
• 第二阶段（正结果模型）：IHS变换后的ihs(y)服从正态分布

library(maxLik)

ihs_hurdle_loglik <- function(params, y, X_participate, X_outcome) {
  # 拆分参数：第一阶段系数、第二阶段系数、标准差
  n_beta1 <- ncol(X_participate)
  beta1 <- params[1:n_beta1]
  beta2 <- params[(n_beta1+1):(n_beta1 + ncol(X_outcome))]
  sigma <- params[length(params)]
  
  # 第一阶段：计算参与概率
  prob_participate <- pnorm(X_participate %*% beta1)
  
  # 第二阶段：计算正结果的条件对数似然
  y_pos <- y[y > 0]
  ihs_y_pos <- ihs(y_pos)
  mu_ihs <- X_outcome[y > 0, ] %*% beta2
  cond_loglik <- dnorm(ihs_y_pos, mean = mu_ihs, sd = sigma, log = TRUE)
  
  # 合并所有观测的对数似然
  total_loglik <- numeric(length(y))
  total_loglik[y == 0] <- log(1 - prob_participate[y == 0])
  total_loglik[y > 0] <- log(prob_participate[y > 0]) + cond_loglik
  
  sum(total_loglik)
}

步骤3：拟合模型 #

# 假设你的数据集为data，y是因变量，var1/var2是参与阶段的自变量，var3/var4是结果阶段的自变量
# 构造自变量矩阵（自动添加截距项）
X_participate <- model.matrix(~ var1 + var2, data = data)
X_outcome <- model.matrix(~ var3 + var4, data = data)

# 设置初始参数（用简单模型的结果初始化，提升收敛效率）
init_beta1 <- coef(glm(y > 0 ~ var1 + var2, data = data, family = binomial(link = "probit")))
init_beta2 <- coef(lm(ihs(y) ~ var3 + var4, data = data[data$y > 0, ]))
init_sigma <- sd(ihs(data$y[data$y > 0]))
init_params <- c(init_beta1, init_beta2, init_sigma)

# 拟合模型
model <- maxLik(ihs_hurdle_loglik, start = init_params, 
                y = data$y, X_participate = X_participate, X_outcome = X_outcome)

# 查看结果
summary(model)

这个手动实现的方法不仅能解决mhurdle的bug，还能让你根据需求调整模型（比如把第一阶段换成logit，或者加入两阶段的相关参数）。

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