# -*- coding: utf-8 -*-
"""
地缘政治冲击下的移民政治态度与社会认同
统计分析脚本

分析步骤：
1. 数据读取与清洗
2. 反向计分与量表合成
3. 虚拟变量生成
4. 信度分析（Cronbach's α）
5. 描述性统计与相关分析
6. 中介效应分析（Model 4，Bootstrap = 2000）
7. 调节效应分析（Model 1）
8. 分组中介分析（按俄语流利度、族群分组）
9. 结果导出

依赖库：pandas, numpy, scipy, statsmodels, pingouin, openpyxl
"""

import pandas as pd
import numpy as np
from scipy import stats
from statsmodels.regression.linear_model import OLS
from statsmodels.tools import add_constant
import pingouin as pg
import warnings
warnings.filterwarnings('ignore')

# ============================================================
# 0. 文件路径与参数设置
# ============================================================
DATA_PATH = "最终数据.xlsx"
OUTPUT_PATH = "统计分析结果.xlsx"
N_BOOTSTRAP = 2000     # 与交付文档“Bootstrap法（2000次重复）”一致
N_BOOTSTRAP_GROUP = 1000  # 分组中介的Bootstrap次数
ALPHA = 0.05

# ============================================================
# 1. 数据读取与清洗
# ============================================================
df_raw = pd.read_excel(DATA_PATH, sheet_name=0, header=0)

var_names = [
    'timestamp', 'consent', 'gender', 'age', 'edu', 'country', 'years_ru',
    'listen', 'speak', 'read', 'write', 'ethnic', 'imm_status', 'unfriendly',
    'q11', 'q12', 'q13', 'q14', 'q15', 'q16', 'q17', 'q18', 'q19', 'q20',
    'q21', 'q22', 'q23', 'q24', 'q25', 'q26', 'q27', 'q28', 'q29', 'q30',
    'q31', 'q32', 'q33', 'q34', 'q35', 'q36', 'q37', 'q38', 'q39', 'q40'
]
df_raw.columns = var_names[:len(df_raw.columns)]

df = df_raw[df_raw['consent'] == 1].copy()
print(f"有效样本量：{len(df)}")

# ============================================================
# 2. 反向计分与量表合成
# ============================================================
df['q15_r'] = 6 - df['q15']

geo_items = ['q11', 'q12', 'q13', 'q14', 'q15_r']
df['geo_shock'] = df[geo_items].mean(axis=1, skipna=True)

disc_items = ['q16', 'q17', 'q18', 'q19', 'q20', 'q21', 'q22']
df['disc_percep'] = df[disc_items].mean(axis=1, skipna=True)

pol_items = ['q23', 'q24', 'q25', 'q26', 'q27', 'q28', 'q29', 'q30']
df['pol_att'] = df[pol_items].mean(axis=1, skipna=True)

id_items = ['q31', 'q32', 'q33', 'q34', 'q35', 'q36', 'q37', 'q38']
df['social_id'] = df[id_items].mean(axis=1, skipna=True)

ru_items = ['listen', 'speak', 'read', 'write']
df['ru_prof'] = df[ru_items].mean(axis=1, skipna=True)

# ============================================================
# 3. 虚拟变量生成
# ============================================================
df['female'] = (df['gender'] == 2).astype(int)
df['male'] = (df['gender'] == 1).astype(int)
df['age_cont'] = df['age']
df['edu_cont'] = df['edu']
df['unfriendly_yes'] = (df['unfriendly'] == 1).astype(int)
df['ethnic_central'] = (df['ethnic'] == 2).astype(int)
df['ethnic_cauc'] = (df['ethnic'] == 3).astype(int)
df['ethnic_other'] = (df['ethnic'] == 4).astype(int)
df['years_ru_cont'] = df['years_ru']

core_vars = ['geo_shock', 'disc_percep', 'pol_att', 'social_id', 'ru_prof',
             'years_ru_cont', 'female', 'age_cont', 'edu_cont', 'unfriendly_yes',
             'ethnic_central', 'ethnic_cauc']
df_clean = df.dropna(subset=core_vars).copy()
print(f"清洗后有效样本量：{len(df_clean)}")

# 协变量：中介与调节分析使用
covariates = ['female', 'age_cont', 'edu_cont', 'unfriendly_yes',
              'ethnic_central', 'ethnic_cauc']

# 协变量：族群作为调节变量时使用（移除族群虚拟变量，避免共线性）
covs_no_ethnic = ['female', 'age_cont', 'edu_cont', 'unfriendly_yes']

# ============================================================
# 4. 信度分析（Cronbach's α）
# ============================================================
def cronbach_alpha(data):
    return pg.cronbach_alpha(data=data)[0]

alpha_results = pd.DataFrame({
    '量表': ['地缘政治冲击感知', '歧视感知', '政治态度', '社会认同'],
    'Cronbach_alpha': [
        cronbach_alpha(df_clean[geo_items]),
        cronbach_alpha(df_clean[disc_items]),
        cronbach_alpha(df_clean[pol_items]),
        cronbach_alpha(df_clean[id_items])
    ],
    '条目数': [len(geo_items), len(disc_items), len(pol_items), len(id_items)]
})
print("\n========== 信度分析 ==========")
print(alpha_results)

