观测数据融合：将实验与观测数据结合的分析框架

第一章：观测数据融合的时代背景与理论基石

1.1 纯实验数据的局限性

随机对照实验（A/B Test）被誉为因果推断的金标准，但在真实商业环境中面临严峻挑战：

这些局限导致实验结果的内部有效性（Internal Validity）虽高，但外部有效性（External Validity）严重不足。我们需要将观测数据（Observational Data）作为补充，构建融合分析框架。

1.2 观测数据的独特价值

观测数据指未经随机干预、自然发生的用户行为数据，其价值体现在：

Ⅰ. 海量样本与长期追踪

• 日志数据：覆盖100%用户，时间跨度可达数年
• 可分析用户生命周期价值（LTV）、长期留存等慢指标

Ⅱ. 自然行为模式

• 无实验干预，反映真实用户习惯
• 可发现实验无法触达的"自然实验"场景

Ⅲ. 成本低廉

• 数据已存在，无需额外开发成本
• 支持探索性分析（Exploratory Analysis）

Ⅳ. 外部有效性高

• 包含全量用户分层，结果可直接指导全量策略

Ⅴ. 因果发现（Causal Discovery）

• 通过结构学习发现隐藏的因果关系链
• 识别实验未考虑的混杂因素（Confounders）

1.3 Mermaid：融合分析的理论框架

第二章：技术架构与数据流水线设计

2.1 整体架构分层

构建融合分析平台需要处理异构数据源、实现因果推断算法、保障结果可复现：

2.2 数据融合ETL核心实现

关键在于统一实验数据与观测数据的schema，并标记数据来源。

# etl/unified_data_pipeline.py
from pyspark.sql import SparkSession
from pyspark.sql.functions import (
    col, lit, when, udf, monotonically_increasing_id
)
from pyspark.sql.types import StringType, BooleanType, StructType
import pyspark.sql.functions as F
from datetime import datetime, timedelta

class UnifiedDataPipeline:
    """实验与观测数据融合ETL管道"""
    
    def __init__(self, spark: SparkSession, config: dict):
        self.spark = spark
        self.config = config
        
    def ingest_treatment_data(self, experiment_id: str, start_date: str, end_date: str):
        """Ⅰ. 摄取实验组数据"""
        
        # 从实验平台API获取数据
        exp_df = self.spark.read.format("jdbc") \
            .option("url", self.config['exp_db_url']) \
            .option("query", f"""
                SELECT 
                    user_id,
                    'experiment' as data_source,
                    '{experiment_id}' as experiment_id,
                    group_name,
                    assignment_time as event_time,
                    variant_params,
                    metadata
                FROM experiments.user_assignments
                WHERE experiment_id = '{experiment_id}'
                    AND assignment_time BETWEEN '{start_date}' AND '{end_date}'
            """) \
            .load()
        
        # 添加数据来源标记
        exp_df = exp_df.withColumn(
            "is_randomized", 
            lit(True).cast(BooleanType())
        ).withColumn(
            "propensity_score",  # 实验组倾向得分为实际分配概率
            when(col("group_name") == "control", 0.5)
            .otherwise(0.5)
        )
        
        return exp_df
    
    def ingest_observational_data(self, start_date: str, end_date: str):
        """Ⅱ. 摄取观测数据"""
        
        # 从行为日志获取，用户主动选择使用新功能
        obs_df = self.spark.read.format("delta") \
            .load(self.config['obs_data_path']) \
            .filter(
                (col("event_time") >= start_date) &
                (col("event_time") <= end_date) &
                (col("feature_usage_new_rec") == True)
            ) \
            .select(
                "user_id",
                lit("observational").alias("data_source"),
                lit(None).cast(StringType()).alias("experiment_id"),
                when(col("feature_usage_new_rec") == True, "treatment")
                    .otherwise("control")
                    .alias("group_name"),
                col("event_time"),
                col("user_features").alias("variant_params"),
                F.struct("device_type", "referrer").alias("metadata")
            )
        
        # 观测数据需要计算倾向得分（基于用户特征）
        obs_df = obs_df.withColumn("is_randomized", lit(False))
        
        return obs_df
    
    def calculate_propensity_scores(self, obs_df, features_col="variant_params"):
        """Ⅲ. 为观测数据计算倾向得分"""
        
        # 将变体参数展开为特征向量
        feature_cols = ["device_type", "account_age", "historical_ctr"]
        
        for col_name in feature_cols:
            obs_df = obs_df.withColumn(
                col_name,
                F.get_json_object(col(features_col), f"$.{col_name}")
            )
        
        # 使用PySpark MLlib训练倾向得分模型
        from pyspark.ml.feature import VectorAssembler, StringIndexer
        from pyspark.ml.classification import LogisticRegression
        
        # 特征向量化
        assembler = VectorAssembler(
            inputCols=feature_cols,
            outputCol="features"
        )
        
        # 训练逻辑回归模型
        lr = LogisticRegression(
            featuresCol="features",
            labelCol="group_name",
            predictionCol="propensity",
            probabilityCol="propensity_score"
        )
        
        # 转换group_name为数值标签
        indexer = StringIndexer(inputCol="group_name", outputCol="label")
        
        # 组装管道
        from pyspark.ml import Pipeline
        pipeline = Pipeline(stages=[assembler, indexer, lr])
        
        # 训练模型（仅使用观测数据）
        model = pipeline.fit(obs_df)
        
        # 预测倾向得分
        scored_df = model.transform(obs_df)
        
        # 提取实验组的倾向得分
        from pyspark.sql.functions import expr
        
        final_df = scored_df.withColumn(
            "propensity_score",
            expr("propensity_score[1]")  # 取treatment概率
        )
        
        return final_df
    
    def merge_and_balance(self, exp_df, obs_df):
        """Ⅳ. 数据融合与平衡性检查"""
        
        # 合并数据
        unified_df = exp_df.unionByName(obs_df, allowMissingColumns=True)
        
        # 添加处理时间戳
        unified_df = unified_df.withColumn(
            "processed_at",
            lit(datetime.utcnow().isoformat())
        )
        
        # 写入Delta Lake（支持ACID事务）
        unified_df.write.format("delta") \
            .mode("overwrite") \
            .partitionBy("data_source", "experiment_id") \
            .option("mergeSchema", "true") \
            .save(self.config['unified_data_path'])
        
        # 生成平衡性报告
        balance_report = self._generate_balance_report(unified_df)
        
        return unified_df, balance_report
    
    def _generate_balance_report(self, df):
        """生成数据平衡性诊断报告"""
        
        report = {}
        
        for source in ["experiment", "observational"]:
            source_df = df.filter(col("data_source") == source)
            
            report[source] = {
                "total_users": source_df.count(),
                "control_users": source_df.filter(col("group_name") == "control").count(),
                "treatment_users": source_df.filter(col("group_name") == "treatment").count(),
                "propensity_score_balance": self._calculate_smd(
                    source_df.filter(col("group_name") == "control"),
                    source_df.filter(col("group_name") == "treatment")
                )
            }
        
        return report
    
    def _calculate_smd(self, control_df, treatment_df):
        """计算标准化均值差异（Standardized Mean Difference）"""
        
        # 倾向得分的SMD应小于0.1才认为平衡
        from pyspark.sql.functions import stddev, mean
        
        control_ps = control_df.select(mean("propensity_score").alias("mean"),
                                     stddev("propensity_score").alias("std"))
        treatment_ps = treatment_df.select(mean("propensity_score").alias("mean"),
                                         stddev("propensity_score").alias("std"))
        
        # 实现略：计算Cohen's d
        # 返回SMD值
        pass

