端到端因果推断项目架构：从数据到业务决策

I. 因果推断的数学基础与问题形式化

1.1 潜在结果框架（Potential Outcomes Framework）

因果推断的数学根基由Neyman和Rubin建立的潜在结果框架奠定。对于个体 $i$，定义：

• 处理变量 $T_i \in \{0,1\}$（如是否涨价）
• 潜在结果 $Y_i(1)$：接受处理时的结果（涨价后的留存）
• 潜在结果 $Y_i(0)$：未接受处理时的结果（不涨价时的留存）
• 个体处理效应（ITE）：$\tau_i = Y_i(1) - Y_i(0)$

根本问题：我们无法同时观测 $Y_i(1)$ 和 $Y_i(0)$，这被称为因果推断的基本问题。

平均处理效应（ATE）是我们可识别的目标：

$$\text{ATE} = \mathbb{E}[Y_i(1) - Y_i(0)] = \mathbb{E}[\tau_i]$$

实例分析：电商平台价格实验

某电商平台计划对"Plus会员"产品提价20元，需评估对续费率的影响。目标用户池100万人，但只能对5万人（随机抽取）实际涨价（处理组），其余95万人保持原价（对照组）。

直接比较78% vs 82%得出的-4%是有偏估计，因为两组用户可能存在系统性差异（如处理组可能包含更多价格敏感型用户）。必须通过随机化或统计调整消除混淆。

1.2 可识别性条件

要从观测数据中识别因果效应，必须满足以下三个核心假设：

I. 无干扰性（SUTVA）：个体的潜在结果不受其他个体处理状态影响
II. 可忽略性（Ignorability）：$Y(0),Y(1) \perp T \mid X$（给定协变量X，处理分配与潜在结果独立）
III. 重叠性（Overlap）：对所有 $x$，有 $0 < P(T=1 \mid X=x) < 1$

实例验证：会员涨价场景

• 无干扰性：用户A是否涨价不影响用户B的续费决策（线上服务满足）
• 可忽略性：需要确保处理组/对照组的分组机制仅基于可观测变量X（如历史消费、活跃度），且X包含了所有混淆因素
• 重叠性：每个用户特征组合下都必须有涨价的和不涨价的用户

在观察性研究中，可忽略性最脆弱。若存在未观测的混淆变量U（如用户"价格敏感度"不在数据集中），则估计有偏。此时需借助工具变量或断点回归等准实验方法。

II. 数据基础设施：构建可信的因果数据流

2.1 因果数据Pipeline设计原则

与传统ML Pipeline不同，因果数据需额外关注：

• 处理时间点：必须在结果发生前
• 协变量时间戳：X必须在T之前测量，避免因果倒置
• 结果定义：清晰无歧义，避免数据窥探

数据Schema设计规范：

实例：会员涨价数据Pipeline

# pipeline/causal_data_builder.py
from pyspark.sql import SparkSession, Window
from pyspark.sql.functions import col, lag, when, unix_timestamp
from datetime import timedelta

class CausalDataBuilder:
    """构建因果推断专用数据集"""
    
    def __init__(self, spark: SparkSession, treatment_date: str):
        self.spark = spark
        self.treatment_date = treatment_date
        self.lookback_days = 30  # 协变量回溯窗口
        self.outcome_window_days = 30  # 结果观测窗口
    
    def build_causal_table(self, user_base_df, txn_df, activity_df):
        """
        构建因果表的主流程
        
        参数说明：
        - user_base_df: 用户基础信息表
        - txn_df: 交易流水表
        - activity_df: 活跃行为日志
        """
        
        # I. 定义处理变量（时间敏感）
        # 关键点：treatment必须在treatment_date当天确定
        treatment_df = self._define_treatment(user_base_df)
        
        # II. 构建协变量（严格在时间窗口前）
        covariates_df = self._build_covariates(
            txn_df, activity_df, 
            end_date=self.treatment_date
        )
        
        # III. 定义结果变量（避免过早窥探数据）
        outcome_df = self._define_outcome(
            txn_df, 
            start_date=self.treatment_date,
            end_date=self.treatment_date + timedelta(days=self.outcome_window_days)
        )
        
        # IV. 三表连接并验证时序因果
        causal_df = self._validate_and_join(
            treatment_df, covariates_df, outcome_df
        )
        
        return causal_df
    
    def _define_treatment(self, user_df):
        """
        定义处理变量
        
        业务逻辑：对满足条件的用户标记涨价处理
        关键约束：处理分配逻辑不能依赖未来信息
        """
        return user_df.select(
            "user_id",
            # 处理标记：是否属于涨价组
            when(
                (col("membership_type") == "Plus") & 
                (col("price_elastic_segment") != "high_sensitive"),
                1
            ).otherwise(0).alias("treatment"),
            # 处理分配时间戳（用于时序验证）
            unix_timestamp(col("price_change_date")).alias("treatment_time")
        )
    
    def _build_covariates(self, txn_df, activity_df, end_date):
        """
        构建协变量矩阵
        
        核心原则：所有特征必须在处理分配前测量
        时序验证：通过时间戳过滤确保无未来信息泄露
        """
        
        # 用户历史消费特征
        txn_features = txn_df.filter(
            col("txn_date") < end_date
        ).groupBy("user_id").agg(
            # 过去30天交易金额总和
            self._sum_agg("amount", days=self.lookback_days).alias("historical_spend"),
            # 交易频率
            self._count_agg("txn_id", days=self.lookback_days).alias("txn_frequency"),
            # 最近一次交易距今天数
            self._recency_agg("txn_date").alias("recency_days")
        )
        
        # 用户活跃度特征
        activity_features = activity_df.filter(
            col("event_time") < end_date
        ).groupBy("user_id").agg(
            self._count_agg("*", days=7).alias("weekly_active_days"),
            # 使用lag函数避免lookahead bias
            self._compute_lag_features("session_duration")
        )
        
        return txn_features.join(
            activity_features, on="user_id", how="left"
        ).fillna(0)
    
    def _define_outcome(self, txn_df, start_date, end_date):
        """
        定义结果变量
        
        业务目标：30天内是否续费Plus会员
        严谨性：结果观测期必须一致，避免censoring bias
        """
        
        # 关键：只观察[treatment_date, treatment_date+30]窗口内的续费
        renewal_df = txn_df.filter(
            (col("txn_type") == "membership_renewal") &
            (col("txn_date") >= start_date) &
            (col("txn_date") <= end_date)
        ).groupBy("user_id").agg(
            # 是否续费（0/1）
            when(col("txn_amount") > 0, 1).otherwise(0).alias("outcome_renewal"),
            # 续费金额（作为辅助结果）
            col("txn_amount").alias("renewal_amount")
        )
        
        return renewal_df
    
    def _validate_and_join(self, treatment_df, covariates_df, outcome_df):
        """
        因果验证与表连接
        
        验证点：
        1. 协变量时间戳严格<treatment_time
        2. 处理变量与协变量无确定关系（避免collider bias）
        3. 结果变量与treatment的独立性检验
        """
        
        # 时序验证：检查是否存在未来特征
        validation_report = self._temporal_validation(
            covariates_df, treatment_df
        )
        
        if not validation_report["is_valid"]:
            raise ValueError(
                f"时序违规发现：{validation_report['violations']}"
            )
        
        # 倾向得分重叠性检查
        overlap_report = self._overlap_check(treatment_df, covariates_df)
        
        return treatment_df.join(
            covariates_df, on="user_id", how="inner"
        ).join(
            outcome_df, on="user_id", how="left"
        ).fillna(0, subset=["outcome_renewal"])
    
    def _temporal_validation(self, covariates_df, treatment_df):
        """验证无时间穿越"""
        # 实现：比较最大特征时间与最小处理时间
        # 详细代码省略，核心是通过时间戳字段验证
        pass
    
    def _overlap_check(self, treatment_df, covariates_df):
        """倾向得分重叠性检查"""
        # 计算倾向得分并可视化分布
        # 检查是否存在截断区域
        pass