# ============================================================
# 5. 描述性统计与相关分析
# ============================================================
desc_vars = ['geo_shock', 'disc_percep', 'pol_att', 'social_id', 'ru_prof', 'years_ru_cont']
desc_table = df_clean[desc_vars].describe().T
desc_table['偏度'] = df_clean[desc_vars].skew()
print("\n========== 描述性统计 ==========")
print(desc_table)

corr_vars = ['geo_shock', 'disc_percep', 'pol_att', 'social_id', 'ru_prof', 'years_ru_cont']
corr_matrix = df_clean[corr_vars].corr()
print("\n========== 相关矩阵 ==========")
print(corr_matrix.round(3))

r, p = stats.pearsonr(df_clean['geo_shock'], df_clean['disc_percep'])
print(f"geo_shock 与 disc_percep：r = {r:.3f}, p = {p:.4f}")

# ============================================================
# 6. 中介效应分析（Model 4）
# ============================================================
def mediation_bootstrap(df, x, m, y, covs, n_bootstrap=N_BOOTSTRAP, alpha=ALPHA):
    """
    简单中介分析（Model 4）
    路径 a：X -> M；路径 b：M -> Y；路径 c'：X -> Y（控制 M）
    返回：间接效应、直接效应、总效应及其 Bootstrap 置信区间
    """
    np.random.seed(42)
    boot_ie, boot_de, boot_te = [], [], []

    a = OLS(df[m], add_constant(df[[x] + covs])).fit().params[x]

    model_b = OLS(df[y], add_constant(df[[x, m] + covs])).fit()
    b = model_b.params[m]
    c_prime = model_b.params[x]

    c = OLS(df[y], add_constant(df[[x] + covs])).fit().params[x]
    ie = a * b

    for _ in range(n_bootstrap):
        sample = df.sample(n=len(df), replace=True)
        try:
            a_boot = OLS(sample[m], add_constant(sample[[x] + covs])).fit().params[x]
            model_b_boot = OLS(sample[y], add_constant(sample[[x, m] + covs])).fit()
            b_boot = model_b_boot.params[m]
            c_prime_boot = model_b_boot.params[x]
            c_boot = OLS(sample[y], add_constant(sample[[x] + covs])).fit().params[x]
            boot_ie.append(a_boot * b_boot)
            boot_de.append(c_prime_boot)
            boot_te.append(c_boot)
        except Exception:
            continue

    ie_ci = np.percentile(boot_ie, [100 * alpha / 2, 100 * (1 - alpha / 2)])
    de_ci = np.percentile(boot_de, [100 * alpha / 2, 100 * (1 - alpha / 2)])
    te_ci = np.percentile(boot_te, [100 * alpha / 2, 100 * (1 - alpha / 2)])

    return {
        'indirect_effect': ie,
        'direct_effect': c_prime,
        'total_effect': c,
        'a': a, 'b': b,
        'ie_ci': ie_ci,
        'de_ci': de_ci,
        'te_ci': te_ci
    }

print("\n========== 中介分析：因变量 = 政治态度 ==========")
med_pol = mediation_bootstrap(df_clean, 'geo_shock', 'disc_percep', 'pol_att', covariates)
print(f"间接效应 = {med_pol['indirect_effect']:.3f}, 95% CI [{med_pol['ie_ci'][0]:.3f}, {med_pol['ie_ci'][1]:.3f}]")
print(f"直接效应 = {med_pol['direct_effect']:.3f}, 95% CI [{med_pol['de_ci'][0]:.3f}, {med_pol['de_ci'][1]:.3f}]")
print(f"总效应   = {med_pol['total_effect']:.3f}, 95% CI [{med_pol['te_ci'][0]:.3f}, {med_pol['te_ci'][1]:.3f}]")

print("\n========== 中介分析：因变量 = 社会认同 ==========")
med_id = mediation_bootstrap(df_clean, 'geo_shock', 'disc_percep', 'social_id', covariates)
print(f"间接效应 = {med_id['indirect_effect']:.3f}, 95% CI [{med_id['ie_ci'][0]:.3f}, {med_id['ie_ci'][1]:.3f}]")
print(f"直接效应 = {med_id['direct_effect']:.3f}, 95% CI [{med_id['de_ci'][0]:.3f}, {med_id['de_ci'][1]:.3f}]")
print(f"总效应   = {med_id['total_effect']:.3f}, 95% CI [{med_id['te_ci'][0]:.3f}, {med_id['te_ci'][1]:.3f}]")

# ============================================================
# 7. 调节效应分析（Model 1）
# ============================================================
def moderation_analysis(df, x, y, moderator, covs):
    """调节效应：y ~ x * moderator + covs"""
    df_temp = df.copy()
    df_temp['inter'] = df_temp[x] * df_temp[moderator]
    X = add_constant(df_temp[[x, moderator, 'inter'] + covs])
    model = OLS(df_temp[y], X).fit()
    return model

print("\n========== 调节效应：因变量 = 政治态度 ==========")