# 使用示例
if __name__ == "__main__":
    spark = SparkSession.builder \
        .appName("UnifiedDataPipeline") \
        .config("spark.jars.packages", "io.delta:delta-core_2.12:2.4.0") \
        .getOrCreate()
    
    config = {
        'exp_db_url': 'jdbc:postgresql://exp-db.company.com:5432/experiments',
        'obs_data_path': 's3://data-lake/user-events/',
        'unified_data_path': 's3://data-lake/unified-experiment-data/'
    }
    
    pipeline = UnifiedDataPipeline(spark, config)
    
    # 运行数据融合
    exp_data = pipeline.ingest_treatment_data(
        experiment_id="rec_algo_v3", 
        start_date="2024-01-01", 
        end_date="2024-01-31"
    )
    
    obs_data = pipeline.ingest_observational_data(
        start_date="2024-01-01", 
        end_date="2024-01-31"
    )
    
    # 为观测数据计算倾向得分
    obs_data_scored = pipeline.calculate_propensity_scores(obs_data)
    
    # 合并并检查平衡性
    unified_df, report = pipeline.merge_and_balance(exp_data, obs_data_scored)
    
    print(f"数据融合完成，平衡性报告: {report}")

代码深度解析：

1. Schema统一策略：强制统一为user_id, data_source, experiment_id, group_name, event_time, variant_params, metadata结构，这是后续分析的基础。

2. 倾向得分计算：对观测数据，我们使用逻辑回归基于用户特征预测接受treatment的概率。关键假设是条件独立假设（CIA）：给定观测特征，treatment分配与潜在结果独立。

3. Delta Lake重要性：mergeSchema=true支持schema演进，partitionBy确保查询性能。实验数据和观测数据物理隔离，但逻辑统一。

4. 平衡性检查：SMD（标准化均值差异）是倾向得分匹配的核心诊断指标，必须小于0.1才认为两组可比较。

2.3 Mermaid：ETL数据流

第三章：因果推断方法体系与代码实现

3.1 双重差分法（Difference-in-Differences, DiD）

当实验数据存在时间序列，观测数据有面板结构时，DiD是首选方法。

# causal_inference/did_analysis.py
import pandas as pd
import numpy as np
from linearmodels import PanelOLS
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan

class DifferenceInDifferencesAnalyzer:
    """双重差分法分析器"""
    
    def __init__(self, data: pd.DataFrame):
        """
        Args:
            data: DataFrame必须包含列：
                - user_id: 用户唯一标识
                - period: 时间周期（0=实验前, 1=实验后）
                - treated: 是否接受处理（0=对照组, 1=实验组）
                - outcome: 结果变量（如GMV、CTR等）
                - covariates: 协变量列表（可选）
        """
        self.data = data.copy()
        self.results = None
        
    def check_parallel_trends(self, pre_periods: int = 4):
        """Ⅰ. 检验平行趋势假设（关键前提）"""
        
        # 需要实验前多期数据
        pre_data = self.data[self.data['period'] < 0].copy()
        pre_data['period_month'] = pre_data['period'].abs()
        
        # 绘制趋势图（实际实现需用matplotlib）
        # 检查实验前两组趋势是否平行
        
        # 统计检验：交互项不显著
        formula = 'outcome ~ period_month * treated'
        model = sm.OLS.from_formula(formula, data=pre_data).fit()
        
        # 获取交互项p值
        interaction_pvalue = model.pvalues['period_month:treated']
        
        parallel_check = {
            "assumption_valid": interaction_pvalue > 0.05,
            "interaction_pvalue": interaction_pvalue,
            "interpretation": "交互项不显著(p>0.05)则平行趋势成立" if interaction_pvalue > 0.05 else "平行趋势违背"
        }
        
        return parallel_check
    
    def estimate_causal_effect(self, covariates: list = None):
        """Ⅱ. 估计因果效应"""
        
        # 创建交互项
        self.data['treated_post'] = self.data['treated'] * self.data['period']
        
        # 基础DiD模型：Y = β0 + β1*period + β2*treated + β3*period*treated + ε
        formula = 'outcome ~ period + treated + treated_post'
        
        if covariates:
            covariate_str = ' + '.join(covariates)
            formula += f' + {covariate_str}'
        
        # 使用面板数据固定效应模型
        model = PanelOLS.from_formula(
            formula,
            data=self.data.set_index(['user_id', 'period']),
            other_effects=self.data['treated']  # 个体固定效应
        )
        
        results = model.fit(cov_type='robust')
        
        # 提取关键系数：交互项系数就是ATT（实验组的平均处理效应）
        att = results.params['treated_post']
        att_se = results.std_errors['treated_post']
        att_pvalue = results.pvalues['treated_post']
        
        self.results = {
            "model_type": "two_way_fixed_effects_did",
            "att": att,  # 平均处理效应
            "att_se": att_se,
            "att_pvalue": att_pvalue,
            "ci_95": [att - 1.96*att_se, att + 1.96*att_se],
            "significant": att_pvalue < 0.05,
            "full_results": results
        }
        
        return self.results
    
    def robustness_checks(self):
        """Ⅲ. 稳健性检验"""
        
        checks = {}
        
        # 检验1：安慰剂检验（Placebo Test）
        # 将处理时间提前，检验效应是否消失
        placebo_data = self.data.copy()
        placebo_data['period'] = placebo_data['period'].replace({1: 0, 0: -1})
        
        placebo_analyzer = DifferenceInDifferencesAnalyzer(placebo_data)
        placebo_analyzer.estimate_causal_effect()
        placebo_att = placebo_analyzer.results['att']
        
        checks['placebo_test'] = {
            "passed": abs(placebo_att) < 0.01,
            "placebo_att": placebo_att,
            "interpretation": "安慰剂效应应接近0，否则存在未观测混杂"
        }
        
        # 检验2：异方差性检验
        if self.results:
            model = self.results['full_results']
            # 提取残差
            residuals = model.resids
            fitted = model.fitted_values
            
            # Breusch-Pagan检验
            bp_stat, bp_pvalue, _, _ = het_breuschpagan(residuals, fitted)
            
            checks['heteroskedasticity'] = {
                "bp_statistic": bp_stat,
                "bp_pvalue": bp_pvalue,
                "passed": bp_pvalue > 0.05,
                "interpretation": "p>0.05说明无异方差问题"
            }
        
        # 检验3：控制个体时间趋势
        # 允许每个个体有自己的时间趋势
        robust_data = self.data.copy()
        robust_data['user_trend'] = robust_data.groupby('user_id').cumcount()
        
        formula_robust = 'outcome ~ period + treated + treated_post + user_trend'
        model_robust = PanelOLS.from_formula(
            formula_robust,
            data=robust_data.set_index(['user_id', 'period'])
        )
        
        robust_results = model_robust.fit(cov_type='robust')
        
        checks['individual_trends'] = {
            "att_robust": robust_results.params['treated_post'],
            "consistent": abs(robust_results.params['treated_post'] - self.results['att']) < 0.5*self.results['att_se'],
            "interpretation": "控制个体趋势后结果应保持一致"
        }
        
        return checks

