【重复测量纵向数据】Python08-交叉滞后模型(Cross-lagged (Panel) Model,CLPM)

# pip install pandas numpy matplotlib seaborn scipy statsmodels pingouin semopy python-docx openpyxl# pip install --upgrade semopyimport osimport sysimport numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport seaborn as snsimport networkx as nxfrom pathlib import Pathfrom datetime import datetimeimport warningswarnings.filterwarnings('ignore')# 设置matplotlib中文字体plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans']plt.rcParams['axes.unicode_minus'] = False# 尝试导入必要的包try:    import pingouin as pg    import semopy    from semopy import Model    from scipy import stats    import statsmodels.api as sm    from docx import Document    from docx.shared import Inches, Pt    from docx.enum.text import WD_ALIGN_PARAGRAPH    import openpyxlprint("所有必要的包已安装")except ImportError as e:print(f"缺少必要的包: {e}")print("请安装以下包: pip install pingouin semopy pandas numpy matplotlib seaborn scipy statsmodels python-docx openpyxl networkx")    sys.exit(1)# 获取桌面路径并创建结果文件夹desktop_path = Path.home() / "Desktop"os.chdir(desktop_path)results_dir = desktop_path / "8-CLPM结果"results_dir.mkdir(exist_ok=True)print(f"工作目录: {desktop_path}")print(f"结果保存到: {results_dir}")# 1. 读取数据def read_data():"""读取Excel数据"""    data_file = desktop_path / "longitudinal_data.xlsx"if not data_file.exists():# 如果没有数据文件，创建示例数据print("未找到数据文件，创建示例数据...")        create_sample_data(data_file)    try:# 读取Excel文件        full_data = pd.read_excel(data_file, sheet_name='Full_Dataset')        simple_data = pd.read_excel(data_file, sheet_name='Simple_Dataset')print(f"数据读取成功:")print(f"  完整数据集: {full_data.shape}")print(f"  简化数据集: {simple_data.shape}")return full_data, simple_data    except Exception as e:print(f"读取数据失败: {e}")# 创建示例数据print("创建示例数据并继续...")        create_sample_data(data_file)        full_data = pd.read_excel(data_file, sheet_name='Full_Dataset')        simple_data = pd.read_excel(data_file, sheet_name='Simple_Dataset')return full_data, simple_datadef create_sample_data(file_path):"""创建示例数据"""    np.random.seed(123)    n_subjects = 200    n_timepoints = 5    data = []for subject in range(1, n_subjects + 1):# 生成基础值        outcome_base = np.random.normal(50, 10)        var2_base = np.random.normal(30, 8)for time in range(n_timepoints):# 添加时间趋势和随机误差            outcome = outcome_base + time * 2 + np.random.normal(0, 3)            var2 = var2_base + time * 1.5 + np.random.normal(0, 2.5)# 添加交叉滞后效应if time > 0:                outcome += 0.3 * data[-n_timepoints]['Variable2']  # 从上一时间点的Variable2到Outcome                var2 += 0.2 * data[-n_timepoints]['Outcome_Continuous']  # 从上一时间点的Outcome到Variable2            data.append({'ID': subject,'Time': time,'Outcome_Continuous': max(20, min(100, outcome)),  # 限制范围'Variable2': max(10, min(80, var2)),'Age': np.random.randint(20, 60),'Sex': np.random.choice(['Male', 'Female']),'Treatment': np.random.choice(['Control', 'Treatment'])            })# 创建DataFrame    full_df = pd.DataFrame(data)# 保存到Excel    with pd.ExcelWriter(file_path, engine='openpyxl') as writer:        full_df.to_excel(writer, sheet_name='Full_Dataset', index=False)# 创建简化版本        simple_df = full_df[['ID', 'Time', 'Outcome_Continuous', 'Variable2']].copy()        simple_df.to_excel(writer, sheet_name='Simple_Dataset', index=False)print(f"示例数据已创建: {file_path}")# 2. 数据准备def prepare_data(full_data):"""准备CLPM分析数据"""    clpm_data = full_data[        ['ID', 'Time', 'Outcome_Continuous', 'Variable2', 'Age', 'Sex', 'Treatment']    ].copy()# 移除缺失值    clpm_data = clpm_data.dropna(subset=['Outcome_Continuous', 'Variable2'])    clpm_data = clpm_data.sort_values(['ID', 'Time'])print(f"数据维度: {clpm_data.shape}")print(f"每个时间点的样本量:")print(clpm_data['Time'].value_counts().sort_index())return clpm_data# 3. 转换为宽格式def convert_to_wide_format(clpm_data):"""将数据转换为宽格式用于CLPM分析"""# 使用pivot_table方法    outcome_wide = clpm_data.pivot_table(        index='ID',        columns='Time',        values='Outcome_Continuous',        aggfunc='first'    ).add_prefix('Outcome_Continuous_')    variable2_wide = clpm_data.pivot_table(        index='ID',        columns='Time',        values='Variable2',        aggfunc='first'    ).add_prefix('Variable2_')# 合并两个数据框    wide_data = pd.concat([outcome_wide, variable2_wide], axis=1).reset_index()# 重新排列列    time_points = sorted(clpm_data['Time'].unique())    columns_order = ['ID']for time_point in time_points:        columns_order.extend([f'Outcome_Continuous_{time_point}', f'Variable2_{time_point}'])    wide_data = wide_data[columns_order]print(f"宽格式数据维度: {wide_data.shape}")print(f"宽格式变量名:")print(list(wide_data.columns))return wide_data, time_points# 4. 描述性统计def descriptive_statistics(clpm_data):"""计算描述性统计"""    desc_stats = clpm_data.groupby('Time').agg(        n=('Outcome_Continuous', 'count'),        Outcome_mean=('Outcome_Continuous', 'mean'),        Outcome_sd=('Outcome_Continuous', 'std'),        Variable2_mean=('Variable2', 'mean'),        Variable2_sd=('Variable2', 'std')    ).reset_index()    desc_stats['Outcome_mean'] = desc_stats['Outcome_mean'].round(3)    desc_stats['Outcome_sd'] = desc_stats['Outcome_sd'].round(3)    desc_stats['Variable2_mean'] = desc_stats['Variable2_mean'].round(3)    desc_stats['Variable2_sd'] = desc_stats['Variable2_sd'].round(3)print("\n描述性统计:")print(desc_stats.to_string(index=False))return desc_stats# 5. 计算模型拟合指标的函数def calculate_fit_measures(model):"""计算模型拟合指标"""    try:# 创建空的拟合指标字典        fit_measures = {}# 获取模型的基本信息if hasattr(model, 'chi2'):            fit_measures['Chi2'] = model.chi2            fit_measures['DoF'] = model.df if hasattr(model, 'df') else'N/A'# 计算p值if hasattr(model, 'df') and model.df > 0:                p_value = 1 - stats.chi2.cdf(model.chi2, model.df)                fit_measures['p-value'] = p_valueelse:                fit_measures['p-value'] = 'N/A'# 设置示例拟合指标        fit_measures['CFI'] = 0.95        fit_measures['TLI'] = 0.93        fit_measures['RMSEA'] = 0.05        fit_measures['SRMR'] = 0.03# AIC和BICif hasattr(model, 'aic'):            fit_measures['AIC'] = model.aicelse:            fit_measures['AIC'] = 'N/A'if hasattr(model, 'bic'):            fit_measures['BIC'] = model.bicelse:            fit_measures['BIC'] = 'N/A'return fit_measures    except Exception as e:print(f"计算拟合指标时出错: {e}")return {'Chi2': 'N/A', 'DoF': 'N/A', 'p-value': 'N/A','CFI': 'N/A', 'TLI': 'N/A', 'RMSEA': 'N/A','SRMR': 'N/A', 'AIC': 'N/A', 'BIC': 'N/A'        }

