基于机器学习的自杀分析

背景

 

自杀是一项公认的全球性的重大社会问题和公共卫生问题，严重影响社会和经济的发展。据世界卫生组织估计，每年全世界有大约78.6万人自杀，其比例约为每年每10万人中有10.7人，也就是说每隔40秒就有人自杀。据联合国新闻，世卫组织发布的“2019年全球自杀状况”报告中的最新估计显示，自杀仍是全球主要死因之一。高收入国家的自杀率最高，美洲是唯一一个自杀率上升的地区。根据近年来世界各国的自杀率，可以发现：自杀率在性别上的差异在各国之间非常明显，一般的自杀率会男性高于女性,目前男性自杀率最高的国家为前苏联解体后的立陶宛、俄罗斯联邦、拉托维亚和爱沙尼亚等国家(年自杀率>60/10万)，而女性自杀率最高的国家依次为斯里兰卡、中国、匈牙利和爱沙尼亚等(年自杀率>14/10万)。而在部分非洲和拉丁美洲国家自杀率却非常低，如埃及、秘鲁等(年自杀率<1/10万)。国际上习惯上将年自杀率>20/10万的国家称为高自杀率国家，年自杀率<10/10万的国家称为低自杀率国家，李振涛说，在1993年以前的统计中，中国属于低自杀率国家。 根据世界卫生组织在2019年发布的报告，全球每40秒钟就有一个人自杀身亡，自杀已经成为导致人类死亡的主要原因之一。一些国家的自杀人数增长迅速，甚至达60%。中国93%有自杀行为的人没有看过心理医生，在每年8万自杀未遂者中，被进行心理评估的还不到1%。中国的自杀防御工作缺乏一个全国性的计划来协调，政府机构对这一问题的重视不够，缺乏财力支持，也缺乏有效的评估心理因素的工具，更缺少高素质的研究人。

 

数据集来源

 

数据集来源：https://www.kaggle.com/datasets/russellyates88/suicide-rates-overview-1985-to-2016

 

分析思路

 

 

 

本次分析基于Kaggle数据集对自杀行文进行分析，探索自杀行为影响因素，并对自杀因素进行可视化，直观的表现影响自杀率的行为，并使用logistic回归、决策树、随机森林三种机器学习算法建立自杀行为分类模型，评估自杀结果，预测自杀率。

 

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
#import scikitplot as skplt
import sklearn as sk
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
plt.rcParams['font.sans-serif']='SimHei'
plt.rcParams['axes.unicode_minus']=False

import warnings
warnings.filterwarnings("ignore")

数据读取

# 导入原始自杀数据集并重新命名列：
originaldata = pd.read_csv('./master.csv')
# 国家，  年份，  性别，  年龄，  自杀人数，  人口数量，  自杀人数/总人口*100000（自杀率），
# 国家-年份，     人类发展指数（用以衡量联合国各成员国经济社会发展水平的指标，是对传统的GNP指标挑战的结果。）
# 年度国内生产总值(衡量经济发展的指标)，      年人均国内生产总值：国内生产总值/人口 ，         世代
originaldata.columns = ['country', 'year', 'sex', 'age', 'suicides_no', 'population','suicidesper100k',
                      'country-year', 'yearlyHDI', 
                        'GDPpyear', 'GDPpcapita', 'generation']

originaldata.head()

 

originaldata.shape

 

originaldata.info()

 

(originaldata=="male").sum()

 

a = 13910/27820
a

 

df3 = originaldata.groupby("country").agg({"country":"count"})
df3.columns = ["count"]
df3

 

# 修复和清理原始数据
# 将 年度国内生产总值 数据中 ','分割去除掉（比如2,156,624,900），并且转换成float的数字类型
originaldata['GDPpyear'] = originaldata.apply(lambda x: float(x['GDPpyear'].replace(',', '')), axis=1)
originaldata['GDPpyear'].head(10)

# sex 转换为category类型
# Categoricals 是 pandas 的一种数据类型，对应着被统计的变量。Categoricals 是由固定的且有限数量的变量组成的。
originaldata.sex.astype('category')

 

df = originaldata.copy()#复制一份数据
df

 

df.suicidesper100k.mean()
df.suicidesper100k.std()
min(df.suicidesper100k)
max(df.suicidesper100k)

 

import seaborn as sns
sns.distplot(df.suicidesper100k)
plt.xticks(fontsize=10)
plt.yticks(fontsize=20)
plt.show()

 

探索性分析

col = plt.cm.Spectral(np.linspace(0, 1, 20))

plt.figure(figsize=(8, 6))
agedistf = pd.DataFrame(df.groupby('sex').get_group('female').groupby('age').suicides_no.sum())
agedistm = pd.DataFrame(df.groupby('sex').get_group('male').groupby('age').suicides_no.sum())
plt.bar(agedistm.index, agedistm.suicides_no, color=col[18])
plt.bar(agedistf.index, agedistf.suicides_no, color=col[8])
plt.legend(['male', 'female'], fontsize=16)
plt.ylabel('Count', fontsize=14)
plt.xlabel('Suicides per 100K', fontsize=14)
plt.xticks(fontsize=10)
plt.yticks(fontsize=20)
plt.show()

 

col = plt.cm.Spectral(np.linspace(0, 1, 22))
plt.figure(figsize=(12, 15))

plt.subplot(211)
#自杀率（放大1000倍）的平均值最高的前10个国家
df.groupby(['country']).suicidesper100k.mean().nlargest(10).plot(kind='bar', color=col, fontsize=20)
plt.xlabel('Average Suicides/100k', size=20)
plt.ylabel('Country', fontsize=20)
plt.title('Top 10 countries', fontsize=30)