# 实战调用示例
if __name__ == "__main__":
    # 模拟面板数据
    np.random.seed(42)
    
    # 1000个用户，每个用户4个周期（-3到0为实验前，1为实验后）
    users = np.repeat(range(1000), 5)
    periods = np.tile([-3, -2, -1, 0, 1], 1000)
    
    # 前500个用户为实验组（treated）
    treated = np.repeat([1]*500 + [0]*500, 5)
    
    # 生成结果变量：实验组在实验后有+0.5的提升
    outcome_base = np.random.normal(10, 2, 5000)
    outcome = outcome_base + treated * (periods == 1) * 0.5
    
    df = pd.DataFrame({
        'user_id': users,
        'period': periods,
        'treated': treated,
        'outcome': outcome
    })
    
    analyzer = DifferenceInDifferencesAnalyzer(df)
    
    # 检查平行趋势
    parallel_check = analyzer.check_parallel_trends()
    print(f"平行趋势检验: {parallel_check}")
    
    # 估计因果效应
    results = analyzer.estimate_causal_effect()
    print(f"ATT: {results['att']:.4f}, p-value: {results['att_pvalue']:.4f}")
    
    # 稳健性检验
    robustness = analyzer.robustness_checks()
    print(f"稳健性检验结果: {robustness}")

3.2 倾向得分匹配（Propensity Score Matching）

# causal_inference/psm_analysis.py
from sklearn.linear_model import LogisticRegression
from sklearn.neighbors import NearestNeighbors
import pandas as pd
import numpy as np
from scipy import stats

class PropensityScoreMatcher:
    """倾向得分匹配分析器"""
    
    def __init__(self, data: pd.DataFrame, treatment_col: str, outcome_col: str, covariates: list):
        self.data = data.copy()
        self.treatment_col = treatment_col
        self.outcome_col = outcome_col
        self.covariates = covariates
        self.matched_data = None
        
    def estimate_propensity_scores(self):
        """Ⅰ. 估计倾向得分"""
        
        X = self.data[self.covariates]
        y = self.data[self.treatment_col]
        
        # 使用逻辑回归
        # 注意：添加正则化防止过拟合
        model = LogisticRegression(
            penalty='l2',
            C=0.1,
            max_iter=1000,
            random_state=42
        )
        
        model.fit(X, y)
        
        # 预测倾向得分（接受处理的概率）
        propensity_scores = model.predict_proba(X)[:, 1]
        
        self.data['propensity_score'] = propensity_scores
        
        # 计算每个协变量的标准化偏差
        bias_reduction = self._calculate_covariate_balance()
        
        return {
            "model": model,
            "propensity_scores": propensity_scores,
            "pre_matching_balance": bias_reduction
        }
    
    def _calculate_covariate_balance(self):
        """计算匹配前的协变量不平衡"""
        
        treated = self.data[self.data[self.treatment_col] == 1]
        control = self.data[self.data[self.treatment_col] == 0]
        
        balance_stats = {}
        
        for covariate in self.covariates:
            treated_mean = treated[covariate].mean()
            control_mean = control[covariate].mean()
            
            # 标准化均值差异
            smd = abs(treated_mean - control_mean) / np.sqrt(
                (treated[covariate].var() + control[covariate].var()) / 2
            )
            
            balance_stats[covariate] = {
                "treated_mean": treated_mean,
                "control_mean": control_mean,
                "smd": smd,
                "imbalanced": smd > 0.1  # 大于0.1认为不平衡
            }
        
        return balance_stats
    
    def match_neighbors(self, caliper: float = 0.05, ratio: int = 1):
        """Ⅱ. 执行最近邻匹配"""
        
        # 分离处理组和对照组
        treated = self.data[self.data[self.treatment_col] == 1].copy()
        control = self.data[self.data[self.treatment_col] == 0].copy()
        
        # 使用caliper限制匹配距离
        # ratio=1表示1:1匹配，ratio=2表示1:2匹配
        
        treated_ps = treated['propensity_score'].values.reshape(-1, 1)
        control_ps = control['propensity_score'].values.reshape(-1, 1)
        
        # 使用KD树加速最近邻搜索
        nn = NearestNeighbors(
            n_neighbors=ratio,
            radius=caliper,
            algorithm='ball_tree'
        )
        nn.fit(control_ps)
        
        # 为每个处理组样本找到匹配
        distances, indices = nn.radius_neighbors(treated_ps, radius=caliper)
        
        matched_treated = []
        matched_control = []
        
        for i, (dist, idx) in enumerate(zip(distances, indices)):
            if len(idx) >= 1:  # 至少找到一个匹配
                matched_treated.append(treated.iloc[i]['user_id'])
                matched_control.extend(control.iloc[idx]['user_id'].tolist())
        
        # 提取匹配后的数据
        matched_users = matched_treated + matched_control
        self.matched_data = self.data[self.data['user_id'].isin(matched_users)].copy()
        
        # 添加匹配权重（1:1匹配时权重相等）
        self.matched_data['matching_weight'] = np.where(
            self.matched_data[self.treatment_col] == 1,
            1.0,
            1.0 / ratio
        )
        
        # 评估匹配后的平衡性
        post_balance = self._calculate_post_matching_balance()
        
        return {
            "matched_sample_size": len(matched_users),
            "treated_matched": len(matched_treated),
            "control_matched": len(matched_control),
            "post_matching_balance": post_balance,
            "balance_improvement": self._compare_balance()
        }
    
    def _calculate_post_matching_balance(self):
        """计算匹配后的协变量平衡性"""
        
        if self.matched_data is None:
            raise ValueError("必须先执行match_neighbors()")
        
        treated = self.matched_data[self.matched_data[self.treatment_col] == 1]
        control = self.matched_data[self.matched_data[self.treatment_col] == 0]
        
        post_balance = {}
        
        for covariate in self.covariates:
            treated_mean = treated[covariate].mean()
            control_mean = control[covariate].mean()
            
            smd = abs(treated_mean - control_mean) / np.sqrt(
                (treated[covariate].var() + control[covariate].var()) / 2
            )
            
            post_balance[covariate] = {
                "smd": smd,
                "balanced": smd < 0.1
            }
        
        return post_balance
    
    def estimate_ate(self):
        """Ⅲ. 估计平均处理效应（ATE）"""
        
        if self.matched_data is None:
            raise ValueError("必须先执行匹配")
        
        # 使用匹配权重计算加权均值
        treated_outcome = self.matched_data[
            self.matched_data[self.treatment_col] == 1
        ][self.outcome_col].mean()
        
        control_outcome = self.matched_data[
            self.matched_data[self.treatment_col] == 0
        ][self.outcome_col].mean()
        
        # 简单ATT
        att = treated_outcome - control_outcome
        
        # 计算标准误（考虑匹配权重）
        # 实现略：使用自助法或解析公式
        
        return {
            "att": att,
            "treated_mean": treated_outcome,
            "control_mean": control_outcome,
            "sample_size": len(self.matched_data)
        }

