【可解释机器学习】-线性回归案例【基础版】（python代码）_机器学习线性回归案例python

可解释机器学习：黑盒模型可解释性理解指南
作者: 【德】 Christoph Molnar
出版社: 电子工业出版社
译者: 朱明超

英文版（作者在持续更新）：Interpretable Machine Learning

原书github：GitHub - christophM/interpretable-ml-book: Book about interpretable machine learning

注意：原书使用R语言，本文的Python代码为我自己开发。内容基本参考原书，我进行了结构调整，并且添加了自己的理解和其他书籍的知识。

本案例模型简述

模型：最基础的多元线性归回

特征选择：无

特征交互：无

本案例数据集简述

目标变量：自行车租赁量，连续型，因此为回归问题

分类特征处理方式： 进行one-hot编码，其中未出现的取值（例如春、晴）作为参照类别 

连续特征处理方式：使用原始值（不进行标准化和归一化）

目录

线性回归介绍

第一节.准备数据集

第二节.验证数据是否满足某些假设

1.判断共线性

2.判断正态性

第三节.建模和解释

第一，将权重表的值组合为一个DataFrame

第二， 如何获取实例的真实值、预测值、置信区间，并画图

第一步：观察模型拟合效果。

第二步：用文本解释特征权重（数值型和分类型有不同的文本模板）。

第三步：解释截距。

第四步：解释特征重要性。

第五步：进一步可视化解释权重。

第六步：可视化效应图。

第七步：通过效应图解释单个实例。

---

线性回归介绍

1.线性：线性回归的最大优势（使其为可解释模型），也是它最大的劣势。线性方程在模块化水平（即权重）上有易于理解的解释，并且权重带有置信区间（权重范围的估计）。如怀疑有特征交互项或特征与目标呈非线性关系，则可考虑添加交互项或使用回归样条。

2.正态性：假设目标结果服从正态分布，如违反，则特征权重的估计置信区间无效。

3.同方差性：假设误差项的方差在整个特征空间内是恒定的。这种假设经常与现象相违背。比如预测房价，预测大别墅时的价格波动空间比普通两居室更大。

4.独立性：每个实例独立于其他实例。不独立的情况举例：对每个患者进行重复测量，此时可使用混合效应模型或广义估计方程。

5.固定特征：指的是输入的特征值没有测量误差。

6.不存在多重共线性：强相关特征会扰乱对权重的估计和解释（但这通常不会影响模型）。举例：房子面积和房间数量是强相关的，回归的结果可能是房子面积和房价的正相关，但是房间数量和房价负相关。

第一节.准备数据集

 数据集地址：UCI Machine Learning Repository: Bike Sharing Dataset Data Set

 原始数据集：

 按照作者的处理方式，对数据集用python进行处理，数据处理要点：

①未使用全部特征。

②分类变量（季节、天气情况）进行one-hot编码，其中未出现的取值（春、晴）作为参照类别。

③原始数据集中的温度、湿度、风速进行了标准化和归一化，本次处理时将取值恢复了原始值。

④链接中的原始数据集修正过季节，但是为了和书上的数据相对应，我使用了旧版的数据集，因此后续对季节的解释可能不符合实际情况。

import numpy as np
import pandas as pd
import time  #统计运行时间用
import copy  #深拷贝的时候用
import _pickle as cPickle
import gc #释放内存使用
import datetime #处理时间数据
import os
from sklearn.preprocessing import OneHotEncoder
 
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
import seaborn as sns
%matplotlib inline
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
warnings.filterwarnings("ignore")
pd.set_option('display.max_columns',100)
pd.set_option('max_colwidth',100)
 
 
path='../data/'
#原始数据集
bike=pd.read_csv(path+'Bike-Sharing-Dataset/day.csv',parse_dates=['dteday'])
#目前下载的原始数据，季节被更新了，无法与作者书中的数字保持一致，因此导入旧数据集
bike_oldseason=pd.read_csv(path+'Bike-Sharing-Dataset/bike_oldseason.csv' )
bike['season']=bike_oldseason['season'].map({'WINTER':1,'SPRING':2,'SUMMER':3,'FALL':4})
bike.head()
 
