基于Python(sklearn)计算PLS中的VIP值
sklearn中PLS回归模型并没有计算VIP值的方法,但VIP又是很重要的筛选变量方法。下附代码思路与完整代码,若有错误,万望指正。
1.首先亮明计算公式:
其中: VIPj:对应于第j个特征的VIP值; p:预测变量的总数; A:PLS成分的总数; R矩阵:A个PLS成分中,每个成分a都对应一套系数wa将X转换为成分得分,系数矩阵写作R,大小为p×A; T矩阵:n个样本中,每个样本会计算出A个成分得分,得分矩阵记做T,大小为n×A,ta代表n个样本的第a个成分的得分列表; 以上有T = XR qa: e是残差,X0和y0代表数据和标签。
2. 按照sklearn中对属性的解释: 该x_rotations_应对应于R矩阵,T矩阵对应于transform后生成的降维数据。 又依据: 有模型中e残差应为0,即y = yhat = X coef_
3. 组织计算步骤 1)计算qa2:qa = yhat.T ta -> Q2 = q2 = np.square(y.T, T) 2) p = X.shape[1] 3) = np.sum(Q2) 4) R = x_rotations_ 4) = R[j, a] / np.linalg.norm(R[:, a])
Ok,万事具备,开搞。 下付完整代码
def compute_VIP(X, y, R, T, A):
"""
计算模型中各预测变量的VIP值
:param X: 数据集X
:param y: 标签y
:param R: A个PLS成分中,每个成分a都对应一套系数wa将X转换为成分得分,系数矩阵写作R,大小为p×A
:param T: 得分矩阵记做T,大小为n×A,ta代表n个样本的第a个成分的得分列表
:param A: PLS成分的总数
:return: VIPs = np.zeros(p)
"""
p = X.shape[1]
Q2 = np.square(np.dot(y.T, T))
VIPs = np.zeros(p)
temp = np.zeros(A)
for j in range(p):
for a in range(A):
temp[a] = Q2[a] * pow(R[j, a] / np.linalg.norm(R[:, a]), 2)
VIPs[j] = np.sqrt(p * np.sum(temp) / np.sum(Q2))
return VIPs
下面为函数调用代码:
X, Y = np.zero()
x_train, x_test, y_train, y_test = model_selection.train_test_split(X, Y, test_size=0.3)
y_train_labels = pd.get_dummies(y_train)
n_component = 3
model = PLSRegression(n_components=n_component)
model.fit(x_train, y_train_labels)
x_test_trans = model.transform(x_test)
VIPs = compute_VIP(x_test, y_test, model.x_rotations_, x_test_trans, n_component)
plt.scatter(np.arange(0, X.shape[1]), VIPs)
plt.show()
|