import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import matplotlib.gridspec as gridspec
%matplotlib inline
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
import statsmodels.formula.api as smf
import statsmodels.api as sm
import copy
import warnings
warnings.filterwarnings("ignore")
from sklearn.linear_model import LinearRegression
from scipy.stats import gaussian_kde
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error

#=======================导入数据============================
df = pd.read_excel('某副省级城市长护险试点数据.xlsx')
df = df[df.disable_date.apply(lambda x:len(str(x))>4)]#删除建档时间为空
df.disable_date = pd.to_datetime(df.disable_date)#建档时间提取为时间对象
df['year'] = df.disable_date.apply(lambda x:x.year)
df = df[df.level.notna()]#保留所有失能存在失能等级且建档的样本
df.level = df.level.apply(lambda x:x[-1])#提取失能等级的级数
df = df.drop_duplicates('id')#剔除ID相同的样本，仅保留首次建档
#----------------------------------------------------------
l = ['2017-12-31','2018-12-31','2019-12-31','2020-12-31'] #设定统计时间区间
l = [pd.to_datetime(i) for i in l ]#将统计时间设定为datetime对象
d1 = df[df.dtime.isna()]#非死亡样本
d2 = df[~df.dtime.isna()]#死亡样本
d2['dtime'] = d2.dtime.apply(lambda x:pd.to_datetime(x))
for i in range(4):#对四个年度区间的建档群体循环执行描述性统计
    dd = d1[d1['disable_date']<l[i]]
    dd = pd.concat([d2[(d2['disable_date']<l[i])&(d2['dtime']>=l[i])],dd])
    if i==0:
        res =dd.groupby('level').describe()['age']
    else:
        res2 =dd.groupby('level').describe()['age']
        res =pd.concat([res,res2])

        
        

#=======================导入数据============================
df = pd.read_excel('某副省级城市长护险试点数据.xlsx')
df = df[df.disable_date.apply(lambda x:len(str(x))>4)]#删除建档时间为空
df.disable_date = pd.to_datetime(df.disable_date)#建档时间提取为时间对象
df['year'] = df.disable_date.apply(lambda x:x.year)
df = df[df.level.notna()]#保留所有失能存在失能等级且建档的样本
df.level = df.level.apply(lambda x:x[-1])#提取失能等级的级数

df2 = pd.read_excel('某副省级城市城镇职工参保结构数据.xlsx')
df2['rate'] = df2['num']/(df2.num.sum())#计算不同年龄参保人数占比，变量num对应人数，变量rate对应比例
#=======================相关函数============================
def get_kde_from_agedata(df,year):#提取不同时间参保人数的kde核密度
    data = df[df['year']==year]['age']
    kde = gaussian_kde(data, bw_method='scott')  #
    x_values = np.linspace(data.min(), data.max(), 1000) 
    kde_values = kde(x_values)  # 计算 KDE 值
    return x_values,kde_values
#=======================可视化代码============================
fig = plt.figure(figsize=(20,8),dpi=600)
for i in range(4):#前四个子图的分布
    ax = fig.add_subplot(2,4,i+1)
    d = df[df.year==i+2017]
    sns.histplot(d.age, bins=40, kde=True, label='{}年的重度失能人口（新评定）'.format(i+2017), color=sns.color_palette()[i], alpha=0.6,stat='density')
    plt.legend(fontsize=10,loc='upper left')
    plt.grid()
ax = fig.add_subplot(2,1,2)
for i in range(4):#前四个子图的kde核密度曲线
    x,y = get_kde_from_agedata(df,2017+i)
    sns.lineplot(x=x-15,y=y)
    ax.fill_between(x-15, y, color=sns.color_palette()[i], alpha=0.5, label='{}年的重度失能人口分布'.format(i+2017))
plt.ylim(0,0.045)
plt.legend(fontsize=12)
ax = ax.twinx()
sns.barplot(x=df2['age'].astype('i'),y=df2['rate'],label='城镇职工基本医疗保险群体参保结构')
plt.xlim(0,100)
ax.set_xticks([i for i in range(0,100)][::5])  
ax.set_xticklabels([i for i in range(15,115)][::5],fontsize=7)
plt.legend(loc='upper left',fontsize=12)
plt.show()
#=======================输出结果============================


