# 多变量线性回归：输入 huyiba_rawdata，结局为有限数值。
outcome_var <- "outcome"
predictor_vars <- c("predictor_continuous", "predictor_group")
reference_levels <- c(predictor_group = "reference")
# 每个无序因子声明一个参考水平；没有因子时用 character(0)。
conf_level <- 0.95 # 系数的 t 区间；提高置信水平通常会使区间变宽。

# 只选所需列；分析样本由研究者先确定，缺失值不由模型静默删除。
model_columns <- c(outcome_var, predictor_vars)
if (!length(predictor_vars) || anyDuplicated(model_columns) ||
    !all(model_columns %in% names(huyiba_rawdata))) {
  stop("请指定一个结局及至少一个不同的预测列；列名必须存在且不可重复。")
}
huyiba_analysis_data <- huyiba_rawdata[model_columns]
is_finite_numeric <- function(x) is.numeric(x) && is.null(dim(x)) && all(is.finite(x))
valid_predictors <- vapply(
  huyiba_analysis_data[predictor_vars],
  function(x) is_finite_numeric(x) || (is.factor(x) && !is.ordered(x)),
  logical(1)
)
if (anyNA(huyiba_analysis_data) ||
    !is_finite_numeric(huyiba_analysis_data[[outcome_var]]) ||
    !all(valid_predictors)) {
  stop("请提供无缺失的有限数值结局；预测列须为有限数值或无序因子。")
}

# 用显式处理对比，使每个分类系数都相对所声明的参考水平。
factor_vars <- predictor_vars[vapply(
  huyiba_analysis_data[predictor_vars], is.factor, logical(1)
)]
if (!setequal(names(reference_levels), factor_vars) ||
    anyDuplicated(names(reference_levels))) {
  stop("请为每个因子声明一个参考水平。")
}
model_contrasts <- NULL
for (variable in factor_vars) {
  huyiba_analysis_data[[variable]] <- stats::relevel(
    droplevels(huyiba_analysis_data[[variable]]), ref = reference_levels[[variable]]
  )
  model_contrasts[[variable]] <- stats::contr.treatment(
    levels(huyiba_analysis_data[[variable]]), base = 1L
  )
}

# 加性模型含截距；每项系数是在其余项固定时的平均结局差。
model_formula <- stats::reformulate(predictor_vars, response = outcome_var)
huyiba_model <- stats::lm(
  model_formula, data = huyiba_analysis_data, contrasts = model_contrasts,
  na.action = stats::na.fail, singular.ok = FALSE
)
if (stats::df.residual(huyiba_model) <= 0L) {
  stop("没有剩余自由度，无法估计系数区间。")
}
coefficient_summary <- summary(huyiba_model)$coefficients
coefficient_ci <- stats::confint(huyiba_model, level = conf_level)

# 保留计算精度：连续预测列按其原单位增加 1，分类列按对比解释。
huyiba_result <- data.frame(
  term = rownames(coefficient_summary),
  estimate = coefficient_summary[, "Estimate"],
  standard_error = coefficient_summary[, "Std. Error"],
  ci_lower = coefficient_ci[, 1],
  ci_upper = coefficient_ci[, 2],
  p_value = coefficient_summary[, "Pr(>|t|)"],
  row.names = NULL
)
print(stats::nobs(huyiba_model))
print(huyiba_result)