Python 根据AIC准则定义向前逐步回归进行变量筛选（二）
 
AIC简介
 
AIC即赤池值，是衡量模型拟合优良性和模型复杂性的一种标准，在建立多元线性回归模型时，变量过多，且有不显著的变量时，可以使用AIC准则结合逐步回归进行变量筛选。AICD数学表达式如下：
 $$AIC=2p+n(log(SSE/n))$$
 其中，$p$是进入模型当中的自变量个数，$n$为样本量，$SSE$是残差平方和，在$n$固定的情况下，$p$越小，$AIC$越小，$SSE$越小，$AIC$越小，而$p$越小代表着模型越简洁，$SSE$越小代表着模型越精准，即拟合度越好，综上所诉，$AIC$越小，即模型就越简洁和精准。
 
逐步回归
 
逐步回归分为三种，分别是向前逐步回归，向后逐步回归，逐步回归。向前逐步回归的特点是将自变量一个一个当如模型中，每当放入一个变量时，都利用相应的检验准则检验，当加入的变量不能使得模型变得更优良时，变量将会被剔除，如此不断迭代，直到没有适合的新变量加入为止。向后逐步回归的特点是，将所有变量都放入模型之后，一个一个的剔除变量，将某一变量拿出模型而使得模型更优良时，将会剔除此变量。如此反复迭代，直到没有合适的变量剔除为止。逐步回归则是结合了以上的向前和向后逐步回归的特点。
 
继上一篇的代码
 
import numpy as np
import pandas as pd
from sklearn.datasets import fetch_california_housing as fch  #加载加利福尼亚房屋价值数据
#加载线性回归需要的模块和库
import statsmodels.api as sm #最小二乘
from statsmodels.formula.api import ols #加载ols模型
data=fch() #导入数据
house_data=pd.DataFrame(data.data) #将自变量转换成dataframe格式，便于查看
house_data.columns=data.feature_names  #命名自变量
house_data.loc[:,"value"]=data.target #合并自变量，因变量数据
house_data.shape #查看数据量
house_data.head(10) #查看前10行数据
#分训练集测试集
import random
random.seed(123) #设立随机数种子
a=random.sample(range(len(house_data)),round(len(house_data)*0.3))
house_test=[]
for i in a:
    house_test.append(house_data.iloc[i])
house_test=pd.DataFrame(house_test)
house_train=house_data.drop(a)
 
定义向前逐步回归函数
 
#定义向前逐步回归函数
def forward_select(data,target):
    variate=set(data.columns)  #将字段名转换成字典类型
    variate.remove(target)  #去掉因变量的字段名
    selected=[]
    current_score,best_new_score=float('inf'),float('inf')  #目前的分数和最好分数初始值都为无穷大（因为AIC越小越好）
    #循环筛选变量
    while variate:
        aic_with_variate=[]
        for candidate in variate:  #逐个遍历自变量
            formula="{}~{}".format(target,"+".join(selected+[candidate]))  #将自变量名连接起来
            aic=ols(formula=formula,data=data).fit().aic  #利用ols训练模型得出aic值
            aic_with_variate.append((aic,candidate))  #将第每一次的aic值放进空列表
        aic_with_variate.sort(reverse=True)  #降序排序aic值
        best_new_score,best_candidate=aic_with_variate.pop()  #最好的aic值等于删除列表的最后一个值，以及最好的自变量等于列表最后一个自变量
        if current_score>best_new_score:  #如果目前的aic值大于最好的aic值
            variate.remove(best_candidate)  #移除加进来的变量名，即第二次循环时，不考虑此自变量了
            selected.append(best_candidate)  #将此自变量作为加进模型中的自变量
            current_score=best_new_score  #最新的分数等于最好的分数
            print("aic is {},continuing!".format(current_score))  #输出最小的aic值
        else:
            print("for selection over!")
            break
    formula="{}~{}".format(target,"+".join(selected))  #最终的模型式子
    print("final formula is {}".format(formula))
    model=ols(formula=formula,data=data).fit()
    return(model)
 
利用向前逐步回归筛选变量
 
forward_select(data=house_train,target="value")
 

 最终保留了7个自变量，其实只是剔除了一个自变量，将这七个自变量放进模型里再运行一遍，查看模型结果
 