# 实战调用
if __name__ == "__main__":
    # 生成模拟数据
    np.random.seed(42)
    n = 10000
    
    # 协变量
    age = np.random.normal(30, 10, n)
    income = np.random.normal(5000, 2000, n)
    account_age = np.random.exponential(2, n)
    
    # 倾向得分（older, richer users more likely to receive treatment）
    ps = 1 / (1 + np.exp(-(0.1*age + 0.0002*income + 0.5*account_age - 8)))
    treatment = np.random.binomial(1, ps, n)
    
    # 结果变量（真实效应为+2）
    outcome = 10 + 0.2*age + 0.0001*income + 2*treatment + np.random.normal(0, 1, n)
    
    df = pd.DataFrame({
        'user_id': range(n),
        'treated': treatment,
        'outcome': outcome,
        'age': age,
        'income': income,
        'account_age': account_age
    })
    
    matcher = PropensityScoreMatcher(
        data=df,
        treatment_col='treated',
        outcome_col='outcome',
        covariates=['age', 'income', 'account_age']
    )
    
    # 估计倾向得分
    ps_result = matcher.estimate_propensity_scores()
    print("匹配前平衡性:", ps_result['pre_matching_balance'])
    
    # 执行1:1匹配，caliper=0.1
    match_result = matcher.match_neighbors(caliper=0.1, ratio=1)
    print(f"匹配样本数: {match_result['matched_sample_size']}")
    print("匹配后平衡性:", match_result['post_matching_balance'])
    
    # 估计ATE
    ate_result = matcher.estimate_ate()
    print(f"ATT估计: {ate_result['att']:.4f}")

3.3 工具变量法（Instrumental Variables）

当存在未观测的混杂因素时，工具变量法提供一致的估计。

# causal_inference/iv_analysis.py
import pandas as pd
import numpy as np
from statsmodels.sandbox.regression.gmm import IV2SLS

class InstrumentalVariableAnalyzer:
    """工具变量法分析器"""
    
    def __init__(self, data: pd.DataFrame, outcome: str, treatment: str, instrument: str, covariates: list):
        self.data = data.copy()
        self.outcome = outcome
        self.treatment = treatment
        self.instrument = instrument
        self.covariates = covariates
        
    def validate_instrument(self):
        """Ⅰ. 验证工具变量的有效性"""
        
        # 工具变量需满足三个条件：
        # 1. 相关性：Z与X相关
        # 2. 外生性：Z与误差项不相关
        # 3. 排他性：Z只通过X影响Y
        
        validation = {}
        
        # 检验1：第一阶段F统计量（相关性）
        X = self.data[[self.instrument] + self.covariates]
        y = self.data[self.treatment]
        
        first_stage = IV2SLS(y, X, None).fit()
        f_stat = first_stage.fvalue
        
        validation['relevance'] = {
            "f_statistic": f_stat,
            "passed": f_stat > 10,  # F>10表示强工具变量
            "interpretation": "F统计量大于10说明工具变量与处理变量强相关"
        }
        
        # 检验2：过度识别检验（如果工具变量多于内生变量）
        # 实现略：Sargan检验或Hansen J检验
        
        return validation
    
    def estimate_late(self):
        """Ⅱ. 估计局部平均处理效应（LATE）"""
        
        # 两阶段最小二乘（2SLS）
        # 阶段1：D = π0 + π1*Z + π2*C + ν
        # 阶段2：Y = β0 + β1*D̂ + β2*C + ε
        
        # 准备数据
        exog = self.data[[self.instrument] + self.covariates]
        endog = self.data[self.treatment]
        dependent = self.data[self.outcome]
        
        # 2SLS估计
        model = IV2SLS(dependent, exog, endog)
        results = model.fit()
        
        # LATE是阶段2中treatment的系数
        late = results.params[self.treatment]
        late_se = results.bse[self.treatment]
        
        return {
            "late": late,
            "late_se": late_se,
            "ci_95": [late - 1.96*late_se, late + 1.96*late_se],
            "first_stage_f_stat": results.first_stage.fvalue,
            "full_results": results
        }

# 实战案例：推荐系统实验
# Z = 功能入口是否可见（随机分配）
# X = 是否实际使用新推荐
# Y = 用户GMV
# C = 用户特征协变量

3.4 Mermaid：因果推断方法选择决策树

第四章：实战案例 - 电商平台推荐算法升级的深度分析（2000+字）

4.1 业务背景与挑战

某跨境电商平台（月活3.2亿）在2024年Q1启动推荐算法升级项目，从GBDT排序模型切换为深度兴趣网络（DIN）。项目面临的核心挑战远超常规A/B测试：

数据获取难题矩阵：

数据可用性分析：

4.2 分析框架设计：三阶段融合策略

阶段一：RCT小样本分析（内部有效性优先）

目标：获取无偏的因果效应估计，作为基准真值

技术实现：

# case_study/rct_baseline.py
import pandas as pd
import numpy as np
from scipy import stats
from analysis.bayesian_analysis import BayesianExperimentAnalyzer

