Ridge-Regression in R mit p-Werten und Anpassungsgüte

Gratregression in RI haben entdeckt

• linearRidgeim ridgePaket - das zu einem Modell passt, Koeffizienten und p-Werte angibt, aber nichts, um die allgemeine Anpassungsgüte zu messen

• lm.ridgeim MASSPaket - das Koeffizienten und GCV, aber keine p-Werte für Parameter meldet

Wie kann ich all diese Dinge (Anpassungsgüte, Koeffizienten und p-Werte) aus derselben Gratregression erhalten? Ich bin neu in R und daher nicht mit Einrichtungen vertraut, die möglicherweise verfügbar sind, z. B. um aus den Daten und angepassten Koeffizienten zu berechnen.$r^2$

Obwohl ich lm.ridgeIhnen als Antwort auf eine frühere Frage empfohlen habe , könnten Sie das glmnetPaket als einen besseren Einstieg in die Ridge-Regression in R betrachten. Es hat den Vorteil, dass Sie dann den Beispielen in Kapitel 6 einer Einführung in folgen können Statistisches Lernen . Das Paket enthält auch Funktionen zur Kreuzvalidierung Ihres Modells. Wenn Sie diesen Prozess durchlaufen, können Sie die Probleme verstehen, die bei der Auswahl des Strafwerts für Ihre Gratregression auftreten. Es ist nicht so einfach wie eine automatische Antwort zu erhalten, aber am Ende haben Sie ein besseres Verständnis dafür, wie das endgültige Modell gebaut wurde, und können es anderen besser erklären (oder verteidigen). Und glmneterlaubt auch Kombinationen von Grat- und LASSO-Regression, die manche bevorzugen.

Wenn Sie daran interessiert sind, Ihr Modell für die Vorhersage zu verwenden, sind Messungen des Vorhersagefehlers durch Kreuzvalidierung oder Bootstrapping viel wichtiger als p- Werte für einzelne Koeffizienten. Für die Vorhersage möchten Sie alle Prädiktoren im Modell behalten (wenn nicht zu viele vorhanden sind), vorausgesetzt, Sie haben die Überanpassung mit einem Ansatz wie der Ridge-Regression minimiert, anstatt diejenigen mit nominell "unbedeutenden" p- Werten wegzulassen . Die im Paket angegebenen p- Werte basieren auf einem Algorithmus in einem Artikel des Autors des Pakets, der anscheinend nicht viel Beachtung gefunden hat, und das Paket ist derzeit verwaist, da die E-Mail des Autors nicht mehr an das R-Repository bei CRAN gemeldet wird funktioniert. Also ich 'linearRidgeridgep- Werte oder das nicht gepflegte Paket.

Wenn Sie aus irgendeinem Grund Koeffizienten- p- Werte benötigen , ist Bootstrapping ein besserer Ansatz. Sobald Sie einen Algorithmus zum Auswählen der Koeffizienten in einer Ridge-Regression haben, erstellen Sie aus Ihren Originaldaten mehrere Bootstrap-Samples (mit Ersetzung) derselben Größe. Wiederholen Sie dann den gesamten Vorgang , um die Regressionskoeffizienten für jedes Bootstrap-Sample mit dem bootPaket zu erhalten. Sie müssen eine Funktion schreiben, um die Regressionskoeffizienten zu melden, aber das ist ziemlich einfach. Die Verteilung dieser Regressionskoeffizienten unter den Bootstrap-Analysen liefert eine Schätzung ihrer Konfidenzgrenzen mit der boot.ciFunktion in derbootPaket. Auf diese Weise integrieren Sie die Variabilität, die sich aus allen Schritten des Modellierungsprozesses ergibt, in die Koeffizientenfehlerschätzungen.

Aber selbst diese Bootstrap- p- Werte könnten irreführend sein, da sie die Kompromisse zwischen den Koeffizienten kollinearer Prädiktoren ignorieren würden, denen Sie offensichtlich aufgrund Ihrer früheren Frage gegenüberstehen . Mit nur 2 kollinearen Prädiktoren variieren ihre individuellen Regressionskoeffizienten wahrscheinlich stark zwischen den Bootstrap-Stichproben, so dass ihre individuellen p- Werte möglicherweise unbedeutend erscheinen. Dennoch könnten ihre gemeinsamen Beiträge zur Regression viel stabiler und damit ihre Kombination in praktischer Hinsicht sehr bedeutsam sein. Überlegen Sie sich also genau, ob Sie wirklich an p- Werten für einzelne Koeffizienten interessiert sind .

Für das Ridge-Paket können Sie leicht entweder AIC oder BIC oder angepasstes R2 als Maß für die Anpassungsgüte berechnen, wenn Sie in diesen Formeln die richtigen effektiven Freiheitsgrade für die Ridge-Regression verwenden, die als Spur der Hutmatrix dienen.

Ridge-Regressionsmodelle sind in der Tat einfach als reguläre lineare Regression geeignet , wobei jedoch die kovariate Matrixzeile mit einer Matrix mitsqrt(lambda) [oder sqrt(lambdas)im Fall einer adaptiven Ridge-Regression] entlang der Diagonale (und pNullen, die der Ergebnisvariablen hinzugefügt werden y) ergänzt wird. Angesichts der Tatsache, dass die Ridge-Regression nur auf eine lineare Regression mit einer erweiterten kovariaten Matrix hinausläuft, können Sie viele der Funktionen regulärer linearer Modellanpassungen weiterhin verwenden. Das Originalpapier, auf dem das Ridge-Paket basiert, ist lesenswert: Signifikanztests bei der Ridge-Regression auf genetische Daten

Die Hauptfrage ist, wie Sie Ihren Lambda-Regularisierungsparameter einstellen. Im Folgenden habe ich das optimale Lambda für die Ridge- oder adaptive Ridge-Regression anhand derselben Daten abgestimmt, die für das endgültige Ridge- oder adaptive Ridge-Regressionsmodell verwendet wurden. In der Praxis ist es möglicherweise sicherer, Ihre Daten in einen Trainings- und Validierungssatz aufzuteilen und den Trainingssatz zum Einstellen von Lambda und den Validierungssatz zum Ableiten von Schlussfolgerungen zu verwenden. Oder machen Sie im Fall der adaptiven Gratregression 3 Teilungen und verwenden Sie 1, um Ihr anfängliches lineares Modell anzupassen, 1, um Lambda für die adaptive Gratregression abzustimmen (unter Verwendung der Koeffizienten der 1. Teilung, um Ihre adaptiven Gewichte zu definieren) und die dritte, um Inferenz zu machen Verwenden Sie das optimierte Lambda und die linearen Modellkoeffizienten, die aus den anderen Abschnitten Ihrer Daten abgeleitet wurden. Es gibt jedoch viele verschiedene Strategien, um Lambda auf Ridge und adaptive Ridge Regression einzustellenreden . Natürlich tendiert die Gratregression dazu, kollineare Variablen beizubehalten und zusammen auszuwählen, im Gegensatz zu z. B. LASSO oder nichtnegativen kleinsten Quadraten. Dies ist natürlich etwas zu beachten.

Die Koeffizienten der regulären Gratregression sind ebenfalls stark verzerrt, so dass dies natürlich auch die p-Werte stark beeinflusst. Dies ist beim adaptiven Grat weniger der Fall.

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  

Normalerweise werden bei der Ridge-Regression die effektiven Freiheitsgrade als Spur der Hutmatrix definiert. Sie müssten jedoch überprüfen, ob die wie oben berechneten BICs diese korrekten Freiheitsgrade verwenden, andernfalls müssten Sie sie selbst berechnen.