尽管我lm.ridge在回答之前的问题时向您推荐过，但您可能会认为该glmnet软件包是在 R 中开始使用岭回归的更好方法。它的优点是您可以遵循An Introduction to第 6 章中的示例统计学习。该包还具有用于交叉验证模型的功能。完成此过程将帮助您了解选择岭回归的惩罚值所涉及的问题。这并不像获得自动答案那么容易，但最终您将更好地了解最终模型是如何构建的，并且能够更好地向他人解释（或捍卫）它。并且glmnet还允许一些人喜欢的岭回归和 LASSO 回归的组合。

如果您有兴趣使用您的模型进行预测，那么来自交叉验证或自举的预测误差的度量要比单个系数的p值重要得多。对于预测，您希望将所有预测变量保留在模型中（如果没有太多），前提是您使用岭回归之类的方法最大限度地减少了过度拟合，而不是省略那些名义上“无关紧要”的p值。包中报告的p值基于包作者的论文中的算法，该算法似乎没有受到太多关注，并且该包目前是孤立的，因为作者的电子邮件报告给 CRAN 的 R 存储库不再作品。所以我'linearRidgeridgep值或未维护的包。

如果出于某种原因确实需要系数p值，则自举将是一种更好的方法。一旦有了用于在岭回归中选择系数的算法，您就可以从原始数据中制作多个相同大小的引导样本（带替换），然后重复整个过程以获得每个带有boot包的引导样本的回归系数。您将不得不编写一个函数来报告回归系数，但这相当简单。这些回归系数在 bootstrap 分析中的分布提供了对其置信限的估计，boot.ci函数在boot包裹。这样，您就可以将建模过程的所有步骤产生的可变性纳入系数误差估计中。

但即使是那些自举的p值也可能会产生误导，因为它们会忽略共线预测变量系数之间的权衡，根据您之前的问题，您显然会面临这种权衡。即使只有 2 个共线预测变量，它们各自的回归系数也可能在 bootstrap 样本之间变化很大，因此它们各自的p值可能显得微不足道。然而，它们对回归的共同贡献可能更加稳定，因此它们的组合在实际中非常重要。因此，请认真考虑您是否真的对单个系数的p值感兴趣。

对于 ridge 包，您可以轻松地计算 AIC 或 BIC 或调整后的 R2 作为拟合优度的度量，如果在这些公式中使用脊回归的正确有效自由度，则作为帽子矩阵的轨迹。

岭回归模型实际上只是作为常规线性回归拟合，但协变量矩阵行增加了一个矩阵，sqrt(lambda) [或sqrt(lambdas)在自适应岭回归的情况下]沿对角线（并且p零添加到结果变量y）。因此，鉴于岭回归只是归结为使用增强协变量矩阵进行线性回归，您可以继续使用常规线性模型拟合的许多特征。ridge 包所基于的原论文值得一读： Significance testing in ridge regression forgenetic data

主要问题是如何调整您的 lambda 正则化参数。下面我在用于拟合最终脊或自适应脊回归模型的相同数据上调整了脊或自适应脊回归的最佳 lambda。在实践中，将数据拆分到训练和验证集中并使用训练集调整 lambda 并使用验证集进行推理可能更安全。或者在自适应岭回归的情况下，进行 3 次分割并使用 1 来拟合您的初始线性模型，使用 1 来调整自适应岭回归的 lambda（使用第一次分割的系数来定义自适应权重），第三次进行推理使用优化的 lambda 和从数据的其他部分得出的线性模型系数。不过，有许多不同的策略可以为岭和自适应岭回归调整 lambda，请参阅此说话。当然，岭回归将倾向于保留共线变量并将它们一起选择，这与 LASSO 或非负最小二乘法不同。当然，这是要牢记的。

常规岭回归的系数也有很大的偏差，所以这当然也会严重影响 p 值。自适应脊的情况较少。