#添加特征days_since_2011，自数据集中第一天（20110101）起的天数），引入此特征为了考虑随时间变化的趋势
bike['days_since_2011']= (bike['dteday']-min(bike['dteday'])).dt.days
#原始数据集中的几个特征进行了标准化，这里还原到原始值
# denormalize weather features:
# temp : Normalized temperature in Celsius. The values are derived via (t-t_min)/(t_max-t_min), t_min=-8, t_max=+39 (only in hourly scale)
bike['temp'] = bike['temp'] * (39 - (-8)) + (-8)
# windspeed: Normalized wind speed. The values are divided to 67 (max)
bike['windspeed'] = 67 * bike['windspeed']
# hum: Normalized humidity. The values are divided to 100 (max)
bike['hum'] = 100 * bike['hum']
 
#保留的特征：租赁数量cnt(目标变量)，季节，是否假期（1是，0不是），是否工作日（1是，0不是），天气情况（晴、雾、雨雪），温度（单位为℃），湿度（相对湿度百分比0~100%），风速（单位为km/h），天数
select_cols=['cnt','season','holiday', 'workingday', 'weathersit', 'temp', 'hum', 'windspeed', 'days_since_2011']
bike=bike[select_cols]
bike.to_csv(path+'处理完数据集/bike.csv',index=False,encoding='utf_8_sig')

# OneHotEncoder：Encode categorical features as a one-hot numeric array(aka 'one-of-K' or 'dummy')
#Reshape your data either using array.reshape(-1, 1) if your data has a single feature or array.reshape(1, -1) if it contains a single sample.
#1.对特征进行独热编码（可以是'one-of-K' or 'dummy'），适用于无序的类别特征。如果对标签进行编码则使用LabelBinarizer
#2.sparse : bool, default=True，默认为稀疏矩阵形式
#3.可以通过参数handle_unknown规定fit中未出现但是transform中出现的值进行'ignore'
#4.可以通过参数categories指定需要保留的取值（指定后可以通过ohe.categories_查看）
#5.灵活使用drop参数，控制如何drop掉多余的列
#更多信息参见ohe?
def f_ohe(df,col,new_names=None,keep_cate='auto'):
    '''
    【功能】对一个特征进行独热编码
    【参数】df：dataframe,数据集
           col：str,需要进行独热编码的特征名
           new_names: dic,如果需要对取值重命名（使特征名更能表达真实意思）,则新建一个字典，默认None则特征名为col_取值
           keep_cate: list,需要保留的取值，如果取值是数值型则需要先排序,例如[[1,3,4]];默认'auto'表示保留所有值
    【返回】dataframe
    【举例】t_season=f_ohe(df=bike,col='season',new_names={1:'冬',3:'夏',4:'秋'},keep_cate=[[1,3,4]])
    
    '''
    ohe=OneHotEncoder(dtype=np.int8,handle_unknown='ignore',sparse=False,categories=keep_cate) 
    ohe.fit(df[col].values.reshape(-1,1))
    tmp=pd.DataFrame(ohe.transform(df[col].values.reshape(-1,1)))
    org_names=ohe.get_feature_names_out([col]).tolist()  #col作为新生成字段的前缀
    if new_names is None: 
        tmp.columns=org_names
    else:
        new_names_keys=list(new_names.keys()) #获取输入的keys
        new_names_keys=[col+'_'+str(item) for item in new_names_keys] #输入的keys加上col前缀
        # print(new_names_keys) 
        new_names=dict(zip(new_names_keys,list(new_names.values()))) #加上前缀的keys和输入的values重新打包为字典
        # print(new_names)        
        a=list(pd.Series(data=org_names).map(new_names).values) #list不能直接map,把原特征名map转为Series后映射为新特征名       
        tmp.columns=[col+'_'+str(item) for item in a]  #加上col前缀
 
    return tmp
 
#对于取值有两类的特征（是否假期holiday、是否工作日workingday），由于本身已是0/1编码，因此不作进一步处理
#对于取值>2类的特征（季节season、天气情况weathersit），进行one-hot编码
#季节，以春（取值2）作为参照类别，参照类别不保留
t_season=f_ohe(df=bike,col='season',new_names={1:'冬',3:'夏',4:'秋'},keep_cate=[[1,3,4]])
t_season
 
#天气情况，以晴天（取值1）作为参照类别，参照类别不保留
t_weathersit=f_ohe(df=bike,col='weathersit',new_names={2:'雾',3:'雨雪'},keep_cate=[[2,3]])
t_weathersit
 
