比值比（Odds Ratio）置信区间计算异常问题排查

我正在做统计分析，计算多个变量在特定时间点的比值比（Odds Ratio, OR）和95%置信区间（CI）。目前OR计算正确，但CI结果异常，怀疑问题出在总系数或标准误的计算逻辑，尤其是交互项求和环节。

简化代码 #

import numpy as np
import pandas as pd
# Define the time points and variables
time_points = [6,12]
variables = ['var1', 'var2']
# Initialize dictionaries for odds ratios and confidence intervals
odds_ratios = {var: [] for var in variables}
conf_intervals = {var: [] for var in variables}
# Create a DataFrame for time points and generate spline terms (function not shown)
all_times = pd.DataFrame({'Time': time_points})
all_spline_terms = create_spline_terms(all_times, column='Time', num_knots=6)
# Loop through time points
for time in time_points:
   base_data = data.copy()
   base_data['Time'] = time
   for col in all_spline_terms.columns:
       base_data[col] = all_spline_terms.loc[all_times['Time'] == time, col].values[0]
   for var in variables:
       # Predict probabilities for var = 1 and var = 0
       base_data_var_1 = base_data.copy()
       base_data_var_1[var] = 1
       proba_1 = result.predict(base_data_var_1)
       odds_1 = proba_1 / (1 - proba_1)
       base_data_var_0 = base_data.copy()
       base_data_var_0[var] = 0
       proba_0 = result.predict(base_data_var_0)
       odds_0 = proba_0 / (1 - proba_0)
       # Calculate the odds ratio
       odds_ratio = odds_1.mean() / odds_0.mean() if odds_0.mean() != 0 else np.nan
       odds_ratios[var].append(odds_ratio)
       # Calculate Wald 95% confidence intervals
       coef_var = result.params[var]
       std_error_var = result.bse[var]
       interaction_coef = 0
       interaction_var = 0
       for col in all_spline_terms.columns:
           interaction_term = f"{var}:{col}"
           if interaction_term in result.params:
               interaction_coef += result.params[interaction_term]
               interaction_var += result.bse[interaction_term] ** 2
       total_coef = coef_var + interaction_coef
       total_std_error = np.sqrt(std_error_var ** 2 + interaction_var)
       lower_ci_log = total_coef - 1.96 * total_std_error
       upper_ci_log = total_coef + 1.96 * total_std_error
       lower_ci = np.exp(lower_ci_log)
       upper_ci = np.exp(upper_ci_log)
       conf_intervals[var].append((lower_ci, upper_ci))
# Output results
for var in variables:
   print(f"\nOdds ratios and 95% CI for {var}:")
   for time, odds_ratio, (lower_ci, upper_ci) in zip(time_points, odds_ratios[var], conf_intervals[var]):
       print(f'  At time {time}: Odds ratio = {odds_ratio:.4f}, 95% CI = ({lower_ci:.4f}, {upper_ci:.4f})')

样条项生成函数 #

def create_spline_terms(df, column='Time', num_knots=6, degree=3):
   knots = np.linspace(df[column].min(), df[column].max(), num_knots + 2)[1:-1]
   spline_terms = patsy.dmatrix(
       f"cr({column}, knots={list(knots)}, constraints='center')",
       df,
       return_type='dataframe'
   )
   return spline_terms.iloc[:, 1:]

具体疑问 #

1. 我的置信区间计算逻辑是否存在错误？
2. 如何确保交互项被正确纳入置信区间计算？
3. 计算Logistic回归中的比值比及置信区间有哪些我忽略的最佳实践？

问题诊断与修正 #

1. 置信区间计算的核心错误 #

你的代码存在两个致命错误：

• 交互项系数未加权：特定时间点的总效应不是主效应系数加所有交互项系数的和，而是主效应系数加上每个交互项系数乘以对应时间点的样条项取值的和——你完全漏掉了样条值的加权步骤，导致所有时间点的总效应计算完全一致，CI自然重复。
• 忽略协方差项：标准误计算不能只加各系数的方差，主效应和交互项、交互项之间存在协方差，这些都会影响总效应的方差。手动拆分方差会严重低估/高估标准误，导致CI异常。

2. 正确纳入交互项的方法 #