#=======================导入数据============================
df = pd.read_excel('某副省级城市长护险试点数据.xlsx')
df = df[df.disable_date.apply(lambda x:len(str(x))>4)]#删除建档时间为空
df.disable_date = pd.to_datetime(df.disable_date)#建档时间提取为时间对象
df['year'] = df.disable_date.apply(lambda x:x.year)
df = df[df.level.notna()]#保留所有失能存在失能等级且建档的样本
df.level = df.level.apply(lambda x:x[-1])#提取失能等级的级数
df2 = df.drop_duplicates('id')
df2.dtime = df2.dtime.apply(lambda x:x if pd.isna(x) else pd.to_datetime(x))
#=======================相关函数============================
def f(df,dtrg,i):#函数注释：计算在档总人数
    t = dtrg[i+1]
    dt = df[df['disable_date']<=t]
    dtime = dt.dtime
    dtime2 = dtime[~dtime.isna()]
    return len(dt)-(dtime2<t).sum()
def f_new(df,dtrg,i):#函数注释：计算新建档人数
    t0 = dtrg[i]
    t1 = dtrg[i+1]
    dt=  df.disable_date
    return ((dt>=t0)&(dt<t1)).sum()
def f_dd(df,dtrg,i):#函数注释：计算死亡人数
    t0 = dtrg[i]
    t1 = dtrg[i+1]
    dt=  df.dtime
    dt = dt[~dt.isna()]
    return ((dt>=t0)&(dt<t1)).sum()

dtrg = pd.date_range(start = pd.to_datetime('2015-12-31'),end = pd.to_datetime('2019-12-31'),freq='M') #按照时间区间生成月份列表
idx = [i for i in range(len(dtrg)-1)]#生成index列表
n = [f(df2,dtrg,i) for i in idx]
n_new = [f_new(df2,dtrg,i) for i in idx]
n_dd = [f_dd(df2,dtrg,i) for i in idx]
sr = pd.Series(dtrg)
sr2 = sr.apply(lambda x:str(x.year)+'-'+str(x.month))
#==================================================================
fig = plt.figure(figsize=(20,6),dpi=800)
ax = fig.add_subplot(2,2,1)#对应第一个子图
ax.bar([i for i in idx],n_new,edgecolor=sns.color_palette()[1],facecolor= sns.color_palette()[0])
ax.set_xticks([i for i in idx]) 
ax.set_xticklabels(sr2.tolist()[1:],fontsize=7)
ax.vlines(18.5,0,3000,colors='black',linestyle='--',lw=1)#对应长护险制度的启动时间
ax.hlines(500,0,48,colors='black',linestyle='--',lw=1)#500为实际月度新增的估计结果
plt.xticks(rotation=45)
plt.xlim(-1,49)
custom_legend = [Patch(edgecolor=sns.color_palette()[1],facecolor= sns.color_palette()[0], label='试点前后政府部门面临的失能群体增量')]
plt.legend(handles=custom_legend,loc='upper right',fontsize=12)
#————————————————————————————————————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(2,2,3)#对应第三个子图，均为理论推测数值
ax.bar([i for i in idx],n_new,edgecolor=sns.color_palette()[0],facecolor=sns.color_palette()[3])
ax.set_xticks([i for i in idx])  # 设置刻度标签位置
ax.set_xticklabels(sr2.tolist()[1:],fontsize=7)
ax.vlines(18.5,0,3000,colors='black',linestyle='--',lw=1)
plt.xticks(rotation=45)
a = 19
ax.bar([i for i in idx[:a]],[500 for i in idx[:a]],edgecolor=sns.color_palette()[0],facecolor=sns.color_palette()[3])
v = n_new[a:]
a2 = 32
ax.bar([i for i in idx[a:a2]],[500 for i in idx[a:a2]],edgecolor=sns.color_palette()[1],facecolor=sns.color_palette()[0])
ax.bar([i for i in idx[a2:]],[i for i in n_new[a2:]],edgecolor=sns.color_palette()[1],facecolor=sns.color_palette()[0])
# plt.title('试点前后实际发生的失能群体增量')
custom_legend = [Patch(edgecolor=sns.color_palette()[0],facecolor=sns.color_palette()[3], label='过去发生的失能新增（理论推测）'),
    Patch(edgecolor=sns.color_palette()[1],facecolor=sns.color_palette()[0], label='长护险试点后发生的失能新增（理论推测）')]