bike_ohe=pd.concat((bike[['cnt','holiday', 'workingday', 'temp', 'hum', 'windspeed', 'days_since_2011']],
                    t_season,t_weathersit),axis=1)
bike_ohe
 
bike_ohe.to_csv(path+'处理完数据集/bike_ohe.csv',index=False,encoding='utf_8_sig')

数据处理结果如下图：

每一行代表某一天的情况。

cnt：自行车租赁数量(目标变量)

holiday：是否假期（1是，0不是)

workingday：是否工作日（1是，0不是），除了假期和双休日的为工作日

temp：温度（单位为℃）

hum：湿度（相对湿度百分比0~100%）

windspeed：风速（单位为km/h）

days_since_2011：自数据集中第一天（20110101）起的天数，引入此特征为了考虑随时间变化的趋势

season_冬、season_夏、season_秋：季节特征（1是，0不是），其中春未出现，可由其他三个季节推知，将春作为参照类别

weathersit_雾、weathersit_雨雪：天气特征（1是，0不是），其中晴未出现，可由其他两个天气情况推知，将晴作为参照类别

第二节.验证数据是否满足某些假设

1.线性：可先不用管。根据基础线性模型建立后的R方判断预测效果。后续如有需要可对基础线性模型进行改造（如添加交互项、使用回归样条等）。

2.正态性：本节对目标变量的正态性进行了检验，检验通过。

3.同方差性：必须建模后再判断，此时无法判断。

4.独立性：除了特殊情况（比如对同一人进行重复测量），一般可认为实例之间相互独立。

5.固定特征：测量无误差，一般都可认为是固定特征。

6.不存在多重共线性：本节进行了检验，检验通过。虽然同时计算了相关系数和VIF，但实践时可以直接以VIF的结论为准。

下面提供检验多重共线性和正态性的代码： 

df=bike_ohe
label='cnt'
 
feas=df.columns.tolist()
feas.remove(label)
print('特征数量:',len(feas))
 
X_train=df[feas]
y_train=df[label]
print('X_train:',X_train.shape)
print('y_train:',y_train.shape)

1.判断共线性

（强相关特征会扰乱对权重的估计和解释（但这通常不会影响模型））
1.1.计算相关系数（只能判断两两的线性关系）  

#1.1.计算相关系数（只能判断两两的线性关系）  
import matplotlib.pyplot as plt
import seaborn as sns
 
corr_res=X_train.corr()
plt.subplots(figsize=(9, 9))
sns.heatmap(corr_res, annot=True, vmax=1, square=True, cmap="coolwarm")
plt.xticks(fontsize=15) #放大横纵坐标刻度线上的特征名字体
plt.yticks(fontsize=15)
plt.show()

如下图，绝对值最大的相关系数是秋和温度（0.68），可接受。 

1.2.计算方差膨胀因子（可以判断三个或更多变量之间的线性关系） 

# 1.2.计算方差膨胀因子（可以判断三个或更多变量之间的线性关系）
from statsmodels.stats.outliers_influence import variance_inflation_factor #计算方差膨胀因子
from statsmodels.tools.tools import add_constant  #添加常量
 
def checkVIF(df):
    '''
    【功能】计算方差膨胀因子
    【参数】df:dataframe,特征集（不含target）
    【返回】dataframe，展示各个特征的VIF
    【参考】当0<VIF<10，不存在多重共线性；当10≤VIF<100，存在较强的多重共线性；当VIF≥100，存在严重多重共线性。
    【来源与介绍】https://blog.csdn.net/nixiang_888/article/details/122342338
    【举例】VIF1=checkVIF(X_train)
    '''
    df = add_constant(df) #添加一列常量const作为截距，全部赋值为1（不会改变原数据集）
    name = df.columns
    x = np.matrix(df)
    VIF_list = [variance_inflation_factor(x,i) for i in range(x.shape[1])]
    VIF = pd.DataFrame({'feature':name,"VIF":VIF_list})
    VIF = VIF.drop(index=0) #删除截距const行
    VIF.sort_values(['VIF'],ascending=False,inplace=True)
    VIF['remark']=np.where(VIF['VIF']>=100,'严重多重共线性',np.where(VIF['VIF']>=10,'较强多重共线性','无多重共线性'))
    
    return VIF
 
VIF1=checkVIF(X_train)
VIF1

结果如下图，不存在多重共线性。

2.判断正态性