class RCTBaseline:
    """RCT小样本分析基线"""
    
    def __init__(self, spark, experiment_id="rec_din_001"):
        self.spark = spark
        self.experiment_id = experiment_id
        
    def extract_rct_data(self):
        """从Delta Lake提取实验数据"""
        
        rct_df = self.spark.sql(f"""
            SELECT 
                user_id,
                group_name,
                assignment_time,
                -- 结果变量
                SUM(CASE WHEN event_type = 'purchase' THEN 1 ELSE 0 END) as conversions,
                SUM(CASE WHEN event_type = 'purchase' THEN revenue ELSE 0 END) as gmv,
                COUNT(DISTINCT session_id) as sessions,
                AVG(CASE WHEN event_type = 'click' THEN 1 ELSE 0 END) as ctr
                
            FROM delta.`s3://data-lake/unified-experiment-data/`
            WHERE experiment_id = '{self.experiment_id}'
                AND data_source = 'experiment'
                AND assignment_time BETWEEN '2024-01-15' AND '2024-02-15'
            GROUP BY user_id, group_name, assignment_time
        """).toPandas()
        
        return rct_df
    
    def analyze_bayesian(self, rct_data):
        """贝叶斯分析获取后验分布"""
        
        # 按组聚合
        grouped = rct_data.groupby('group_name').agg({
            'conversions': 'sum',
            'gmv': 'sum',
            'user_id': 'count',
            'ctr': 'mean'
        }).rename(columns={'user_id': 'users'})
        
        # CTR分析
        analyzer = BayesianExperimentAnalyzer()
        ctr_result = analyzer.analyze_conversion_rate(
            control_conversions=int(grouped.loc['control', 'conversions']),
            control_users=int(grouped.loc['control', 'users']),
            treatment_conversions=int(grouped.loc['treatment', 'conversions']),
            treatment_users=int(grouped.loc['treatment', 'users'])
        )
        
        # GMV分析（连续变量）
        # 使用PyMC建模
        gmv_result = self._bayesian_gmv_analysis(
            control_gmv=rct_data[rct_data['group_name'] == 'control']['gmv'].values,
            treatment_gmv=rct_data[rct_data['group_name'] == 'treatment']['gmv'].values
        )
        
        return {
            "experiment_id": self.experiment_id,
            "sample_size": len(rct_data),
            "ctr_effect": ctr_result,
            "gmv_effect": gmv_result,
            "internal_validity": "HIGH"
        }
    
    def _bayesian_gmv_analysis(self, control_gmv, treatment_gmv, 
                               num_samples=5000):
        """连续变量的贝叶斯分析"""
        
        import pymc as pm
        
        with pm.Model() as model:
            # 先验：正态分布
            mu_control = pm.Normal('mu_control', mu=0, sigma=100)
            mu_treatment = pm.Normal('mu_treatment', mu=0, sigma=100)
            
            sigma_control = pm.HalfNormal('sigma_control', sigma=50)
            sigma_treatment = pm.HalfNormal('sigma_treatment', sigma=50)
            
            # 似然
            pm.Normal('control', mu=mu_control, sigma=sigma_control, 
                     observed=control_gmv)
            pm.Normal('treatment', mu=mu_treatment, sigma=sigma_treatment, 
                     observed=treatment_gmv)
            
            # 效应量
            effect = pm.Deterministic('effect', mu_treatment - mu_control)
            
            # 采样
            trace = pm.sample(num_samples, chains=4, cores=2, return_inferencedata=True)
        
        # 后验统计
        effect_mean = float(trace.posterior['effect'].mean())
        effect_hdi = pm.hdi(trace, var_names=['effect'], hdi_prob=0.95)
        
        return {
            "mean_effect": effect_mean,
            "hdi_95_lower": float(effect_hdi['effect'][0]),
            "hdi_95_upper": float(effect_hdi['effect'][1]),
            "probability_positive": float((trace.posterior['effect'] > 0).mean())
        }

# 执行基线分析
if __name__ == "__main__":
    from pyspark.sql import SparkSession
    
    spark = SparkSession.builder.appName("RCTBaseline").getOrCreate()
    
    analyzer = RCTBaseline(spark)
    rct_data = analyzer.extract_rct_data()
    baseline_results = analyzer.analyze_bayesian(rct_data)
    
    print(f"RCT基线结果: {baseline_results}")

RCT基线结果解读（30天数据）：

关键洞察：RCT显示DIN算法在核心指标上有显著正向提升，但由于仅20%流量参与实验，样本量不足以分析长尾效应和小众品类表现。

阶段二：观测数据偏差校正（扩大代表性）

目标：利用80%观测数据，校正选择偏差后验证RCT结论

技术挑战深度解析：

观测数据中用户主动使用"新推荐模式"存在严重自选择偏差。使用新功能的用户特征如下：

技术实现：

# case_study/observational_correction.py
from pyspark.sql import SparkSession
from pyspark.ml.feature import VectorAssembler, StandardScaler
from pyspark.ml.classification import LogisticRegression
from causal_inference.psm_analysis import PropensityScoreMatcher

class ObservationalBiasCorrection:
    """观测数据偏差校正"""
    
    def __init__(self, spark: SparkSession):
        self.spark = spark
        
    def extract_observational_cohort(self):
        """提取观测对照组与实验组"""
        
        # 观测"实验组"：主动使用新推荐功能的用户
        obs_treatment = self.spark.sql("""
            SELECT 
                user_id,
                'observational' as data_source,
                first_use_time as assignment_time,
                device_type,
                account_age_days,
                historical_gmv,
                weekly_active_days,
                exploration_score
            FROM user_features.features_daily
            WHERE feature_new_rec_active = true
                AND first_use_time BETWEEN '2024-01-15' AND '2024-02-15'
        """)
        
        # 观测"对照组"：未使用新功能的相似用户
        # 使用倾向得分匹配前的候选池
        obs_control_pool = self.spark.sql("""
            SELECT 
                user_id,
                'observational' as data_source,
                '2024-01-15' as assignment_time,  -- 对齐时间
                device_type,
                account_age_days,
                historical_gmv,
                weekly_active_days,
                exploration_score
            FROM user_features.features_daily
            WHERE feature_new_rec_active = false
                AND last_active_time >= '2024-01-15'
        """)
        
        return obs_treatment, obs_control_pool
    
    def build_propensity_model(self, obs_treatment, obs_control_pool):
        """构建倾向得分模型"""
        
        # 合并数据并标记处理状态
        obs_treatment = obs_treatment.withColumn("treated", F.lit(1))
        obs_control_pool = obs_control_pool.withColumn("treated", F.lit(0))
        
        # 采样控制组（避免数据不平衡）
        obs_control = obs_control_pool.sample(
            fraction=0.3,
            seed=42
        )
        
        # 合并
        obs_data = obs_treatment.union(obs_control)
        
        # 特征工程
        feature_cols = ['account_age_days', 'historical_gmv', 
                       'weekly_active_days', 'exploration_score']
        
        # 类别变量编码
        obs_data = obs_data.withColumn(
            "device_mobile", 
            (F.col("device_type") == "mobile").cast("int")
        )
        feature_cols.append("device_mobile")
        
        # 向量化
        assembler = VectorAssembler(
            inputCols=feature_cols,
            outputCol="features"
        )
        
        obs_vectorized = assembler.transform(obs_data)
        
        # 标准化
        scaler = StandardScaler(
            inputCol="features",
            outputCol="scaled_features"
        )
        scaler_model = scaler.fit(obs_vectorized)
        obs_scaled = scaler_model.transform(obs_vectorized)
        
        # 逻辑回归
        lr = LogisticRegression(
            featuresCol="scaled_features",
            labelCol="treated",
            probabilityCol="propensity_score"
        )
        
        propensity_model = lr.fit(obs_scaled)
        
        # 预测倾向得分
        obs_with_ps = propensity_model.transform(obs_scaled)
        
        # 模型评估
        summary = propensity_model.summary
        
        # 提取平衡性指标
        balance_check = self._evaluate_balance(obs_with_ps, feature_cols)
        
        return obs_with_ps, balance_check
    
    def _evaluate_balance(self, df, feature_cols):
        """评估匹配前后协变量平衡性"""
        
        # 转换为Pandas进行统计分析（Spark DataFrame操作复杂度高）
        pd_df = df.select(
            "treated", "propensity_score", *feature_cols
        ).toPandas()
        
        matcher = PropensityScoreMatcher(
            data=pd_df,
            treatment_col="treated",
            outcome_col="propensity_score",  # 先平衡性检查
            covariates=feature_cols
        )
        
        # 匹配前SMD
        pre_balance = matcher._calculate_covariate_balance()
        
        # 执行匹配
        match_result = matcher.match_neighbors(caliper=0.2, ratio=1)
        
        # 检查平衡性改善
        improvement = {}
        for cov in feature_cols:
            pre_smd = pre_balance[cov]['smd']
            post_smd = match_result['post_matching_balance'][cov]['smd']
            
            improvement[cov] = {
                "pre_smd": pre_smd,
                "post_smd": post_smd,
                "improvement_pct": ((pre_smd - post_smd) / pre_smd) * 100,
                "balanced": post_smd < 0.1
            }
        
        return {
            "pre_matching_balance": pre_balance,
            "post_matching_balance": match_result['post_matching_balance'],
            "improvement": improvement,
            "match_rate": match_result['treated_matched'] / pd_df['treated'].sum()
        }

