version 0.4-2

DESCRIPTION

@@ -1,13 +1,13 @@
 Package: grplasso
 Type: Package
 Title: Fitting user specified models with Group Lasso penalty
-Version: 0.4-1
-Date: 2009-03-30
+Version: 0.4-2
+Date: 2009-04-01
 Author: Lukas Meier
 Maintainer: Lukas Meier <lukas.meier@gmx.net>
 Description: Fits user specified (GLM-) models with Group Lasso penalty
 Depends: methods
 License: GPL
-Packaged: Mon Mar 30 15:52:36 2009; meier
+Packaged: Wed Apr 1 11:36:10 2009; meier
 Repository: CRAN
-Date/Publication: 2009-03-30 18:25:43
+Date/Publication: 2009-04-01 11:43:14

R/grplasso.R

@@ -114,6 +114,9 @@ grplasso.default <- function(x, y, index, weights = rep(1, length(y)),
   if(!is.numeric(index))
     stop("index has to be of type 'numeric'!")
 
+  if(length(index) != ncol(x))
+    stop("length(index) not equal ncol(x)!")
+
   if(is.unsorted(rev(lambda)))
     warning("lambda values should be sorted in decreasing order")
 
@@ -169,18 +172,27 @@ grplasso.default <- function(x, y, index, weights = rep(1, length(y)),
   intercept.which <- which(apply(x == 1, 2, all))
   has.intercept   <- length(intercept.which)
 
+  if(!has.intercept & center){
+    message("Couldn't find intercept. Setting center = FALSE.")
+    center <- FALSE
+  }
+
   if(length(intercept.which) > 1)
     stop("Multiple intercepts!")
 
-  has.intercept.notpen <- is.na(index[intercept.which])
-
+  if(has.intercept)
+    has.intercept.notpen <- is.na(index[intercept.which])
+  else
+    has.intercept.notpen <- FALSE
+
   others.notpen   <- nrnotpen - has.intercept.notpen
   notpen.int.only <- has.intercept.notpen & !others.notpen
 
-  if(notpen.int.only & !center)
-    warning("Are you sure you want uncentered predictors in a model with intercept?")
+  if(has.intercept & !center & standardize)
+    warning("Are you sure that you don't want to perform centering in a model with intercept and standardized predictors?")
 
-  if(center & others.notpen)
+  ##if(center & others.notpen)
+  if(others.notpen)
     warning("Penalization not adjusted to non-penalized predictors.")
 
   ## Index vector of the penalized parameter groups
@@ -194,22 +206,19 @@ grplasso.default <- function(x, y, index, weights = rep(1, length(y)),
     ipen.which <- split((1:ncolx), ipen)
   }
 
-  nrpen      <- length(ipen.which)
-  dict.pen   <- sort(unique(ipen))
+  nrpen    <- length(ipen.which)
+  dict.pen <- sort(unique(ipen))
 
   ## Table of degrees of freedom
   ipen.tab   <- table(ipen)[as.character(dict.pen)]
 
   x.old <- x
 
   if(center){
-    ##if(length(intercept.which) == 0)
-    ##  stop("Need intercept term")
-
-    ##if(length(intercept.which) == 1 & !is.na(index[intercept.which]))
-    ##  stop("Need *unpenalized* intercept")
+    if(!has.intercept) ## could be removed; already handled above
+      stop("Need intercept term when using center = TRUE")
 
-    mu.x <- apply(x[,-intercept.which], 2, mean)
+    mu.x                 <- apply(x[,-intercept.which], 2, mean)
     x[,-intercept.which] <- sweep(x[,-intercept.which], 2, mu.x)
   }
 

R/helpers.R

@@ -215,21 +215,28 @@ lambdamax.default <- function(x, y, index, weights = rep(1, length(y)),
   ## Indices of parameter groups
   ipen.which <- split((1:ncol(x))[!is.na(index)], ipen)
 
-  intercept.which  <- which(apply(x == 1, 2, all))
+  intercept.which <- which(apply(x == 1, 2, all))
+  has.intercept   <- length(intercept.which)
 
-  if(center){
-    if(length(intercept.which) == 0)
-      stop("Need intercept term")
+  if(!has.intercept & center){
+    message("Couldn't find intercept. Setting center = FALSE.")
+    center <- FALSE
+  }
 
-    if(length(intercept.which) > 1)
-      stop("Multiple intercepts!")
+  if(length(intercept.which) > 1)
+    stop("Multiple intercepts!")
 
-    if(length(intercept.which) == 1 & !is.na(index[intercept.which]))
-      stop("Need unpenalized intercept")
+  if(center){
+    if(!has.intercept) ## Could be removed; already handled above
+      stop("Need intercept term when using center = TRUE")
+
+    ##if(length(intercept.which) == 1 & !is.na(index[intercept.which]))
+    ##  stop("Need unpenalized intercept")
 
     mu.x <- apply(x[,-intercept.which], 2, mean)
     x[,-intercept.which] <- sweep(x[,-intercept.which], 2, mu.x)
   }
+
   if(standardize){
     stand        <- blockstand(x, ipen.which, inotpen.which)
     x            <- stand$x

man/grplasso.Rd

@@ -16,8 +16,8 @@ grplasso(x, ...)
 
 \method{grplasso}{default}(x, y, index, weights = rep(1, length(y)), offset = rep(0,
          length(y)), lambda, coef.init = rep(0, ncol(x)),
-         penscale = sqrt, model = LogReg(), center = TRUE, standardize = TRUE,
-         control = grpl.control(), ...)
+         penscale = sqrt, model = LogReg(), center = TRUE,
+         standardize = TRUE, control = grpl.control(), ...)
 }
 \arguments{
   \item{x}{design matrix (including intercept)}
@@ -52,7 +52,8 @@ grplasso(x, ...)
   \item{center}{logical. If true, the columns of the design matrix will be
     centered (except a possible intercept column).}
   \item{standardize}{logical. If true, the design matrix will be
-    blockwise orthonormalized such that for each block \eqn{X^TX = n 1}.}
+    blockwise orthonormalized such that for each block \eqn{X^TX = n 1}
+         (*after* possible centering).}
   \item{control}{options for the fitting algorithm, see
     \code{\link{grpl.control}}.}
   \item{contrasts}{an optional list. See the 'contrasts.arg' of

man/lambdamax.Rd

@@ -17,8 +17,8 @@ lambdamax(x, ...)
 
 \method{lambdamax}{default}(x, y, index, weights = rep(1, length(y)),
           offset = rep(0, length(y)), coef.init = rep(0, ncol(x)),
-          penscale = sqrt, model = LogReg(), center = TRUE, standardize = TRUE,
-          nlminb.opt = list(), ...)
+          penscale = sqrt, model = LogReg(), center = TRUE,
+          standardize = TRUE, nlminb.opt = list(), ...)
 }
 \arguments{
   \item{x}{design matrix (including intercept)}
@@ -49,7 +49,8 @@ lambdamax(x, ...)
   \item{center}{logical. If true, the columns of the design matrix will be
     centered (except a possible intercept column).}
   \item{standardize}{logical. If true, the design matrix will be blockwise
-          orthonormalized, such that for each block \eqn{X^TX = n 1}}
+          orthonormalized, such that for each block \eqn{X^TX = n 1}
+          (*after* possible centering).}
   \item{contrasts}{an (optional) list with the contrasts for the factors
     in the model.}
   \item{nlminb.opt}{arguments to be supplied to \code{\link{nlminb}}.}

tests/test_centerscaling.R

@@ -10,9 +10,6 @@ y <- 4 * sin(x[,1]) + x[,2]^2 + rnorm(n)
 
 x.new <- matrix(runif(n * p, min = -2.5, max = 2.5), nrow = n, ncol = p)
 
-center  <- TRUE
-scaling <- TRUE
-
 shift <- 10
 scale <- 10
 
@@ -44,3 +41,64 @@ stopifnot(all.equal(int.resc, coef(fit.resc)[1,], tol = 10^-7))
 stopifnot(all.equal(predict(fit, newdata = cbind(1,x.new)),
                     predict(fit.resc,
                             newdata = cbind(1, shift + scale * x.new))))
+
+## Check whether every case is running, including function lambda.max
+
+## center = TRUE & unpenalized intercept
+
+x.use <- cbind(1, x)
+index <- c(NA, 1:10)
+
+lambda.max <- lambdamax(x.use, y, index, model = LinReg())
+lambda     <- lambda.max * c(1, 0.1)
+
+fit1 <- grplasso(x.use, y, index, model = LinReg(), lambda = lambda,
+                 center = TRUE)
+
+## center = TRUE & penalized intercept
+
+x.use <- cbind(1, x)
+index <- c(99, 1:10)
+
+lambda.max <- lambdamax(x.use, y, index, model = LinReg())
+lambda     <- lambda.max * c(1, 0.1)
+
+fit2 <- grplasso(x.use, y, index, model = LinReg(), lambda = lambda,
+                 center = TRUE)
+
+## center = TRUE & *no* intercept
+x.use <- cbind(x)
+index <- c(1:10)
+
+lambda.max <- lambdamax(x.use, y, index, model = LinReg())
+lambda     <- lambda.max * c(1, 0.1)
+
+fit3 <- grplasso(x.use, y, index, model = LinReg(), lambda = lambda,
+                 center = TRUE)
+
+## center = FALSE & *no* intercept
+x.use <- cbind(x)
+index <- c(1:10)
+
+lambda.max <- lambdamax(x.use, y, index, model = LinReg(),
+                        center = FALSE)
+lambda     <- lambda.max * c(1, 0.1)
+
+fit4 <- grplasso(x.use, y, index, model = LinReg(), lambda = lambda,
+                 center = FALSE)
+
+## Check whether fit3 and fit4 are the same
+stopifnot(all.equal(coef(fit3), coef(fit4)))
+
+
+## center = FALSE & standardize = TRUE
+x.use <- cbind(1,x)
+index <- c(NA, 1:10)
+
+lambda.max <- lambdamax(x.use, y, index, model = LinReg(),
+                        center = FALSE)
+lambda     <- lambda.max * c(1, 0.1)
+
+fit5 <- grplasso(x.use, y, index, model = LinReg(), lambda = lambda,
+                 center = FALSE)
+