（假设目标结果服从正态分布，如违反，则特征权重的估计置信区间无效 ）

def checkNORM(se,p=0.05,alt='two-sided',if_plot=True):
    '''
    【功能】检验一组数据是否符合正态分布
    【参数】se:Series
           p：float,p值，默认0.05
           alt：str,默认双侧检验'two-sided'，可选'less', 'greater'
           if_plot,是否画图，默认True
    【返回】dataframe，展示各个特征的VIF
    【参考】结果返回两个值：statistic → D值，pvalue → P值
    【备注】import matplotlib.pyplot as plt
           %matplotlib inline 
    【举例】 res=checkNORM(y_train)
    '''
    from scipy import stats
 
    print('数据量：',len(se))
    
    u = se.mean()  # 计算均值
    std = se.std()  # 计算标准差
    res=stats.kstest(rvs=se, cdf='norm',args= (u, std), alternative=alt)
    print('p值为:',res[1])
    if res[1]>p:
        print('p值>',p,'符合正态分布')
    else:
         print('p值<=',p,'不符合正态分布')
   
    if if_plot==True:
        fig = plt.figure(figsize = (10,6))         
        se.hist(bins=30,alpha = 0.5) #直方图 alpha表示透明度
        se.plot(kind = 'kde', secondary_y=True) #核密度估计KDE
        plt.show()
 
 
    return res
 
res=checkNORM(y_train)

结果如下图，符合正态分布。

第三节.建模和解释

本案例的基础线性模型可以选择使用1）statsmodels库，2）sklearn库。其中statsmodels库的信息更多、更全。

本文先用statsmodels库进行详细解释，文末提供sklearn库的代码。

from statsmodels.regression.linear_model import OLS,GLS #Ordinary least squares普通最小二乘法
import statsmodels.formula.api as smf
import statsmodels.api as sm
 
#建模方式1：使用smf.ols，自己编写formula，会自动添加常数列
#cnt为目标变量，分类特征可使用C(season)进行编码，由于本数据集的分类特征都已事先编码，因此不需要添加c()
model=smf.ols(formula='cnt ~  season_夏 +  season_秋 + season_冬 + holiday + workingday +weathersit_雾 + weathersit_雨雪 + temp + hum + windspeed +days_since_2011 ',data=df)
results=model.fit()
results.summary()

 

#建模方式2：使用sm.OLS，需要先sm.add_constant添加常数列以计算截距
X = sm.add_constant(X_train)
model=sm.OLS(endog=y_train,exog=X)
results=model.fit()
results.summary()

本文采用#建模方式1：使用smf.ols。

在进行模型解释前，进行两个操作。

第一，将权重表的值组合为一个DataFrame

df_coef=pd.DataFrame(results.params ) #权重
df_coef.reset_index(inplace=True)
df_coef.columns=['feature','coef']
df_coef['lw']=results.conf_int(alpha=0.05)[0].values #获取权重的置信区间下限
df_coef['up']=results.conf_int(alpha=0.05)[1].values #获取权重的置信区间上限
df_coef['SE']=results.bse.values #权重的标准误std err
df_coef['t']=results.tvalues.values #权重的t统计量，等于权重/标准误
df_coef['p']=results.pvalues.values #参数的t统计的双尾 p 值
df_coef['t_abs']=abs(df_coef['t']) #求绝对值
#根据已有的权重和置信区间计算上下误差，计算完毕后发现上下误差相同
df_coef['lw_err']=df_coef['coef']-df_coef['lw'] 
df_coef['up_err']=df_coef['up']-df_coef['coef']
df_coef

第二， 获取实例的真实值、预测值、置信区间，并画图

#获取置信区间的上下限
pred_ols = results.get_prediction()
iv_l = pred_ols.summary_frame()["obs_ci_lower"]
iv_u = pred_ols.summary_frame()["obs_ci_upper"]
 
#results.fittedvalues为模型预测值
target_df=pd.concat((y_train,results.fittedvalues,iv_l,iv_u),axis=1)
target_df.columns=['true','predict','ci_lower','ci_upper']
target_df['resid']=results.resid #残差
target_df
 
#按实际租赁量排序，reset_index是必须的
plot_df=target_df.sort_values(['true']).reset_index(drop=True)
 
fig, ax = plt.subplots(figsize=(20, 9))
 