# 执行观测数据偏差校正
if __name__ == "__main__":
    spark = SparkSession.builder.appName("ObservationalCorrection").getOrCreate()
    
    corrector = ObservationalBiasCorrection(spark)
    obs_treatment, obs_control = corrector.extract_observational_cohort()
    obs_data, balance_report = corrector.build_propensity_model(obs_treatment, obs_control)
    
    print(f"倾向得分模型完成，平衡性改善: {balance_report['improvement']}")

偏差校正结果深度解读：

匹配后，观测数据的平衡性显著改善：

关键发现：虽然并非所有协变量完全达标（GMV的SMD=0.12），但已极大缓解偏差。匹配后样本量从2.56亿缩减到420万对匹配样本，统计功效与RCT相当但代表性更广。

观测数据校正后的效应估计：

# 对匹配后的观测数据计算效应
matched_obs_data = obs_data.filter("propensity_score IS NOT NULL")

obs_results = analyzer.analyze_bayesian(matched_obs_data.toPandas())

print(f"观测数据校正后CTR效应: {obs_results['ctr_effect']['mean_effect']:.4f}")

结果：观测数据校正后CTR提升+4.89%，与RCT的+5.22%接近，偏差仅6.3%，验证了RCT结论的外部有效性。

阶段三：数据融合与因果增强（1+1>2）

目标：结合实验与观测数据，提升估计精度，分析长尾效应

核心挑战：如何加权融合两种数据源？

技术方案：逆方差加权融合

# case_study/data_fusion.py
import numpy as np
from typing import Dict

class DataFusionEngine:
    """实验与观测数据融合引擎"""
    
    def __init__(self, rct_results: Dict, obs_results: Dict):
        self.rct = rct_results
        self.obs = obs_results
        
    def inverse_variance_weighting(self):
        """逆方差加权融合"""
        
        # 提取点估计和标准误
        rct_est = self.rct['mean_effect']
        rct_se = self.rct['se']
        
        obs_est = self.obs['mean_effect']
        obs_se = self.obs['se']
        
        # 计算权重（精度权重）
        # 精度 = 1/方差
        rct_precision = 1 / (rct_se ** 2)
        obs_precision = 1 / (obs_se ** 2)
        
        total_precision = rct_precision + obs_precision
        
        w_rct = rct_precision / total_precision
        w_obs = obs_precision / total_precision
        
        # 融合估计
        fused_estimate = w_rct * rct_est + w_obs * obs_est
        
        # 融合标准误（方差加权平均）
        fused_se = np.sqrt(1 / total_precision)
        
        # 95%置信区间
        ci_lower = fused_estimate - 1.96 * fused_se
        ci_upper = fused_estimate + 1.96 * fused_se
        
        return {
            "fused_estimate": fused_estimate,
            "fused_se": fused_se,
            "ci_95": [ci_lower, ci_upper],
            "weights": {
                "rct_weight": w_rct,
                "obs_weight": w_obs,
                "rct_precision": rct_precision,
                "obs_precision": obs_precision
            },
            "efficiency_gain": (1 / fused_se ** 2) / (1 / rct_se ** 2) - 1
        }
    
    def bayesian_synthesis(self, prior_strength: float = 0.5):
        """贝叶斯合成：RCT作为观测数据的先验"""
        
        # RCT结果作为强先验
        prior_mean = self.rct['mean_effect']
        prior_sd = self.rct['se'] / prior_strength  # 先验强度可调
        
        # 观测数据作为似然
        likelihood_mean = self.obs['mean_effect']
        likelihood_sd = self.obs['se']
        
        # 后验分布（共轭先验：正态-正态）
        posterior_precision = (1 / prior_sd ** 2) + (1 / likelihood_sd ** 2)
        posterior_sd = np.sqrt(1 / posterior_precision)
        posterior_mean = (
            (prior_mean / prior_sd ** 2) + 
            (likelihood_mean / likelihood_sd ** 2)
        ) / posterior_precision
        
        return {
            "posterior_mean": posterior_mean,
            "posterior_sd": posterior_sd,
            "ci_95": [
                posterior_mean - 1.96 * posterior_sd,
                posterior_mean + 1.96 * posterior_sd
            ],
            "prior_strength": prior_strength
        }
    
    def analyze_heterogeneous_effects(self):
        """异质性效应分析：哪些用户受益更多？"""
        
        # 使用观测数据的大样本量进行亚组分析
        # 1. 新用户 vs 老用户
        # 2. 不同品类偏好
        # 3. 不同设备类型
        
        # 实现：在观测数据上按特征分层，每层计算效应
        # 再用RCT数据验证关键亚组
        
        # 返回CATE（条件平均处理效应）分布
        pass

# 融合分析执行
if __name__ == "__main__":
    # 假设已有RCT和观测分析结果
    rct_results = {
        "mean_effect": 0.0522,  # 5.22%
        "se": 0.0085
    }
    
    obs_results = {
        "mean_effect": 0.0489,  # 4.89%
        "se": 0.0092
    }
    
    fusion = DataFusionEngine(rct_results, obs_results)
    
    # 方法1：逆方差加权
    ivw_result = fusion.inverse_variance_weighting()
    print(f"逆方差加权融合效应: {ivw_result['fused_estimate']:.4f} (效率提升: {ivw_result['efficiency_gain']:.2%})")
    
    # 方法2：贝叶斯合成
    bayes_result = fusion.bayesian_synthesis(prior_strength=1.0)
    print(f"贝叶斯合成后验均值: {bayes_result['posterior_mean']:.4f}")

融合结果深度解读：

核心洞察：

1. 效率提升88%：融合后标准误从0.85%降低到0.62%，等价于样本量扩大3.5倍，这是观测数据带来的最大价值。

2. "数据增效"原理 ：观测数据虽存在偏差，但经过倾向得分匹配后，其残余偏差较小，而方差降低显著，因此加权融合能提升总体精度。

3. 权重分配的智慧：RCT权重54.2%，观测权重45.8%，说明系统自动平衡了内部有效性与样本量的关系，无需人工设定。

阶段四：长期效应与业务落地

目标：分析3个月长期效应，指导全量上线决策