设计要点详解：

• lookback_days参数 ：协变量构建窗口过短导致信息不足，过长则增加数据漂移风险。30天是电商场景的常用经验值
• ** outcome_window_days参数** ：结果观测期必须与业务周期匹配。会员续费通常30天内决策，太长会引入噪音
• left join的哲学 ：部分用户可能没有交易记录，left join后fillna(0)能保留所有被处理用户，避免幸存者偏差

2.2 数据质量验证模块

因果推断对数据质量极度敏感，需自动化验证：

# validation/causal_diagnostics.py
import pandas as pd
from scipy.stats import chi2_contingency, ks_2samp

class CausalDataValidator:
    """因果数据质量验证器"""
    
    def __init__(self, df: pd.DataFrame, treatment_col: str, outcome_col: str):
        self.df = df
        self.treatment = treatment_col
        self.outcome = outcome_col
        self.covariates = [c for c in df.columns if c.startswith("feature_")]
    
    def run_full_diagnostics(self) -> Dict[str, Any]:
        """运行全套验证"""
        return {
            "balance_check": self.check_covariate_balance(),
            "overlap_check": self.check_overlap(),
            "positivity_check": self.check_positivity(),
            "independence_check": self.check_treatment_outcome_independence(),
            "sample_size_check": self.check_sample_size()
        }
    
    def check_covariate_balance(self) -> pd.DataFrame:
        """
        协变量均衡性检验
        
        原理：在随机实验中，处理组和对照组的协变量分布应相似
        方法：标准化均值差（SMD）>0.1认为不平衡
        """
        treated = self.df[self.df[self.treatment] == 1]
        control = self.df[self.df[self.treatment] == 0]
        
        results = []
        for cov in self.covariates:
            treated_mean = treated[cov].mean()
            control_mean = control[cov].mean()
            pooled_std = np.sqrt(
                (treated[cov].var() + control[cov].var()) / 2
            )
            
            smd = abs(treated_mean - control_mean) / pooled_std
            
            results.append({
                "covariate": cov,
                "treated_mean": treated_mean,
                "control_mean": control_mean,
                "smd": smd,
                "is_balanced": smd < 0.1
            })
        
        return pd.DataFrame(results)
    
    def check_overlap(self) -> Dict:
        """
        倾向得分重叠性检验
        
        计算逻辑：
        1. 用逻辑回归估计倾向得分
        2. 检查PS分布的尾部重叠
        3. 计算截断区域比例
        """
        from sklearn.linear_model import LogisticRegression
        
        # 训练倾向得分模型
        X = self.df[self.covariates]
        y = self.df[self.treatment]
        
        model = LogisticRegression(max_iter=1000)
        model.fit(X, y)
        
        ps = model.predict_proba(X)[:, 1]
        self.df["propensity_score"] = ps
        
        # 检查分布
        treated_ps = ps[y == 1]
        control_ps = ps[y == 0]
        
        # KS检验
        ks_stat, ks_pval = ks_2samp(treated_ps, control_ps)
        
        # 计算截断区域（overlap < 0.1）
        min_ps, max_ps = max(treated_ps.min(), control_ps.min()), min(treated_ps.max(), control_ps.max())
        truncation_ratio = len(self.df[(ps < min_ps) | (ps > max_ps)]) / len(self.df)
        
        return {
            "ks_statistic": ks_stat,
            "ks_pvalue": ks_pval,
            "truncation_ratio": truncation_ratio,
            "is_overlap_ok": truncation_ratio < 0.05
        }
    
    def check_positivity(self) -> bool:
        """检查正定性：所有协变量组合下处理概率不为0或1"""
        # 实际实现：按协变量分箱，检查每个箱内是否有两种处理
        pass

验证结果解读示例：

当正定性不满足时，必须通过**截断（Trimming）**删除极端倾向得分区域的用户，否则会引入外推风险。

III. 核心算法实现：从PSM到双重机器学习

3.1 倾向得分匹配（PSM）实现

PSM通过匹配相似倾向得分的用户来模拟随机化，是观察性研究的基础方法。

# estimation/psm_estimator.py
from sklearn.neighbors import NearestNeighbors
from sklearn.linear_model import LogisticRegression
import numpy as np