library(MASS)
data=longley
data=data.frame(apply(data,2,function(col) scale(col))) # we standardize all columns
# UNPENALIZED REGRESSION MODEL
lmfit = lm(y~.,data)
summary(lmfit)
# Coefficients:
#                  Estimate Std. Error t value Pr(>|t|)  
# (Intercept)     9.769e-16  2.767e-02   0.000   1.0000  
# GNP             2.427e+00  9.961e-01   2.437   0.0376 *
# Unemployed      3.159e-01  2.619e-01   1.206   0.2585  
# Armed.Forces    7.197e-02  9.966e-02   0.722   0.4885  
# Population     -1.120e+00  4.343e-01  -2.578   0.0298 *
# Year           -6.259e-01  1.299e+00  -0.482   0.6414  
# Employed        7.527e-02  4.244e-01   0.177   0.8631  
lmcoefs = coef(lmfit)

# RIDGE REGRESSION MODEL
# function to augment covariate matrix with matrix with sqrt(lambda) along diagonal to fit ridge penalized regression
dataaug=function (lambda, data) { p=ncol(data)-1 # nr of covariates; data contains y in first column
                                  data.frame(rbind(as.matrix(data[,-1]),diag(sqrt(lambda),p)),yaugm=c(data$y,rep(0,p))) } 
# function to calculate optimal penalization factor lambda of ridge regression on basis of BIC value of regression model
BICval_ridge = function (lambda, data) BIC(lm(yaugm~.,data=dataaug(lambda, data)))
lambda_ridge = optimize(BICval_ridge, interval=c(0,10), data=data)$minimum
lambda_ridge # ridge lambda optimized based on BIC value = 5.575865e-05
ridgefit = lm(yaugm~.,data=dataaug(lambda_ridge, data))
summary(ridgefit)
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)   
# (Intercept)  -0.000388   0.018316  -0.021  0.98338   
# GNP           2.411782   0.769161   3.136  0.00680 **
# Unemployed    0.312543   0.202167   1.546  0.14295   
# Armed.Forces  0.071125   0.077050   0.923  0.37057   
# Population   -1.115754   0.336651  -3.314  0.00472 **
# Year         -0.608894   1.002436  -0.607  0.55266   
# Employed      0.072204   0.328229   0.220  0.82885   

# ADAPTIVE RIDGE REGRESSION MODEL
# function to calculate optimal penalization factor lambda of adaptive ridge regression 
# (with gamma=2) on basis of BIC value of regression model
BICval_adridge = function (lambda, data, init.coefs) { dat = dataaug(lambda*(1/(abs(init.coefs)+1E-5)^2), data)
                                                       BIC(lm(yaugm~.,data=dat)) }
lamvals = 10^seq(-12,-1,length.out=100)
BICvals = sapply(lamvals,function (lam) BICval_adridge (lam, data, lmcoefs))
firstderivBICvals = function (lambda) splinefun(x=lamvals, y=BICvals)(lambda, deriv=1)
plot(lamvals, BICvals, type="l", ylab="BIC", xlab="Adaptive ridge lambda", log="x")
lambda_adridge = lamvals[which.min(lamvals*firstderivBICvals(lamvals))] # we place optimal lambda at middle flat part of BIC as coefficients are most stable there
lambda_adridge # 4.641589e-08
abline(v=lambda_adridge, col="red")

adridgefit = lm(yaugm~.,data=dataaug(lambda_adridge, data))
summary(adridgefit)
# Coefficients:
#                 Estimate Std. Error t value Pr(>|t|)   
# (Intercept)   -1.121e-05  1.828e-02  -0.001  0.99952   
# GNP            2.427e+00  7.716e-01   3.146  0.00666 **
# Unemployed     3.159e-01  2.029e-01   1.557  0.14026   
# Armed.Forces   7.197e-02  7.719e-02   0.932  0.36590   
# Population    -1.120e+00  3.364e-01  -3.328  0.00459 **
# Year          -6.259e-01  1.006e+00  -0.622  0.54326   
# Employed       7.527e-02  3.287e-01   0.229  0.82197  

通常在岭回归中，有效自由度被定义为帽子矩阵的轨迹 - 你必须检查如上计算的 BIC 是否使用这些正确的自由度，否则你必须自己计算它们。