ax.plot(plot_df['true'], "b-", label="True")
ax.plot(plot_df['predict'], "r", label="Pred")
# ax.plot(plot_df['ci_lower'], "r--",alpha=0.5) #置信区间虚线
# ax.plot(plot_df['ci_upper'], "r--",alpha=0.5) #置信区间虚线
plt.fill_between(plot_df.index,plot_df['ci_lower'],plot_df['ci_upper'],color='blue',alpha=0.15)
ax.legend(loc="best")
 
 
plt.ylabel('自行车租赁量',fontsize=18)
plt.title('真实值与预测值对比',fontsize=20)
plt.show()

上图中，蓝线为真实值，红色实线为预测值，蓝紫色为置信区间。

由上图可知，左侧租赁量较小时，部分预测值远高于真实值且波动较大；右侧租赁量较大时，预测值整体偏低。

由此图也可以看出，前文中提到的线性回归模型的【同方差性】在现实中是很难满足的。

本案例数据量较小，如果数据量较大，可以随机抽样后再画图。 

#残差——同方差性
#1.应该为均值是0的正态分布
sns.set(style="whitegrid",font_scale=1.2)#设置主题，文本大小
plt.hist(target_df['resid'])
plt.show()

 

#2.残差与predict之间应该不相关
#regplot默认参数线性回归图
plt.figure(figsize=(8, 8))
sns.set(style="whitegrid",font_scale=1.2)#设置主题，文本大小
g=sns.regplot(x='resid', y='predict', data=target_df,
             color='#000000',#设置marker及线的颜色
             # marker='*',#设置marker形状
             )

 

下面进行模型解释。

第一步：观察模型拟合效果。

观察调整后R方（解释一个该值很低的模型是没有意义的，对权重的任何解释都没有意义）。

该模型的调整后R方是0.790，表示该模型解释了目标结果79%的总方差，拟合度较优，可解释程度较高。

第二步：用文本解释特征权重（数值型和分类型有不同的文本模板）。

数值特征文本模板：当所有其他特征保持不变时，特征x增加一个单位，则预测结果y增加β。

分类特征文本模板：当所有其他特征保持不变时，将特征x从参照类别改变为其他类别时，预测结果y会增加β。

观察权重（图中coef列）（由于本数据集使用了原始值，即未进行标准化和归一化，因此可以直接进行表述）：

温度（数值特征）：当所有其他特征保持不变时，将温度升高1℃，自行车的预测数量增加110辆。

风速（数值特征）：当所有其他特征保持不变时，风速增加1km/h，自行车的预测数量减少42辆。

天气情况（分类特征）：所有其他特征保持不变，雨雪天气的自行车的预测数量比晴天减少了1901辆；天气为雾时，自行车的预测数量比晴天减少了379辆。

季节（分类特征）：所有其他特征保持不变，夏天自行车的预测数量比春天增加了899辆，秋天自行车的预测数量比春天增加了138辆，冬天自行车的预测数量比春天增加了425辆。

是否假期（分类特征）：所有其他特征保持不变，假期的自行车预测数量比非假期减少了425辆。

第三步：解释截距。

截距权重：对所有数值特征为0和分类特征为参照类别的实例，模型预测值即为截距的权重。上述解释通常没有意义（特征全部为0 的实例通常无实际含义）。但是，当特征标准化（均值为0，标准差为1）时，这种解释将会有实际含义，此时截距反应所有特征都处于均值时实例的预测结果。

本例中，截距的权重为2399，表示当处于春天、晴天、假期、工作日，温度、湿度、风速都为0，且为2011年1月1号时，预测的自行车数量为2399辆。以上解释无实际意义。

第四步：解释特征重要性。

使用t-统计量的绝对值解释特征重要性，（t=权重/SE，其中SE是标准误）。特征的重要性随着权重的增加而增加，随着方差的增加而减小（方差越大表明对正确值的把握越小）。

本例中，t统计量已经被计算出来了，上图中的 t=coef/std err。

进行数据处理（求绝对值、排序）后画图。

plot_df=df_coef.drop(index='Intercept') #删除截距行
plot_df=plot_df.sort_values(['t_abs']) #排序
 
fig = plt.figure(figsize = (9,5))
plt.barh(plot_df.index,plot_df['t_abs']) #画水平条形图
 
