# K-fold or hold-out cross validation for ridge regression using R

I am working on cross-validation of prediction of my data with 200 subjects and 1000 variables. I am interested ridge regression as number of variables (I want to use) is greater than number of sample. So I want to use shrinkage estimators. The following is made up example data:

 #random population of 200 subjects with 1000 variables
M <- matrix(rep(0,200*100),200,1000)
for (i in 1:200) {
set.seed(i)
M[i,] <- ifelse(runif(1000)<0.5,-1,1)
}
rownames(M) <- 1:200

#random yvars
set.seed(1234)
u <- rnorm(1000)
g <- as.vector(crossprod(t(M),u))
h2 <- 0.5
set.seed(234)
y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))

myd <- data.frame(y=y, M)
myd[1:10,1:10]

y X1 X2 X3 X4 X5 X6 X7 X8 X9
1   -7.443403 -1 -1  1  1 -1  1  1  1  1
2  -63.731438 -1  1  1 -1  1  1 -1  1 -1
3  -48.705165 -1  1 -1 -1  1  1 -1 -1  1
4   15.883502  1 -1 -1 -1  1 -1  1  1  1
5   19.087484 -1  1  1 -1 -1  1  1  1  1
6   44.066119  1  1 -1 -1  1  1  1  1  1
7  -26.871182  1 -1 -1 -1 -1  1 -1  1 -1
8  -63.120595 -1 -1  1  1 -1  1 -1  1  1
9   48.330940 -1 -1 -1 -1 -1 -1 -1 -1  1
10 -18.433047  1 -1 -1  1 -1 -1 -1 -1  1


I would like to do following for cross validation -

(1) split data into two halts - use first half as training and second half as test

(2) K-fold cross validation (say 10 fold or suggestion on any other appropriate fold for my case are welcome)

I can simply sample the data into two (gaining and test) and use them:

# using holdout (50% of the data) cross validation
training.id <- sample(1:nrow(myd), round(nrow(myd)/2,0), replace = FALSE)
test.id <- setdiff(1:nrow(myd), training.id)

myd_train <- myd[training.id,]
myd_test  <- myd[test.id,]


I am using lm.ridge from MASS R package.

library(MASS)
out.ridge=lm.ridge(y~., data=myd_train, lambda=seq(0, 100,0.001))
plot(out.ridge)
select(out.ridge)

lam=0.001
abline(v=lam)

out.ridge1 =lm.ridge(y~., data=myd_train, lambda=lam)
hist(out.ridge1$coef) out.ridge1$ym
hist(out.ridge1$xm)  I have two questions - (1) How can I predict the test set and calculate accuracy (as correlation of predicted vs actual)? (2) How can I perform K-fold validation? say 10-fold? • this question is helpful, partially - stats.stackexchange.com/questions/23548/… Jun 26, 2014 at 17:16 • You might look at the R rms package ols, calibrate, and validate function with quadratic penalization (ridge regression). Jun 26, 2014 at 17:41 • @FrankHarrell I tried to extend your suggestion as answer for benefit of all. Please have a look ! Jul 1, 2014 at 18:10 ## 3 Answers You can use caret package (vignnettes, paper ) for this type of things, which can wrap a number of machine learning models or you can use your own customized models . As you are interested in ridge regression here is just custom codes for ridge regression, you might want to adopt to your situation more precisely. For simple split in data: set.seed(107) # stratiﬁed random split of the data inTrain <- createDataPartition(y = myd$y, p = .5,list = FALSE)
training <- myd[ inTrain,]
testing <- myd[-inTrain,]


For K-fold validation and other type of CV including default boot

ridgeFit1 <- train(y ~ ., data = training,method = 'ridge',
preProc = c("center", "scale"), metric = "ROC")
plot(ridgeFit1)


Here is discussion on how to use train function. Note the ridge method depends upon the package elasticnet functions ( and its dependency on lars, should or need to be installed). If not installed in the system will ask if you want to do so.

the type of resampling used, The simple bootstrap is used by default.To modify the resampling method, a trainControl function is used

The option method controls the type of resampling and defaults to "boot". Another method, "repeatedcv", is used to specify repeated K–fold cross–validation (and the argument repeats controls the number of repetitions). K is controlled by the number argument and defaults to 10.

 ctrl <- trainControl(method = "repeatedcv", repeats = 5)

