如何在R的mlogit包中实现部分退化嵌套Logit的参数归一化？

一、mlogit包实现（推荐） #

1. 数据格式转换 #

先将宽格式数据转换为mlogit要求的长格式，同时明确嵌套结构：

library(mlogit)

# 假设原始数据为宽格式，包含个体id、各选项属性（如cost、yield）、选择变量choice
data_long <- mlogit.data(
  your_data,
  choice = "choice",          # 选择变量
  shape = "wide",            # 原始数据为宽格式
  chid.var = "id",           # 个体标识
  alt.levels = c("crop1", "crop2", "forest"), # 所有选项
  nests = list(
    clear = c("crop1", "crop2"), # 清理土地的多选项嵌套
    keep = c("forest")           # 保留森林的单选项退化嵌套
  )
)

2. 模型拟合（含约束设置） #

根据Hunt(2000)的要求，需将退化嵌套（keep）的log-sum系数固定为1，仅估计clear嵌套的log-sum系数。通过un.nest.el = FALSE开启独立log-sum系数估计，再用constraints参数固定目标系数：

# 构建约束矩阵：固定keep嵌套的log-sum系数为1
temp_model <- mlogit(choice ~ cost + yield | 1, data = data_long, 
                     nests = list(clear = c("crop1", "crop2"), keep = c("forest")), 
                     un.nest.el = FALSE)
constraints <- matrix(0, nrow = 1, ncol = length(coefnames(temp_model)))
colnames(constraints) <- coefnames(temp_model)
constraints[1, "logsum(keep)"] <- 1

# 拟合最终模型
degen_nested_model <- mlogit(
  choice ~ cost + yield | 1,  # 效用函数：属性变量 + 截距项
  data = data_long,
  nests = list(clear = c("crop1", "crop2"), keep = c("forest")),
  un.nest.el = FALSE,         # 启用独立嵌套的log-sum系数
  constraints = list(fixed = constraints, values = c(1)) # 固定退化嵌套的log-sum为1
)

# 查看模型结果
summary(degen_nested_model)

• 结果中logsum(clear)为估计得到的参数，logsum(keep)会显示为固定值1。

二、nnet包替代方案（自定义似然） #

若需用nnet实现，需手动构建嵌套Logit的似然函数，适合熟悉模型原理的场景：

library(nnet)

# 自定义退化嵌套Logit的负对数似然函数（供optim最小化）
degen_nll <- function(params, data) {
  # 拆分参数：属性系数（beta_cost, beta_yield）、clear嵌套的log-sum系数（lambda_clear）
  beta_cost <- params[1]
  beta_yield <- params[2]
  lambda_clear <- params[3]
  
  # 计算各选项的效用
  u_crop1 <- beta_cost * data$cost_crop1 + beta_yield * data$yield_crop1
  u_crop2 <- beta_cost * data$cost_crop2 + beta_yield * data$yield_crop2
  u_forest <- beta_cost * data$cost_forest + beta_yield * data$yield_forest
  
  # 计算下层嵌套的inclusive value（退化嵌套log-sum固定为1）
  iv_clear <- lambda_clear * log(exp(u_crop1 / lambda_clear) + exp(u_crop2 / lambda_clear))
  iv_keep <- u_forest  # 等价于1 * log(exp(u_forest/1))
  
  # 上层选择概率
  prob_clear <- exp(iv_clear) / (exp(iv_clear) + exp(iv_keep))
  prob_keep <- 1 - prob_clear
  
  # 嵌套内的条件选择概率
  prob_crop1_given_clear <- exp(u_crop1 / lambda_clear) / (exp(u_crop1 / lambda_clear) + exp(u_crop2 / lambda_clear))
  prob_crop2_given_clear <- 1 - prob_crop1_given_clear
  prob_forest_given_keep <- 1
  
  # 计算个体实际选择的概率
  individual_prob <- ifelse(data$choice == "crop1", prob_clear * prob_crop1_given_clear,
                            ifelse(data$choice == "crop2", prob_clear * prob_crop2_given_clear,
                                   prob_keep * prob_forest_given_keep))
  
  # 返回负对数似然之和
  -sum(log(individual_prob))
}

# 初始参数值
start_params <- c(0, 0, 1)

# 优化拟合
fit <- optim(
  par = start_params,
  fn = degen_nll,
  data = your_data,
  method = "BFGS"
)

# 查看估计参数
cat("Estimated parameters:\n")
cat("beta_cost:", fit$par[1], "\n")
cat("beta_yield:", fit$par[2], "\n")
cat("lambda_clear:", fit$par[3], "\n")

关键注意事项 #

• 必须固定退化嵌套的log-sum系数为1，否则模型存在识别问题（单选项嵌套的log-sum无法从数据中估计）。
• mlogit包的实现更简洁且自带模型诊断（如拟合优度、参数显著性检验），优先推荐使用。

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