要准确计算特定时间点的总效应和CI，必须利用模型的协方差矩阵，自动处理所有系数间的协方差：

• 针对当前时间点，获取对应的样条项取值
• 构造包含主效应和交互项的参数向量
• 构造对应时间点的设计向量（主效应取1，交互项取对应样条值）
• 用设计向量和参数向量计算总效应，用设计向量和协方差矩阵计算总方差
• 基于总效应和标准误计算CI

3. Logistic回归OR与CI计算的最佳实践 #

• 避免手动计算方差：直接用模型的cov_params()获取协方差矩阵，不要手动拆分系数和方差，减少人为错误。
• 明确OR类型：你当前计算的是边际OR（对所有样本的odds取平均再求比值），如果模型包含复杂交互项，要明确你需要的是边际OR还是条件OR。
• 用工具函数验证：可以用statsmodels的get_margeff()方法计算边际效应，再转换为OR，验证自己的计算结果。
• CI方法选择：Wald CI在样本量小或效应极端时表现不佳，可考虑似然比CI作为替代，但大样本下Wald CI足够可靠。

修正后的代码 #

import numpy as np
import pandas as pd
# Define the time points and variables
time_points = [6,12]
variables = ['var1', 'var2']
# Initialize dictionaries for odds ratios and confidence intervals
odds_ratios = {var: [] for var in variables}
conf_intervals = {var: [] for var in variables}
# Create a DataFrame for time points and generate spline terms
all_times = pd.DataFrame({'Time': time_points})
all_spline_terms = create_spline_terms(all_times, column='Time', num_knots=6)
# Loop through time points
for time in time_points:
    base_data = data.copy()
    base_data['Time'] = time
    # Get spline values for current time point
    spline_vals = all_spline_terms.loc[all_times['Time'] == time].values.flatten()
    for col, val in zip(all_spline_terms.columns, spline_vals):
        base_data[col] = val
    for var in variables:
        # Original OR calculation (remains correct)
        base_data_var_1 = base_data.copy()
        base_data_var_1[var] = 1
        proba_1 = result.predict(base_data_var_1)
        odds_1 = proba_1 / (1 - proba_1)
        base_data_var_0 = base_data.copy()
        base_data_var_0[var] = 0
        proba_0 = result.predict(base_data_var_0)
        odds_0 = proba_0 / (1 - proba_0)
        odds_ratio = odds_1.mean() / odds_0.mean() if odds_0.mean() != 0 else np.nan
        odds_ratios[var].append(odds_ratio)
        # Corrected CI calculation
        # Get all relevant parameter names
        param_names = [var] + [f"{var}:{col}" for col in all_spline_terms.columns]
        # Filter parameters that exist in the model
        valid_params = [p for p in param_names if p in result.params.index]
        if not valid_params:
            conf_intervals[var].append((np.nan, np.nan))
            continue
        # Get coefficient vector and design vector
        beta = result.params[valid_params].values
        X = []
        for p in valid_params:
            if p == var:
                X.append(1)
            else:
                # Get corresponding spline column name
                spline_col = p.split(f"{var}:")[-1]
                idx = all_spline_terms.columns.get_loc(spline_col)
                X.append(spline_vals[idx])
        X = np.array(X)
        # Get covariance matrix subset
        cov_matrix = result.cov_params().loc[valid_params, valid_params].values
        # Calculate total effect and variance
        total_coef = X @ beta
        total_variance = X @ cov_matrix @ X.T
        total_std_error = np.sqrt(total_variance)
        # Compute 95% CI
        lower_ci_log = total_coef - 1.96 * total_std_error
        upper_ci_log = total_coef + 1.96 * total_std_error
        lower_ci = np.exp(lower_ci_log)
        upper_ci = np.exp(upper_ci_log)
        conf_intervals[var].append((lower_ci, upper_ci))
# Output results
for var in variables:
    print(f"\nOdds ratios and 95% CI for {var}:")
    for time, odds_ratio, (lower_ci, upper_ci) in zip(time_points, odds_ratios[var], conf_intervals[var]):
        print(f'  At time {time}: Odds ratio = {odds_ratio:.4f}, 95% CI = ({lower_ci:.4f}, {upper_ci:.4f})')

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