# 2.3.1 Logistic回归：拟合二分类模型，提取OR和拟合样本概率。
# 输入：数据框 huyiba_rawdata；结局为数值0/1，1是阳性事件。
# 连续预测变量为有限数值，分类预测变量为无序factor；所用列无缺失。

# 1. 集中修改变量角色、分类参考水平与置信水平。
outcome_var <- "outcome"
predictor_vars <- c("predictor_continuous", "predictor_group")
# 每个分类预测变量都要声明参考水平；无分类变量时用character(0)。
reference_levels <- c(predictor_group = "reference")
conf_level <- 0.95                    # OR的Wald置信水平；提高时区间通常更宽。
# 连续变量的OR按原列每增加1单位解释；换单位应在提供输入前完成。

# 2. 仅取结局和入模变量，拒绝缺失及不符合输入契约的编码。
huyiba_analysis_data <- huyiba_rawdata[c(outcome_var, predictor_vars)]
if (anyNA(huyiba_analysis_data) ||
    !is.numeric(huyiba_analysis_data[[outcome_var]]) ||
    !setequal(huyiba_analysis_data[[outcome_var]], c(0, 1))) {
  stop("所用列须无缺失，结局须为同时含0和1的数值列。")
}
valid_predictors <- vapply(
  huyiba_analysis_data[predictor_vars],
  function(x) {
    (is.numeric(x) && all(is.finite(x))) || (is.factor(x) && !is.ordered(x))
  },
  logical(1)
)
if (!all(valid_predictors)) {
  stop("连续预测变量须为有限数值，分类预测变量须为无序factor。")
}

# 3. 明确所有分类变量的参考组，并固定使用处理对比。
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("reference_levels须为每个分类预测变量恰好指定一个参考水平。")
}
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
  )
}

# 4. 按角色生成加性模型公式；不在本模板中筛选变量。
model_formula <- stats::reformulate(predictor_vars, response = outcome_var)
huyiba_model <- stats::glm(
  model_formula, data = huyiba_analysis_data,
  family = stats::binomial(link = "logit"),
  contrasts = model_contrasts,         # 在拟合时固定分类对比，不修改输入因子的属性。
  # binomial指定二分类分布；logit把事件概率的对数优势与预测变量相连。
  # probit、cloglog是其他连接函数，改用后exp(系数)不能按本页OR解释。
  na.action = stats::na.fail,           # 遇到缺失即停止，避免样本量悄然变化。
  singular.ok = FALSE                 # 设计矩阵不满秩时停止，不报告不可估计项。
)
coefficient_summary <- summary(huyiba_model)$coefficients
if (!isTRUE(huyiba_model$converged) || isTRUE(huyiba_model$boundary) ||
    any(!is.finite(coefficient_summary))) {
  stop("模型未得到可用的有限估计，请检查收敛、分离及共线性问题。")
}

# 5. exp(beta)为OR；confint.default在系数尺度构造Wald近似区间。
# confint(model)通常给出轮廓似然区间，不能与这里的Wald区间混称。
or_intervals <- exp(stats::confint.default(huyiba_model, level = conf_level))
report_terms <- setdiff(rownames(coefficient_summary), "(Intercept)")
huyiba_result <- data.frame(
  term = report_terms,
  beta = coefficient_summary[report_terms, "Estimate"],
  standard_error = coefficient_summary[report_terms, "Std. Error"],
  odds_ratio = exp(coefficient_summary[report_terms, "Estimate"]),
  ci_lower = or_intervals[report_terms, 1],
  ci_upper = or_intervals[report_terms, 2],
  p_value = coefficient_summary[report_terms, "Pr(>|z|)"],
  row.names = NULL
)

# 6. response返回结局=1的概率；link返回对数优势，不能当作概率。
# 此处每条记录均参与了拟合，结果仅用于查看模型计算，不评价预测性能。
huyiba_fitted_probability <- stats::predict(huyiba_model, type = "response")
print(table(huyiba_analysis_data[[outcome_var]]))
print(coefficient_summary)
print(huyiba_result, row.names = FALSE)
print(utils::head(huyiba_fitted_probability))