lm_1=ols("value~MedInc+HouseAge+Latitude+Longitude+AveBedrms+AveRooms+AveOccup",data=house_train).fit()
lm_1.summary()
 

 由以上结果可以看出，其实模型的R方几乎没有变化。
 本篇文章主要是想讲述如何利用statsmodels和AIC准则定义向前逐步回归函数筛选自变量，在日后遇到比较多自变量的时候，方便进行自变量筛选。
 
参考文献
 
常国珍，赵仁乾，张秋剑.Python数据科学技术详解于商业实战[M]. 北京：中国人民大学出版社，2018.

 
 

  
通用文件处理：

  
import numpy as np
//文件名和文件中每行的分隔符
def loadDataSet(fileName,dotSplit):
    numFeat = len(open(fileName).readline().split(dotSplit))
    dataMat = []; labelMat = []
    fr = open(fileName)
//该数据集默认是最后一列是因变量
    for line in fr.readlines():
        lineArr=[]
        curline = line.split(dotSplit)
        for i in range(0,numFeat-1):
            lineArr.append(float(curline[i]))
        dataMat.append(lineArr)
        labelMat.append(float(curline[numFeat-1]))

    xMat = np.mat(dataMat)
    yMat = np.mat(labelMat).T
    return xMat,yMat
  
这里是处理岭回归的实现：

  
import numpy as np

def ridgeRegres(xMat,yMat,lam=0.2):
    xTx = xMat.T*xMat
    denom = xTx + np.eye(np.shape(xMat)[1])*lam
    print np.shape(xMat)[0]
    if np.linalg.det(denom) == 0.0:
        print "wrong"
        return
    ws = denom.I*(xMat.T*yMat)
    return ws

def normalizing(xMat,yMat):
    yMean = np.mean(yMat,0)
    y = yMat-yMean
    xMeans = np.mean(xMat,0)
    xVar = np.var(xMat,0)
    x = (xMat-xMeans)/xVar
    return x,y

def ridgeTest(xM,yM):
    xMat,yMat = normalizing(xM,yM)
    numTestPts = 30
    wMat = np.zeros((numTestPts,np.shape(xMat)[1]))
    print wMat
    for i in range(numTestPts):
        ws = ridgeRegres(xMat,yMat,np.exp(i-10))
        wMat[i,:] = ws.T
    return wMat
  
向前逐步回归：
  
import numpy as np

def rssError(yArr,yHatArr):
    return ((yArr-yHatArr)**2).sum()

def stageWise(xM,yM,eps=0.01,numIt=100):
    m,n = np.shape(xM)
    returnMat = np.zeros((numIt,n))
    ws = np.zeros((n,1));wsTest = ws.copy();wsMax = ws.copy()
    lowestError = 0
    for i in range(numIt):
        print ws.T
        for j in range(n):
            for sign in [-1,1]:
                wsTest = ws.copy()
                wsTest[j] += eps*sign
                yTest = xM*wsTest
                rssE = rssError(yM.A,yTest.A)
                if i == 0:
                    lowestError = rssE
                if rssE < lowestError:
                    lowestError = rssE
                    print lowestError
                    wsMax = wsTest
        ws = wsMax.copy()
        returnMat[i,:] = ws.T
    return returnMat
  


 
 
 
转载于:https://blog.51cto.com/wukong0716/1693328

逐步回归算法
 
监督式学习算法的任务是拟合多个特征与标签的关系，在实际应用中，所收集的数据中，所有的特征可能并不是都与标签有关联。模型中如果包含与标签无关联的特征，不仅会增加数据规模和计算量，还会影响模型对标签的预测效果。因此特征选择是监督式学习算法的一个重要组成部分。
 逐步回归是一个贪心算法。它的运行效率较高，但也因为其是贪心算法，有时会做出次优的最优选择。
 
向前逐步回归算法
 
向前逐步回归算法是最简单的一种特征选择方法，大概做法为：在向前逐步回归算法的初始阶段，先选定第一个特征，然后重复执行以下几个步骤----首先计算只使用当前选定的特征的线性回归的均方误差，然后逐一引入尚未选取的特征，选择能最大程度降低均方误差的一个特征，判断该特征是否在统计意义上显著地降低均方误差，如果是，就将该特征加入模型。重复循环上述过程，直至没有能够被继续选中的特征为止。
 