plt.legend(handles=custom_legend,loc='upper right',fontsize=12)
#===================================================================================================================================
dtrg2 = pd.date_range(start = pd.to_datetime('2017-8-31'),end = pd.to_datetime('2020-8-31'),freq='M')#按照时间区间生成月份列表
idx2 = [i for i in range(len(dtrg2)-1)]
n2 = [f(df2,dtrg2,i) for i in idx2]
n_new2 = [f_new(df2,dtrg2,i) for i in idx2]
n_dd2 = [f_dd(df2,dtrg2,i) for i in idx2]
sr = pd.Series(dtrg2)
sr3 = sr.apply(lambda x:str(x.year)+'-'+str(x.month))
#==================================================================================================
ax = fig.add_subplot(2,2,2)#对应第二个子图
ax.bar([i for i in idx2],n2,edgecolor=sns.color_palette()[1],facecolor= sns.color_palette()[2])
ax.set_xticks([i for i in idx2])  # 设置刻度标签位置
ax.set_xticklabels(sr3.tolist()[1:],fontsize=7)
plt.xticks(rotation=45)
plt.ylim(0,20000)
custom_legend = [Patch(edgecolor=sns.color_palette()[1],facecolor= sns.color_palette()[2], label='长护险试点后政策口径的失能群体总量')]
plt.legend(handles=custom_legend,loc='upper right',fontsize=12)
#==================================================================================================
ax = fig.add_subplot(2,2,4)#对应第四个子图
sns.lineplot(x=idx2,y=[n_dd2[i]/n2[i] for i in idx2],marker='o')
ax.set_xticks([i for i in idx2])  # 设置刻度标签位置
ax.set_xticklabels(sr3.tolist()[1:],fontsize=7)
plt.xticks(rotation=45)
custom_legend = [Patch(edgecolor=sns.color_palette()[0],facecolor= sns.color_palette()[0], label='长护险试点后失能群体的平均死亡率变化')]
plt.legend(handles=custom_legend,loc='upper right',fontsize=12)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.1, hspace=0.2)
plt.show()



#=======================导入数据============================
df = pd.read_excel('某副省级城市长护险试点数据.xlsx')
df = df[df.disable_date.apply(lambda x:len(str(x))>4)]#删除建档时间为空
df.disable_date = pd.to_datetime(df.disable_date)#建档时间提取为时间对象
df['year'] = df.disable_date.apply(lambda x:x.year)
df = df[df.level.notna()]#保留所有失能存在失能等级且建档的样本
df.level = df.level.apply(lambda x:x[-1])#提取失能等级的级数
df['birth_year'] = df['year'].astype('int')-df['age'].astype('int')#birth_year对应出生时间
#=========================================================
fig = plt.figure(figsize=(30,12),dpi=800)
df2 = df[df.level.isin(['1','2'])]
df2.level =df2.level.apply(lambda x:'重度失能1级' if x=='1' else '重度失能2级')
ax = fig.add_subplot(2,3,1)#第一个子图，不同年度失能群体的ADL评分分布（adl变量即ADL评分）
sns.boxplot(x='year',y='adl',hue='level',data=df2, showfliers=False)
plt.legend(fontsize=15,loc='upper right')
plt.xlabel('')
plt.ylabel('ADL评分',fontsize=15)
#=========================================================
ax = fig.add_subplot(2,3,2)#第二个子图，不同年度失能群体的年龄分布
sns.boxplot(x='year',y='age',hue='level',data=df2, showfliers=False)
plt.legend(fontsize=15,loc='upper right')
plt.xlabel('')
plt.ylabel('年龄',fontsize=15)
#=========================================================
ax = fig.add_subplot(2,3,3)#第三个子图，不同年度失能群体的出生时间分布
sns.boxplot(x='year',y='birth_year',hue='level',data=df2, showfliers=False)
plt.legend(fontsize=15,loc='upper right')
plt.xlabel('')
plt.ylabel('出生时间',fontsize=15)
#===================================================================
df2=df[df.year.isin([2020,2017])]
df2 = df2[df2.level.isin(['1','2'])]
ax = fig.add_subplot(2,3,4)#第四个子图，2017年与2020年度失能群体的ADL评分分布差异
sns.violinplot(data=df2, x="level", y="adl", hue="year",
               split=True, inner="quart", fill=False)
plt.xlabel('')
plt.ylabel('ADL评分',fontsize=15)
plt.legend(fontsize=15)
ax.set_xticklabels(['重度失能1级','重度失能2级'],fontsize=15)
#===================================================================
ax = fig.add_subplot(2,3,5)#第五个子图，2017年与2020年度失能群体的年龄分布差异
sns.violinplot(data=df2, x="level", y="age", hue="year",
               split=True, inner="quart", fill=False)
ax.set_xticklabels(['重度失能1级','重度失能2级'],fontsize=15)
plt.xlabel('')
plt.ylabel('年龄',fontsize=15)
plt.legend(fontsize=15)
#===================================================================
ax = fig.add_subplot(2,3,6)#第六个子图，2017年与2020年度失能群体的出生时间分布差异
sns.violinplot(data=df2, x="level", y="birth_year", hue="year",
               split=True, inner="quart", fill=False)
ax.set_xticklabels(['重度失能1级','重度失能2级'],fontsize=15)
plt.xlabel('')
plt.ylabel('出生时间',fontsize=15)
plt.legend(fontsize=15)