plt.figure(figsize=(12, 15))
plt.subplot(212)
#自杀人数的平均值最高的前10个国家
df.groupby(['country']).suicides_no.mean().nlargest(10).plot(kind='bar', color=col, fontsize=20)
plt.xlabel('Average Suicides_no', size=20)
plt.ylabel('Country', fontsize=20);
plt.title('Top 10 countries', fontsize=30)
plt.show()

 

 

plt.figure(figsize=(10, 16))

#总人口的各个年龄段的性别分布
plt.subplot(311)
sns.barplot(x='sex', y='population', hue='age', data=df, palette="Greens")  #hue按年龄分组
plt.xticks(ha='right', fontsize=20)
plt.ylabel('Population', fontsize=20)
plt.xlabel('Sex', fontsize=20)
plt.title("不同年龄段的男、女总人口数")
plt.legend(fontsize=14, loc='best')

#自杀人数的各个年龄段的性别分布
plt.subplot(312)
sns.barplot(x='sex', y='suicides_no', hue='age', data=df, palette="Greens")
plt.xticks(ha='right', fontsize=20)
plt.ylabel('suicides incidences', fontsize=20)
plt.xlabel('Sex', fontsize=20)
plt.title("不同年龄段的男、女自杀人口数")
plt.legend(fontsize=14)

#自杀率（放大1000倍）的各个年龄段的性别分布
plt.subplot(313)
sns.barplot(x='sex', y='suicidesper100k', hue='age', data=df,palette="Greens")
plt.xticks(ha='right', fontsize=20);
plt.ylabel('suicidesper100k',fontsize=20);
plt.xlabel('Sex',fontsize=20);
plt.title("不同年龄段的男、女自杀率")
plt.legend(fontsize=14);

plt.subplots_adjust(top=1.2)
plt.show()

 

plt.figure(figsize=(12, 16))

# 具体按性别和年份 来看 男女的自杀率
plt.subplot(311)
sns.lineplot(x='year', y='suicidesper100k', hue='sex', data=df, palette="hot")  #hue按年龄分组
plt.xticks(ha='right', fontsize=20)
plt.ylabel('suicidesper100k', fontsize=20)
plt.xlabel('year', fontsize=20)
plt.legend(fontsize=14, loc='best')  
plt.title("性别年份与自杀率关系图")
plt.show()

 

plt.figure(figsize=(8, 6))

#选出year列并且去除重复的年份
year = originaldata.groupby('year').year.unique()

#各个年份的自杀人数汇总
#使用seaborn进行可视化,输入的数据必须为dataframe
totalpyear = pd.DataFrame(originaldata.groupby('year').suicides_no.sum())   

plt.plot(year.index[0:31], totalpyear[0:31], color=col[18])  #选取范围为[0:31]  1985年到2015年
plt.xlabel('year', fontsize=15)
plt.ylabel('Total number of suicides in the world', fontsize=15)
plt.show()

 

plt.figure(figsize=(20, 8))

# 自杀率（放大1000倍）的分布情况，y轴为个数
plt.subplot(121)
plt.hist(df.suicidesper100k, bins=30, color=col[18])   #bins:条形数
plt.xlabel('Suicides per 100K of population', fontsize=25)
plt.xticks(rotation = 0,                 fontsize = 20)  
plt.ylabel('count', fontsize=25)
plt.yticks(                              fontsize = 20)   

# 年人均国内生产总值的分布情况，y轴为个数
plt.subplot(122)
plt.hist(df.GDPpcapita, bins=30, color=col[7])
plt.xlabel('GDP', fontsize=25)
plt.xticks(rotation = 0,fontsize = 20)  
plt.ylabel('count', fontsize=25)
plt.yticks(fontsize = 20)   
plt.show()

 

corr = df.corr() #相关系数矩阵，即给出了任意两个变量之间的相关系数