判断选择的特征是否在统计意义上显著降低均方误差
 
在向前逐步回归算法中，采用F检验来判断均方误差的减小是否具有统计显著性。
 给定两个均方误差mse1 和mse2，设mse1>mse2。用F检验计算mse1>mse2的置信度p。
 置信度p是在重新采样训练数据并对其重复向前逐步选择算法时再次出现mse1>mse2的概率。
 如果置信度p>95%，则认为mse1>mse2这一结论具有统计显著性。
 
import numpy as np
from scipy.stats import f
# 向前逐步选择
class StepwiseRegression():
    def fit(self,X,y):
        return np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)

    # 计算均方误差
    def compute_mse(self,X,y):
        w = self.fit(X,y)
        r = y - X.dot(w)
        return r.T.dot(r)

    # 判断引入新特征是否有效降低均方误差
    def f_tesy(self,mse_A,mse_min,m):
        if mse_min > mse_A:
            return False
        F = mse_A/mse_min
        p_value = f.cdf(F,m,m)
        return p_value>0.95


    # 特征选择
    def forward_selection(self,X,y):
        m, n = X.shape
        # 记录特征列数，A为已经选择的特征，C为备选特征
        A, C =[0],[i for i in range(1,n)]
        while len(C)>0:
            # 取训练数据中的A包含的特征
            MSE_A = self.compute_mse(X[:,A],y)
            MSE_min = float("inf")
            j_min = -1
            #遍历未选择的特征，保留均方误差最小值以及代表的特征列
            for j in C:
                MSE_j = self.compute_mse(X[:,A+[j]],y)
                if MSE_j < MSE_min:
                    MSE_min ,j_min = MSE_j,j
            if self.f_tesy(MSE_A,MSE_min,m):
                A.append(j_min)
                C.remove(j_min)
            else:
                break
        self.w = self.fit(X[:,A],y)
        self.A = A


    # 计算标签值
    def predict(self,X):
        return X[:,self.A].dot(self.w)
 
向前逐步回归算法的应用
 
# 应用向前逐步回归算法进行特征选择

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
from machine_learning_lib_stepwise_regression import StepwiseRegression
import matplotlib.pyplot as plt
from machine_learning_lib_stepwise_regressiontwo import Stepwiseregression as Step

# 构造数据
def generate_sample(m):
    X = 2 * (np.random.rand(m,1)-0.5)
    y = X + np.random.normal(0,0.3,(m,1))
    return X,y

np.random.seed(100)
# 生成10个点
X,y = generate_sample(10)
plt.scatter(X,y)
poly = PolynomialFeatures(degree=10)
X_poly = poly.fit_transform(X)
model = StepwiseRegression()
model.forward_selection(X_poly,y)
y_pred = X_poly[:,[0,1]].dot(model.w)
plt.plot(X,y_pred)
print(model.A,model.w)
plt.show()
 
为了直观的观察特征选择的效果，训练数据是通过随机函数产生的，并对其多项式化，产生11个特征，标签以y=x，赋予一定程度的浮动人为产生。
 
该训练数据如果不进行特征的选择，直接对全部特征进行线性回归，将产生过度拟合的线性回归模型。
 
未使用特征选择：
 
 
使用特征选择：
 
 
总结：
 
特征选择能有效的剔除与标签无关联的特征，避免无关特征对模型的影响，导致过度的拟合，但向前逐步回归法是一个贪心算法，在算法的开始，选择了第一个特征为起始，但第一个特征与标签的关系并未进行判断，若第一个特征与标签无关联，则将会产生次优解。
 
思考:
 
向后逐步回归算法可以有效解决向前逐步回归算法的影响，向后逐步回归在起始的时候选取全部特征，之后逐一 剔除不能有效降低均方误差的特征。但是向后逐步回归算法也是一个贪心算法，依然可能产生次优解，并且当前选定的特征已经过度拟合，从而具有较小的均方误差，向后逐步回归算法将无法继续从模型中剔除与标签无关的特征。
 
处理方案：
 
