【重复测量纵向数据】Python05广义线性混合模型(Generalized Linear Mixed Models,GLMM)

# ============================================# 广义线性混合模型 (GLMM) 分析 - Python实现版# ============================================import osimport sysimport warningsimport numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport seaborn as snsfrom pathlib import Pathimport statsmodels.api as smfrom statsmodels.formula.api import mixedlm, gee, poisson, negativebinomialimport statsmodels.formula.api as smffrom statsmodels.genmod.families import Binomial, Poisson, NegativeBinomialfrom statsmodels.genmod.generalized_linear_model import GLMfrom statsmodels.tools.tools import add_constantfrom sklearn.preprocessing import LabelEncoder, StandardScalerfrom scipy import statsimport pingouin as pgimport jsonfrom datetime import datetimeimport matplotlibmatplotlib.use('Agg')  # 使用非交互式后端warnings.filterwarnings('ignore')# ============================================# 0. 环境设置与包安装# ============================================print("=" * 60)print("广义线性混合模型 (GLMM) 分析 - Python实现")print("=" * 60)# 尝试导入所需包，如果失败则安装required_packages = ['pandas', 'numpy', 'scipy', 'statsmodels','matplotlib', 'seaborn', 'pingouin', 'openpyxl']import subprocessimport importlibfor package in required_packages:    try:        importlib.import_module(package)    except ImportError:print(f"正在安装 {package}...")        subprocess.check_call([sys.executable, "-m", "pip", "install", package])# 重新导入import pandas as pdimport numpy as npimport statsmodels.api as smimport matplotlib.pyplot as pltimport seaborn as snsimport pingouin as pg# 设置中文字体（如果系统支持）try:    plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS', 'DejaVu Sans']    plt.rcParams['axes.unicode_minus'] = Falseexcept:    pass# ============================================# 1. 路径设置与文件夹创建# ============================================# 自动查找桌面路径if os.name == 'nt':  # Windows    desktop_path = os.path.join(os.path.join(os.environ['USERPROFILE']), 'Desktop')else:  # Mac/Linux    desktop_path = os.path.join(os.path.join(os.path.expanduser('~')), 'Desktop')# 设置工作路径work_dir = os.path.join(desktop_path, '5-GLMM结果')if not os.path.exists(work_dir):    os.makedirs(work_dir)print(f"创建结果文件夹: {work_dir}")# 数据文件路径data_file = os.path.join(desktop_path, 'longitudinal_data.xlsx')if not os.path.exists(data_file):print(f"警告: 数据文件未找到: {data_file}")print("请确保 'longitudinal_data.xlsx' 文件在桌面上")# 创建示例数据用于测试print("创建示例数据用于演示...")    create_sample_data = Trueelse:    create_sample_data = False# ============================================# 2. 数据读取与预处理# ============================================def load_and_preprocess_data():"""加载和预处理数据"""print("\n1. 数据读取与预处理")print("-" * 40)if create_sample_data:# 创建示例数据        np.random.seed(123)        n_subjects = 100        n_timepoints = 5        n_total = n_subjects * n_timepoints        data = pd.DataFrame({'ID': np.repeat(range(1, n_subjects + 1), n_timepoints),'Time': np.tile(range(1, n_timepoints + 1), n_subjects),'Outcome_Continuous': np.random.normal(50, 10, n_total),'Outcome_Binary': np.random.binomial(1, 0.8, n_total),'Outcome_Count': np.random.poisson(5, n_total),'Stress_Level': np.random.normal(50, 15, n_total),'Quality_of_Life': np.random.normal(70, 10, n_total),'Medication_Adherence': np.random.normal(80, 10, n_total),'Baseline_Score': np.random.normal(60, 12, n_total),'School_ID': np.random.choice(['A', 'B', 'C'], n_total),'Genetic_Marker': np.random.choice([0, 1], n_total),'SES': np.random.choice(['Low', 'Medium', 'High'], n_total),'Latent_Class': np.random.choice([1, 2, 3], n_total),'Sex': np.random.choice(['M', 'F'], n_total),'Treatment': np.random.choice(['Control', 'Treatment'], n_total),'Age': np.random.normal(45, 10, n_total),'Variable2': np.random.normal(30, 5, n_total),'Survival_Time': np.random.exponential(100, n_total),'Event_Status': np.random.choice([0, 1], n_total, p=[0.7, 0.3])        })# 添加一些缺失值for col in ['Outcome_Continuous', 'Stress_Level', 'Quality_of_Life']:            idx = np.random.choice(n_total, size=int(n_total * 0.05), replace=False)            data.loc[idx, col] = np.nanprint("使用示例数据 (模拟数据)")else:# 读取真实数据        try:            data = pd.read_excel(data_file, sheet_name='Full_Dataset')print(f"成功读取数据: {data_file}")        except Exception as e:print(f"读取数据失败: {e}")print("创建示例数据用于演示...")return load_and_preprocess_data()  # 递归调用创建示例数据# 查看数据结构print(f"数据形状: {data.shape}")print(f"变量名: {list(data.columns)}")print("\n数据前5行:")print(data.head())# 检查二分类结局的分布if'Outcome_Binary'in data.columns:        binary_dist = data['Outcome_Binary'].value_counts(dropna=False)print(f"\n二分类结局分布:\n{binary_dist}")if binary_dist.get(0, 0) < 20 or binary_dist.get(1, 0) < 20:print("注意: Outcome_Binary 变量变异较少，但可以尝试分析。")# 检查计数结局的分布if'Outcome_Count'in data.columns:        count_stats = data['Outcome_Count'].describe()print(f"\n计数结局描述统计:\n{count_stats}")        variance = data['Outcome_Count'].var()        mean = data['Outcome_Count'].mean()if variance / mean > 1.5:print(f"注意: Outcome_Count 变量过度离散，方差/均值 = {variance / mean:.2f}")# 变量重命名（简化）    rename_dict = {'ID': 'ID','Time': 'Time','Outcome_Continuous': 'Outcome_Cont','Outcome_Binary': 'Outcome_Bin','Outcome_Count': 'Outcome_Count','Stress_Level': 'Stress','Quality_of_Life': 'QoL','Medication_Adherence': 'Med_Adhere','Baseline_Score': 'Base_Score','School_ID': 'School','Genetic_Marker': 'Genetic','SES': 'SES','Latent_Class': 'Latent','Sex': 'Sex','Treatment': 'Treatment','Age': 'Age','Variable2': 'Variable2','Survival_Time': 'Survival_Time','Event_Status': 'Event_Status'    }# 只重命名存在的列    rename_dict = {k: v for k, v in rename_dict.items() if k in data.columns}    data = data.rename(columns=rename_dict)# 转换分类变量    categorical_vars = ['ID', 'Time', 'Outcome_Bin', 'Sex', 'Treatment','School', 'Genetic', 'SES', 'Latent', 'Event_Status']for var in categorical_vars:if var in data.columns:            data[var] = data[var].astype('category')# 检查缺失值    missing_summary = data.isnull().sum()    missing_summary = missing_summary[missing_summary > 0]if len(missing_summary) > 0:print(f"\n缺失值统计:\n{missing_summary}")else:print("\n没有缺失值")# 保存处理后的数据    processed_data_path = os.path.join(work_dir, 'processed_data.csv')    data.to_csv(processed_data_path, index=False)print(f"\n处理后的数据已保存: {processed_data_path}")return data# ============================================# 3. 描述性统计# ============================================def descriptive_statistics(data):"""计算描述性统计"""print("\n2. 描述性统计")print("-" * 40)# 连续变量描述    cont_vars = ['Outcome_Cont', 'Outcome_Count', 'Variable2', 'Med_Adhere','Stress', 'QoL', 'Age', 'Base_Score', 'Survival_Time']# 只选择存在的变量    cont_vars = [var for var in cont_vars if var in data.columns]    desc_cont = pd.DataFrame()for var in cont_vars:if var in data.columns:            var_data = data[var].dropna()            desc_cont = pd.concat([desc_cont, pd.DataFrame({'Variable': [var],'Mean': [var_data.mean()],'SD': [var_data.std()],'Min': [var_data.min()],'Max': [var_data.max()],'N': [len(var_data)],'Missing': [data[var].isnull().sum()]            })], ignore_index=True)# 四舍五入到3位小数    desc_cont = desc_cont.round(3)print("连续变量描述统计:")print(desc_cont)# 分类变量频数    cat_vars = ['Outcome_Bin', 'Sex', 'Treatment', 'Genetic', 'SES', 'Latent', 'Event_Status']    cat_vars = [var for var in cat_vars if var in data.columns]    cat_tables = {}for var in cat_vars:        cat_tables[var] = data[var].value_counts(dropna=False)print("\n分类变量频数统计:")for var, table in cat_tables.items():print(f"\n{var}:")print(table)# 保存描述性统计    desc_cont_path = os.path.join(work_dir, 'descriptive_stats.csv')    desc_cont.to_csv(desc_cont_path, index=False)# 保存分类变量频数    cat_freq = pd.DataFrame()for var, table in cat_tables.items():        temp_df = pd.DataFrame({'Variable': var,'Category': table.index.astype(str),'Frequency': table.values        })        cat_freq = pd.concat([cat_freq, temp_df], ignore_index=True)    cat_freq_path = os.path.join(work_dir, 'categorical_frequencies.csv')    cat_freq.to_csv(cat_freq_path, index=False)print(f"\n描述性统计已保存到: {work_dir}")return desc_cont, cat_tables# ============================================# 4. 广义线性混合模型 (GLMM) 分析# ============================================def glmm_analysis(data):"""执行GLMM分析"""print("\n3. 广义线性混合模型 (GLMM) 分析")print("-" * 40)# 创建新二分类变量（基于Outcome_Count的中位数）if'Outcome_Count'in data.columns:        median_count = data['Outcome_Count'].median()        data['High_Count'] = (data['Outcome_Count'] > median_count).astype(int)        data['High_Count'] = data['High_Count'].astype('category')print(f"创建新二分类变量 High_Count，中位数阈值: {median_count}")print(f"High_Count 分布:\n{data['High_Count'].value_counts()}")    results = {}