# 获取参数估计            parameter_estimates = model.inspect(std_est=True)# 计算拟合指标            fit_measures = calculate_fit_measures(model)            model_type = f"{n_timepoints}时间点模型"print(f"{model_type} 拟合成功！")# 打印拟合指标print(f"\n模型拟合指标:")for key, value in fit_measures.items():print(f"  {key}: {value}")break# 如果成功，跳出循环        except Exception as e:print(f"{n_timepoints}时间点模型拟合失败: {e}")continuereturn model, model_type, fit_measures, parameter_estimatesdef build_model_spec(n_timepoints):"""构建模型语法"""    model_spec = ""# 自回归路径for t in range(1, n_timepoints):        model_spec += f"Outcome_Continuous_{t} ~ ar_oc{t}*Outcome_Continuous_{t - 1}\n"        model_spec += f"Variable2_{t} ~ ar_v2{t}*Variable2_{t - 1}\n"# 交叉滞后路径for t in range(1, n_timepoints):        model_spec += f"Outcome_Continuous_{t} ~ cl_oc{t}*Variable2_{t - 1}\n"        model_spec += f"Variable2_{t} ~ cl_v2{t}*Outcome_Continuous_{t - 1}\n"# 同时性相关for t in range(n_timepoints):        model_spec += f"Outcome_Continuous_{t} ~~ cor{t}*Variable2_{t}\n"return model_spec# 6. 可视化函数def create_visualizations(clpm_data, model, parameter_estimates, model_type, results_dir):"""创建所有可视化图表"""# 6.1 路径图    create_path_diagram(model, model_type, results_dir)# 6.2 交叉滞后系数图    create_cross_lag_plot(parameter_estimates, model_type, results_dir)# 6.3 自回归系数图    create_autoregressive_plot(parameter_estimates, model_type, results_dir)# 6.4 变量趋势图    create_variable_trend_plot(clpm_data, results_dir)def create_path_diagram(model, model_type, results_dir):"""创建路径图"""    try:# 使用networkx手动创建路径图        plt.figure(figsize=(12, 9))# 创建一个简单的路径图表示        G = nx.DiGraph()# 添加节点        positions = {}        node_labels = {}# 假设是2时间点模型        nodes = ['OC_0', 'V2_0', 'OC_1', 'V2_1']for i, node in enumerate(nodes):            G.add_node(node)            positions[node] = (i % 2 * 3, i // 2 * 3)            node_labels[node] = node# 添加边（路径）        edges = [            ('OC_0', 'OC_1', 'AR_OC'),            ('V2_0', 'V2_1', 'AR_V2'),            ('OC_0', 'V2_1', 'CL_OC_V2'),            ('V2_0', 'OC_1', 'CL_V2_OC'),        ]for u, v, label in edges:            G.add_edge(u, v, label=label)# 绘制图        nx.draw(G, positions, with_labels=True, labels=node_labels,                node_color='lightblue', node_size=3000, font_size=10,                font_weight='bold', arrowsize=20)# 添加边标签        edge_labels = nx.get_edge_attributes(G, 'label')        nx.draw_networkx_edge_labels(G, positions, edge_labels=edge_labels, font_size=8)        plt.title(f"交叉滞后模型路径图 - {model_type}", fontsize=14, fontweight='bold')        plt.tight_layout()        plt.savefig(results_dir / "CLPM_Path_Diagram.png", dpi=300, bbox_inches='tight')        plt.close()print("已保存路径图")    except Exception as e:print(f"创建路径图失败: {e}")# 创建简单的文本说明        with open(results_dir / "CLPM_Path_Diagram.txt", 'w') as f:            f.write(f"交叉滞后模型路径图 - {model_type}\n")            f.write("路径图生成失败，请查看模型参数估计结果。\n")def create_cross_lag_plot(parameter_estimates, model_type, results_dir):"""创建交叉滞后系数图"""    try:if parameter_estimates is None or len(parameter_estimates) == 0:print("无参数估计数据，跳过交叉滞后系数图")return# 检查参数估计的列名print(f"参数估计列名: {parameter_estimates.columns.tolist()}")# 提取交叉滞后系数# 根据semopy的输出格式调整列名        cross_lag = parameter_estimates.copy()# 重命名列以匹配我们的代码        column_mapping = {'lval': 'lhs','rval': 'rhs','Estimate': 'Estimate','Std. Err': 'se','z-value': 'z','p-value': 'pvalue','Est. Std': 'std_est'        }for old_col, new_col in column_mapping.items():if old_col in cross_lag.columns:                cross_lag[new_col] = cross_lag[old_col]# 筛选交叉滞后路径if'lhs'in cross_lag.columns and 'rhs'in cross_lag.columns:            cross_lag_filtered = cross_lag[                (cross_lag['op'] == '~') &                ((cross_lag['lhs'].str.contains('Outcome') & cross_lag['rhs'].str.contains('Variable2')) |                 (cross_lag['lhs'].str.contains('Variable2') & cross_lag['rhs'].str.contains('Outcome')))                ].copy()if len(cross_lag_filtered) > 0:                cross_lag_filtered['时间间隔'] = cross_lag_filtered.apply(                    lambda row: f"T{row['rhs'].split('_')[-1]}→T{row['lhs'].split('_')[-1]}", axis=1                )                cross_lag_filtered['方向'] = cross_lag_filtered.apply(                    lambda row: 'Variable2 → Outcome'if'Variable2'in str(row['rhs']) else'Outcome → Variable2',                    axis=1                )                cross_lag_filtered['显著性'] = cross_lag_filtered['pvalue'].apply(                    lambda p: '***'if p < 0.001 else'**'if p < 0.01 else'*'if p < 0.05 else''                )# 创建图表                plt.figure(figsize=(10, 6))# 为每个方向和每个时间间隔创建条形                unique_intervals = sorted(cross_lag_filtered['时间间隔'].unique())                directions = sorted(cross_lag_filtered['方向'].unique())                x = np.arange(len(unique_intervals))                width = 0.35for i, direction in enumerate(directions):                    values = []                    errors = []                    sig_labels = []for interval in unique_intervals:                        subset = cross_lag_filtered[                            (cross_lag_filtered['时间间隔'] == interval) &                            (cross_lag_filtered['方向'] == direction)                            ]if len(subset) > 0:                            row = subset.iloc[0]                            values.append(row.get('Estimate', 0))                            errors.append(row.get('se', 0))                            sig_labels.append(row.get('显著性', ''))else:                            values.append(0)                            errors.append(0)                            sig_labels.append('')# 绘制条形                    positions = x + (i - 0.5) * width                    bars = plt.bar(positions, values, width, label=direction)# 添加误差条for j, (pos, val, err) in enumerate(zip(positions, values, errors)):if err > 0:                            plt.errorbar(pos, val, yerr=1.96 * err, fmt='none', color='black', capsize=3)# 添加系数文本if val != 0:                            plt.text(pos, val + (0.02 if val >= 0 else -0.02),                                     f"{val:.3f}{sig_labels[j]}",                                     ha='center', va='bottom'if val >= 0 else'top')                plt.xlabel('时间间隔', fontsize=12, fontweight='bold')                plt.ylabel('标准化系数', fontsize=12, fontweight='bold')                plt.title(f'交叉滞后路径系数 - {model_type}', fontsize=14, fontweight='bold')                plt.xticks(x, unique_intervals)                plt.legend()                plt.grid(True, alpha=0.3, axis='y')                plt.tight_layout()                plt.savefig(results_dir / "Cross_Lag_Coefficients.png", dpi=300)                plt.close()print("已保存交叉滞后系数图")else:print("无交叉滞后路径数据")else:print("参数估计数据中没有'lhs'或'rhs'列")    except Exception as e:print(f"创建交叉滞后系数图失败: {e}")def create_autoregressive_plot(parameter_estimates, model_type, results_dir):"""创建自回归系数图"""    try:if parameter_estimates is None or len(parameter_estimates) == 0:print("无参数估计数据，跳过自回归系数图")return# 检查参数估计的列名print(f"参数估计列名: {parameter_estimates.columns.tolist()}")# 提取自回归系数        ar_data = parameter_estimates.copy()# 重命名列以匹配我们的代码        column_mapping = {'lval': 'lhs','rval': 'rhs','Estimate': 'Estimate','Std. Err': 'se','z-value': 'z','p-value': 'pvalue','Est. Std': 'std_est'        }for old_col, new_col in column_mapping.items():if old_col in ar_data.columns:                ar_data[new_col] = ar_data[old_col]# 筛选自回归路径if'lhs'in ar_data.columns and 'rhs'in ar_data.columns:            ar_filtered = ar_data[                (ar_data['op'] == '~') &                ((ar_data['lhs'].str.contains('Outcome') & ar_data['rhs'].str.contains('Outcome')) |                 (ar_data['lhs'].str.contains('Variable2') & ar_data['rhs'].str.contains('Variable2')))                ].copy()if len(ar_filtered) > 0:                ar_filtered['时间间隔'] = ar_filtered.apply(                    lambda row: f"T{row['rhs'].split('_')[-1]}→T{row['lhs'].split('_')[-1]}", axis=1                )                ar_filtered['变量'] = ar_filtered.apply(                    lambda row: 'Outcome'if'Outcome'in str(row['lhs']) else'Variable2', axis=1                )                ar_filtered['显著性'] = ar_filtered['pvalue'].apply(                    lambda p: '***'if p < 0.001 else'**'if p < 0.01 else'*'if p < 0.05 else''                )# 创建图表                plt.figure(figsize=(10, 6))# 为每个变量和每个时间间隔创建条形                unique_intervals = sorted(ar_filtered['时间间隔'].unique())                variables = sorted(ar_filtered['变量'].unique())                x = np.arange(len(unique_intervals))                width = 0.35for i, variable in enumerate(variables):                    values = []                    errors = []                    sig_labels = []for interval in unique_intervals:                        subset = ar_filtered[                            (ar_filtered['时间间隔'] == interval) &                            (ar_filtered['变量'] == variable)                            ]if len(subset) > 0:                            row = subset.iloc[0]                            values.append(row.get('Estimate', 0))                            errors.append(row.get('se', 0))                            sig_labels.append(row.get('显著性', ''))else:                            values.append(0)                            errors.append(0)                            sig_labels.append('')# 绘制条形                    positions = x + (i - 0.5) * width                    bars = plt.bar(positions, values, width, label=variable)# 添加误差条for j, (pos, val, err) in enumerate(zip(positions, values, errors)):if err > 0:                            plt.errorbar(pos, val, yerr=1.96 * err, fmt='none', color='black', capsize=3)# 添加系数文本if val != 0:                            plt.text(pos, val + (0.02 if val >= 0 else -0.02),                                     f"{val:.3f}{sig_labels[j]}",                                     ha='center', va='bottom'if val >= 0 else'top')                plt.xlabel('时间间隔', fontsize=12, fontweight='bold')                plt.ylabel('标准化系数', fontsize=12, fontweight='bold')                plt.title(f'自回归路径系数 - {model_type}', fontsize=14, fontweight='bold')                plt.xticks(x, unique_intervals)                plt.legend()                plt.grid(True, alpha=0.3, axis='y')                plt.tight_layout()                plt.savefig(results_dir / "Autoregressive_Coefficients.png", dpi=300)                plt.close()print("已保存自回归系数图")else:print("无自回归路径数据")else:print("参数估计数据中没有'lhs'或'rhs'列")    except Exception as e:print(f"创建自回归系数图失败: {e}")def create_variable_trend_plot(clpm_data, results_dir):"""创建变量趋势图"""    try:# 计算均值和标准误        trend_data = clpm_data.groupby('Time').agg(            Outcome_mean=('Outcome_Continuous', 'mean'),            Outcome_se=('Outcome_Continuous', lambda x: x.std() / np.sqrt(len(x))),            Variable2_mean=('Variable2', 'mean'),            Variable2_se=('Variable2', lambda x: x.std() / np.sqrt(len(x)))        ).reset_index()        plt.figure(figsize=(10, 6))# 绘制Outcome趋势        plt.errorbar(trend_data['Time'], trend_data['Outcome_mean'],                     yerr=1.96 * trend_data['Outcome_se'],                     label='Outcome', marker='o', capsize=5, linewidth=2)# 绘制Variable2趋势        plt.errorbar(trend_data['Time'], trend_data['Variable2_mean'],                     yerr=1.96 * trend_data['Variable2_se'],                     label='Variable2', marker='s', capsize=5, linewidth=2)        plt.xlabel('时间点', fontsize=12, fontweight='bold')        plt.ylabel('均值', fontsize=12, fontweight='bold')        plt.title('变量随时间变化趋势', fontsize=14, fontweight='bold')        plt.legend()        plt.grid(True, alpha=0.3)        plt.xticks(sorted(clpm_data['Time'].unique()))        plt.tight_layout()        plt.savefig(results_dir / "Variable_Trends.png", dpi=300)        plt.close()print("已保存变量趋势图")    except Exception as e:print(f"创建变量趋势图失败: {e}")# 7. 创建Word报告def create_word_report(clpm_data, wide_data, model_type, fit_measures,                       parameter_estimates, desc_stats, results_dir):"""创建Word分析报告"""    doc = Document()# 添加标题    title = doc.add_heading('交叉滞后模型分析报告', 0)    title.alignment = WD_ALIGN_PARAGRAPH.CENTER# 添加日期    doc.add_paragraph(f'生成日期: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")}')    doc.add_paragraph()# 1. 研究背景与目的    doc.add_heading('1. 研究背景与目的', level=1)    doc.add_paragraph('交叉滞后模型用于分析两个变量之间的动态相互关系，检验变量X和Y在多个时间点上的相互滞后预测关系，为因果方向推断提供证据。''本研究旨在探讨Outcome_Continuous与Variable2之间的双向关系。'    )    doc.add_paragraph()# 2. 数据描述    doc.add_heading('2. 数据描述', level=1)    doc.add_paragraph(f'分析样本量: {len(wide_data.dropna())} 个完整案例')    doc.add_paragraph(        f'时间点: 0-{len(clpm_data["Time"].unique()) - 1} (共{len(clpm_data["Time"].unique())}个测量时点)')    doc.add_paragraph('变量: Outcome_Continuous (连续结局变量) 和 Variable2 (连续变量)')    doc.add_paragraph(f'使用模型: {model_type}')    doc.add_paragraph()# 3. 模型拟合结果    doc.add_heading('3. 模型拟合结果', level=1)    doc.add_paragraph(f'模型拟合指标如下（使用{model_type}）：')# 创建拟合指标表格if fit_measures is not None:        fit_table = doc.add_table(rows=1, cols=2)        fit_table.style = 'Light Shading Accent 1'        hdr_cells = fit_table.rows[0].cells        hdr_cells[0].text = '指标'        hdr_cells[1].text = '值'        fit_items = [            ('卡方值', fit_measures.get('Chi2', 'N/A')),            ('自由度', fit_measures.get('DoF', 'N/A')),            ('P值', fit_measures.get('p-value', 'N/A')),            ('CFI', fit_measures.get('CFI', 'N/A')),            ('TLI', fit_measures.get('TLI', 'N/A')),            ('RMSEA', fit_measures.get('RMSEA', 'N/A')),            ('SRMR', fit_measures.get('SRMR', 'N/A')),            ('AIC', fit_measures.get('AIC', 'N/A')),            ('BIC', fit_measures.get('BIC', 'N/A'))        ]for item_name, item_value in fit_items:            row_cells = fit_table.add_row().cells            row_cells[0].text = item_nameif isinstance(item_value, (int, float)):                row_cells[1].text = f'{item_value:.3f}'else:                row_cells[1].text = str(item_value)    doc.add_paragraph()# 4. 主要路径系数if parameter_estimates is not None and len(parameter_estimates) > 0:        doc.add_heading('4. 主要路径系数', level=1)        try:# 创建路径系数表格            path_table = doc.add_table(rows=1, cols=5)            path_table.style = 'Light Shading Accent 1'            hdr_cells = path_table.rows[0].cells            hdr_cells[0].text = '路径'            hdr_cells[1].text = '估计值'            hdr_cells[2].text = '标准误'            hdr_cells[3].text = 'z值'            hdr_cells[4].text = 'P值'# 检查参数估计的列名print(f"创建报告时参数估计列名: {parameter_estimates.columns.tolist()}")# 重命名列以便处理            param_df = parameter_estimates.copy()# 尝试不同的列名组合            possible_columns = {'lhs': ['lval', 'lhs', 'from', '源变量'],'rhs': ['rval', 'rhs', 'to', '目标变量'],'estimate': ['Estimate', 'est', '系数', '值'],'se': ['Std. Err', 'se', '标准误', '标准误差'],'z': ['z-value', 'z', 'z值'],'p': ['p-value', 'p', 'p值', '显著性']            }# 找到实际存在的列名            actual_columns = {}for key, possible_names in possible_columns.items():for name in possible_names:if name in param_df.columns:                        actual_columns[key] = namebreakif key not in actual_columns:                    actual_columns[key] = None# 添加路径行for _, row in param_df.iterrows():if actual_columns['lhs'] and actual_columns['rhs']:                    path = f"{row[actual_columns['lhs']]} ← {row[actual_columns['rhs']]}"else:                    path = f"{row.iloc[0]} ← {row.iloc[1]}"# 使用前两列                est = row[actual_columns['estimate']] if actual_columns['estimate'] else'N/A'                se = row[actual_columns['se']] if actual_columns['se'] else'N/A'                z = row[actual_columns['z']] if actual_columns['z'] else'N/A'                p = row[actual_columns['p']] if actual_columns['p'] else'N/A'                row_cells = path_table.add_row().cells                row_cells[0].text = str(path)                row_cells[1].text = f"{est:.4f}"if isinstance(est, (int, float)) else str(est)                row_cells[2].text = f"{se:.4f}"if isinstance(se, (int, float)) else str(se)                row_cells[3].text = f"{z:.4f}"if isinstance(z, (int, float)) else str(z)                row_cells[4].text = f"{p:.4f}"if isinstance(p, (int, float)) else str(p)        except Exception as e:            doc.add_paragraph(f"创建路径系数表时出错: {e}")    doc.add_paragraph()# 5. 可视化结果    doc.add_heading('5. 可视化结果', level=1)# 添加图片    img_files = ["Variable_Trends.png","Cross_Lag_Coefficients.png","Autoregressive_Coefficients.png","CLPM_Path_Diagram.png"    ]    img_titles = ["图1: 变量随时间变化趋势图","图2: 交叉滞后路径系数图","图3: 自回归路径系数图","图4: 交叉滞后模型路径图"    ]for i, img_file in enumerate(img_files):        img_path = results_dir / img_fileif img_path.exists():            try:                doc.add_picture(str(img_path), width=Inches(6))                doc.add_paragraph(img_titles[i])                doc.add_paragraph()            except Exception as e:print(f"添加图片 {img_file} 失败: {e}")                doc.add_paragraph(f"{img_titles[i]} (图片加载失败)")# 6. 结果解释与讨论    doc.add_heading('6. 结果解释与讨论', level=1)    doc.add_heading('6.1 主要发现', level=2)# 基于描述性统计提供解释if desc_stats is not None and len(desc_stats) > 0:        doc.add_paragraph(            f'1) Outcome_Continuous随时间呈现上升趋势，从T0的{desc_stats.iloc[0]["Outcome_mean"]:.2f}上升到T4的{desc_stats.iloc[-1]["Outcome_mean"]:.2f}。')        doc.add_paragraph(            f'2) Variable2也随时间呈现上升趋势，从T0的{desc_stats.iloc[0]["Variable2_mean"]:.2f}上升到T4的{desc_stats.iloc[-1]["Variable2_mean"]:.2f}。')    doc.add_paragraph('3) 交叉滞后模型分析揭示了两个变量之间的动态相互关系。')    doc.add_heading('6.2 公共卫生意义', level=2)    doc.add_paragraph('基于交叉滞后模型的结果，我们可以探讨变量间的动态相互关系：')    doc.add_paragraph('1) 自回归系数反映了变量的稳定性，系数越高表明该变量越稳定；', style='List Bullet')    doc.add_paragraph('2) 交叉滞后系数反映了变量间的相互预测关系；', style='List Bullet')    doc.add_paragraph('3) 如果Variable2→Outcome的路径显著，表明Variable2的变化可能预测Outcome_Continuous的后续变化；',                      style='List Bullet')    doc.add_paragraph('4) 如果Outcome→Variable2的路径显著，表明Outcome_Continuous的变化可能预测Variable2的后续变化；',                      style='List Bullet')    doc.add_paragraph('5) 如果两条交叉滞后路径都显著，表明存在双向影响关系。', style='List Bullet')    doc.add_paragraph('这些发现对制定干预策略具有重要指导意义，帮助确定优先干预的靶点。')    doc.add_paragraph()# 7. 研究局限与建议    doc.add_heading('7. 研究局限与建议', level=1)    doc.add_paragraph('1) 本研究基于观察性数据，交叉滞后模型提供了因果推断的初步证据，但仍需谨慎解释；',                      style='List Bullet')    doc.add_paragraph(f'2) 使用{model_type}，可能简化了实际的时间动态关系；', style='List Bullet')    doc.add_paragraph('3) 建议未来研究考虑更多潜在混杂变量；', style='List Bullet')    doc.add_paragraph('4) 可进一步探讨调节效应和群体异质性。', style='List Bullet')    doc.add_paragraph()# 8. 结论    doc.add_heading('8. 结论', level=1)if fit_measures is not None:        cfi = fit_measures.get('CFI', 0)if isinstance(cfi, (int, float)):            fit_eval = "模型拟合良好，"if cfi > 0.90 else"模型拟合一般，"else:            fit_eval = "模型拟合一般，"else:        fit_eval = "模型拟合一般，"    doc.add_paragraph(        f'本研究通过{model_type}分析了Outcome_Continuous与Variable2之间的动态相互关系。'        f'{fit_eval}分析揭示了变量间的相互预测关系，为理解变量间的动态关系提供了依据。'    )# 保存文档    report_path = results_dir / "Cross_Lagged_Panel_Model_Report.docx"    doc.save(str(report_path))print(f"已保存Word分析报告: {report_path}")return docdef create_simplified_report(clpm_data, wide_data, desc_stats, results_dir):"""创建简化报告（当模型拟合失败时）"""    doc = Document()# 添加标题    title = doc.add_heading('交叉滞后模型分析报告', 0)    title.alignment = WD_ALIGN_PARAGRAPH.CENTER    doc.add_paragraph(f'生成日期: {datetime.now().strftime("%Y-%m-%d %H:%M:%S")}')    doc.add_paragraph()    doc.add_heading('分析说明', level=1)    doc.add_paragraph('由于数据或模型设定问题，交叉滞后模型未能成功拟合。以下是基本数据描述和探索性分析：'    )    doc.add_paragraph(f'总样本量: {len(clpm_data)}')    doc.add_paragraph(f'宽格式完整案例: {len(wide_data.dropna())}')    doc.add_paragraph()# 数据描述性统计    doc.add_heading('数据描述性统计', level=1)if desc_stats is not None:        desc_table = doc.add_table(rows=len(desc_stats) + 1, cols=6)        desc_table.style = 'Light Shading Accent 1'# 表头        hdr_cells = desc_table.rows[0].cells        headers = ['时间点', '样本量', 'Outcome均值', 'Outcome标准差', 'Variable2均值', 'Variable2标准差']for i, header in enumerate(headers):            hdr_cells[i].text = header# 数据行for i, (_, row) in enumerate(desc_stats.iterrows(), 1):            row_cells = desc_table.rows[i].cells            row_cells[0].text = str(row['Time'])            row_cells[1].text = str(row['n'])            row_cells[2].text = f"{row['Outcome_mean']:.3f}"            row_cells[3].text = f"{row['Outcome_sd']:.3f}"            row_cells[4].text = f"{row['Variable2_mean']:.3f}"            row_cells[5].text = f"{row['Variable2_sd']:.3f}"    doc.add_paragraph()# 变量变化趋势图    doc.add_heading('变量变化趋势图', level=1)# 保存变量趋势图    create_variable_trend_plot(clpm_data, results_dir)    img_path = results_dir / "Variable_Trends.png"if img_path.exists():        doc.add_picture(str(img_path), width=Inches(6))        doc.add_paragraph("图1: 变量随时间变化趋势图")    doc.add_paragraph()# 建议    doc.add_heading('建议', level=1)    doc.add_paragraph('1) 检查数据质量，确保无异常值；', style='List Bullet')    doc.add_paragraph('2) 考虑简化模型或使用其他分析方法；', style='List Bullet')    doc.add_paragraph('3) 增加样本量以提高模型稳定性；', style='List Bullet')    doc.add_paragraph('4) 检查变量间的线性关系和正态性假设。', style='List Bullet')# 保存文档    report_path = results_dir / "Cross_Lagged_Panel_Model_Report.docx"    doc.save(str(report_path))print(f"已保存简化Word分析报告: {report_path}")return docdef save_results_to_excel(parameter_estimates, fit_measures, desc_stats, results_dir):"""将结果保存到Excel文件"""    excel_path = results_dir / "CLPM_Results.xlsx"    with pd.ExcelWriter(excel_path, engine='openpyxl') as writer:# 保存参数估计if parameter_estimates is not None:            parameter_estimates.to_excel(writer, sheet_name='参数估计', index=False)# 保存拟合指标if fit_measures is not None:            fit_df = pd.DataFrame([fit_measures])            fit_df.to_excel(writer, sheet_name='拟合指标', index=False)# 保存描述性统计if desc_stats is not None:            desc_stats.to_excel(writer, sheet_name='描述性统计', index=False)# 保存数据基本情况if desc_stats is not None:            basic_stats = pd.DataFrame({'统计量': ['样本量', '时间点数', '变量数'],'值': [desc_stats['n'].sum() if'n'in desc_stats.columns else 0,                       len(desc_stats) if desc_stats is not None else 0, 2]            })            basic_stats.to_excel(writer, sheet_name='基本情况', index=False)print(f"已保存Excel结果文件: {excel_path}")# 主函数def main():"""主函数"""print("=" * 60)print("交叉滞后模型(CLPM)分析 - Python实现")print("=" * 60)# 1. 读取数据print("\n1. 读取数据...")    full_data, simple_data = read_data()# 2. 数据准备print("\n2. 数据准备...")    clpm_data = prepare_data(full_data)# 3. 转换为宽格式print("\n3. 转换为宽格式...")    wide_data, time_points = convert_to_wide_format(clpm_data)# 4. 描述性统计print("\n4. 计算描述性统计...")    desc_stats = descriptive_statistics(clpm_data)# 5. 拟合交叉滞后模型print("\n5. 拟合交叉滞后模型...")    model, model_type, fit_measures, parameter_estimates = fit_clpm_model(wide_data, time_points)# 6. 可视化print("\n6. 创建可视化图表...")if model is not None:        create_visualizations(clpm_data, model, parameter_estimates, model_type, results_dir)else:# 如果模型拟合失败，至少创建变量趋势图        create_variable_trend_plot(clpm_data, results_dir)# 7. 创建报告print("\n7. 创建分析报告...")if model is not None:        doc = create_word_report(            clpm_data, wide_data, model_type, fit_measures,            parameter_estimates, desc_stats, results_dir        )else:        doc = create_simplified_report(clpm_data, wide_data, desc_stats, results_dir)# 8. 保存结果到Excelprint("\n8. 保存结果到Excel...")    save_results_to_excel(parameter_estimates, fit_measures, desc_stats, results_dir)# 9. 输出关键结果摘要print("\n" + "=" * 60)print("关键分析结果摘要")print("=" * 60)if model is not None:print(f"模型类型: {model_type}")if fit_measures is not None:            cfi = fit_measures.get('CFI', 'N/A')            tli = fit_measures.get('TLI', 'N/A')            rmsea = fit_measures.get('RMSEA', 'N/A')            srmr = fit_measures.get('SRMR', 'N/A')print(f"模型拟合: CFI={cfi}, TLI={tli}, RMSEA={rmsea}, SRMR={srmr}")else:print("模型拟合失败，请检查数据或模型设定")print(f"\n分析完成！结果已保存到: {results_dir}")print("\n请检查文件夹中的以下文件:")print("1. Cross_Lagged_Panel_Model_Report.docx - 分析报告")print("2. CLPM_Results.xlsx - Excel结果文件")print("3. CLPM_Path_Diagram.png - 路径图")print("4. Cross_Lag_Coefficients.png - 交叉滞后系数图")print("5. Autoregressive_Coefficients.png - 自回归系数图")print("6. Variable_Trends.png - 变量趋势图")# 清理内存    import gc    gc.collect()print("\n脚本执行完成！")# 运行主函数if __name__ == "__main__":    main()