技术挑战：实验仅运行30天，如何预测90天效应？

解决方案：合成控制法 + 结构模型

# case_study/long_term_projection.py
import pandas as pd
import numpy as np
from causal_inference.synthetic_control import SyntheticControl
from sklearn.ensemble import RandomForestRegressor

class LongTermEffectProjector:
    """长期效应预测器"""
    
    def __init__(self, short_term_data: pd.DataFrame, observational_long_term: pd.DataFrame):
        """
        Args:
            short_term_data: 30天实验数据
            observational_long_term: 90天观测数据（包含处理组与对照组）
        """
        self.short_term = short_term_data
        self.long_term_obs = observational_long_term
        
    def fit_synthetic_control(self):
        """Ⅰ. 构建合成控制（基于观测对照组）"""
        
        # 使用观测数据中未接受处理的用户
        control_pool = self.long_term_obs[
            self.long_term_obs['treated'] == 0
        ]
        
        # 选择预测变量：前30天的行为特征
        predictors = [
            'daily_gmv_mean', 'daily_sessions_mean', 'category_diversity',
            'search_frequency', 'cart_add_rate'
        ]
        
        # 实验组（短期）作为目标
        treatment_unit = self.short_term[
            self.short_term['group_name'] == 'treatment'
        ].groupby('day')[predictors].mean()
        
        # 对照组作为捐赠池
        control_units = control_pool.groupby(['user_id', 'day'])[predictors].mean().unstack(level=0)
        
        sc = SyntheticControl()
        sc.fit(
            treatment_unit=treatment_unit,
            control_units=control_units,
            optimization_period=30,
            predictor_periods=range(0, 30)
        )
        
        # 预测90天反事实（假设未接受处理）
        counterfactual_90d = sc.predict(periods=90)
        
        return counterfactual_90d
    
    def estimate_long_term_lift(self):
        """Ⅱ. 估计长期提升效应"""
        
        # 获取合成控制预测
        counterfactual = self.fit_synthetic_control()
        
        # 实际观测的90天实验组表现
        actual_treatment = self.long_term_obs[
            self.long_term_obs['treated'] == 1
        ].groupby('day')['gmv'].mean()
        
        # 计算每日提升
        daily_lift = (actual_treatment - counterfactual) / counterfactual
        
        # 关键洞察：效应随时间衰减
        # 前30天：+5.08%（与RCT一致）
        # 31-60天：+4.12%（习惯效应递减）
        # 61-90天：+3.89%（稳定期）
        
        return {
            "month1_lift": daily_lift[0:30].mean(),
            "month2_lift": daily_lift[30:60].mean(),
            "month3_lift": daily_lift[60:90].mean(),
            "stability": "效应衰减15%后稳定",
            "annual_projected_lift": daily_lift.mean() * 0.85  # 保守估计
        }
    
    def build_causal_forest(self):
        """Ⅲ. 异质性长期效应：因果森林"""
        
        from econml.grf import CausalForest
        
        # 准备数据
        features = self.long_term_obs[
            ['account_age', 'historical_gmv', 'category_preference']
        ]
        T = self.long_term_obs['treated']
        Y = self.long_term_obs['gmv_90d_sum']
        
        # 训练因果森林
        cf = CausalForest(
            n_estimators=1000,
            min_samples_leaf=500,
            max_depth=10,
            random_state=42
        )
        
        cf.fit(Y, T, X=features)
        
        # 预测CATE
        cate_pred = cf.effect(features)
        
        # 识别受益最大的用户群
        high_value_segments = self._identify_segments(cate_pred, features)
        
        return high_value_segments

# 长期效应执行
projector = LongTermEffectProjector(
    short_term_data=rct_data,
    observational_long_term=obs_90d_data
)

long_term_results = projector.estimate_long_term_lift()
print(f"90天长期效应预测: {long_term_results}")

# 输出：90天效应衰减至+3.89%，但仍显著为正

4.3 业务决策与ROI计算

最终融合分析结论：

ROI计算与上线决策：

年度增量GMV计算：
- 短期效应：¥47.8 - ¥45.2 = ¥2.6/人/月（5.75%提升）
- 长期衰减调整后：¥2.6 * 0.85 = ¥2.21/人/月
- 全量月活：3.2亿
- 影响用户比例：85%（15%低活用户不受影响）
- 年度增量GMV = ¥2.21 * 320M * 0.85 * 12 = **¥72.3亿**（约$10.3B）

成本分析：
- 算法研发：¥8M
- 工程部署：¥15M
- 算力增加：¥3M/月 = ¥36M/年
- 总成本：¥59M

ROI = (¥7.23B - ¥59M) / ¥59M = **122.6倍**

决策：✅ **立即全量上线**

4.4 业务影响复盘（上线后6个月真实数据）

预测 vs 真实对比：

关键成功要素总结：

第五章：生产级部署与MLOps集成

5.1 Airflow DAG生产化

# dags/fusion_analysis_dag.py
from airflow import DAG
from airflow.operators.python import PythonOperator
from airflow.providers.postgres.operators.postgres import PostgresOperator
from datetime import datetime, timedelta
import json

default_args = {
    'owner': 'causal_inference_team',
    'depends_on_past': True,
    'email_on_failure': True,
    'email_on_retry': False,
    'retries': 2,
    'retry_delay': timedelta(minutes=15),
    'execution_timeout': timedelta(hours=3)
}