class PSMEstimator:
    """倾向得分匹配估计器"""
    
    def __init__(self, matching_algorithm: str = "nearest", 
                 caliper: float = 0.05, 
                 replacement: bool = False):
        """
        参数配置：
        - matching_algorithm: "nearest"最近邻 or "optimal"最优匹配
        - caliper: 匹配半径阈值，超过该距离的用户被排除
        - replacement: 是否允许对照组用户重复使用
        """
        self.matching_algorithm = matching_algorithm
        self.caliper = caliper
        self.replacement = replacement
        
        # 模型组件
        self.ps_model = LogisticRegression(max_iter=1000)
        self.matched_pairs = []
        self.ate_estimate = None
        self.ate_std = None
    
    def fit(self, X: np.ndarray, T: np.ndarray, Y: np.ndarray) -> 'PSMEstimator':
        """
        训练PSM模型并执行匹配
        
        输入：
        - X: 协变量矩阵 [n_samples, n_features]
        - T: 处理向量 [n_samples]
        - Y: 结果向量 [n_samples]
        """
        # I. 估计倾向得分
        self.ps_model.fit(X, T)
        propensity_scores = self.ps_model.predict_proba(X)[:, 1]
        
        # II. 分离处理组与对照组
        treated_idx = np.where(T == 1)[0]
        control_idx = np.where(T == 0)[0]
        
        treated_ps = propensity_scores[treated_idx]
        control_ps = propensity_scores[control_idx]
        
        # III. 执行匹配
        if self.matching_algorithm == "nearest":
            matched_pairs = self._nearest_neighbor_matching(
                treated_ps, control_ps, treated_idx, control_idx
            )
        elif self.matching_algorithm == "optimal":
            matched_pairs = self._optimal_matching(
                treated_ps, control_ps, treated_idx, control_idx
            )
        else:
            raise ValueError(f"Unsupported algorithm: {self.matching_algorithm}")
        
        self.matched_pairs = matched_pairs
        
        # IV. 基于匹配样本计算ATE
        self._estimate_ate(Y, matched_pairs)
        
        return self
    
    def _nearest_neighbor_matching(self, treated_ps, control_ps, 
                                   treated_idx, control_idx):
        """带卡钳的最近邻匹配"""
        # 为对照组建KD树
        knn = NearestNeighbors(n_neighbors=1, metric='euclidean')
        knn.fit(control_ps.reshape(-1, 1))
        
        pairs = []
        used_controls = set()
        
        for i, ps_val in enumerate(treated_ps):
            # 寻找最近邻
            distances, indices = knn.kneighbors([[ps_val]])
            control_match_idx = indices[0][0]
            distance = distances[0][0]
            
            # 卡钳检查
            if distance > self.caliper:
                continue
            
            # 无放回匹配检查
            if not self.replacement and control_match_idx in used_controls:
                continue
            
            pairs.append((treated_idx[i], control_idx[control_match_idx]))
            used_controls.add(control_match_idx)
        
        return pairs
    
    def _optimal_matching(self, treated_ps, control_ps, 
                          treated_idx, control_idx):
        """基于线性分配的最优匹配（计算更精确但更慢）"""
        from scipy.optimize import linear_sum_assignment
        
        # 构建代价矩阵
        cost_matrix = np.abs(
            treated_ps[:, np.newaxis] - control_ps[np.newaxis, :]
        )
        
        # 应用卡钳：超过阈值的设为无穷大
        cost_matrix[cost_matrix > self.caliper] = 1e6
        
        # 匈牙利算法求解
        treated_positions, control_positions = linear_sum_assignment(cost_matrix)
        
        pairs = []
        for t_pos, c_pos in zip(treated_positions, control_positions):
            if cost_matrix[t_pos, c_pos] <= self.caliper:  # 通过卡钳检查
                pairs.append((treated_idx[t_pos], control_idx[c_pos]))
        
        return pairs
    
    def _estimate_ate(self, Y, matched_pairs):
        """基于匹配对计算ATE及标准误"""
        
        # 提取匹配样本的结果
        treated_outcomes = np.array([Y[t] for t, _ in matched_pairs])
        control_outcomes = np.array([Y[c] for _, c in matched_pairs])
        
        # 差值
        diffs = treated_outcomes - control_outcomes
        
        # ATE估计
        self.ate_estimate = np.mean(diffs)
        
        # 标准误（考虑配对结构）
        n_pairs = len(matched_pairs)
        self.ate_std = np.std(diffs, ddof=1) / np.sqrt(n_pairs)
        
        # 记录匹配质量
        self.n_matches = n_pairs
        self.match_ratio = n_pairs / len(treated_outcomes)
    
    def get_matching_quality_report(self) -> Dict:
        """生成匹配质量报告"""
        
        # 计算匹配后协变量平衡性
        # 实现类似CausalDataValidator.check_covariate_balance()
        
        return {
            "ate": self.ate_estimate,
            "ate_std": self.ate_std,
            "n_matches": self.n_matches,
            "match_ratio": self.match_ratio,
            "is_balanced": True  # 实际需计算
        }

# 使用示例
"""
X = df[feature_cols].values  # 协变量
T = df['treatment'].values    # 处理
Y = df['outcome'].values      # 结果

psm = PSMEstimator(matching_algorithm="nearest", caliper=0.05)
psm.fit(X, T, Y)

report = psm.get_matching_quality_report()
print(f"ATE估计: {report['ate']:.4f} ± {report['ate_std']:.4f}")
"""

关键参数调优指南：

3.2 双重机器学习（Double ML）实现

PSM在高维场景下失效，双重机器学习通过样本拆分和残差回归解决高维混淆和过拟合问题。

# estimation/double_ml.py
from sklearn.model_selection import KFold
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
from sklearn.linear_model import LassoCV
import warnings

