Python中岭回归中的方差膨胀因子

Question

Python中岭回归中的方差膨胀因子

9

我正在对一些高度共线的数据进行岭回归分析。其中一个用于确定稳定拟合的方法是利用岭迹图。感谢在Scikit-learn上提供的优秀示例，我能够做到这一点。另一种方法是计算每个变量的方差膨胀因子(VIF)，随着k的增加而增加。当VIF降至<5时，表明拟合效果令人满意。Statsmodels有关于VIF的代码，但是它是针对OLS回归的。我尝试修改它以处理岭回归。

我正在检查我的结果是否与《Regression Analysis by Example》第五版第10章所述的一致。我的代码可以生成k=0.000时的正确结果，但之后就不行了。虽然有SAS工作代码可供使用，但我不是SAS用户，也不知道其实现与scikit-learn和/或statsmodels的区别。

我已经卡住了几天，任何帮助都将不胜感激。

#http://www.ats.ucla.edu/stat/sas/examples/chp/chp_ch10.htm

from __future__ import division
import numpy as np
import pandas as pd
example = pd.read_csv('by_example_import.csv')
example.dropna(inplace=True)

from sklearn import preprocessing
scaler = preprocessing.StandardScaler().fit(example)
scaler.transform(example)

X = example.drop(['year', 'import'], axis=1)
#c_matrix = X.corr()
y = example['import']
#w, v = np.linalg.eig(c_matrix)

import pylab as pl
from sklearn import linear_model

###############################################################################
# Compute paths

alphas = [0.000, 0.001, 0.003, 0.005, 0.007, 0.009, 0.010, 0.012, 0.014, 0.016, 0.018,
          0.020, 0.022, 0.024, 0.026, 0.028, 0.030, 0.040, 0.050, 0.060, 0.070, 0.080,
          0.090, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0]
clf = linear_model.Ridge(fit_intercept=False)
clf2 = linear_model.Ridge(fit_intercept=False)
coefs = []
vif_list = [[] for x in range(X.shape[1])]
for a in alphas:
    clf.set_params(alpha=a)
    clf.fit(X, y)
    coefs.append(clf.coef_)

    for j, data in enumerate(X.columns):
        cols = [col for col in X.columns if col not in [data]]
        Z = X[cols]
        yy = X.iloc[:,j]
        clf2.set_params(alpha=a)
        clf2.fit(Z, yy)

        r_squared_j = clf2.score(Z, yy)
        vif = 1. / (1. - r_squared_j)
        print r_squared_j
        vif_list[j].append(vif)

pd.DataFrame(vif_list, columns = alphas).T
pd.DataFrame(coefs, index=alphas)

###############################################################################
# Display results

ax = pl.gca()
ax.set_color_cycle(['b', 'r', 'g', 'c', 'k', 'y', 'm'])

ax.plot(alphas, coefs)
pl.vlines(ridge_cv.alpha_, np.min(coefs), np.max(coefs), linestyle='dashdot')
pl.xlabel('alpha')
pl.ylabel('weights')
pl.title('Ridge coefficients as a function of the regularization')
pl.axis('tight')
pl.show()

- zerovector

1

我刚读了一篇1970年的文章，提到了Ridge的VIF。 Ridge回归的参数估计的协方差矩阵具有三明治形式，我认为您不能直接使用与OLS相同的模式。如果您没有得到更快的答案，我应该会在几天内提供针对statsmodels的代码。 - Josef

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Josef · Accepted Answer

岭回归的方差膨胀因子只需要三行代码。我在UCLA统计页面上的示例中进行了验证。

这个变体将会在下一个statsmodels发布版本中出现。这是我的当前函数:

def vif_ridge(corr_x, pen_factors, is_corr=True):
    """variance inflation factor for Ridge regression

    assumes penalization is on standardized variables
    data should not include a constant

    Parameters
    ----------
    corr_x : array_like
        correlation matrix if is_corr=True or original data if is_corr is False.
    pen_factors : iterable
        iterable of Ridge penalization factors
    is_corr : bool
        Boolean to indicate how corr_x is interpreted, see corr_x

    Returns
    -------
    vif : ndarray
        variance inflation factors for parameters in columns and ridge
        penalization factors in rows

    could be optimized for repeated calculations
    """
    corr_x = np.asarray(corr_x)
    if not is_corr:
        corr = np.corrcoef(corr_x, rowvar=0, bias=True)
    else:
        corr = corr_x

    eye = np.eye(corr.shape[1])
    res = []
    for k in pen_factors:
        minv = np.linalg.inv(corr + k * eye)
        vif = minv.dot(corr).dot(minv)
        res.append(np.diag(vif))
    return np.asarray(res)