with DAG(
    'observational_fusion_pipeline',
    default_args=default_args,
    description='实验与观测数据融合分析流水线',
    schedule_interval='0 3 * * *',  # 每天凌晨3点
    start_date=datetime(2024, 1, 1),
    catchup=False,
    tags=['causal_inference', 'experiment', 'production'],
    max_active_runs=1
) as dag:
    
    # Ⅰ. 数据摄取与融合
    ingest_unified_data = PythonOperator(
        task_id='ingest_and_merge_data',
        python_callable=run_unified_pipeline,
        op_kwargs={
            'experiment_id': 'rec_din_001',
            'date_range': '{{ ds }}'
        }
    )
    
    # Ⅱ. 倾向得分计算与匹配（PySpark任务）
    calculate_propensity = PythonOperator(
        task_id='calculate_propensity_scores',
        python_callable=run_propensity_matching,
        op_kwargs={
            'output_table': 'matched_data_{{ ds_nodash }}'
        }
    )
    
    # Ⅲ. 因果推断分析（三方法并行）
    with TaskGroup('causal_estimation') as causal_tg:
        
        rct_analysis = PythonOperator(
            task_id='rct_baseline',
            python_callable=run_rct_analysis,
            op_kwargs={'method': 'bayesian'}
        )
        
        did_analysis = PythonOperator(
            task_id='did_analysis',
            python_callable=run_did_analysis,
            op_kwargs={
                'panel_data_table': 'panel_data_{{ ds_nodash }}'
            }
        )
        
        psm_analysis = PythonOperator(
            task_id='psm_analysis',
            python_callable=run_psm_analysis,
            op_kwargs={
                'matched_data_table': 'matched_data_{{ ds_nodash }}'
            }
        )
        
        [rct_analysis, did_analysis, psm_analysis]
    
    # Ⅳ. 数据融合与效应合成
    fusion_task = PythonOperator(
        task_id='fuse_and_synthesize',
        python_callable=run_fusion_synthesis,
        provide_context=True
    )
    
    # Ⅴ. 质量检查与稳健性检验
    quality_check = PythonOperator(
        task_id='robustness_checks',
        python_callable=run_robustness_checks,
        op_kwargs={
            'checks': ['placebo_test', 'heteroskedasticity', 'sensitivity_analysis']
        }
    )
    
    # Ⅵ. 报告生成与发布
    generate_report = PythonOperator(
        task_id='generate_decision_report',
        python_callable=generate_html_report,
        op_kwargs={
            'template': 'causal_inference_report.html',
            'output_bucket': 's3://experiment-reports/'
        }
    )
    
    # Ⅶ. 元数据记录
    log_metadata = PostgresOperator(
        task_id='log_analysis_metadata',
        postgres_conn_id='metadata_db',
        sql="""
            INSERT INTO causal_inference_metadata 
            (execution_date, experiment_id, methods_used, fused_effect, ci_lower, ci_upper, decision_recommendation)
            VALUES (
                '{{ ds }}',
                '{{ task_instance.xcom_pull(task_ids="fuse_and_synthesize", key="experiment_id") }}',
                '{{ task_instance.xcom_pull(task_ids="fuse_and_synthesize", key="methods") }}',
                {{ task_instance.xcom_pull(task_ids="fuse_and_synthesize", key="fused_effect") }},
                {{ task_instance.xcom_pull(task_ids="fuse_and_synthesize", key="ci_lower") }},
                {{ task_instance.xcom_pull(task_ids="fuse_and_synthesize", key="ci_upper") }},
                '{{ task_instance.xcom_pull(task_ids="robustness_checks", key="recommendation") }}'
            );
        """
    )
    
    # 依赖关系
    ingest_unified_data >> calculate_propensity >> causal_tg >> fusion_task
    fusion_task >> quality_check >> generate_report >> log_metadata

5.2 Mermaid：生产级数据流

5.3 数据质量监控与异常检测

# monitoring/data_quality_monitor.py
from great_expectations import DataContext
from great_expectations.checkpoint import Checkpoint
import pandas as pd
from datetime import datetime

class CausalDataQualityMonitor:
    """因果分析数据质量监控"""
    
    def __init__(self, ge_context_root: str):
        self.context = DataContext(ge_context_root)
        
    def validate_unified_data(self, dataset_name: str):
        """验证融合数据质量"""
        
        # Ⅰ. 基础完整性
        expectation_suite = self.context.create_expectation_suite(
            f"{dataset_name}_fusion_validation",
            overwrite_existing=True
        )
        
        # 检查1：倾向得分的合理性
        expectation_suite.add_expectation(
            {
                "expectation_type": "expect_column_values_to_be_between",
                "kwargs": {
                    "column": "propensity_score",
                    "min_value": 0.0,
                    "max_value": 1.0
                }
            }
        )
        
        # 检查2：匹配后的SMD必须小于0.1
        expectation_suite.add_expectation(
            {
                "expectation_type": "expect_column_values_to_be_between",
                "kwargs": {
                    "column": "smd_after_matching",
                    "min_value": 0.0,
                    "max_value": 0.1
                }
            }
        )
        
        # 检查3：实验组与观测组数据比例合理性
        expectation_suite.add_expectation(
            {
                "expectation_type": "expect_column_proportion_of_unique_values_to_be_between",
                "kwargs": {
                    "column": "data_source",
                    "min_proportion": 0.45,
                    "max_proportion": 0.55
                }
            }
        )
        
        # Ⅱ. 业务逻辑验证
        # 检查4：结果变量不存在异常值（使用IQR规则）
        expectation_suite.add_expectation(
            {
                "expectation_type": "expect_column_values_to_not_be_null",
                "kwargs": {"column": "outcome_gmv"}
            }
        )
        
        # 创建checkpoint并运行
        checkpoint_config = {
            "name": f"{dataset_name}_checkpoint",
            "config_version": 1.0,
            "class_name": "SimpleCheckpoint",
            "validations": [
                {
                    "batch_request": {
                        "datasource_name": "fusion_data",
                        "data_connector_name": "default_inferred_data_connector_name",
                        "data_asset_name": dataset_name,
                    },
                    "expectation_suite_name": f"{dataset_name}_fusion_validation",
                }
            ],
        }
        
        checkpoint = Checkpoint(**checkpoint_config, data_context=self.context)
        results = checkpoint.run()
        
        return results.success

# 在Airflow中集成
def run_data_quality_checks(**context):
    monitor = CausalDataQualityMonitor("/opt/great_expectations")
    
    execution_date = context['execution_date'].strftime("%Y%m%d")
    dataset_name = f"unified_data_{execution_date}"
    
    success = monitor.validate_unified_data(dataset_name)
    
    if not success:
        raise ValueError("数据质量检查失败，请查看Great Expectations报告")
    
    return True

5.4 模型版本控制与实验追踪

# utils/mlflow_integration.py
import mlflow
import mlflow.sklearn
from datetime import datetime

def track_causal_model(experiment_name: str, model_params: dict, metrics: dict, artifacts: dict):
    """使用MLflow追踪因果推断模型"""
    
    mlflow.set_experiment(experiment_name)
    
    with mlflow.start_run(run_name=f"causal_analysis_{datetime.utcnow().isoformat()}"):
        
        # Ⅰ. 记录参数
        mlflow.log_params(model_params)
        
        # Ⅱ. 记录指标
        mlflow.log_metrics(metrics)
        
        # Ⅲ. 记录模型
        if 'propensity_model' in artifacts:
            mlflow.sklearn.log_model(artifacts['propensity_model'], "propensity_model")
        
        # Ⅳ. 记录融合效应估计
        mlflow.log_dict(artifacts['fused_results'], "fused_results.json")
        
        # Ⅴ. 记录稳健性检验结果
        mlflow.log_dict(artifacts['robustness_checks'], "robustness.json")
        
        # Ⅵ. 记录报告
        if 'report_html' in artifacts:
            mlflow.log_artifact(artifacts['report_html'], "reports")
        
        # 注册模型版本
        model_uri = f"runs:/{mlflow.active_run().info.run_id}/propensity_model"
        mlflow.register_model(model_uri, "propensity_score_model")
        
        return mlflow.active_run().info.run_id

# 使用示例
track_causal_model(
    experiment_name="recommendation_din_fusion",
    model_params={
        "propensity_model": "LogisticRegression",
        "matching_method": "nearest_neighbor",
        "fusion_method": "inverse_variance_weighting"
    },
    metrics={
        "fused_effect_size": 0.0508,
        "fused_se": 0.0062,
        "rct_weight": 0.542,
        "efficiency_gain": 0.88
    },
    artifacts={
        "propensity_model": lr_model,
        "fused_results": fusion_result,
        "robustness_checks": robustness_report
    }
)