9
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ this question is helpful, partially - stats.stackexchange.com/questions/23548/… $\endgroup$ – Ram Sharma Jun 26 '14 at 17:16
  • 4
    $\begingroup$ You might look at the R rms package ols, calibrate, and validate function with quadratic penalization (ridge regression). $\endgroup$ – Frank Harrell Jun 26 '14 at 17:41
  • $\begingroup$ @FrankHarrell I tried to extend your suggestion as answer for benefit of all. Please have a look ! $\endgroup$ – Ram Sharma Jul 1 '14 at 18:10
2
$\begingroup$

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)
# stratified 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)
$\endgroup$
4
+50
$\begingroup$

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

Your 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)

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.

$\endgroup$
  • $\begingroup$ Looks good. Also use pentrace function. $\endgroup$ – Frank Harrell Jul 1 '14 at 19:29
  • $\begingroup$ @FrankHarrell thanks for looking at. Please have a look at my current version, I hit few issues including error while executing penetrance function $\endgroup$ – Ram Sharma Jul 2 '14 at 15:10
  • $\begingroup$ 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$. $\endgroup$ – Frank Harrell Jul 2 '14 at 15:57
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ how can we do general linear regression in this package ? $\endgroup$ – rdorlearn Jun 29 '14 at 18:13
  • $\begingroup$ 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. $\endgroup$ – shadowtalker Jun 29 '14 at 20:58
  • $\begingroup$ glmnet seems interesting package, thanks for the information $\endgroup$ – rdorlearn Jun 29 '14 at 22:32
  • 1
    $\begingroup$ @rdorlearn Linear regression is just a GLM with an identity link function. $\endgroup$ – Joe Jul 1 '14 at 18:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.