#=====================导入城市统计年鉴数据====================================
d = pd.read_excel('中国城市统计年鉴数据.xlsx')
l =['年份', '省份', '城市','人均地区生产总值(元)','医院、卫生院床位数(张)'] 
d0 = d[l]
d0.columns = ['year','c1','c2','pgdp','hos']
d0 = d0[d0.year==2019]
d0 = d0[['c2','pgdp','hos']]#pgdp：人均GDP数据，量化城市经济
#=====================导入城市统计年鉴数据====================================
df = pd.read_excel('第七次人口普查分县数据.xlsx',header=2)
l2 = ['Unnamed: 0','人口数（人）','65岁及以上人口比重（%）']
df0 = df[l2]
df0.columns = ['c2','p_num','old_rate']#old_rate：65岁以上老龄人口比重，量化人口结构
#=====================合并数据====================================
df0 = df0[df0.c2.isin(d0.c2)]
res = pd.merge(d0,df0)
res['phos'] = res['hos']/res['p_num']*1000 #phos：每千人床位数，量化医疗资源
#=====================合并数据====================================
fig = plt.figure(figsize=(30,8),dpi=600)
ax = fig.add_subplot(1,3,1)
sns.histplot(res.pgdp, bins=20, kde=True, label='人均GDP', color=sns.color_palette()[0], alpha=0.6,stat='density')
ax.axvline(x=103386, color='red', linestyle='--',linewidth=10, label='样本城市')
plt.xlabel('经济发展',fontsize=25)
plt.grid()
plt.legend(fontsize=20,loc='upper right')
#==============================================================
ax = fig.add_subplot(1,3,2)
sns.histplot(res.old_rate, bins=20, kde=True, label='65岁以上人口占比', color=sns.color_palette()[1], alpha=0.6,stat='density')
ax.axvline(x=13.62, color='red', linestyle='--',linewidth=10, label='样本城市')
plt.grid()
plt.xlabel('人口结构',fontsize=25)
plt.legend(fontsize=20,loc='upper left')
#==============================================================
ax = fig.add_subplot(1,3,3)
sns.histplot(res.phos, bins=20, kde=True, label='每千人口床位数', color=sns.color_palette()[2], alpha=0.6,stat='density')
ax.axvline(x=5.910948, color='red', linestyle='--',linewidth=10, label='样本城市')
plt.grid()
plt.legend(fontsize=20,loc='upper right')
plt.xlabel('医疗资源',fontsize=25)


gm = pd.read_excel('德国长护险数据.xlsx')
gm2 = pd.read_excel('德国老龄人口数据.xlsx')
jp = pd.read_excel('日本长护险数据.xlsx')
idx =[i for i in range(len(d))]
fig = plt.figure(figsize=(10,8),dpi=800)
ax = fig.add_subplot(2,1,1)
plt.bar(gm2.year,gm2['age2_65upper']/10,color=sns.color_palette()[0],label='德国老龄人口数量(65岁以上)')#age2_65upper：德国老龄人口数量(65岁以上)
handles, labels = ax.get_legend_handles_labels()
plt.ylim(1100,1800)
plt.ylabel('德国老龄人口数量（万人）')
#========================================
ax =ax.twinx()
sns.lineplot(x=gm.year,y=gm['all']/10,marker='o',label='德国长护险政策支持人数',color=sns.color_palette()[1])#all：德国长护险政策支持人数
ax.set_xticks(gm2.year)  # 设置刻度标签位置
ax.set_xticklabels(gm2.year.to_list(),fontsize=10)
handles2, labels2 = ax.get_legend_handles_labels()
plt.legend(handles+handles2,labels+labels2,fontsize=10,loc='upper left')
plt.ylabel('德国长护险政策支持人数（万人）')
#========================================
ax = fig.add_subplot(2,1,2)
plt.bar(jp.year,jp['n65'],color=sns.color_palette()[0],label='日本老龄人口数量(65岁以上),第1号被保险者')#jp['n65']：日本老龄人口数量(65岁以上),第1号被保险者
handles, labels = ax.get_legend_handles_labels()
plt.ylim(1500,3700)
plt.ylabel('日本老龄人口数量（万人）')
#========================================
ax =ax.twinx()
sns.lineplot(x=jp.year,y=jp['r'],marker='o',label='日本长护险政策支持人数(照护等级3-5)',color=sns.color_palette()[1])#jp['r']：日本长护险政策支持人数(照护等级3-5)
ax.set_xticks(jp.year)  # 设置刻度标签位置
ax.set_xticklabels(jp.year.to_list(),fontsize=10)
handles2, labels2 = ax.get_legend_handles_labels()
plt.legend(handles+handles2,labels+labels2,fontsize=10,loc='upper left')
plt.ylabel('日本长护险政策支持人数（万人）')
plt.show()