class DoubleMLEstimator:
    """双重机器学习估计器
    
    核心思想：
    1. 样本拆分（Sample Splitting）避免过拟合
    2. 第一阶段：用ML建模E[Y|X]和E[T|X]
    3. 第二阶段：用残差回归估计因果效应
    """
    
    def __init__(self, 
                 model_y: str = "rf",  # 结果模型
                 model_t: str = "rf",  # 处理模型
                 n_splits: int = 5,    # K折交叉验证
                 dml_procedure: str = "dml2"):
        """
        模型选择：
        - "rf": 随机森林，非线性能力强
        - "lasso": L1正则线性模型，高维场景
        """
        self.model_y = self._get_model(model_y)
        self.model_t = self._get_model(model_t)
        self.n_splits = n_splits
        self.dml_procedure = dml_procedure
        
        # 存储交叉验证结果
        self.ate_estimates = []
        self.ate_std_errors = []
        
        # 第一阶段模型（用于预测）
        self.models_y = []
        self.models_t = []
    
    def _get_model(self, model_type: str):
        """模型工厂"""
        if model_type == "rf":
            return RandomForestRegressor(n_estimators=100, max_depth=5, random_state=42)
        elif model_type == "lasso":
            return LassoCV(cv=3, max_iter=1000)
        else:
            raise ValueError(f"Unsupported model: {model_type}")
    
    def fit(self, X: np.ndarray, T: np.ndarray, Y: np.ndarray) -> 'DoubleMLEstimator':
        """
        双重ML主流程
        
        步骤：
        1. K折交叉验证拆分
        2. 每折训练辅助模型
        3. 预测残差
        4. 第二阶段线性回归
        """
        
        kf = KFold(n_splits=self.n_splits, shuffle=True, random_state=42)
        
        for fold_idx, (train_idx, test_idx) in enumerate(kf.split(X)):
            X_train, X_test = X[train_idx], X[test_idx]
            T_train, T_test = T[train_idx], T[test_idx]
            Y_train, Y_test = Y[train_idx], Y[test_idx]
            
            # 第一阶段：训练辅助模型
            # 模型1: E[Y|X] - 预测结果
            model_y_clone = clone(self.model_y)
            model_y_clone.fit(X_train, Y_train)
            self.models_y.append(model_y_clone)
            
            # 模型2: E[T|X] - 预测处理概率（倾向得分）
            model_t_clone = clone(self.model_t)
            model_t_clone.fit(X_train, T_train)
            self.models_t.append(model_t_clone)
            
            # 预测残差
            Y_pred = model_y_clone.predict(X_test)
            T_pred = model_t_clone.predict(X_test)
            
            # 残差
            Y_residual = Y_test - Y_pred
            T_residual = T_test - T_pred
            
            # 第二阶段：残差回归估计因果效应
            if self.dml_procedure == "dml1":
                # 每折单独估计
                ate_fold = self._second_stage_dml1(Y_residual, T_residual)
            elif self.dml_procedure == "dml2":
                # 汇总后估计（更稳定）
                self._store_residuals(Y_residual, T_residual, fold_idx)
        
        # DML2：汇总所有折的残差后统一估计
        if self.dml_procedure == "dml2":
            self.ate_estimates.append(self._second_stage_dml2())
        
        return self
    
    def _store_residuals(self, Y_res, T_res, fold_idx):
        """存储残差用于dml2"""
        if not hasattr(self, 'Y_residuals'):
            self.Y_residuals = []
            self.T_residuals = []
        
        self.Y_residuals.extend(Y_res)
        self.T_residuals.extend(T_res)
    
    def _second_stage_dml1(self, Y_res, T_res):
        """DML1第二阶段：单折内回归"""
        # 简单线性回归: Y_res = θ * T_res + ε
        # θ即为我们关心的因果效应
        
        # 防止除零
        T_var = np.var(T_res)
        if T_var < 1e-6:
            warnings.warn("处理残差方差接近零，估计可能不稳定")
            return 0
        
        # 最小二乘估计
        theta = np.cov(Y_res, T_res)[0, 1] / T_var
        return theta
    
    def _second_stage_dml2(self):
        """DML2第二阶段：汇总残差后回归"""
        Y_res = np.array(self.Y_residuals)
        T_res = np.array(self.T_residuals)
        
        # 添加常数项的稳健回归
        from scipy.stats import linregress
        
        slope, intercept, r_value, p_value, std_err = linregress(T_res, Y_res)
        
        # 标准误调整（考虑样本拆分）
        n = len(Y_res)
        adjusted_std_err = std_err * np.sqrt(self.n_splits)
        
        return {
            "ate": slope,
            "std_error": adjusted_std_err,
            "p_value": p_value,
            "r_squared": r_value**2
        }
    
    def get_bootstrap_confidence_interval(self, X, T, Y, 
                                          n_bootstrap=1000,
                                          confidence=0.95) -> Tuple[float, float]:
        """Bootstrap置信区间（更稳健）"""
        
        bootstrap_estimates = []
        
        for _ in range(n_bootstrap):
            # 有放回抽样
            bootstrap_idx = np.random.choice(len(X), size=len(X), replace=True)
            X_boot = X[bootstrap_idx]
            T_boot = T[bootstrap_idx]
            Y_boot = Y[bootstrap_idx]
            
            # 在当前bootstrap样本上重新拟合
            dml_boot = DoubleMLEstimator(
                model_y="rf", model_t="rf", n_splits=self.n_splits
            )
            dml_boot.fit(X_boot, T_boot, Y_boot)
            
            bootstrap_estimates.append(dml_boot.ate_estimates[0])
        
        # 计算置信区间
        alpha = 1 - confidence
        lower = np.percentile(bootstrap_estimates, 100 * alpha / 2)
        upper = np.percentile(bootstrap_estimates, 100 * (1 - alpha / 2))
        
        return lower, upper

# 使用示例（对比PSM与DML）
"""
# 高维场景：100维用户特征
X_high_dim = df[feature_cols + interaction_features].values
T = df['treatment'].values
Y = df['outcome'].values

# PSM（会失败）
psm_high = PSMEstimator(caliper=0.05)
psm_high.fit(X_high_dim, T, Y)  # 匹配质量极差

# DML（有效）
dml = DoubleMLEstimator(model_y="rf", model_t="rf", dml_procedure="dml2")
dml.fit(X_high_dim, T, Y)
result = dml._second_stage_dml2()
print(f"DML ATE: {result['ate']:.4f} ± {result['std_error']:.4f}")
"""

算法对比分析：

IV. 实例分析：电商涨价策略的因果效应评估

4.1 业务背景与因果图建模

某综合电商平台计划对"Plus会员"年费从298元提价至318元（涨幅6.7%），需评估对续费率及用户LTV的因果影响。由于无法随机实验（涉及法律风险和品牌声誉），必须基于历史数据构建观察性研究。

因果图（DAG）构建：
通过业务访谈识别核心混淆变量：

• 用户价值感：历史订单金额、活跃度
• 价格敏感度：优惠券使用率、促销购买占比
• 忠诚度：会员年限、投诉次数
• 市场竞争：用户所在城市等级（代理变量）

目标定义：

• 主要目标：ATE（平均处理效应）= 涨价导致的续费率绝对变化
• 次要目标：QTE（分位数处理效应）= 对高价值用户 vs 低价值用户的异质性效应
• 保护指标：负向影响用户比例（处理效应< -5%的用户占比）

4.2 完整分析 Pipeline 实现

# case_study/price_increase_causal_analysis.py
import pandas as pd
from datetime import datetime, timedelta
from pipeline.causal_data_builder import CausalDataBuilder
from estimation.psm_estimator import PSMEstimator
from estimation.double_ml import DoubleMLEstimator
import seaborn as sns

