基于逐步选择的迁移占比建模：替代模型与实施步骤咨询

一、替代建模方法（支持AIC与stepAIC()） #

准二项式模型基于拟似然框架，没有严格的AIC值，无法直接用stepAIC()做变量筛选。推荐以下带完整似然的模型，完美适配你的需求：

1. Beta回归（优先推荐） #

专门适配0-1区间的比例数据（你的迁移占比正好符合），有完整似然函数，能计算AIC，直接支持向前选择。

• R代码示例：

  # 加载依赖包
  library(betareg)
  library(MASS)
  
  # 预处理：Beta回归要求响应值严格在(0,1)之间，处理0/1边界值
  df$fraction <- ifelse(df$fraction == 0, 1e-6, 
                        ifelse(df$fraction == 1, 1 - 1e-6, df$fraction))
  
  # 拟合包含所有预测变量的全模型
  full_model <- betareg(fraction ~ family_count + distance + purchasing_power + other_vars, 
                        data = df)
  
  # 从空模型开始执行向前变量选择
  null_model <- betareg(fraction ~ 1, data = df)
  step_model <- stepAIC(null_model, 
                        scope = list(lower = null_model, upper = full_model),
                        direction = "forward")
  

2. 标准二项式回归（若有原始计数数据） #

如果你能拿到每个邮编对的总迁移人数和该方向迁移人数（而非直接的占比），可以用标准二项式模型，天然支持AIC筛选：

• R代码示例：

  library(MASS)
  
  # 构造二项式响应：绑定迁移人数与非迁移人数
  df$non_migrate <- df$total_migrate - df$migrate_count
  
  # 拟合全模型
  full_model <- glm(cbind(migrate_count, non_migrate) ~ family_count + distance + purchasing_power + other_vars,
                    family = binomial, data = df)
  
  # 执行向前变量选择
  step_model <- stepAIC(full_model, direction = "forward")
  

3. 广义加性模型（GAM） #

如果预测变量和迁移占比存在非线性关系（比如距离的影响可能先降后升），GAM既能捕捉非线性效应，也支持AIC筛选：

• R代码示例：

  library(mgcv)
  library(MASS)
  
  # 拟合包含非线性项的全GAM模型
  full_model <- gam(fraction ~ s(family_count) + s(distance) + s(purchasing_power) + other_vars,
                    family = betar(link = "logit"), data = df)
  
  # 执行向前变量选择
  step_model <- stepAIC(full_model, direction = "forward")
  

二、若坚持使用准二项式模型的变量筛选方案 #

如果必须保留准二项式模型，可跳过stepAIC()，用以下替代方法：

1. 逐步F检验 #

每次添加一个变量，用方差分析比较模型显著性，手动完成向前选择：

# 初始空模型
model <- glm(fraction ~ 1, family = quasibinomial, data = df)

# 候选预测变量列表
vars <- c("family_count", "distance", "purchasing_power", "other_vars")

selected_vars <- c()
while(length(vars) > 0){
  # 测试每个候选变量的显著性
  p_values <- sapply(vars, function(var){
    temp_model <- update(model, formula = paste(". ~ . +", var))
    anova(model, temp_model, test = "F")$Pr[2]
  })
  # 选择p值最小的变量
  best_var <- names(which.min(p_values))
  if(min(p_values) < 0.05){ # 设定显著性阈值
    model <- update(model, formula = paste(". ~ . +", best_var))
    selected_vars <- c(selected_vars, best_var)
    vars <- setdiff(vars, best_var)
  } else {
    break # 无显著变量可添加，停止循环
  }
}

2. Lasso正则化（高维数据首选） #

用glmnet包的Lasso回归自动压缩不重要变量的系数到0，实现变量筛选：

library(glmnet)

# 构造预测矩阵和响应变量
X <- model.matrix(fraction ~ family_count + distance + purchasing_power + other_vars - 1, data = df)
y <- df$fraction

# 拟合Lasso模型
lasso_model <- glmnet(X, y, family = "quasibinomial", alpha = 1)

# 交叉验证选择最优正则化参数
cv_lasso <- cv.glmnet(X, y, family = "quasibinomial", alpha = 1)
best_lambda <- cv_lasso$lambda.min

# 提取筛选后的变量（系数非零的变量）
coef(lasso_model, s = best_lambda)

三、关键注意事项 #

• 空间自相关性：邮编属于空间数据，迁移行为可能存在空间依赖，可加入空间滞后项（比如用spdep包）或空间回归模型优化效果。
• 变量共线性：筛选前用car::vif()检查方差膨胀因子，排除VIF>5的变量，避免筛选结果不稳定。
• 邮编对对称性：从A到B和B到A的迁移可能存在关联，可加入交互项或控制反向迁移变量。

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