导入数据¶
In [4]:
#导入15年数据
import os
import pandas as pd
import glob
import numpy as np
#导入15年的CHARLS调查数据
# 获取当前工作目录下的所有dta文件的路径
dta_files = glob.glob('D:\\论文\\数据\\2015年全国追踪调查\\*.dta')
# 创建一个空字典，用于存放dataframe对象和名称
dataframes = {}
# 遍历每个dta文件，读取数据并存入字典
for file in dta_files:
# 去掉文件名后缀，作为dataframe的名称，并在后面加上15，表示这是15年的
name = os.path.basename(file.split('.')[0])+ '15'
# 读取dta文件，返回一个dataframe对象
df = pd.read_stata(file)
# 将dataframe对象和名称存入字典
dataframes[name] = df
# 使用globals()函数将字典中的键和值分别赋值给全局变量和dataframe对象
globals().update(dataframes)
In [5]:
#导入18年数据
import os
import pandas as pd
import glob
import numpy as np#用同样的方法导入18年的CHARLS数据
# 获取当前工作目录下的所有dta文件的路径
dta_filess = glob.glob('D:\\论文\\数据\\2018年全国追踪调查\\*.dta')
# 创建一个空字典，用于存放dataframe对象和名称
dataframess = {}
# 遍历每个dta文件，读取数据并存入字典
for file in dta_filess:
# 去掉文件名后缀，作为dataframe的名称，并在后面加上15，表示这是15年的
name = os.path.basename(file.split('.')[0])+ '18'
# 读取dta文件，返回一个dataframe对象
dff = pd.read_stata(file)
# 将dataframe对象和名称存入字典
dataframess[name] = dff
# 使用globals()函数将字典中的键和值分别赋值给全局变量和dataframe对象
globals().update(dataframess)
导入消费和医疗数据¶
In [62]:
# 导入15年的 消费c0 以及 医疗支出数据m0
# 拿出消费包含的那些子项目
cols = ["ge006", "ge007", "ge008", "ge009_1", "ge009_2", "ge009_3", "ge009_4", "ge009_5", "ge009_6", "ge009_7", "ge010_1", "ge010_2", "ge010_3", "ge010_4", "ge010_5", "ge010_7", "ge010_8", "ge010_9", "ge010_10", "ge010_11", "ge010_12", "ge010_13"]
# 把空值都用0给填上嘻嘻嘻
Household_Income15[cols] = Household_Income15[cols].fillna(0)
# 创建新列c0喽，表示2015年的消费集合
Household_Income15["c0"] = (Household_Income15["ge006"] + Household_Income15["ge007"] + Household_Income15["ge008"]) * 52 + (Household_Income15["ge009_1"] + Household_Income15["ge009_2"] + Household_Income15["ge009_3"] + Household_Income15["ge009_4"] + Household_Income15["ge009_5"] + Household_Income15["ge009_6"] + Household_Income15["ge009_7"])*12 + Household_Income15["ge010_1"] + Household_Income15["ge010_2"] + Household_Income15["ge010_3"] + Household_Income15["ge010_4"] + Household_Income15["ge010_5"] + Household_Income15["ge010_7"] + Household_Income15["ge010_8"] + Household_Income15["ge010_9"] + Household_Income15["ge010_10"] + Household_Income15["ge010_11"] + Household_Income15["ge010_12"] + Household_Income15["ge010_13"]
# 把这些数据都整理到一个新的表c0m0里，其中m0是ge010_6这一列重命名的
c0m0 = Household_Income15 [ ["ID", "householdID", "communityID","c0", "ge010_6"]].copy()
c0m0.rename(columns={'ge010_6': 'm0'}, inplace=True)
# 使用平均值填充 c0列中的空值和0的行
c0m0['c0'] = c0m0['c0'].replace(0, np.nan)
mean_c0 = c0m0['c0'].mean()
c0m0['c0'] = c0m0['c0'].fillna(mean_c0)
#使用0替换m0中的空值
c0m0['m0'] = c0m0['m0'].fillna(0)
# 计算c0列的1%分位数和99%分位数
q1_15 = c0m0["c0"].quantile(0.01) # 计算hc15列的1%分位数
q99_15 = c0m0["c0"].quantile(0.99) # 计算hc15列的99%分位数
# 过滤掉前1%和后1%的行
c0m0_filtered = c0m0[(c0m0["c0"] > q1_15) & (c0m0["c0"] < q99_15)] # 保留c0列在(q1, q99)区间内的行
# 把c0m0_filtered改名为c0m0，也就是用过滤后的dataframe替换原来的dataframe
c0m0 = c0m0_filtered # 把c0m0_filtered赋值给c0
del c0m0_filtered # 删除c0m0_filtered变量，释放内存空间
c0m0
Out[62]:
| ID | householdID | communityID | c0 | m0 | |
|---|---|---|---|---|---|
| 0 | 094004126002 | 0940041260 | 0940041 | 6566.960000 | 0.0 |
| 1 | 094004103001 | 0940041030 | 0940041 | 2713.347277 | 0.0 |
| 2 | 094004111001 | 0940041110 | 0940041 | 3286.800000 | 20000.0 |
| 3 | 094004110001 | 0940041100 | 0940041 | 452.184000 | 170.0 |
| 4 | 094004108001 | 0940041080 | 0940041 | 195.880000 | 0.0 |
| ... | ... | ... | ... | ... | ... |
| 12056 | 294099314002 | 2940993140 | 2940993 | 13578.800000 | 0.0 |
| 12057 | 294099320001 | 2940993200 | 2940993 | 437.244000 | 0.0 |
| 12058 | 294099304001 | 2940993040 | 2940993 | 1574.510000 | 25000.0 |
| 12059 | 294099303001 | 2940993030 | 2940993 | 1884.100000 | 10000.0 |
| 12060 | 294099311002 | 2940993110 | 2940993 | 2392.989600 | 8000.0 |
11817 rows × 5 columns
计算健康的熵权¶
In [66]:
import pandas as pd
import numpy as np
import re
# ================= 工具函数 =================
def lead_int(x):
"""提取字符串中的首个整数；无则 NaN"""
if pd.isna(x): return np.nan
m = re.search(r'-?\d+', str(x))
return int(m.group()) if m else np.nan
def df_elemwise_map(df_, func):
"""兼容老版本 pandas 的 DataFrame 元素级映射"""
return df_.apply(lambda s: s.map(func))
def to_float(s):
"""稳健转为 float（包含 categorical/混合类型）"""
return pd.to_numeric(s, errors='coerce').astype(float)
def yn_to01(x): # 1 Yes -> 1, 2 No/空 -> 0
v = lead_int(x)
if np.isnan(v): v = 2
return 1 if v == 1 else 0
def likert_1_4_or_mean(x): # 抑郁：1–4，8/9/空 -> 2.5
v = lead_int(x)
if np.isnan(v) or v in (8, 9):
return 2.5
return float(v)
def is_correct_conservative(val, correct): # 认知：空白=错误
if pd.isna(val) or str(val).strip()=='':
return 0
v = lead_int(val)
if np.isnan(v):
return 0
return 1 if v == correct else 0
# ================= 1) 2015 指标（固定量程 & 保守缺失） =================
df15 = Health_Status_and_Functioning15.copy()
id15 = 'ID' if 'ID' in df15.columns else df15.columns[0]
# 自评健康：空=3(Fair)，映射到0..4，再/4 → [0,1]
srh15_raw = df15['da002'].apply(lead_int).fillna(3)
srh15_map = srh15_raw.map({1:4, 2:3, 3:2, 4:1, 5:0})
srh15 = (to_float(srh15_map).fillna(2.0)) / 4.0 # 显式转 float 后再除
# 残疾 5项：空=No → 0；健康分=1 - count/5
da005_15 = [c for c in [f'da005_{i}_' for i in range(1,6)] if c in df15.columns]
if da005_15:
dis15_cnt = df_elemwise_map(df15[da005_15], yn_to01).sum(axis=1)
else:
dis15_cnt = pd.Series(0.0, index=df15.index)
dis15_cnt = to_float(dis15_cnt)
dis15_health = 1 - dis15_cnt/5.0
# 慢性病 14项（w2 命名），非危重: 1,2,3,6,10,13,14；危重: 4,5,7,8,9,12（#11 不计入）
da007_15_all = [c for c in df15.columns if c.startswith('da007_w2_2_')]
def id_from_15(c):
m = re.search(r'da007_w2_2_(\d+)_?$', c)
return int(m.group(1)) if m else None
col2id15 = {c: id_from_15(c) for c in da007_15_all}
nonsev_ids = {1,2,3,6,10,13,14}
sev_ids = {4,5,7,8,9,12}
non15 = [c for c,i in col2id15.items() if i in nonsev_ids]
sev15 = [c for c,i in col2id15.items() if i in sev_ids]
if non15:
nonsev15_cnt = df_elemwise_map(df15[non15], yn_to01).sum(axis=1)
else:
nonsev15_cnt = pd.Series(0.0, index=df15.index)
if sev15:
severe15_cnt = df_elemwise_map(df15[sev15], yn_to01).sum(axis=1)
else:
severe15_cnt = pd.Series(0.0, index=df15.index)
nonsev15_cnt = to_float(nonsev15_cnt)
severe15_cnt = to_float(severe15_cnt)
nonsev15_health = 1 - nonsev15_cnt/7.0
severe15_health = 1 - severe15_cnt/6.0
# 认知 9题：空白=错；健康分=正确率
cog15_key = {
'dc031_w3_1':9, 'dc031_w3_2':10, 'dc031_w3_3':4,
'dc033_w3_1':4, 'dc033_w3_2':13, 'dc033_w3_3':2,
'dc034_w3_1':9, 'dc034_w3_2':5, 'dc034_w3_3':15,
}
cog15_cols = [c for c in cog15_key if c in df15.columns]
if cog15_cols:
correct15 = pd.DataFrame({c: df15[c].apply(is_correct_conservative, correct=cog15_key[c]) for c in cog15_cols})
cog15_health = to_float(correct15.sum(axis=1)) / float(len(cog15_cols))
else:
cog15_health = pd.Series(0.0, index=df15.index)
# 抑郁 10题：8/9/空=2.5；dc013/dc016 反向；题均分 → 健康分 (4-mean)/3
dep_cols = [f'dc{str(i).zfill(3)}' for i in range(9,19)]
dep15_cols = [c for c in dep_cols if c in df15.columns]
if dep15_cols:
dep15_num = df_elemwise_map(df15[dep15_cols], likert_1_4_or_mean)
for pos in ('dc013','dc016'):
if pos in dep15_num.columns:
dep15_num[pos] = 5 - dep15_num[pos]
dep15_mean = to_float(dep15_num.mean(axis=1))
dep15_health = (4.0 - dep15_mean)/3.0
else:
dep15_health = pd.Series(0.0, index=df15.index)
X15 = pd.DataFrame({
'srh' : srh15,
'disability': dis15_health,
'nonsevere' : nonsev15_health,
'severe' : severe15_health,
'cognition' : cog15_health,
'depression': dep15_health,
}, index=df15.index).clip(0,1)
# ================= 2) 2018 指标（固定量程 & 保守缺失） =================
df18 = Health_Status_and_Functioning18.copy()
cog18 = Cognition18.copy()
id18 = 'ID' if 'ID' in df18.columns else df18.columns[0]
need_cog = ['dc029_w3_1','dc030_w3_1','dc031_w3_1','dc031_w3_2','dc031_w3_3']
need_dep = [f'dc{str(i).zfill(3)}' for i in range(9,19)]
keep_cog = [id18] + [c for c in need_cog if c in cog18.columns] + [c for c in need_dep if c in cog18.columns]
cog18_sub = cog18[keep_cog].copy()
df18m = df18.merge(cog18_sub, on=id18, how='left')
# 自评健康：同口径
srh18_raw = df18m['da002'].apply(lead_int).fillna(3)
srh18_map = srh18_raw.map({1:4, 2:3, 3:2, 4:1, 5:0})
srh18 = (to_float(srh18_map).fillna(2.0)) / 4.0
# 残疾 5项
da005_18 = [c for c in [f'da005_{i}_' for i in range(1,6)] if c in df18m.columns]
if da005_18:
dis18_cnt = df_elemwise_map(df18m[da005_18], yn_to01).sum(axis=1)
else:
dis18_cnt = pd.Series(0.0, index=df18m.index)
dis18_cnt = to_float(dis18_cnt)
dis18_health = 1 - dis18_cnt/5.0
# 慢性病 14项（2018 命名 da007_i_）
da007_18_all = [c for c in [f'da007_{i}_' for i in range(1,15)] if c in df18m.columns]
def id_from_18(c):
m = re.match(r'^da007_(\d+)_$', c)
return int(m.group(1)) if m else None
col2id18 = {c: id_from_18(c) for c in da007_18_all}
non18 = [c for c,i in col2id18.items() if i in nonsev_ids]
sev18 = [c for c,i in col2id18.items() if i in sev_ids]
if non18:
nonsev18_cnt = df_elemwise_map(df18m[non18], yn_to01).sum(axis=1)
else:
nonsev18_cnt = pd.Series(0.0, index=df18m.index)
if sev18:
severe18_cnt = df_elemwise_map(df18m[sev18], yn_to01).sum(axis=1)
else:
severe18_cnt = pd.Series(0.0, index=df18m.index)
nonsev18_cnt = to_float(nonsev18_cnt)
severe18_cnt = to_float(severe18_cnt)
nonsev18_health = 1 - nonsev18_cnt/7.0
severe18_health = 1 - severe18_cnt/6.0
# 认知 5题 → 正确率（空白=错）
cog18_key = {'dc029_w3_1':6,'dc030_w3_1':5,'dc031_w3_1':9,'dc031_w3_2':10,'dc031_w3_3':4}
cog18_cols = [c for c in cog18_key if c in df18m.columns]
if cog18_cols:
correct18 = pd.DataFrame({c: df18m[c].apply(is_correct_conservative, correct=cog18_key[c]) for c in cog18_cols})
cog18_health = to_float(correct18.sum(axis=1)) / float(len(cog18_cols))
else:
cog18_health = pd.Series(0.0, index=df18m.index)
# 抑郁 10题：题均分→健康分
dep18_cols = [c for c in need_dep if c in df18m.columns]
if dep18_cols:
dep18_num = df_elemwise_map(df18m[dep18_cols], likert_1_4_or_mean)
for pos in ('dc013','dc016'):
if pos in dep18_num.columns:
dep18_num[pos] = 5 - dep18_num[pos]
dep18_mean = to_float(dep18_num.mean(axis=1))
dep18_health = (4.0 - dep18_mean)/3.0
else:
dep18_health = pd.Series(0.0, index=df18m.index)
X18 = pd.DataFrame({
'srh' : srh18,
'disability': dis18_health,
'nonsevere' : nonsev18_health,
'severe' : severe18_health,
'cognition' : cog18_health,
'depression': dep18_health,
}, index=df18m.index).clip(0,1)
# ================= 3) 合并计算统一熵权（pooled） =================
X_pool = pd.concat([X15, X18], axis=0, ignore_index=True)
# 熵权（列归一、熵、差异系数、权重）
eps = 1e-12
n = max(2, len(X_pool)) # 防止 log(1)=0
col_sums = X_pool.sum(axis=0)
P = X_pool.div(col_sums.replace(0, np.nan), axis=1).fillna(1.0/n)
E = -(P.replace(0, eps) * np.log(P.replace(0, eps))).sum(axis=0) / np.log(n)
D = 1 - E
W = (D / D.sum()).replace([np.inf, -np.inf], 0).fillna(0)
W.name = 'weight'
print("Unified (pooled) entropy weights:")
print(pd.DataFrame({'indicator': W.index, 'weight': W.values}).to_string(index=False))
# ================= 4) 生成新指数并评估改善比例 =================
hsf15_new = (X15 @ W) * 100
hsf18_new = (X18 @ W) * 100
hsf15 = pd.DataFrame({ 'ID': df15[id15],'householdID': df15['householdID'], 'hsf15': hsf15_new })
hsf18 = pd.DataFrame({ 'ID': df18m[id18], 'householdID': df18m['householdID'],'hsf18': hsf18_new })
print("\n2015（前10行）:")
print(hsf15.head(10).to_string(index=False))
print("\n2018（前10行）:")
print(hsf18.head(10).to_string(index=False))
Unified (pooled) entropy weights:
indicator weight
srh 0.177299
disability 0.009149
nonsevere 0.009182
severe 0.005538
cognition 0.705520
depression 0.093311
2015（前10行）:
ID householdID hsf15
094004126001 0940041260 15.917436
094004126002 0940041260 43.789299
094004103001 0940041030 43.571006
094004111002 0940041110 15.917436
094004111001 0940041110 36.868519
094004110001 0940041100 41.923078
094004108001 0940041080 23.756552
094004112001 0940041120 19.518710
094004112002 0940041120 15.917436
094004106001 0940041060 42.413983
2018（前10行）:
ID householdID hsf18
094004113002 0940041130 71.648015
094004111002 0940041110 83.145229
094004111001 0940041110 69.389463
094004112001 0940041120 27.581738
094004118001 0940041180 79.470728
094004118002 0940041180 70.103841
094004121001 0940041210 15.917436
094004119002 0940041190 95.436347
094004122001 0940041220 75.649306
094004122002 0940041220 80.392828
计算健康的边际效用phi¶
In [33]:
#找到15年效用
import pandas as pd
# 提取 ID 和 dc028 列
utility_selected15 = Health_Status_and_Functioning15[['ID', 'dc028']].copy()
# 删除 dc028 为空的行
utility_selected15 = utility_selected15.dropna(subset=['dc028'])
# 仅保留 dc028 列的第一个字符
utility_selected15['dc028'] = utility_selected15['dc028'].astype(str).str[0]
# 将 dc028 转换为数值类型（例如整数）
utility_selected15['dc028'] = pd.to_numeric(utility_selected15['dc028'], errors='coerce')
# 重命名 dc028 列为 utility
utility_selected15.rename(columns={'dc028': 'utility15'}, inplace=True)
# 最终结果存储在 df_selected 中
print(utility_selected15)
ID utility15 1 094004126002 3 2 094004103001 2 4 094004111001 3 5 094004110001 3 7 094004112001 2 ... ... ... 20962 294099304001 2 20963 294099303002 2 20964 294099303001 3 20965 294099311001 2 20966 294099311002 3 [19468 rows x 2 columns]
In [5]:
#找到18年效用
import pandas as pd
# 提取 ID 和 dc028 列
utility_selected18 = Cognition18[['ID', 'dc028']].copy()
# 删除 dc028 为空的行
utility_selected18 = utility_selected18.dropna(subset=['dc028'])
# 仅保留 dc028 列的第一个字符
utility_selected18['dc028'] = utility_selected18['dc028'].astype(str).str[0]
# 将 dc028 转换为数值类型（例如整数）
utility_selected18['dc028'] = pd.to_numeric(utility_selected18['dc028'], errors='coerce')
# 重命名 dc028 列为 utility
utility_selected18.rename(columns={'dc028': 'utility18'}, inplace=True)
# 最终结果存储在 df_selected 中
print(utility_selected18)
ID utility18 0 094004113002 3 1 094004111002 3 2 094004111001 3 3 094004112001 1 4 094004118001 4 ... ... ... 19739 294099304001 2 19740 294099303002 3 19741 294099303001 2 19742 294099311001 1 19743 294099311002 3 [18152 rows x 2 columns]
In [37]:
z151 = pd.merge(utility_selected15,c0m0,on="ID",how="inner")
z152 = pd.merge(z151,hsf15,on="ID",how="inner")
phi15 = z152[ ["ID", "utility15","c0","hsf15"]].copy()
phi15
Out[37]:
| ID | utility15 | c0 | hsf15 | |
|---|---|---|---|---|
| 0 | 094004126002 | 3 | 79120.000000 | 81.448436 |
| 1 | 094004103001 | 2 | 32690.931043 | 77.549631 |
| 2 | 094004111001 | 3 | 39600.000000 | 65.595125 |
| 3 | 094004110001 | 3 | 5448.000000 | 76.186413 |
| 4 | 094004112001 | 2 | 15700.000000 | 79.362932 |
| ... | ... | ... | ... | ... |
| 11293 | 294099314002 | 2 | 163600.000000 | 72.156289 |
| 11294 | 294099320001 | 3 | 5268.000000 | 64.829910 |
| 11295 | 294099304001 | 2 | 18970.000000 | 45.094159 |
| 11296 | 294099303001 | 3 | 22700.000000 | 56.367055 |
| 11297 | 294099311002 | 3 | 28831.200000 | 46.321275 |
11298 rows × 4 columns
In [43]:
# =========================
# 估计健康的边际效用（phi）
# - u: 生活满意度（1=最好, 5=最差）→ 反向为 u_good（越大=越满意/效用越高）
# - h: hsf15（已确认越大=越健康）→ 直接使用 h_good = hsf15
# - c: c0 → x(c) = (c^(1-sigma))/(1-sigma)，此处 sigma=3
# 主结果：有序 Probit；稳健性：OLS(HC1)
# =========================
import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.miscmodels.ordinal_model import OrderedModel
from scipy.stats import spearmanr
# ---------- 0) 参数 ----------
SIGMA = 3
COL_U = "utility15" # 1(最好)~5(最差) 的满意度等级
COL_C = "c0"
COL_H = "hsf15"
# ---------- 1) 取数与清理 ----------
df = phi15[[COL_U, COL_C, COL_H]].copy()
# 严格保留三列都非空的观测
df = df.dropna(subset=[COL_U, COL_C, COL_H]).reset_index(drop=True)
# ---------- 2) 构造变量 ----------
# 2.1 效用反向：u_good 越大=越满意（更高的效用）
df["u_good"] = 6 - df[COL_U].astype(int)
# 2.2 健康：已确认 hsf15 越大=越健康
df["h_good"] = df[COL_H].astype(float)
# 2.3 经济项：x(c) = (c^(1-sigma))/(1-sigma)；此处 sigma=3 → -0.5 * c^{-2}
df["x"] = (df[COL_C].astype(float) ** (1 - SIGMA)) / (1 - SIGMA)
# （可选）标准化，有助于数值稳定与系数可比性（不影响有序模型识别方向）
# 注：你可以把下两行注释/启用来比较标准化与否的差异
# df["x_z"] = (df["x"] - df["x"].mean()) / df["x"].std(ddof=0)
# df["h_z"] = (df["h_good"] - df["h_good"].mean()) / df["h_good"].std(ddof=0)
# ---------- 3) 快速体检（诊断，不影响估计） ----------
print("样本量：", len(df))
print("u_good 描述：\n", df["u_good"].describe())
print("h_good 描述：\n", df["h_good"].describe())
print("x(c) 描述：\n", df["x"].describe())
# u_good 与 h_good 的秩相关（只看排序，方向应为正更符合直觉）
rho_uh = spearmanr(df["u_good"], df["h_good"]).correlation
print("\nSpearman rho(u_good, h_good) =", rho_uh)
# ---------- 4) 有序 Probit（主结果） ----------
# 注：OrderedModel 不要手动加常数项，阈值会吸收截距
X_ord = df[["x", "h_good"]]
y_ord = df["u_good"].astype(int)
mod = OrderedModel(endog=y_ord, exog=X_ord, distr="probit")
res = mod.fit(method="bfgs", disp=False)
print("\n========== Ordered Probit (主结果) ==========")
print(res.summary())
phi_hat_ordered = float(res.params["h_good"])
print("\nphi_hat (Ordered Probit, latent scale) =", phi_hat_ordered)
# 边际效应：h_good 对各等级概率的边际影响（单位：h_good 增加 1 的效应）
# 注：若报 NotImplemented，可能是 statsmodels 版本差异，可忽略或升级
try:
me = res.get_margeff(at="overall", method="dydx", dummy=False)
print("\n------ 边际效应（dP/df）------")
print(me.summary())
except Exception as e:
print("\n[提示] 边际效应计算未提供/失败：", e)
# ---------- 5) OLS 稳健性（把 u_good 当作近似连续评分；HC1 异方差稳健） ----------
X_ols = sm.add_constant(df[["x", "h_good"]])
y_ols = df["u_good"].astype(float)
ols = sm.OLS(y_ols, X_ols).fit(cov_type="HC1")
print("\n========== OLS 稳健性 (HC1) ==========")
print(ols.summary())
phi_hat_ols = float(ols.params["h_good"])
print("\nphi_hat (OLS, HC1) =", phi_hat_ols)
# ---------- 6) 可选：仅基于排序的“秩回归”（极快，作进一步稳健性） ----------
# 把三列都转成分位秩（0-1），回归 rank(u) ~ rank(x) + rank(h)
try:
ru = df["u_good"].rank(pct=True)
rx = df["x"].rank(pct=True)
rh = df["h_good"].rank(pct=True)
X_rank = sm.add_constant(pd.DataFrame({"rx": rx, "rh": rh}))
rank_ols = sm.OLS(ru, X_rank).fit(cov_type="HC1")
print("\n========== Rank-based 回归（仅排序信息） ==========")
print(rank_ols.summary())
print("\nrank-based gamma on h =", float(rank_ols.params["rh"]))
except Exception as e:
print("\n[提示] 秩回归未执行：", e)
# ---------- 7) 结果小结 ----------
print("\n====== 小结 ======")
print(f"Ordered Probit: phi_hat (h_good) = {phi_hat_ordered:.6f} [潜在效用尺度]")
print(f"OLS (HC1): phi_hat (h_good) = {phi_hat_ols:.6f} [把 u_good 视作连续]")
print("方向应一致：h_good 系数 > 0 表示“健康更好 → 效用更高”。")
样本量： 11298
u_good 描述：
count 11298.000000
mean 3.352540
std 0.799599
min 1.000000
25% 3.000000
50% 3.000000
75% 4.000000
max 5.000000
Name: u_good, dtype: float64
h_good 描述：
count 11298.000000
mean 66.059231
std 12.906474
min 17.346899
25% 57.525920
50% 66.276734
75% 75.374650
max 100.000000
Name: h_good, dtype: float64
x(c) 描述：
count 1.129800e+04
mean -2.625157e-08
std 1.280439e-07
min -2.143264e-06
25% -6.716999e-09
50% -1.490095e-09
75% -4.392177e-10
max -7.466585e-12
Name: x, dtype: float64
Spearman rho(u_good, h_good) = 0.20007606676913314
========== Ordered Probit (主结果) ==========
OrderedModel Results
==============================================================================
Dep. Variable: u_good Log-Likelihood: -12929.
Model: OrderedModel AIC: 2.587e+04
Method: Maximum Likelihood BIC: 2.591e+04
Date: Wed, 08 Oct 2025
Time: 23:18:40
No. Observations: 11298
Df Residuals: 11292
Df Model: 2
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
x 2.139e-07 nan nan nan nan nan
h_good 0.0197 nan nan nan nan nan
1/2 -0.8176 nan nan nan nan nan
2/3 -0.2473 nan nan nan nan nan
3/4 0.4662 nan nan nan nan nan
4/5 0.2561 nan nan nan nan nan
==============================================================================
phi_hat (Ordered Probit, latent scale) = 0.019743378970743634
[提示] 边际效应计算未提供/失败： 'OrderedResults' object has no attribute 'get_margeff'
========== OLS 稳健性 (HC1) ==========
OLS Regression Results
==============================================================================
Dep. Variable: u_good R-squared: 0.054
Model: OLS Adj. R-squared: 0.054
Method: Least Squares F-statistic: 3.836
Date: Wed, 08 Oct 2025 Prob (F-statistic): 0.0502
Time: 23:18:40 Log-Likelihood: -13188.
No. Observations: 11298 AIC: 2.638e+04
Df Residuals: 11295 BIC: 2.640e+04
Df Model: 2
Covariance Type: HC1
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 2.4078 0.042 57.947 0.000 2.326 2.489
x 1.245e+05 6.36e+04 1.959 0.050 -86.065 2.49e+05
h_good 0.0144 0.001 23.942 0.000 0.013 0.016
==============================================================================
Omnibus: 62.128 Durbin-Watson: 1.935
Prob(Omnibus): 0.000 Jarque-Bera (JB): 83.703
Skew: -0.070 Prob(JB): 6.67e-19
Kurtosis: 3.398 Cond. No. 5.26e+08
==============================================================================
Notes:
[1] Standard Errors are heteroscedasticity robust (HC1)
[2] The condition number is large, 5.26e+08. This might indicate that there are
strong multicollinearity or other numerical problems.
phi_hat (OLS, HC1) = 0.014351112236961797
========== Rank-based 回归（仅排序信息） ==========
OLS Regression Results
==============================================================================
Dep. Variable: u_good R-squared: 0.041
Model: OLS Adj. R-squared: 0.041
Method: Least Squares F-statistic: 244.8
Date: Wed, 08 Oct 2025 Prob (F-statistic): 8.18e-105
Time: 23:18:40 Log-Likelihood: -749.38
No. Observations: 11298 AIC: 1505.
Df Residuals: 11295 BIC: 1527.
Df Model: 2
Covariance Type: HC1
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 0.3936 0.006 63.809 0.000 0.382 0.406
rx 0.0355 0.009 4.131 0.000 0.019 0.052
rh 0.1774 0.009 20.694 0.000 0.161 0.194
==============================================================================
Omnibus: 4653.553 Durbin-Watson: 1.926
Prob(Omnibus): 0.000 Jarque-Bera (JB): 625.063
Skew: 0.173 Prob(JB): 1.86e-136
Kurtosis: 1.901 Cond. No. 4.93
==============================================================================
Notes:
[1] Standard Errors are heteroscedasticity robust (HC1)
rank-based gamma on h = 0.17737240799015952
====== 小结 ======
Ordered Probit: phi_hat (h_good) = 0.019743 [潜在效用尺度]
OLS (HC1): phi_hat (h_good) = 0.014351 [把 u_good 视作连续]
方向应一致：h_good 系数 > 0 表示“健康更好 → 效用更高”。
D:\softwares\Anaconda\Lib\site-packages\statsmodels\base\model.py:595: HessianInversionWarning: Inverting hessian failed, no bse or cov_params available
warnings.warn('Inverting hessian failed, no bse or cov_params '
D:\softwares\Anaconda\Lib\site-packages\statsmodels\base\model.py:1894: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 2, but rank is 1
warnings.warn('covariance of constraints does not have full '