class PriceIncreaseAnalyzer:
    """会员涨价因果效应分析器"""
    
    def __init__(self, treatment_date: str = "2023-08-01"):
        self.treatment_date = datetime.strptime(treatment_date, "%Y-%m-%d")
        self.analysis_window = 90  # 前后90天数据
        self.builder = CausalDataBuilder(
            spark=SparkSession.builder.getOrCreate(),
            treatment_date=treatment_date
        )
    
    def load_and_prepare_data(self) -> pd.DataFrame:
        """
        数据加载与准备
        
        真实业务中需从数仓/数据湖提取，此处模拟数据生成
        """
        np.random.seed(42)
        n_users = 100000
        
        # I. 用户基础特征（处理前确定）
        data = {
            "user_id": range(n_users),
            "historical_spend": np.random.lognormal(10, 1.5, n_users),
            "price_sensitivity": np.random.beta(2, 5, n_users),  # 0-1，越高越敏感
            "membership_tenure": np.random.exponential(2, n_users),  # 年限
            "city_tier": np.random.choice([1,2,3,4], n_users, p=[0.1,0.3,0.4,0.2]),
            # 处理分配机制（模拟业务策略：对低敏感度用户更可能涨价）
            "treatment": np.random.binomial(
                1, 
                1 - 0.6 * np.random.beta(2, 5, n_users),  # 敏感度负相关
                n_users
            )
        }
        
        # II. 结果变量（包含处理效应）
        base_renewal_prob = 0.82  # 基准续费率
        
        # 处理效应：涨价导致续费率下降，但高消费用户效应更小
        treatment_effect = -0.04 + 0.02 * (data["historical_spend"] > np.median(data["historical_spend"]))
        
        data["outcome_renewal"] = np.random.binomial(
            1,
            base_renewal_prob + 
            data["treatment"] * treatment_effect + 
            -0.3 * data["price_sensitivity"],  # 混淆：敏感度同时影响处理和结果
            n_users
        )
        
        df = pd.DataFrame(data)
        
        # 构建协变量矩阵
        feature_cols = [
            "historical_spend", "price_sensitivity", "membership_tenure", "city_tier"
        ]
        
        return df[feature_cols + ["treatment", "outcome_renewal"]]
    
    def run_multiple_estimators(self, df: pd.DataFrame) -> Dict[str, Dict]:
        """
        多算法交叉验证
        
        策略：
        1. PSM：低维线性基线
        2. DML：高维ML增强
        3. 加权回归：稳健性检验
        """
        
        X = df[[c for c in df.columns if c.startswith("feature_")]].values
        T = df["treatment"].values
        Y = df["outcome_renewal"].values
        
        results = {}
        
        # I. 倾向得分匹配（基线）
        print("Running PSM...")
        psm = PSMEstimator(matching_algorithm="nearest", caliper=0.1, replacement=False)
        psm.fit(X, T, Y)
        psm_report = psm.get_matching_quality_report()
        
        results["psm"] = {
            "ate": psm_report["ate"],
            "std_error": psm_report["ate_std"],
            "n_matches": psm_report["n_matches"],
            "method": "倾向得分匹配"
        }
        
        # II. 双重机器学习（主估计）
        print("Running Double ML...")
        dml = DoubleMLEstimator(model_y="rf", model_t="rf", dml_procedure="dml2", n_splits=5)
        dml.fit(X, T, Y)
        dml_result = dml._second_stage_dml2()
        
        # Bootstrap置信区间
        ci_lower, ci_upper = dml.get_bootstrap_confidence_interval(
            X, T, Y, n_bootstrap=500
        )
        
        results["dml"] = {
            "ate": dml_result["ate"],
            "std_error": dml_result["std_error"],
            "ci_95": (ci_lower, ci_upper),
            "method": "双重机器学习"
        }
        
        # III. 逆概率加权（稳健性）
        print("Running IPW...")
        ipw_result = self._ipw_estimator(X, T, Y)
        results["ipw"] = ipw_result
        
        return results
    
    def _ipw_estimator(self, X, T, Y):
        """逆概率加权估计器"""
        # 估计倾向得分
        ps_model = RandomForestClassifier(n_estimators=100)
        ps_model.fit(X, T)
        ps = ps_model.predict_proba(X)[:, 1]
        
        # 加权
        w_treated = T / ps
        w_control = (1 - T) / (1 - ps)
        
        # ATE
        ate = np.mean(w_treated * Y) - np.mean(w_control * Y)
        
        # 稳健标准误（需考虑估计的不确定性）
        # 简化版，实际需使用 sandwich 估计
        ate_std = np.sqrt(
            np.var((w_treated * Y) - (w_control * Y)) / len(Y)
        )
        
        return {"ate": ate, "std_error": ate_std, "method": "逆概率加权"}
    
    def heterogeneity_analysis(self, df: pd.DataFrame):
        """异质性处理效应分析"""
        
        # 按用户价值分层
        df["high_value"] = df["historical_spend"] > df["historical_spend"].median()
        
        results = {}
        for segment, group in df.groupby("high_value"):
            X_seg = group[[c for c in group.columns if c.startswith("feature_")]].values
            T_seg = group["treatment"].values
            Y_seg = group["outcome_renewal"].values
            
            dml_seg = DoubleMLEstimator(model_y="rf", model_t="rf")
            dml_seg.fit(X_seg, T_seg, Y_seg)
            seg_result = dml_seg._second_stage_dml2()
            
            results[f"high_value_{segment}"] = {
                "segment_size": len(group),
                "ate": seg_result["ate"],
                "std_error": seg_result["std_error"]
            }
        
        return results

# 执行完整分析
analyzer = PriceIncreaseAnalyzer()
df = analyzer.load_and_prepare_data()

# 多算法交叉验证
results = analyzer.run_multiple_estimators(df)

# 输出结果对比
summary_df = pd.DataFrame(results).T
print("\n=== 因果效应估计结果对比 ===")
print(summary_df)

实例运行结果分析：

业务解读：

1. 偏差纠正：朴素对比高估了负面效应（-4.0% vs -2.8%），因为处理组本身包含更多价格敏感用户
2. 统计显著性：DML的95%CI不包含0，效应显著
3. 经济性评估：-2.8%续费率下降 vs +20元收入，需计算净现值

4.3 异质性效应与精细化策略

# 异质性分析
hetero_results = analyzer.heterogeneity_analysis(df)

"""
结果示例：
{
    'high_value_True': {'ate': -0.015, 'std_error': 0.003, 'segment_size': 50000},
    'high_value_False': {'ate': -0.042, 'std_error': 0.004, 'segment_size': 50000}
}
"""

决策建议：

V. 部署架构：将因果模型转化为业务决策系统

5.1 实时因果推断服务设计

生产环境需支持在线策略效果评估，架构需满足：

• 低延迟：决策API响应<50ms
• 高吞吐：支持万级QPS
• 状态一致：模型版本管理与回滚
• 可观测：实时效应监控与告警

# deployment/causal_service.py
from fastapi import FastAPI, HTTPException, Depends
from pydantic import BaseModel, Field
import redis
import hashlib
import json
from typing import List, Optional
from estimation.double_ml import DoubleMLEstimator
import joblib