#设置x轴y轴
plt.xlabel('t-value绝对值',fontsize=18)
plt.ylabel('特征',fontsize=18)
plt.xticks(fontsize=12) #放大横纵坐标刻度线上的特征名字体
plt.yticks(fontsize=12)
plt.title('特征重要性',fontsize=20)
plt.show()

由上图可知，最重要的特征为距离2011年第一天的天数，排名第二的为温度，排名第三的为是否是雨雪天气。

问题来了，为什么不能直接使用权重（coef）代表特征重要性呢？

其中一个原因：如果改变特征的量纲，权重就会发生变化。

为进行验证，将风速的单位从km/h转换为km/分钟，即将风速的原始值除以60，得出以下结果：

变化：其他不变，仅将风速除以60

结论：

其他特征的估计均不变。风速的部分估计发生了变化

权重（coef）：扩大了60倍

权重标准误（std err）：扩大了60倍

t统计量：不变（t=权重/标准误）

p值：不变

置信区间：扩大了60倍

所以在风速的本质并未发生改变的情况下，如果采用权重（coef）作为特征重要性的度量依据，就会发现其值会随着量纲的变化而变化，但t统计量却可以保持一致。

第五步：进一步可视化解释权重。

第二步已通过文本解释了权重（coef）的实际含义，这一步根据权重和置信区间画权重图。

plot_df=df_coef.drop(index='Intercept') #删除截距行
fig = plt.figure(figsize = (8,8))
#由于上下误差相同，因此直接用 xerr=plot_df['lw_err']，否则可以使用xerr=plot_df[['lw_err','up_err']].T.values来分别规定上下限
plt.errorbar(x=plot_df['coef'], y=plot_df.index,xerr=plot_df['lw_err'], color="black", capsize=3,
             linestyle="None",
             marker="s", markersize=7, mfc="black", mec="black")
 
plt.grid(True) #显示网格线
 
plt.xlabel('权重估计',fontsize=18)
plt.ylabel('特征',fontsize=18)
plt.xticks(fontsize=12) #放大横纵坐标刻度线上的特征名字体
plt.yticks(fontsize=12)
plt.title('权重估计图',fontsize=20)
 
plt.axvline(c="c",ls="--",lw=2) #原点竖线
plt.show()

由上图可知：

1.雨雪天气对自行车租赁量有很大的负效应。

2.是否工作日的权重接近于0，并且95%的置信区间包含0，这表明该效应在统计上不显著。

3.温度的置信区间很短，估计值接近于0，但特征效应在统计上是显著的。

权重图的问题：

各个特征的量纲不一样，比如天气情况反映了晴天和雨雪天的差异，但是温度只反映了1℃的变化情况。

因此可以通过在建模前对特征进行标准化（均值为0，标准差为1），使估计的权重更具有可比性。

第六步：可视化效应图。

效应图（effect plot）帮助了解权重和特征的组合对数据预测的贡献程度。特征效应为每个特征的权重乘以实例的特征值。如改变特征的量纲，则权重会发生变化，但特征效应不会改变。

通过画箱线图（注意，分类特征总结为一个箱线图），可以观察下面几个方面：1）特征效应的正负性；2）特征效应的绝对值大小；3）特征效应的方差（如果方差小，则意味着这个特征几乎在所有实例中都有类似的贡献）。

#求特征效应——每个特征的权重乘以实例的特征值
w=df_coef['coef'].values
w_order=[] #将特征权重与实例中的顺序一一对应
my_dict={0:4,1:5,2:8,3:9,4:10,5:11,6:3,7:1,8:2,9:6,10:7} #权重表与数据集中特征的对应顺序
for i in range(11):
    w_order.insert(i,w[my_dict[i]])
    
#计算特征效应    
effect=X_train*w_order 
 
#分类特征合并
effect['season']=np.sum(effect[['season_冬','season_夏','season_秋']],axis=1)
effect['weathersit']=np.sum(effect[['weathersit_雾','weathersit_雨雪']],axis=1)
effect

plt.subplots(figsize=(9, 9))
 
cols=['holiday', 'workingday', 'temp', 'hum', 'windspeed', 'days_since_2011',
       'season', 'weathersit']
sns.boxplot(data=effect[cols],orient="h",width=0.5,whis=0.5, palette="Set2")
 
plt.grid(True) #显示网格线
 
plt.xlabel('特征效应',fontsize=18)
plt.ylabel('特征',fontsize=18)
plt.xticks(fontsize=12) #放大横纵坐标刻度线上的特征名字体
plt.yticks(fontsize=12)
plt.title('特征效应图',fontsize=20)
 
