复现sampleSelection包结果失败：代码报错求修正

错误原因分析 #

1. 初始参数逻辑错误：选择方程的初始参数误用了y对x1、z的OLS回归，正确做法应该用selection作为因变量做probit回归，否则初始参数严重偏离合理范围，导致似然函数无法求值。
2. 似然函数公式错误：核心的inside项计算错误，原代码分子误写为rho*sigma_u*(y_in-xb)，正确公式应为rho*(y_in-xb)/sigma_u（即rho*z_adj），该错误直接导致似然计算出现NaN或极端值。
3. 未处理边界值问题：当概率趋近于0或1时，log()会返回无穷大，需通过pmax添加极小值避免计算报错。
4. 官方包代码笔误：原官方包代码中data = df应为data = df_sel，否则会因数据集不存在报错。

修正后的完整代码 #

数据集与官方包验证代码 #

library(mvtnorm)
set.seed(123)
n = 1000
sigmau <- 1; sigmav <- 1; rho <- 0.5
Sigma <- matrix(c(sigmau^2,rho*sigmau*sigmav,rho*sigmau*sigmav,sigmav^2),2,2)

x1 <- rnorm(n,0,1); x2 <- rnorm(n,0,1); z <- rnorm(n,0,1)
u = rmvnorm(n,c(0,0),Sigma)[,1]; v = rmvnorm(n,c(0,0),Sigma)[,2]

selection_latent <- 1 + x1 + z + v
selection <- as.numeric(selection_latent >= 0)
y <- ifelse(selection == 1, 1 + x1 + x2 + u, NA)  # 未选择样本设为NA更符合逻辑

# 创建数据集
df_sel <- data.frame(id = 1:n,x1,x2,z,u = u,v,selection_latent,selection,y)

## 官方包估计
library(sampleSelection)
library(texreg)
sample_selection1 <- selection(selection = selection ~ 1 + x1 + z,
                               outcome = y ~ 1 + x1 + x2, 
                               data = df_sel)
screenreg(sample_selection1)

修正后的自定义函数代码 #

## 自定义样本选择模型函数
Selection_Reg <- function(y_in, x_in, z_in, selection, df) {
  # 构造带截距的自变量矩阵
  x <- cbind(1, as.matrix(x_in))
  z <- cbind(1, as.matrix(z_in))
  
  # 初始参数估计：beta用被选择样本OLS，theta用probit（符合选择模型逻辑）
  lm1 <- lm(y ~ 1 + x1 + x2, data = df, subset = selection == 1)
  glm2 <- glm(selection ~ 1 + x1 + z, data = df, family = binomial(link = "probit"))
  # 参数顺序：sigma_u, rho, beta(截距, x1, x2), theta(截距, x1, z)
  parameters <- c(1, 0.5, as.numeric(lm1$coef), as.numeric(glm2$coef))
  
  f_selection <- function(parameters, selection, y_in, x, z) {
    sigma_u <- parameters[1]
    rho <- parameters[2]
    sigma_v <- 1  # 选择方程误差方差固定为1（识别条件）
    beta <- parameters[3:(3 + ncol(x) - 1)]
    theta <- parameters[(3 + ncol(x)):(length(parameters))]
    
    # 计算线性预测
    xb <- x %*% beta
    zt <- z %*% theta
    
    # 计算被选择样本的似然部分
    z_adj <- (y_in - xb) / sigma_u
    # 修正后的inside计算公式
    inside <- (zt + rho * z_adj) / sqrt(1 - rho^2)
    
    # 计算对数似然，用pmax避免log(0)
    log_lik_selected <- -log(sigma_u) + log(dnorm(z_adj)) + log(pmax(pnorm(inside), 1e-10))
    log_lik_unselected <- log(pmax(1 - pnorm(zt), 1e-10))
    
    # 仅对对应样本计算似然
    log_lik <- ifelse(selection == 1, log_lik_selected, log_lik_unselected)
    return(-sum(log_lik, na.rm = TRUE))  # optim默认最小化，返回负对数似然和
  }
  
  # 使用BFGS优化方法（更适合光滑似然函数）
  result <- optim(par = parameters, fn = f_selection, 
                  selection = selection, y_in = y_in, x = x, z = z,
                  method = "BFGS", control = list(maxit = 1000))
  
  # 为参数命名，方便查看
  names(result$par) <- c("sigma_u", "rho", 
                         paste0("beta_", colnames(x)), 
                         paste0("theta_", colnames(z)))
  return(result)
}

# 运行自定义函数
fit_custom <- Selection_Reg(y_in = df_sel$y, 
                            x_in = df_sel[, c("x1", "x2")], 
                            z_in = df_sel[, c("x1", "z")],
                            selection = df_sel$selection, 
                            df = df_sel)

# 查看估计结果
print("自定义函数估计参数：")
fit_custom$par

# 对比官方包结果
print("官方包估计参数：")
coef(sample_selection1)

关键修正点说明 #

• 初始参数：选择方程用probit回归估计theta，beta用被选择样本的OLS回归，更符合模型逻辑。
• 似然公式：修正inside项的计算，匹配样本选择模型的联合似然公式。
• 边界处理：用pmax限制概率范围，避免log(0)导致的无穷大错误。
• 数据合理性：将未选择样本的y设为NA，并在似然计算中跳过，更符合样本选择模型的定义。

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