# 7.1 俄语流利度
mod_ru = moderation_analysis(df_clean, 'geo_shock', 'pol_att', 'ru_prof', covariates)
print("\n调节变量：俄语流利度")
print(mod_ru.summary().tables[1])

# 7.2 在俄时长
mod_years = moderation_analysis(df_clean, 'geo_shock', 'pol_att', 'years_ru_cont', covariates)
print("\n调节变量：在俄时长")
print(mod_years.summary().tables[1])

# 7.3 族群身份（中亚裔 vs 俄族），协变量移除族群虚拟变量
mod_eth = moderation_analysis(df_clean, 'geo_shock', 'pol_att', 'ethnic_central', covs_no_ethnic)
print("\n调节变量：族群身份（中亚裔 vs 俄族）")
print(mod_eth.summary().tables[1])

# ============================================================
# 8. 分组中介分析
# ============================================================
# 说明：本研究的有调节的中介采用分组中介法，
# 即在调节变量的不同水平下分别估计间接效应，
# 通过比较各组的间接效应判断中介路径是否受调节。

# 8.1 按俄语流利度高低分组（均值 ± 0.5 SD）
ru_mean = df_clean['ru_prof'].mean()
ru_sd = df_clean['ru_prof'].std()
low_group = df_clean[df_clean['ru_prof'] < ru_mean - 0.5 * ru_sd]
high_group = df_clean[df_clean['ru_prof'] > ru_mean + 0.5 * ru_sd]

print("\n========== 分组中介：俄语流利度 ==========")
for label, sub in zip(['低流利度', '高流利度'], [low_group, high_group]):
    if len(sub) > 20:
        res = mediation_bootstrap(sub, 'geo_shock', 'disc_percep', 'pol_att',
                                  covariates, n_bootstrap=N_BOOTSTRAP_GROUP)
        print(f"{label} (n={len(sub)})：间接效应 = {res['indirect_effect']:.3f}, "
              f"95% CI [{res['ie_ci'][0]:.3f}, {res['ie_ci'][1]:.3f}]")

# 8.2 按族群分组（中亚裔、高加索裔、俄族）
# 注意：族群分组时协变量须移除族群虚拟变量，避免组内完全共线性
print("\n========== 分组中介：族群 ==========")
for ethnic_code, label in [(2, '中亚裔'), (3, '高加索裔'), (1, '俄罗斯族')]:
    sub = df_clean[df_clean['ethnic'] == ethnic_code]
    if len(sub) > 20:
        res = mediation_bootstrap(sub, 'geo_shock', 'disc_percep', 'pol_att',
                                  covs_no_ethnic, n_bootstrap=N_BOOTSTRAP_GROUP)
        print(f"{label} (n={len(sub)})：间接效应 = {res['indirect_effect']:.3f}, "
              f"95% CI [{res['ie_ci'][0]:.3f}, {res['ie_ci'][1]:.3f}]")
    else:
        print(f"{label} 样本量不足（n={len(sub)}），跳过")

# ============================================================
# 9. 结果导出到 Excel
# ============================================================
with pd.ExcelWriter(OUTPUT_PATH, engine='openpyxl') as writer:
    alpha_results.to_excel(writer, sheet_name='信度分析', index=False)
    desc_table.to_excel(writer, sheet_name='描述统计')
    corr_matrix.to_excel(writer, sheet_name='相关矩阵')

    med_results = pd.DataFrame({
        '因变量': ['政治态度', '政治态度', '政治态度', '社会认同', '社会认同', '社会认同'],
        '效应类型': ['间接效应', '直接效应', '总效应', '间接效应', '直接效应', '总效应'],
        '效应值': [med_pol['indirect_effect'], med_pol['direct_effect'], med_pol['total_effect'],
                  med_id['indirect_effect'], med_id['direct_effect'], med_id['total_effect']],
        '95%CI下限': [med_pol['ie_ci'][0], med_pol['de_ci'][0], med_pol['te_ci'][0],
                    med_id['ie_ci'][0], med_id['de_ci'][0], med_id['te_ci'][0]],
        '95%CI上限': [med_pol['ie_ci'][1], med_pol['de_ci'][1], med_pol['te_ci'][1],
                    med_id['ie_ci'][1], med_id['de_ci'][1], med_id['te_ci'][1]]
    })
    med_results.to_excel(writer, sheet_name='中介效应', index=False)

    mod_rows = []
    for name, model in [('俄语流利度', mod_ru), ('在俄时长', mod_years), ('族群(中亚裔)', mod_eth)]:
        inter_name = [n for n in model.params.index if 'inter' in n][0]
        mod_rows.append({
            '调节变量': name,
            '交互项系数': model.params[inter_name],
            '标准误': model.bse[inter_name],
            't值': model.tvalues[inter_name],
            'p值': model.pvalues[inter_name],
            '95%CI下限': model.conf_int().loc[inter_name, 0],
            '95%CI上限': model.conf_int().loc[inter_name, 1]
        })
    pd.DataFrame(mod_rows).to_excel(writer, sheet_name='调节效应', index=False)

print(f"\n分析完成，结果已保存至：{OUTPUT_PATH}")