# 相关矩阵的上三角部分与下三角对称。因此，热图不需要显示整个矩阵。在下一步隐藏上三角形。
# 设置mask隐藏上三角
# np.zeros_like() 返回一个零数组，其形状和类型与给定的数组相同。
# 该 dtype=np.bool 参数会覆盖数据类型，因此我们的数组是一个布尔数组。
# np.triu_indices_from(mask) 返回数组上三角形的索引。
# 现在，我们将上三角形设置为True。 mask[np.triu_indices_from(mask)]= True
mask = np.zeros_like(corr, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True

# 在Seaborn中创建热图
f, ax = plt.subplots(figsize=(10, 8))
# 生成自定义发散颜色图
cmap = sns.diverging_palette(220, 10, as_cmap=True)

# 绘制热图
# 数据为 corr
# vmax,vmin:分别是热力图的颜色取值最大和最小范围
# center:数据表取值有差异时，设置热力图的色彩中心对齐值；通过设置center值，可以调整生成的图像颜色的整体深浅
# square:设置热力图矩阵小块形状，默认值是False
# linewidths(矩阵小块的间隔),
# cbar_kws:热力图侧边绘制颜色刻度条时，相关字体设置，默认值是None
sns.heatmap(corr, mask=mask, cmap=cmap, vmax=1, vmin=-1, center=0,  
            square=True, linewidths=0.2, cbar_kws={"shrink": 0.8},annot=True)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.show()

 

数据预处理

#复制数据
total = df.copy()
# 对总人口用均值填充
total.population.fillna(total.population.mean(), inplace=True)
# 小于平均值的为低风险，大于平均值的为高风险
total['risk'] = np.where(total.suicidesper100k < total.suicidesper100k.mean(), 0, 1)

plt.figure(figsize=(16, 5))
plt.subplot(121)
plt.hist(total.risk, color=col[8])
plt.ylabel('counts', fontsize=20)
plt.xlabel('Suicide risk', fontsize=20)
plt.title("自杀高低风险数量")

plt.subplot(122)
sns.distplot(total.suicidesper100k[total.risk == 0], bins=10)
sns.distplot(total.suicidesper100k[total.risk == 1], bins=20)  
plt.xlabel('Suicides', fontsize=20)
plt.title("自杀高低风险分布情况")
plt.show()

 

(total['risk']==1).sum()
(total['risk']==0).sum()
total.head(5)
total['yearlyHDI'].isnull().sum()
a = 19456/27820
a

 

total['yearlyHDI'].describe()
total['yearlyHDI'].isnull().sum()
total

 

fill_1 = total['yearlyHDI']
fill_1 = pd.DataFrame(fill_1)
fill_1

 

from scipy.interpolate import lagrange #导入拉格朗日函数
def ployinterp_column(s,n,k=2,result_type='broadcast'): # k=2表示用空值的前后两个数值来拟合曲线，从而预测空值
    y = s[list(range(n+2-k,n+2)) + list(range(n+3,n+3-k))] # 取值，range函数返回一个左闭右开（[left,right)）的序列数
    y = y[y.notnull()] # 取上一行中取出数值列表中的非空值，保证y的每行都有数值，便于拟合函数
    return lagrange(y.index,list(y))(n) # 调用拉格朗日函数，并添加索引
for i in fill_1.columns: # 如果i在data的列名中，data.columns生成的是data的全部列名
    for j in range(len(fill_1)): # len(data)返回了data的长度，若此长度为5，则range(5)会产生从0开始计数的整数列表
        if (fill_1[i].isnull())[j]: # 如果data[i][j]为空，则调用函数ployinterp_column为其插值
            fill_1[i][j] = ployinterp_column(fill_1[i],j,result_type='broadcast')  
fill_1

 

# 对国家进行 标准化标签，将标签值统一转换成range(标签值个数-1)范围内
# 相当于fit(X).transform(X),意思就是先进行fit(),进行数据拟合，然后在进行transform() 进行标准化处理

from sklearn.preprocessing import LabelEncoder
le = LabelEncoder()
total.country = le.fit_transform(total.country)
total.country.unique()

# 为建立模型准备数据
X = total.drop(['suicidesper100k','risk','country-year'],axis = 1)
#用0进行填补
X['yearlyHDI'] = fill_1['yearlyHDI']
y = total['risk']
#对X进行独热编码处理
XS = pd.get_dummies(X)

 

#GDP和人口的分布
ax1 = total[total['risk'] == 1][0:200].plot(kind='scatter', x='GDPpcapita', y='population', color='DarkRed',
                                            label='high risk', figsize=(6, 5), fontsize=12)
total[total['risk'] == 0][0:200].plot(kind='scatter', x='GDPpcapita', y='population', color='DarkBlue',
                                            label='low risk', ax=ax1)

plt.ylabel('population', fontsize=16)
plt.xlabel('GDP per capita', fontsize=16)
plt.legend(fontsize=14)

 

#GDP和人类发展指数的分布
ax1 = total[total['risk'] == 1][0:200].plot(kind='scatter', x='GDPpcapita', y='yearlyHDI', color='DarkRed',
                                            label='high risk', figsize=(6, 5), fontsize=12)
total[total['risk'] == 0][0:200].plot(kind='scatter', x='GDPpcapita', y='yearlyHDI', color='DarkBlue',
                                            label='low risk', ax=ax1)

plt.ylabel('yearlyHDI', fontsize=16)
plt.xlabel('GDP per capita', fontsize=16)
plt.legend(fontsize=14)

 

# 关于年份的高低自杀风险
fig = plt.figure(figsize=(30, 30))

plt.subplot(4, 3, 1)
sns.distplot(total[total.columns[0]][total.risk == 0], label='low risk')
sns.distplot(total[total.columns[0]][total.risk == 1], label='high risk')       
plt.legend(loc='best', fontsize=18)  
plt.xlabel(total.columns[0], fontsize=18)
plt.show()

 

建模

from sklearn.model_selection import train_test_split
# 将数据集拆分为训练集和测试集
# 首先将该数据3/4作为训练集，1/4作为测试集；再由训练集划分1/5作为验证集
X_train, X_test, y_train, y_test = train_test_split(XS, y, test_size=0.25, random_state=4)
X_train,X_valid,y_train,y_valid=train_test_split(X_train,y_train,test_size=0.2,random_state=4)
print('Train set:', X_train.shape, y_train.shape)
print('Test set:', X_test.shape, y_test.shape)
print('Test set:', X_valid.shape, y_valid.shape)

 

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import recall_score,accuracy_score
from sklearn.metrics import precision_recall_fscore_support
from sklearn.metrics import confusion_matrix, classification_report

LR = LogisticRegression(C=0.001, solver='liblinear')
LR.fit(X_train, y_train)

# 预测类别:0还是1
yLRhat = LR.predict(X_valid)
 
# 预测 0或1的概率（例如 [0.54689436, 0.45310564] 预测出来为0）
yLRhat_prob = LR.predict_proba(X_valid)

cm = confusion_matrix(y_valid, yLRhat)
print('\n confusion matrix：混淆矩阵 \n', cm)
print('此时模型准确率为:',accuracy_score(y_valid, yLRhat))
print('********************************************************')
print('\n')
print('逻辑回归模型的分类报告\n', classification_report(y_valid, yLRhat))
#skplt.metrics.plot_confusion_matrix(y_valid,yLRhat,normalize=True)

 

X_train.head()

 

from sklearn.model_selection import GridSearchCV

param={"tol":[1e-4, 1e-3,1e-2], "C":[0.4, 0.6, 0.8]}
grid = GridSearchCV(LogisticRegression(),param_grid=param, cv=5)#这里是5折交叉验证
grid.fit(X_train,y_train)
print(grid.best_params_)
print(grid.best_score_)

#得到最好的逻辑回归分类器
best_LR=grid.best_estimator_ 

#进行训练、预测
best_LR.fit(X_train,y_train)
pred=best_LR.predict(X_valid)

cm = confusion_matrix(y_valid, pred)
print('\n confusion matrix：混淆矩阵 \n', cm)
print('最好的逻辑回归模型准确率为:',accuracy_score(y_valid, pred))
print('********************************************************')
print('\n')
print('最好的逻辑回归模型的分类报告\n', classification_report(y_valid, pred))

 

from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix, classification_report

# 决策树学习
# 函数为创建一个决策树模型
# criterion：gini或者entropy,前者是基尼系数，后者是信息熵。
# max_depth：  int or None, optional (default=None) 设置决策随机森林中的决策树的最大深度，深度越大，越容易过拟合，推荐树的深度为：5-20之间。
# max_leaf_nodes： 通过限制最大叶子节点数，可以防止过拟合，默认是"None”，即不限制最大的叶子节点数。
DT = DecisionTreeClassifier(criterion="entropy", max_depth=7, max_leaf_nodes=30)
DT.fit(X_train, y_train)
ydthat = DT.predict(X_valid)

print('***决策树模型***')

#决策树模型性能评估
print('验证集上的准确率:', DT.score(X_valid, y_valid))
print('训练集上的准确率:', DT.score(X_train, y_train))

# 混淆矩阵
print('CM\n', confusion_matrix(y_valid, ydthat))
print('********************************************************')

print('决策树模型的分类报告\n', classification_report(y_valid, ydthat))

DTfeat_importance = DT.feature_importances_
DTfeat_importance = pd.DataFrame([X_train.columns, DT.feature_importances_]).T

# 特征重要性排序
print("特征重要性排名如下:")
print(DTfeat_importance.sort_values(by=1, ascending=False))

 

# 通过cv_results观察过程并做图
max_depth=np.arange(0,20,3)
max_leaf_nodes=np.arange(0,30,3)
min_samples_leaf =np.arange(0,20,2)
param1= {'criterion':['entropy'],'max_depth':max_depth}
param2= {'criterion':['entropy'],'min_samples_leaf':min_samples_leaf}
param3= {'criterion':['entropy'],'max_leaf_nodes':max_leaf_nodes}


clf1 = GridSearchCV(DecisionTreeClassifier(),param_grid=param1,cv=6)
clf1.fit(X_train,y_train)

clf2 = GridSearchCV(DecisionTreeClassifier(),param_grid=param2,cv=6)
clf2.fit(X_train,y_train)


clf3 = GridSearchCV(DecisionTreeClassifier(),param_grid=param3,cv=6)
clf3.fit(X_train,y_train)

 

fig = plt.figure(figsize=(20,10), dpi=100)
ax = fig.add_subplot(311)
ax.plot(min_samples_leaf,clf2.cv_results_['mean_test_score'],'g*-',)
plt.title('决策树模型网格搜索的训练过程',fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)

ax = fig.add_subplot(312)
ax.plot(max_depth,clf1.cv_results_['mean_test_score'],'g*-')
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
ax = fig.add_subplot(313)
ax.plot(max_leaf_nodes,clf3.cv_results_['mean_test_score'],'g*-')
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.show()

 

best_DT = DecisionTreeClassifier(criterion="entropy", max_depth=15, max_leaf_nodes=25,min_samples_leaf=2)
best_DT.fit(X_train, y_train)
ydthat =best_DT.predict(X_valid)

print('***决策树模型***')

#决策树模型性能评估
print('验证集上的准确率:',best_DT.score(X_valid, y_valid))
print('训练集上的准确率:',best_DT.score(X_train, y_train))

# 混淆矩阵
print('CM\n', confusion_matrix(y_valid, ydthat))
print('********************************************************')

print('决策树模型的分类报告\n', classification_report(y_valid, ydthat))

DTfeat_importance = best_DT.feature_importances_
DTfeat_importance = pd.DataFrame([X_train.columns, DT.feature_importances_]).T

# 特征重要性排序
print("特征重要性排名如下:")
print(DTfeat_importance.sort_values(by=1, ascending=False))

 

# 随机森林可以视为多颗决策树的集成，鲁棒性更强，泛化能力更好，不易产生过拟合现象。但是噪声比较大的情况下会过拟合。

from sklearn.ensemble import RandomForestClassifier

random_forest = RandomForestClassifier(n_estimators=20, max_depth=10, min_samples_split=2, min_samples_leaf=5,
                                       max_leaf_nodes=20, max_features=len(X_train.columns)) 

random_forest.fit(X_train, y_train)

yrfhat = random_forest.predict(X_valid)
feat_importance = random_forest.feature_importances_
rffeat_importance = pd.DataFrame([X_train.columns, random_forest.feature_importances_]).T

print('******************Random forest classifier**************')
print('训练集上的准确率', random_forest.score(X_train, y_train))
print('验证集上的准确率', random_forest.score(X_valid,y_valid))
print('混淆矩阵\n', confusion_matrix(y_valid, yrfhat))
print('********************************************************')
print('随机森林的分类报告\n', classification_report(y_valid, yrfhat))
print(rffeat_importance.sort_values(by=1, ascending=False))

 

#我们首先对n_estimators进行网格搜索：
param_test1 = {'n_estimators':[50,120,160,200,250]}
gsearch1 = GridSearchCV(estimator = RandomForestClassifier(), param_grid = param_test1,cv=5)
gsearch1.fit(X_train,y_train)
print( gsearch1.best_params_, gsearch1.best_score_)
#接着我们对决策树最大深度max_depth和内部节点再划分所需最小样本数min_samples_split进行网格搜索。
param_test2 = {'max_depth':[1,2,3,5,7,9,11,13,20],'min_samples_split':[5,10,25,50,100,120,150]}
gsearch2 = GridSearchCV(estimator = RandomForestClassifier(n_estimators=160),param_grid = param_test2, cv=5)
gsearch2.fit(X_train,y_train)
print( gsearch2.best_params_, gsearch2.best_score_)

 

#最好的随机森林模型
best_rf=RandomForestClassifier(n_estimators=160, max_depth=20, min_samples_split=5) 

best_rf .fit(X_train, y_train)

yrfhat = best_rf .predict(X_valid)
feat_importance = best_rf.feature_importances_
rffeat_importance = pd.DataFrame([X_train.columns, random_forest.feature_importances_]).T

print('******************随机森林模型**************')
print('训练集上的准确率',best_rf.score(X_train, y_train))
print('验证集上的准确率',best_rf.score(X_valid,y_valid))
print('混淆矩阵\n', confusion_matrix(y_valid, yrfhat))
print('********************************************************')
print('最好的随机森林模型的分类报告\n', classification_report(y_valid, yrfhat))
print(rffeat_importance.sort_values(by=1, ascending=False))

 

 

全代码附上

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
#import scikitplot as skplt
import sklearn as sk
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
plt.rcParams['font.sans-serif']='SimHei'
plt.rcParams['axes.unicode_minus']=False

import warnings
warnings.filterwarnings("ignore")
# 导入原始自杀数据集并重新命名列：
originaldata = pd.read_csv('./master.csv')
# 国家，  年份，  性别，  年龄，  自杀人数，  人口数量，  自杀人数/总人口*100000（自杀率），
# 国家-年份，     人类发展指数（用以衡量联合国各成员国经济社会发展水平的指标，是对传统的GNP指标挑战的结果。）
# 年度国内生产总值(衡量经济发展的指标)，      年人均国内生产总值：国内生产总值/人口 ，         世代
originaldata.columns = ['country', 'year', 'sex', 'age', 'suicides_no', 'population','suicidesper100k',
                      'country-year', 'yearlyHDI', 
                        'GDPpyear', 'GDPpcapita', 'generation']

originaldata.head()
originaldata.shape
originaldata.info()
(originaldata=="male").sum()
a = 13910/27820
a
df3 = originaldata.groupby("country").agg({"country":"count"})
df3.columns = ["count"]
df3
# 修复和清理原始数据
# 将 年度国内生产总值 数据中 ','分割去除掉（比如2,156,624,900），并且转换成float的数字类型
originaldata['GDPpyear'] = originaldata.apply(lambda x: float(x['GDPpyear'].replace(',', '')), axis=1)
originaldata['GDPpyear'].head(10)

# sex 转换为category类型
# Categoricals 是 pandas 的一种数据类型，对应着被统计的变量。Categoricals 是由固定的且有限数量的变量组成的。
originaldata.sex.astype('category')
df = originaldata.copy()#复制一份数据
df
df.suicidesper100k.mean()
df.suicidesper100k.std()
min(df.suicidesper100k)
max(df.suicidesper100k)
import seaborn as sns
sns.distplot(df.suicidesper100k)
plt.xticks(fontsize=10)
plt.yticks(fontsize=20)
plt.show()
df.suicidesper100k.isnull().sum()
col = plt.cm.Spectral(np.linspace(0, 1, 20))

plt.figure(figsize=(8, 6))
agedistf = pd.DataFrame(df.groupby('sex').get_group('female').groupby('age').suicides_no.sum())
agedistm = pd.DataFrame(df.groupby('sex').get_group('male').groupby('age').suicides_no.sum())
plt.bar(agedistm.index, agedistm.suicides_no, color=col[18])
plt.bar(agedistf.index, agedistf.suicides_no, color=col[8])
plt.legend(['male', 'female'], fontsize=16)
plt.ylabel('Count', fontsize=14)
plt.xlabel('Suicides per 100K', fontsize=14)
plt.xticks(fontsize=10)
plt.yticks(fontsize=20)
plt.show()
col = plt.cm.Spectral(np.linspace(0, 1, 22))
plt.figure(figsize=(12, 15))

plt.subplot(211)
#自杀率（放大1000倍）的平均值最高的前10个国家
df.groupby(['country']).suicidesper100k.mean().nlargest(10).plot(kind='bar', color=col, fontsize=20)
plt.xlabel('Average Suicides/100k', size=20)
plt.ylabel('Country', fontsize=20)
plt.title('Top 10 countries', fontsize=30)

plt.figure(figsize=(12, 15))
plt.subplot(212)
#自杀人数的平均值最高的前10个国家
df.groupby(['country']).suicides_no.mean().nlargest(10).plot(kind='bar', color=col, fontsize=20)
plt.xlabel('Average Suicides_no', size=20)
plt.ylabel('Country', fontsize=20);
plt.title('Top 10 countries', fontsize=30)
plt.show()
plt.figure(figsize=(10, 16))

#总人口的各个年龄段的性别分布
plt.subplot(311)
sns.barplot(x='sex', y='population', hue='age', data=df, palette="Greens")  #hue按年龄分组
plt.xticks(ha='right', fontsize=20)
plt.ylabel('Population', fontsize=20)
plt.xlabel('Sex', fontsize=20)
plt.title("不同年龄段的男、女总人口数")
plt.legend(fontsize=14, loc='best')

#自杀人数的各个年龄段的性别分布
plt.subplot(312)
sns.barplot(x='sex', y='suicides_no', hue='age', data=df, palette="Greens")
plt.xticks(ha='right', fontsize=20)
plt.ylabel('suicides incidences', fontsize=20)
plt.xlabel('Sex', fontsize=20)
plt.title("不同年龄段的男、女自杀人口数")
plt.legend(fontsize=14)

#自杀率（放大1000倍）的各个年龄段的性别分布
plt.subplot(313)
sns.barplot(x='sex', y='suicidesper100k', hue='age', data=df,palette="Greens")
plt.xticks(ha='right', fontsize=20);
plt.ylabel('suicidesper100k',fontsize=20);
plt.xlabel('Sex',fontsize=20);
plt.title("不同年龄段的男、女自杀率")
plt.legend(fontsize=14);

plt.subplots_adjust(top=1.2)
plt.show()
plt.figure(figsize=(12, 16))

# 具体按性别和年份 来看 男女的自杀率
plt.subplot(311)
sns.lineplot(x='year', y='suicidesper100k', hue='sex', data=df, palette="hot")  #hue按年龄分组
plt.xticks(ha='right', fontsize=20)
plt.ylabel('suicidesper100k', fontsize=20)
plt.xlabel('year', fontsize=20)
plt.legend(fontsize=14, loc='best')  
plt.title("性别年份与自杀率关系图")
plt.show()
plt.figure(figsize=(8, 6))

#选出year列并且去除重复的年份
year = originaldata.groupby('year').year.unique()

#各个年份的自杀人数汇总
#使用seaborn进行可视化,输入的数据必须为dataframe
totalpyear = pd.DataFrame(originaldata.groupby('year').suicides_no.sum())   

plt.plot(year.index[0:31], totalpyear[0:31], color=col[18])  #选取范围为[0:31]  1985年到2015年
plt.xlabel('year', fontsize=15)
plt.ylabel('Total number of suicides in the world', fontsize=15)
plt.show()
plt.figure(figsize=(20, 8))

# 自杀率（放大1000倍）的分布情况，y轴为个数
plt.subplot(121)
plt.hist(df.suicidesper100k, bins=30, color=col[18])   #bins:条形数
plt.xlabel('Suicides per 100K of population', fontsize=25)
plt.xticks(rotation = 0,                 fontsize = 20)  
plt.ylabel('count', fontsize=25)
plt.yticks(                              fontsize = 20)   

# 年人均国内生产总值的分布情况，y轴为个数
plt.subplot(122)
plt.hist(df.GDPpcapita, bins=30, color=col[7])
plt.xlabel('GDP', fontsize=25)
plt.xticks(rotation = 0,fontsize = 20)  
plt.ylabel('count', fontsize=25)
plt.yticks(fontsize = 20)   
plt.show()
corr = df.corr() #相关系数矩阵，即给出了任意两个变量之间的相关系数

# 相关矩阵的上三角部分与下三角对称。因此，热图不需要显示整个矩阵。在下一步隐藏上三角形。
# 设置mask隐藏上三角
# np.zeros_like() 返回一个零数组，其形状和类型与给定的数组相同。
# 该 dtype=np.bool 参数会覆盖数据类型，因此我们的数组是一个布尔数组。
# np.triu_indices_from(mask) 返回数组上三角形的索引。
# 现在，我们将上三角形设置为True。 mask[np.triu_indices_from(mask)]= True
mask = np.zeros_like(corr, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True

# 在Seaborn中创建热图
f, ax = plt.subplots(figsize=(10, 8))
# 生成自定义发散颜色图
cmap = sns.diverging_palette(220, 10, as_cmap=True)

# 绘制热图
# 数据为 corr
# vmax,vmin:分别是热力图的颜色取值最大和最小范围
# center:数据表取值有差异时，设置热力图的色彩中心对齐值；通过设置center值，可以调整生成的图像颜色的整体深浅
# square:设置热力图矩阵小块形状，默认值是False
# linewidths(矩阵小块的间隔),
# cbar_kws:热力图侧边绘制颜色刻度条时，相关字体设置，默认值是None
sns.heatmap(corr, mask=mask, cmap=cmap, vmax=1, vmin=-1, center=0,  
            square=True, linewidths=0.2, cbar_kws={"shrink": 0.8},annot=True)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.show()
#复制数据
total = df.copy()
# 对总人口用均值填充
total.population.fillna(total.population.mean(), inplace=True)
# 小于平均值的为低风险，大于平均值的为高风险
total['risk'] = np.where(total.suicidesper100k < total.suicidesper100k.mean(), 0, 1)

plt.figure(figsize=(16, 5))
plt.subplot(121)
plt.hist(total.risk, color=col[8])
plt.ylabel('counts', fontsize=20)
plt.xlabel('Suicide risk', fontsize=20)
plt.title("自杀高低风险数量")

plt.subplot(122)
sns.distplot(total.suicidesper100k[total.risk == 0], bins=10)
sns.distplot(total.suicidesper100k[total.risk == 1], bins=20)  
plt.xlabel('Suicides', fontsize=20)
plt.title("自杀高低风险分布情况")
plt.show()
(total['risk']==1).sum()
(total['risk']==0).sum()
total.head(5)
total['yearlyHDI'].isnull().sum()
a = 19456/27820
a
total['yearlyHDI'].describe()
total['yearlyHDI'].isnull().sum()
total
fill_1 = total['yearlyHDI']
fill_1 = pd.DataFrame(fill_1)
fill_1
from scipy.interpolate import lagrange #导入拉格朗日函数
def ployinterp_column(s,n,k=2,result_type='broadcast'): # k=2表示用空值的前后两个数值来拟合曲线，从而预测空值
    y = s[list(range(n+2-k,n+2)) + list(range(n+3,n+3-k))] # 取值，range函数返回一个左闭右开（[left,right)）的序列数
    y = y[y.notnull()] # 取上一行中取出数值列表中的非空值，保证y的每行都有数值，便于拟合函数
    return lagrange(y.index,list(y))(n) # 调用拉格朗日函数，并添加索引
for i in fill_1.columns: # 如果i在data的列名中，data.columns生成的是data的全部列名
    for j in range(len(fill_1)): # len(data)返回了data的长度，若此长度为5，则range(5)会产生从0开始计数的整数列表
        if (fill_1[i].isnull())[j]: # 如果data[i][j]为空，则调用函数ployinterp_column为其插值
            fill_1[i][j] = ployinterp_column(fill_1[i],j,result_type='broadcast')  
fill_1
# 对国家进行 标准化标签，将标签值统一转换成range(标签值个数-1)范围内
# 相当于fit(X).transform(X),意思就是先进行fit(),进行数据拟合，然后在进行transform() 进行标准化处理

from sklearn.preprocessing import LabelEncoder
le = LabelEncoder()
total.country = le.fit_transform(total.country)
total.country.unique()

# 为建立模型准备数据
X = total.drop(['suicidesper100k','risk','country-year'],axis = 1)
#用0进行填补
X['yearlyHDI'] = fill_1['yearlyHDI']
y = total['risk']
#对X进行独热编码处理
XS = pd.get_dummies(X)
#GDP和人口的分布
ax1 = total[total['risk'] == 1][0:200].plot(kind='scatter', x='GDPpcapita', y='population', color='DarkRed',
                                            label='high risk', figsize=(6, 5), fontsize=12)
total[total['risk'] == 0][0:200].plot(kind='scatter', x='GDPpcapita', y='population', color='DarkBlue',
                                            label='low risk', ax=ax1)

plt.ylabel('population', fontsize=16)
plt.xlabel('GDP per capita', fontsize=16)
plt.legend(fontsize=14)
#GDP和人类发展指数的分布
ax1 = total[total['risk'] == 1][0:200].plot(kind='scatter', x='GDPpcapita', y='yearlyHDI', color='DarkRed',
                                            label='high risk', figsize=(6, 5), fontsize=12)
total[total['risk'] == 0][0:200].plot(kind='scatter', x='GDPpcapita', y='yearlyHDI', color='DarkBlue',
                                            label='low risk', ax=ax1)

plt.ylabel('yearlyHDI', fontsize=16)
plt.xlabel('GDP per capita', fontsize=16)
plt.legend(fontsize=14)
# 关于年份的高低自杀风险
fig = plt.figure(figsize=(30, 30))

plt.subplot(4, 3, 1)
sns.distplot(total[total.columns[0]][total.risk == 0], label='low risk')
sns.distplot(total[total.columns[0]][total.risk == 1], label='high risk')       
plt.legend(loc='best', fontsize=18)  
plt.xlabel(total.columns[0], fontsize=18)
plt.show()
from sklearn.model_selection import train_test_split
# 将数据集拆分为训练集和测试集
# 首先将该数据3/4作为训练集，1/4作为测试集；再由训练集划分1/5作为验证集
X_train, X_test, y_train, y_test = train_test_split(XS, y, test_size=0.25, random_state=4)
X_train,X_valid,y_train,y_valid=train_test_split(X_train,y_train,test_size=0.2,random_state=4)
print('Train set:', X_train.shape, y_train.shape)
print('Test set:', X_test.shape, y_test.shape)
print('Test set:', X_valid.shape, y_valid.shape)
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import recall_score,accuracy_score
from sklearn.metrics import precision_recall_fscore_support
from sklearn.metrics import confusion_matrix, classification_report

LR = LogisticRegression(C=0.001, solver='liblinear')
LR.fit(X_train, y_train)

# 预测类别:0还是1
yLRhat = LR.predict(X_valid)
 
# 预测 0或1的概率（例如 [0.54689436, 0.45310564] 预测出来为0）
yLRhat_prob = LR.predict_proba(X_valid)

cm = confusion_matrix(y_valid, yLRhat)
print('\n confusion matrix：混淆矩阵 \n', cm)
print('此时模型准确率为:',accuracy_score(y_valid, yLRhat))
print('********************************************************')
print('\n')
print('逻辑回归模型的分类报告\n', classification_report(y_valid, yLRhat))
#skplt.metrics.plot_confusion_matrix(y_valid,yLRhat,normalize=True)
X_train.head()
from sklearn.model_selection import GridSearchCV

param={"tol":[1e-4, 1e-3,1e-2], "C":[0.4, 0.6, 0.8]}
grid = GridSearchCV(LogisticRegression(),param_grid=param, cv=5)#这里是5折交叉验证
grid.fit(X_train,y_train)
print(grid.best_params_)
print(grid.best_score_)

#得到最好的逻辑回归分类器
best_LR=grid.best_estimator_ 

#进行训练、预测
best_LR.fit(X_train,y_train)
pred=best_LR.predict(X_valid)

cm = confusion_matrix(y_valid, pred)
print('\n confusion matrix：混淆矩阵 \n', cm)
print('最好的逻辑回归模型准确率为:',accuracy_score(y_valid, pred))
print('********************************************************')
print('\n')
print('最好的逻辑回归模型的分类报告\n', classification_report(y_valid, pred))
from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix, classification_report

# 决策树学习
# 函数为创建一个决策树模型
# criterion：gini或者entropy,前者是基尼系数，后者是信息熵。
# max_depth：  int or None, optional (default=None) 设置决策随机森林中的决策树的最大深度，深度越大，越容易过拟合，推荐树的深度为：5-20之间。
# max_leaf_nodes： 通过限制最大叶子节点数，可以防止过拟合，默认是"None”，即不限制最大的叶子节点数。
DT = DecisionTreeClassifier(criterion="entropy", max_depth=7, max_leaf_nodes=30)
DT.fit(X_train, y_train)
ydthat = DT.predict(X_valid)

print('***决策树模型***')

#决策树模型性能评估
print('验证集上的准确率:', DT.score(X_valid, y_valid))
print('训练集上的准确率:', DT.score(X_train, y_train))

# 混淆矩阵
print('CM\n', confusion_matrix(y_valid, ydthat))
print('********************************************************')

print('决策树模型的分类报告\n', classification_report(y_valid, ydthat))

DTfeat_importance = DT.feature_importances_
DTfeat_importance = pd.DataFrame([X_train.columns, DT.feature_importances_]).T

# 特征重要性排序
print("特征重要性排名如下:")
print(DTfeat_importance.sort_values(by=1, ascending=False))
# 通过cv_results观察过程并做图
max_depth=np.arange(0,20,3)
max_leaf_nodes=np.arange(0,30,3)
min_samples_leaf =np.arange(0,20,2)
param1= {'criterion':['entropy'],'max_depth':max_depth}
param2= {'criterion':['entropy'],'min_samples_leaf':min_samples_leaf}
param3= {'criterion':['entropy'],'max_leaf_nodes':max_leaf_nodes}


clf1 = GridSearchCV(DecisionTreeClassifier(),param_grid=param1,cv=6)
clf1.fit(X_train,y_train)

clf2 = GridSearchCV(DecisionTreeClassifier(),param_grid=param2,cv=6)
clf2.fit(X_train,y_train)


clf3 = GridSearchCV(DecisionTreeClassifier(),param_grid=param3,cv=6)
clf3.fit(X_train,y_train)

fig = plt.figure(figsize=(20,10), dpi=100)
ax = fig.add_subplot(311)
ax.plot(min_samples_leaf,clf2.cv_results_['mean_test_score'],'g*-',)
plt.title('决策树模型网格搜索的训练过程',fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)

ax = fig.add_subplot(312)
ax.plot(max_depth,clf1.cv_results_['mean_test_score'],'g*-')
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
ax = fig.add_subplot(313)
ax.plot(max_leaf_nodes,clf3.cv_results_['mean_test_score'],'g*-')
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.show()
best_DT = DecisionTreeClassifier(criterion="entropy", max_depth=15, max_leaf_nodes=25,min_samples_leaf=2)
best_DT.fit(X_train, y_train)
ydthat =best_DT.predict(X_valid)

print('***决策树模型***')

#决策树模型性能评估
print('验证集上的准确率:',best_DT.score(X_valid, y_valid))
print('训练集上的准确率:',best_DT.score(X_train, y_train))

# 混淆矩阵
print('CM\n', confusion_matrix(y_valid, ydthat))
print('********************************************************')

print('决策树模型的分类报告\n', classification_report(y_valid, ydthat))

DTfeat_importance = best_DT.feature_importances_
DTfeat_importance = pd.DataFrame([X_train.columns, DT.feature_importances_]).T

# 特征重要性排序
print("特征重要性排名如下:")
print(DTfeat_importance.sort_values(by=1, ascending=False))
# 随机森林可以视为多颗决策树的集成，鲁棒性更强，泛化能力更好，不易产生过拟合现象。但是噪声比较大的情况下会过拟合。

from sklearn.ensemble import RandomForestClassifier

random_forest = RandomForestClassifier(n_estimators=20, max_depth=10, min_samples_split=2, min_samples_leaf=5,
                                       max_leaf_nodes=20, max_features=len(X_train.columns)) 

random_forest.fit(X_train, y_train)

yrfhat = random_forest.predict(X_valid)
feat_importance = random_forest.feature_importances_
rffeat_importance = pd.DataFrame([X_train.columns, random_forest.feature_importances_]).T

print('******************Random forest classifier**************')
print('训练集上的准确率', random_forest.score(X_train, y_train))
print('验证集上的准确率', random_forest.score(X_valid,y_valid))
print('混淆矩阵\n', confusion_matrix(y_valid, yrfhat))
print('********************************************************')
print('随机森林的分类报告\n', classification_report(y_valid, yrfhat))
print(rffeat_importance.sort_values(by=1, ascending=False))
#我们首先对n_estimators进行网格搜索：
param_test1 = {'n_estimators':[50,120,160,200,250]}
gsearch1 = GridSearchCV(estimator = RandomForestClassifier(), param_grid = param_test1,cv=5)
gsearch1.fit(X_train,y_train)
print( gsearch1.best_params_, gsearch1.best_score_)
#接着我们对决策树最大深度max_depth和内部节点再划分所需最小样本数min_samples_split进行网格搜索。
param_test2 = {'max_depth':[1,2,3,5,7,9,11,13,20],'min_samples_split':[5,10,25,50,100,120,150]}
gsearch2 = GridSearchCV(estimator = RandomForestClassifier(n_estimators=160),param_grid = param_test2, cv=5)
gsearch2.fit(X_train,y_train)
print( gsearch2.best_params_, gsearch2.best_score_)
#最好的随机森林模型
best_rf=RandomForestClassifier(n_estimators=160, max_depth=20, min_samples_split=5) 

best_rf .fit(X_train, y_train)

yrfhat = best_rf .predict(X_valid)
feat_importance = best_rf.feature_importances_
rffeat_importance = pd.DataFrame([X_train.columns, random_forest.feature_importances_]).T

print('******************随机森林模型**************')
print('训练集上的准确率',best_rf.score(X_train, y_train))
print('验证集上的准确率',best_rf.score(X_valid,y_valid))
print('混淆矩阵\n', confusion_matrix(y_valid, yrfhat))
print('********************************************************')
print('最好的随机森林模型的分类报告\n', classification_report(y_valid, yrfhat))
print(rffeat_importance.sort_values(by=1, ascending=False))