The model
We simulate from a linear model \[
Y=X\beta + \epsilon
\] where \(Y\in\mathbb R^n\),
\(\beta\in\mathbb R^p\), and \(\epsilon\sim N(0,\sigma^2I_n)\). We have
taken \(X\) to have iid \(N(0,1)\) entries and treat it as fixed in
this simulation. We define a Model object, which specifies
the parameters and, most importantly, describes how to simulate
data.
make_sparse_linear_model <- function(n, p, k, snr) {
x <- matrix(rnorm(n * p), n, p)
beta <- rep(c(1, 0), c(k, p - k))
mu <- as.numeric(x %*% beta)
sigma <- sqrt(sum(mu^2) / (n * snr)) # taking snr = ||mu||^2 / (n * sigma^2)
new_model(name = "slm", label = sprintf("n = %s, p = %s, k = %s", n, p, k),
params = list(x = x, beta = beta, mu = mu, sigma = sigma, n = n,
p = p, k = k),
simulate = function(mu, sigma, nsim) {
return(lapply(1:nsim, function(i) mu + sigma * rnorm(n)))
})
}
We will typically put the code above in a file named
model_functions.R.
The methods
We compare the lasso and ridge. Both of these methods depend on
tuning parameters, so we compute a sequence of solutions.
library(glmnet)
lasso <- new_method("lasso", "Lasso",
method = function(model, draw, lambda = NULL) {
if (is.null(lambda))
fit <- glmnet(x = model$x, y = draw, nlambda = 50,
intercept = FALSE)
else {
fit <- glmnet(x = model$x, y = draw, lambda = lambda,
intercept = FALSE)
}
list(beta = fit$beta, yhat = model$x %*% fit$beta,
lambda = fit$lambda, df = fit$df)
})
ridge <- new_method("ridge", "Ridge",
method = function(model, draw, lambda = NULL) {
sv <- svd(model$x)
df_fun <- function(lam) {
# degrees of freedom when tuning param is lam
sum(sv$d^2 / (sv$d^2 + lam))
}
if (is.null(lambda)) {
nlambda <- 50
get_lam <- function(target_df) {
f <- function(lam) df_fun(lam) - target_df
uniroot(f, c(0, 100 * max(sv$d^2)))$root
}
lambda <- sapply(seq(1, nrow(model$x),
length = nlambda), get_lam)
}
df <- sapply(lambda, df_fun)
beta <- sapply(lambda, function(r) {
d <- sv$d / (sv$d^2 + r)
return(sv$v %*% (d * crossprod(sv$u, draw)))
})
list(beta = beta, yhat = model$x %*% beta,
lambda = lambda, df = df)
})
Methods can return different items. However, aspects of the method
that will be used downstream in the simulation and compared across
methods should be in a common format. Thus beta,
yhat, and df in each case are in the same
format. These will be the items used when evaluating the methods’
performances.
We will typically put the code above in a file named
method_functions.R.
The metrics
When we compare methods in plots and tables, there are usually a
number of “metrics” we use. An object of class Metric
specifies how to go from a model’s parameters and the output of a method
and return some quantity of interest.
sqrerr <- new_metric("sqrerr", "squared error",
metric = function(model, out) {
colMeans(as.matrix(out$beta - model$beta)^2)
})
best_sqrerr <- new_metric("best_sqrerr", "best squared error",
metric = function(model, out) {
min(colMeans(as.matrix(out$beta - model$beta)^2))
})
nnz <- new_metric("nnz", "number of nonzeros",
metric = function(model, out) {
colSums(as.matrix(out$beta) != 0)
})
df <- new_metric("df", "degrees of freedom",
metric = function(model, out) out$df)
Observe that out refers to the list returned by our
methods and model refers to the Model object
that is generated by make_sparse_linear_model. The
$ operator can be used to get parameters that are stored in
the Model object.
We will typically put the code above in a file named
eval_functions.R.
Extensions of the methods
In the simulations above, we saw two examples in which we studied
methods that were extended versions of other methods. While the relaxed
lasso could have been treated as a completely separate method from the
lasso, this would have been wasteful in a comparison to the lasso since
it would have required us to fit the lasso twice (once when calling the
lasso and once within the relaxed lasso). Instead, we created what is
called an ExtendedMethod object, which behaves like a
Method object, except it allows one to use the output of
another method.
If we only care about relaxing the lasso, we might directly create an
ExtendedMethod. However, in some situations, we may wish to
extend multiple methods in the same fashion. For example, perhaps we may
later wish to “relax” both the lasso and another variable selection
procedure, such as SCAD or MCP. In the spirit of code reusability, we
therefore created a MethodExtension object called
refit. A MethodExtension object when added to
a Method object generates an ExtendedMethod.
In our example of the relaxed lasso, we created a
MethodExtension called refit that takes any
method’s output and performs least squares on the nonzeros.
refit <- new_method_extension(name = "refit", label = "refitted",
method_extension = function(model, draw, out,
base_method) {
beta <- apply(out$beta, 2, function(b) {
ii <- b != 0
if (sum(ii) == 0)
return(b)
xtx <- crossprod(model$x[, ii])
# add small ridge in case solution has more
# than n nonzeros:
diag(xtx) <- diag(xtx) + 1e-4
bb <- solve(xtx,
crossprod(model$x[, ii], draw))
b[ii] <- bb
return(b)
})
return(list(beta = beta,
yhat = model$x %*% beta,
lambda = out$lambda,
df = rep(NA, ncol(beta))))
})
Having created the MethodExtension, we can get the
relaxed lasso ExtendedMethod through the simple addition of
objects:
## ExtendedMethod Component
## name: lasso_refit
## label: Lasso refitted
If we had defined a different variable selection procedure, say
mcp, we could just as easily generated its relaxed version,
writing no new code:
relaxed_mcp <- mcp + refit
The other example of a method extension in the simulations above was
cross-validation. We wrote this method extension so that it could be
used to create extended versions of multiple methods.
#' Make folds for cross validation
#'
#' Divides the indices \code{1:n} into \code{nfolds} random folds of about the same size.
#'
#' @param n sample size
#' @param nfolds number of folds
make_folds <- function(n, nfolds) {
nn <- round(n / nfolds)
sizes <- rep(nn, nfolds)
sizes[nfolds] <- sizes[nfolds] + n - nn * nfolds
b <- c(0, cumsum(sizes))
ii <- sample(n)
folds <- list()
for (i in seq(nfolds))
folds[[i]] <- ii[seq(b[i] + 1, b[i + 1])]
folds
}
cv <- new_method_extension("cv", "cross validated",
method_extension = function(model, draw, out,
base_method) {
nfolds <- 5
err <- matrix(NA, ncol(out$beta), nfolds)
ii <- make_folds(model$n, nfolds)
for (i in seq_along(ii)) {
train <- model
train@params$x <- model@params$x[-ii[[i]], ]
train@params$n <- nrow(train@params$x)
train_draw <- draw[-ii[[i]]]
test <- model
test@params$x <- model@params$x[ii[[i]], ]
test@params$n <- nrow(test@params$x)
test_draw <- draw[ii[[i]]]
fit <- base_method@method(model = train,
draw = train_draw,
lambda = out$lambda)
yhat <- test$x %*% fit$beta
ll <- seq(ncol(yhat))
err[ll, i] <- colMeans((yhat - test_draw)^2)
}
m <- rowMeans(err)
se <- apply(err, 1, sd) / sqrt(nfolds)
imin <- which.min(m)
ioneserule <- max(which(m <= m[imin] + se[imin]))
list(err = err, m = m, se = se, imin = imin,
ioneserule = ioneserule,
beta = out$beta[, imin],
yhat = model$x %*% out$beta[, imin])
})
We will typically put the code above in the file named
method_functions.R, though if that file is large, we may
prefer to have a separate file named, perhaps,
method_extensions_functions.R.