向前逐步回归算法与向后逐步回归算法的合并使用。
 向前逐步回归算法起始选择随机，多次运行，选择最优情况。

Python的statsmodels包含了一些R风格的统计模型和工具。在内部实现上，statsmodels使用patsy包将数据转化为矩阵并建立线性模型，具体信息参见pasty主页http://patsy.readthedocs.io/en/latest/overview.html。但是，Python的statsmodels工具中没有向前逐步回归算法。逐步回归的基本思想是将变量逐个引入模型，每引入一个解释变量后都要进行F检验，并对已经选入的解释变量逐个进行t检验，当原来引入的解释变量由于后面解释变量的引入变得不再显著时，则将其删除。以确保每次引入新的变量之前回归方程中只包含显著性变量。这是一个反复的过程，直到既没有显著的解释变量选入回归方程，也没有不显著的解释变量从回归方程中剔除为止。以保证最后所得到的解释变量集是最优的（https://baike.baidu.com/item/%E9%80%90%E6%AD%A5%E5%9B%9E%E5%BD%92/585832?fr=aladdin）。
网上有人用statsmodels写了一个向前逐步回归的工具，具体网址见https://planspace.org/20150423-forward_selection_with_statsmodels/。 我试了一下，速度还不错，比我用sklearn写的要好。具体代码如下：
import statsmodels.formula.api as smf
import pandas as pd

def forward_selected(data, response):
    """前向逐步回归算法，源代码来自https://planspace.org/20150423-forward_selection_with_statsmodels/
    使用Adjusted R-squared来评判新加的参数是否提高回归中的统计显著性
    Linear model designed by forward selection.

    Parameters:
    -----------
    data : pandas DataFrame with all possible predictors and response

    response: string, name of response column in data

    Returns:
    --------
    model: an "optimal" fitted statsmodels linear model
           with an intercept
           selected by forward selection
           evaluated by adjusted R-squared
    """
    remaining = set(data.columns)
    remaining.remove(response)
    selected = []
    current_score, best_new_score = 0.0, 0.0
    while remaining and current_score == best_new_score:
        scores_with_candidates = []
        for candidate in remaining:
            formula = "{} ~ {} + 1".format(response,
                                           ' + '.join(selected + [candidate]))
            score = smf.ols(formula, data).fit().rsquared_adj
            scores_with_candidates.append((score, candidate))
        scores_with_candidates.sort()
        best_new_score, best_candidate = scores_with_candidates.pop()
        if current_score < best_new_score:
            remaining.remove(best_candidate)
            selected.append(best_candidate)
            current_score = best_new_score
    formula = "{} ~ {} + 1".format(response,
                                   ' + '.join(selected))
    model = smf.ols(formula, data).fit()

    return model


def main():

    '''
    首先从网上读取普林斯顿大学52位员工的工资信息，工资文件一共有6列，各列意义解释如下：
    sx = Sex, coded 1 for female and 0 for male
    rk = Rank, coded
        1 for assistant professor,
        2 for associate professor, and
        3 for full professor
    yr = Number of years in current rank
    dg = Highest degree, coded 1 if doctorate, 0 if masters
    yd = Number of years since highest degree was earned
    sl = Academic year salary, in dollars.

    '''
    url = "http://data.princeton.edu/wws509/datasets/salary.dat"
    data = pd.read_csv(url, sep='\\s+')

    #将sl（年收入）设为目标变量
    model = forward_selected(data, 'sl')

    #打印出最后的回归模型
    print(model.model.formula)
    # sl ~ rk + yr + 1
    print(model.params)
    # Intercept          16203.268154
    # rk[T.associate]     4262.284707
    # rk[T.full]          9454.523248
    # yr                   375.695643
    # dtype: float64
    # 0.835190760538

    print(model.rsquared_adj)
    # 0.835190760538

if __name__ == '__main__':
    main()


最好，该方法也有缺点。随着模型中参数增多，Adjusted R-squared的增加不太明显，所以后来加的很多参数对模型的影响不大。下面是在我自己的数据集（包含98个参数）上做的测试图，我们看到，当参数个数多于7个时，再增加参数，模型的R2并不能显著增加，所以最终模型取前7个参数即可。