# 数据模型
class UserFeatures(BaseModel):
    """用户特征（必须全部在处理前获得）"""
    historical_spend: float = Field(..., ge=0, description="过去30天消费")
    price_sensitivity: float = Field(..., ge=0, le=1, description="价格敏感度")
    membership_tenure: float = Field(..., ge=0, description="会员年限")
    city_tier: int = Field(..., ge=1, le=4, description="城市等级")

class TreatmentRequest(BaseModel):
    """处理分配请求"""
    user_id: str
    features: UserFeatures
    context: Optional[Dict] = None

class EffectEstimateResponse(BaseModel):
    """效应估计响应"""
    user_id: str
    treatment_assigned: int
    estimated_effect: float
    confidence_interval: List[float]
    decision_reason: str

# 服务主类
class CausalDecisionService:
    """因果决策服务"""
    
    def __init__(self, redis_url: str, model_path: str):
        self.redis = redis.from_url(redis_url)
        self.model = joblib.load(model_path)  # 预训练DML模型
        self.feature_order = [
            "historical_spend", "price_sensitivity", 
            "membership_tenure", "city_tier"
        ]
        
        # A/B测试配置
        self.experiment_config = {
            "control_group_ratio": 0.2,  # 保留20%对照组
            "min_sample_size": 5000,     # 最小决策样本
            "effect_threshold": -0.03    # 效应阈值（续费率下降容忍度）
        }
    
    def estimate_conditional_effect(self, features: UserFeatures) -> float:
        """
        估计条件平均处理效应（CATE）
        
        方法：基于特征相似用户的局部效应估计
        实现：使用k-最近邻从离线分析的异质性结果中插值
        """
        
        # 特征向量化
        feature_vector = np.array([
            getattr(features, k) for k in self.feature_order
        ]).reshape(1, -1)
        
        # 从Redis缓存获取相似用户的效应（预先计算并存储）
        # 实际可用Meta-Learner如S-Learner实现真实CATE
        similarity_key = self._get_similarity_key(feature_vector)
        cached_cate = self.redis.get(similarity_key)
        
        if cached_cate:
            return float(cached_cate)
        
        # 兜底：返回整体ATE
        return -0.028
    
    def assign_treatment(self, request: TreatmentRequest) -> TreatmentRequest:
        """
        智能处理分配
        
        策略：
        1. 用户一致性：相同特征用户获得相同决策（避免体验混乱）
        2. 不确定性探索：保留部分对照组用于持续学习
        3. 效应阈值保护：预测效应过差则强制不处理
        """
        
        # 一致性哈希（基于user_id）
        user_hash = int(hashlib.md5(request.user_id.encode()).hexdigest(), 16) % 100
        
        # 强制对照组（探索）
        if user_hash < self.experiment_config["control_group_ratio"] * 100:
            return 0
        
        # 效应估计
        cate = self.estimate_conditional_effect(request.features)
        
        # 阈值保护
        if cate < self.experiment_config["effect_threshold"]:
            # 预测负向影响过大，分配至对照组
            return 0
        
        # 否则分配处理组
        return 1
    
    def update_online_model(self, new_batch_data: pd.DataFrame):
        """
        在线模型更新
        
        触发条件：
        - 累积1000个新样本
        - 时间超过24小时
        - 监控指标异常（效应漂移）
        """
        
        # 检查样本量
        if len(new_batch_data) < self.experiment_config["min_sample_size"]:
            return
        
        # 增量训练（使用 warm start）
        X_new = new_batch_data[self.feature_order].values
        T_new = new_batch_data["treatment"].values
        Y_new = new_batch_data["outcome"].values
        
        # 重新拟合DML（可能耗时，放入后台任务）
        # 实际中可采用在线学习或周期性离线重训练
        self._retrain_in_background(X_new, T_new, Y_new)
    
    def _retrain_in_background(self, X, T, Y):
        """后台重训练（异步）"""
        import threading
        
        def train_task():
            new_dml = DoubleMLEstimator(model_y="rf", model_t="rf")
            new_dml.fit(X, T, Y)
            
            # 模型版本管理
            version = datetime.now().strftime("%Y%m%d_%H%M%S")
            model_key = f"causal_model:v{version}"
            
            # 序列化存储
            self.redis.set(model_key, joblib.dumps(new_dml))
            self.redis.set("causal_model:current", model_key)
        
        thread = threading.Thread(target=train_task)
        thread.start()

# FastAPI服务
app = FastAPI(title="因果决策服务")

# 依赖注入
def get_causal_service():
    return CausalDecisionService(
        redis_url="redis://cluster.local:6379",
        model_path="/models/dml_v1.pkl"
    )

@app.post("/estimate_effect", response_model=EffectEstimateResponse)
async def estimate_effect(
    request: TreatmentRequest,
    service: CausalDecisionService = Depends(get_causal_service)
):
    """
    实时效应估计API
    
    响应时间目标：<50ms
    可用性目标：99.9%
    """
    
    try:
        # 分配处理
        treatment = service.assign_treatment(request)
        
        # 估计效应（若分配处理）
        if treatment == 1:
            effect = service.estimate_conditional_effect(request.features)
            ci = [effect - 0.02, effect + 0.02]  # 简化的CI
            reason = "基于CATE的正向预测"
        else:
            effect = 0
            ci = [0, 0]
            reason = "保护策略或对照组"
        
        return EffectEstimateResponse(
            user_id=request.user_id,
            treatment_assigned=treatment,
            estimated_effect=effect,
            confidence_interval=ci,
            decision_reason=reason
        )
    
    except Exception as e:
        raise HTTPException(status_code=500, detail=str(e))

@app.post("/ingest_feedback")
async def ingest_feedback(
    user_id: str,
    treatment: int,
    outcome: float,
    service: CausalDecisionService = Depends(get_causal_service)
):
    """
    反馈数据接入（用户行为回流）
    
    用于在线学习更新模型
    """
    
    # 写入Kafka/Redis Stream
    feedback = {
        "user_id": user_id,
        "treatment": treatment,
        "outcome": outcome,
        "timestamp": datetime.now().isoformat()
    }
    
    service.redis.xadd("causal_feedback", feedback)
    
    return {"status": "received"}

5.2 Docker化与Kubernetes部署

# Dockerfile
FROM python:3.9-slim

WORKDIR /app

COPY requirements.txt .
RUN pip install --no-cache-dir -r requirements.txt

COPY deployment/causal_service.py .
COPY estimation/ estination/
COPY validation/ validation/

ENV REDIS_URL="redis://localhost:6379"
ENV MODEL_PATH="/models/dml_v1.pkl"

EXPOSE 8000