ridgeFit <- train(y ~ ., data = training,method = 'ridge',
preProc = c("center", "scale"),trControl = ctrl, metric = "ROC")

plot(ridgefit)


For predictions:

plsClasses <- predict(ridgeFit, newdata = testing)


This is extension of the suggestion by Frank in the comments. Dr. Harrel please correct if I am wrong (appreciate corrections).

#random population of 200 subjects with 1000 variables
M <- matrix(rep(0,200*100),200,1000)
for (i in 1:200) {
set.seed(i)
M[i,] <- ifelse(runif(1000)<0.5,-1,1)
}
rownames(M) <- 1:200

#random yvars
set.seed(1234)
u <- rnorm(1000)
g <- as.vector(crossprod(t(M),u))
h2 <- 0.5
set.seed(234)
y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))

myd <- data.frame(y=y, M)


Install rms package and load it.

require(rms)


ols function is used for Linear Model Estimation Using Ordinary Least Squares where can specify penalty term.

As suggested below in comments I added petrace function. This function trace AIC and BIC vs Penalty.

# using holdout (50% of the data) cross validation
training.id <- sample(1:nrow(myd), round(nrow(myd)/2,0), replace = FALSE)
test.id <- setdiff(1:nrow(myd), training.id)

myd_train <- myd[training.id,]
myd_test  <- myd[test.id,]

frm <- as.formula(paste("y~",paste(names(myd_train)[2:100],collapse="+")))


Important note I could not use all 1000 of the variables as the program complains if number of variable exceeds 100. Also y~. type formula designation did not work. So see above way of doing same creating formula object frm

f <- ols(frm, data = myd_train, method="qr", x=TRUE, y=TRUE)

p <- pentrace(f, seq(.2,1,by=.05))

Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x),  :
'data' must be of a vector type, was 'NULL'

plot(p)


"For an ordinary unpenalized fit from lrm or ols and for a vector or list of penalties, fits a series of logistic or linear models using penalized maximum likelihood estimation, and saves the effective degrees of freedom, Akaike Information Criterion (AIC), Schwarz Bayesian Information Criterion (BIC), and Hurvich and Tsai's corrected AIC (AIC_c). Optionally pentrace can use the nlminb function to solve for the optimum penalty factor or combination of factors penalizing different kinds of terms in the model." from rms package manual.

calibrate function is for Resampling Model Calibration and Uses bootstrapping or cross-validation to get bias-corrected (overfitting- corrected) estimates of predicted vs. observed values based on subsetting predictions into intervals. The validate function does resampling validation of a regression model, with or without backward step-down variable deletion. B = number of repetitions. For method="crossvalidation", is the number of groups of omitted observations

cal <- calibrate(f, method = "cross validation", B=20)
plot(cal)


You can use Predict function to compute predicted values and confidence limits. I am not sure I this works in test situation.

• Looks good. Also use pentrace function. Jul 1, 2014 at 19:29
• @FrankHarrell thanks for looking at. Please have a look at my current version, I hit few issues including error while executing penetrance function Jul 2, 2014 at 15:10
• You didn't specify x=TRUE, y=TRUE to ols. But there is a problem with pentrace when the model is completely overfit (error d.f. of zero) in that pentrace tries to examine an unpenalized model, which has $R^{2} = 1.0$. For the next release of rms I've added a new argument to pentrace: noaddzero=TRUE to not add zero to the list of penalties to try. Note that your example isn't the best one, as the optimum penalty is $\infty$. Jul 2, 2014 at 15:57

R package glmnet (vignette) has a wrapper function that does exactly what you want, called cv.glmnet (doc). I just used it yesterday, it works like a dream.

• how can we do general linear regression in this package ? Jun 29, 2014 at 18:13
• For linear regression, there's cv.lm in package:DAAG, and for a GLM there's cv.glm in package:boot. But, I just realized Frank Harrell suggested rms. Basically you should do whatever he tells you. It also seems like it's a more general framework than the piecemeal one I'm suggesting anyway. Jun 29, 2014 at 20:58
• glmnet seems interesting package, thanks for the information Jun 29, 2014 at 22:32
• @rdorlearn Linear regression is just a GLM with an identity link function.
– Joe
Jul 1, 2014 at 18:56