plt.axvline(c="c",ls="--",lw=2) #原点竖线
plt.show()

 

 由上图可知：

1.对预测自行车租赁数量正向贡献最大的来自温度和天数。

2.天气的情况参照类别为晴天，图中说明除了晴天外的天气（雾、雨雪）都会对自行车租赁量产生负向影响。

第七步：通过效应图解释单个实例。

通过计算单个实例的特征效应，可以得到这个实例的各个特征对预测有多大的贡献。

步骤1：得到这个实例的预测值、所有实例的平均预测值、这个实例的实际值。将这个实例的预测值与所有实例的平均预测值进行对比，如果差距很大，则进一步解释原因。

步骤2：计算这个实例的特征效应，然后加入第六步的特征效应图中。即将这个实例的特征效应与所有实例的特征效应分布进行比较，得出结论。

single_idx=5 #第6个实例
print(bike_ohe.loc[single_idx])
 
target_predict=target_df.loc[single_idx,'predict']
target_predict_mean=np.mean(target_df['predict'])
target_true=target_df.loc[single_idx,'true']
print('该实例预测值',target_predict)
print('所有实例平均预测值',target_predict_mean)
print('该实例实际值',target_true)

plt.subplots(figsize=(9, 9))
 
cols=['holiday', 'workingday', 'temp', 'hum', 'windspeed', 'days_since_2011',
       'season', 'weathersit']
sns.boxplot(data=effect[cols],orient="h",width=0.5,whis=0.5, palette="Set2")
 
plt.grid(True) #显示网格线
 
plt.xlabel('特征效应',fontsize=18)
plt.ylabel('特征',fontsize=18)
plt.xticks(fontsize=12) #放大横纵坐标刻度线上的特征名字体
plt.yticks(fontsize=12)
plt.title('单个实例的特征效应图',fontsize=20)
 
plt.axvline(c="c",ls="--",lw=2) #原点竖线
 
 
 
#画单个实例中每个特征的效应
for col in cols:
    col_index=cols.index(col) #获取某个特征在特征名列表的索引位置
    plt.plot(effect.loc[single_idx,col], col_index,'rx', ms=10)  #rx 红色叉号，ms控制大小
    
plt.show()

 

以数据集中第6个实例为例：
相较于所有实例的平均预测值4504辆，该实例的预测值很小，只有1571辆自行车被租赁。

效应图揭示了原因：

1.该实例温度的特征效应较小，这一天温度仅为1.6℃，与其他大多数日期的温度相比较低（温度权重为正）。

2.该实例天数的特征效应也较小，该实例自第一天起仅过了5天（天数权重为正）。

也可使用sklearn库进行建模，但是该库提供的信息十分有限，如下图：

---

完毕。

【可解释机器学习】-线性回归案例【lasso特征选择版】（python）

基础版线性模型的参照实例为构造的一个数据点，其中所数值型特征为0，分类型特征为参照类别，这通常是一个毫无意义的实例，不可能存在于现实中。

如果将所有数值特征以均值为中心（特征原始值减去该特征平均值），并且所有分类特征都采用效应编码，那么参照实例就是所有特征取平均特征值的数据点。这个数据点可能也并不存在，但可能会更有意义。

### 分类型：one-hot编码 →  效应编码
### 连续型：原值 →  原值-均值
### 结果有如下变化： 

相关系数：连续型+二分类型：不变；多分类型：变化  
VIF值：连续型+二分类型：不变；多分类型：变化  

目标变量正态性：未处理目标变量，因此不变
模型预测效果：R方、实例预测情况不变  
连续型：权重等均不变  
分类型：权重、t值等均发生变化（且排序也会发生变化）。由于效应编码中，总体均值为各个类别（目标）均值的均值，因此解释起来会比较困难，比如雾的系数是正的（表示大于均值），因为雨雪天很少，直接计算总体的目标均值和计算各个类别目标均值的均值结果非常不同。  
建议还是使用one-hot编码。  

效应编码参考：https://stats.oarc.ucla.edu/other/mult-pkg/faq/general/faqwhat-is-effect-coding/   
参照类别全赋值为-1，其余与one-hot相同。参照类别可由其他类别推得。

附图（one-hot编码 →  效应编码；原值 →  原值-均值）