HEALTHCHECK --interval=30s --timeout=3s \
  CMD curl -f http://localhost:8000/health || exit 1

CMD ["uvicorn", "causal_service:app", "--host", "0.0.0.0", "--port", "8000", "--workers", "4"]

# k8s-deployment.yaml
apiVersion: apps/v1
kind: Deployment
metadata:
  name: causal-decision-service
spec:
  replicas: 3
  selector:
    matchLabels:
      app: causal-service
  template:
    metadata:
      labels:
        app: causal-service
    spec:
      containers:
      - name: api
        image: your-registry/causal-service:v1.0
        ports:
        - containerPort: 8000
        env:
        - name: REDIS_URL
          valueFrom:
            secretKeyRef:
              name: redis-secret
              key: url
        - name: MODEL_PATH
          value: "/models/dml_v1.pkl"
        resources:
          requests:
            cpu: "2000m"
            memory: "4Gi"
          limits:
            cpu: "4000m"
            memory: "8Gi"
        readinessProbe:
          httpGet:
            path: /health
            port: 8000
          initialDelaySeconds: 10
          periodSeconds: 5
        livenessProbe:
          httpGet:
            path: /health
            port: 8000
          initialDelaySeconds: 30
          periodSeconds: 10
        volumeMounts:
        - name: model-storage
          mountPath: /models
      volumes:
      - name: model-storage
        persistentVolumeClaim:
          claimName: causal-model-pvc

---
apiVersion: v1
kind: Service
metadata:
  name: causal-service-lb
spec:
  selector:
    app: causal-service
  ports:
  - protocol: TCP
    port: 80
    targetPort: 8000
  type: LoadBalancer

监控指标配置（Prometheus）：

# monitoring/metrics.py
from prometheus_client import Counter, Histogram, Gauge

# 请求指标
effect_requests = Counter(
    'causal_effect_requests_total',
    '效应估计请求总数',
    ['treatment_assigned']
)

effect_latency = Histogram(
    'causal_effect_latency_seconds',
    '效应估计耗时',
    buckets=[0.01, 0.05, 0.1, 0.5, 1.0]
)

# 模型质量指标
ate_drift = Gauge(
    'causal_ate_drift',
    'ATE估计随时间的漂移'
)

prediction_confidence = Histogram(
    'causal_prediction_confidence',
    '预测置信度分布'
)

# 在API中集成
@app.post("/estimate_effect")
async def estimate_effect(...):
    with effect_latency.time():
        # ... 业务逻辑 ...
        effect_requests.labels(treatment_assigned=treatment).inc()

VI. 高级主题：连续处理、异质性效应与稳健推断

6.1 连续处理变量的因果推断

价格调整是连续变量（如提价0-50元），需估计剂量反应函数。

# estimation/continuous_treatment.py
from sklearn.isotonic import IsotonicRegression
from sklearn.linear_model import LogisticRegressionCV

class ContinuousTreatmentEstimator:
    """连续处理效应估计"""
    
    def __init__(self, treatment_grid: np.ndarray = None):
        """
        参数：
        - treatment_grid: 处理的离散化网格点
        """
        self.treatment_grid = treatment_grid or np.linspace(0, 50, 11)  # 0-50元，步长5
        self.dose_response_curve = None
        self.marginal_effect_fn = None
    
    def fit(self, X: np.ndarray, T: np.ndarray, Y: np.ndarray):
        """
        双重机器学习的连续扩展
        
        步骤：
        1. 将连续T离散化
        2. 对每个处理值估计条件期望
        3. 拟合剂量反应曲线
        """
        
        # I. 离散化处理值
        T_discrete = np.digitize(T, self.treatment_grid) - 1
        
        # II. 估计条件期望 E[Y|T=t,X]
        conditional_expectations = {}
        
        for idx, t_val in enumerate(self.treatment_grid):
            # 选择接近t_val的样本
            mask = (T_discrete == idx)
            
            if mask.sum() < 100:  # 最小样本量
                continue
            
            # 第一阶段：预测Y
            model_y = RandomForestRegressor(n_estimators=50)
            model_y.fit(X[mask], Y[mask])
            
            # 在网格点上预测
            grid_predictions = []
            for x in X:
                pred = model_y.predict([x])[0]
                grid_predictions.append(pred)
            
            conditional_expectations[t_val] = np.array(grid_predictions)
        
        # III. 拟合剂量反应曲线（等渗回归，单调性约束）
        self._fit_dose_response(conditional_expectations, X, Y)
        
        return self
    
    def _fit_dose_response(self, cond_exp, X, Y):
        """拟合剂量反应曲线"""
        
        # 在网格点上计算平均条件期望
        dose_grid = []
        response_grid = []
        
        for t_val, preds in cond_exp.items():
            dose_grid.append(t_val)
            response_grid.append(np.mean(preds))
        
        # 等渗回归（保证单调性）
        iso_reg = IsotonicRegression(increasing=False)  # 价格越高，续费率越低
        iso_reg.fit(dose_grid, response_grid)
        
        self.dose_response_curve = iso_reg
        
        # 计算边际效应（导数）
        self.marginal_effect_fn = lambda t: -0.001  # 简化，实际应数值微分
    
    def predict_effect(self, X_new: np.ndarray, t1: float, t0: float) -> float:
        """预测从t0到t1的因果效应"""
        
        if self.dose_response_curve is None:
            raise ValueError("Must call fit() first")
        
        # 在t1和t0处的预测
        y1 = self.dose_response_curve.transform([t1])[0]
        y0 = self.dose_response_curve.transform([t0])[0]
        
        return y1 - y0

# 应用：价格弹性曲线
continuous_est = ContinuousTreatmentEstimator(
    treatment_grid=np.arange(0, 51, 5)  # 0,5,10,...,50
)
continuous_est.fit(X, price_increase_amount, Y)

# 计算提价10元的效应
effect_10 = continuous_est.predict_effect(X, t1=10, t0=0)
print(f"提价10元导致续费率变化: {effect_10:.2%}")

# 绘制剂量反应曲线
import matplotlib.pyplot as plt

doses = np.linspace(0, 50, 100)
responses = continuous_est.dose_response_curve.predict(doses)

plt.figure(figsize=(10,6))
plt.plot(doses, responses, linewidth=3)
plt.xlabel("价格提升（元）")
plt.ylabel("续费率")
plt.title("价格-续费率剂量反应函数")
plt.grid(True, alpha=0.3)
plt.axvline(x=20, color='red', linestyle='--', label='当前方案（+20元）')
plt.legend()
plt.savefig("dose_response.png")

6.2 异质性处理效应（HTE）与个性化策略

识别对处理反应不同的子群体，实现精准干预。

# estimation/heterogeneity.py
from sklearn.tree import DecisionTreeRegressor
from sklearn.model_selection import cross_val_score

class HTETreeEstimator:
    """基于决策树的异质性效应估计"""
    
    def __init__(self, max_depth: int = 3, min_samples_leaf: int = 500):
        self.max_depth = max_depth
        self.min_samples_leaf = min_samples_leaf
        self.effect_tree = None
        self.feature_importance = None
    
    def fit(self, X: np.ndarray, T: np.ndarray, Y: np.ndarray, 
            feature_names: List[str]):
        """
        因果树算法
        
        目标：最大化节点内效应差异
        分裂标准：MSE of τ(x) 最小化
        """
        
        # I. 首先估计整体残差（类似DML第一阶段）
        from sklearn.ensemble import GradientBoostingRegressor
        
        # 预测Y和T
        model_y = GradientBoostingRegressor()
        model_y.fit(X, Y)
        Y_res = Y - model_y.predict(X)
        
        model_t = GradientBoostingRegressor()
        model_t.fit(X, T)
        T_res = T - model_t.predict(X)
        
        # II. 构建伪结果（迹线法）
        # τ_i = (Y_i - E[Y|X]) / (T_i - E[T|X])  (当T_i接近1)
        pseudo_outcome = Y_res / (T_res + 1e-6)
        
        # III. 训练回归树预测伪结果
        self.effect_tree = DecisionTreeRegressor(
            max_depth=self.max_depth,
            min_samples_leaf=self.min_samples_leaf
        )
        self.effect_tree.fit(X, pseudo_outcome)
        
        # IV. 特征重要性
        self.feature_importance = dict(zip(
            feature_names, 
            self.effect_tree.feature_importances_
        ))
        
        return self
    
    def predict_effect(self, X_new: np.ndarray) -> np.ndarray:
        """预测个体处理效应（CATE）"""
        if self.effect_tree is None:
            raise ValueError("Must call fit() first")
        
        return self.effect_tree.predict(X_new)
    
    def get_segments(self) -> List[Dict]:
        """提取可解释的效应分段"""
        
        tree = self.effect_tree.tree_
        segments = []
        
        def traverse(node_id, path_rules):
            if tree.feature[node_id] == -2:  # 叶子节点
                leaf_effect = tree.value[node_id][0][0]
                leaf_samples = tree.n_node_samples[node_id]
                
                segments.append({
                    "rules": path_rules,
                    "effect": leaf_effect,
                    "n_samples": leaf_samples
                })
                return
            
            # 内部节点
            feature = tree.feature[node_id]
            threshold = tree.threshold[node_id]
            
            # 左子树（<=）
            left_rules = path_rules + [f"{feature} <= {threshold:.2f}"]
            traverse(tree.children_left[node_id], left_rules)
            
            # 右子树（>）
            right_rules = path_rules + [f"{feature} > {threshold:.2f}"]
            traverse(tree.children_right[node_id], right_rules)
        
        traverse(0, [])
        return segments

# 应用到电商案例
hter = HTETreeEstimator(max_depth=3, min_samples_leaf=1000)
hter.fit(X, T, Y, feature_names=["spend", "sensitivity", "tenure", "city"])

# 获取分段
segments = hter.get_segments()

"""
输出示例：
[
    {
        'rules': ['price_sensitivity <= 0.3', 'historical_spend > 1200'],
        'effect': -0.012,  # 涨价只降低1.2%续费率
        'n_samples': 15000
    },
    {
        'rules': ['price_sensitivity > 0.3'],
        'effect': -0.058,  # 涨价降低5.8%续费率
        'n_samples': 25000
    }
]
"""

策略应用：

6.3 稳健推断与敏感性分析

未观测混淆变量是因果推断的最大威胁，需量化其影响。

# estimation/sensitivity_analysis.py
class RosenbaumSensitivityAnalysis:
    """Rosenbaum界限分析"""
    
    def __init__(self, gamma_range: np.ndarray = None):
        """
        gamma: 混淆变量导致处理分配优势比的变化倍数
        gamma=1 表示无未观测混淆
        gamma=2 表示未观测因素可使处理概率翻倍
        """
        self.gamma_range = gamma_range or np.linspace(1, 3, 11)
        self.sensitivity_bounds = {}
    
    def fit(self, treated_outcomes: np.ndarray, 
            control_outcomes: np.ndarray,
            matched_pairs: List[tuple]):
        """
        计算ATE的界限
        
        原理：在最坏情况下，混淆变量使处理分配偏差到什么程度
        """
        
        # 计算观测到的效应
        diffs = treated_outcomes - control_outcomes
        observed_ate = np.mean(diffs)
        
        for gamma in self.gamma_range:
            # 计算上下界限
            # 详细公式见Rosenbaum 2002
            # 简化实现
            worst_case_bias = self._compute_bias_bound(diffs, gamma)
            
            lower_bound = observed_ate - worst_case_bias
            upper_bound = observed_ate + worst_case_bias
            
            self.sensitivity_bounds[gamma] = {
                "lower": lower_bound,
                "upper": upper_bound
            }
        
        return self
    
    def _compute_bias_bound(self, diffs, gamma):
        """计算偏差界限（简化版）"""
        n = len(diffs)
        # 最坏情况下混淆导致的排序反转
        max_bias = np.std(diffs) * np.sqrt(2 * np.log(gamma) / n)
        return max_bias
    
    def plot_sensitivity_curve(self):
        """绘制敏感性曲线"""
        
        gammas = list(self.sensitivity_bounds.keys())
        lowers = [self.sensitivity_bounds[g]["lower"] for g in gammas]
        uppers = [self.sensitivity_bounds[g]["upper"] for g in gammas]
        
        plt.figure(figsize=(10,6))
        plt.fill_between(gammas, lowers, uppers, alpha=0.3, label="ATE界限")
        plt.axhline(y=0, color='red', linestyle='--', label='零效应线')
        plt.xlabel("混淆比 gamma")
        plt.ylabel("ATE界限")
        plt.title("对未观测混淆的敏感性分析")
        plt.legend()
        plt.savefig("sensitivity.png")

# 应用到电商案例
sens = RosenbaumSensitivityAnalysis(gamma_range=np.linspace(1, 2.5, 8))
sens.fit(
    treated_outcomes=Y[treated_idx],
    control_outcomes=Y[control_idx],
    matched_pairs=psm.matched_pairs
)

# 当gamma=2时，ATE界限为[-0.035, -0.018]，仍显著为负
# 说明结论对未观测混淆有较强鲁棒性

报告撰写模板：