## PALM_STLS Class' methods
PALM_STLS$methods(
initialize = function(X, y, ...) {
if(nrow(X) != length(y)) {
stop("nrow(X) should be equal to size of y.")
}
.self$X = as.matrix(X)
.self$y = as.numeric(y)
.self$lambda = numeric(0)
.self$nlambda = 100L
.self$lmin_ratio = 1e-5
.self$alpha = 1.0
.self$maxit = 10000L
.self$eps_abs = 1e-6
.self$eps_rel = 1e-6
.self$center = TRUE
.self$warm_start = FALSE
}
)
## print PALM_STLS object
PALM_STLS$methods(
show_common = function() {
cat(sprintf("X := %d x %d Matrix", nrow(.self$X), ncol(.self$X)))
cat(sprintf("y := %d Vector"), length(.self$y))
fields = setdiff(names(.refClassDef@fieldClasses), c("X", "y"))
for(field in fields) {
cat(field, " : ", paste(.self$field(field), collapse = " "), "\n", sep = "")
}
},
show = function()
{
cat("Sparse total least square model in PALM \n \n")
show_common()
}
)
## penalty parameter set-up
PALM_STLS$methods(
penalty = function(lambda = NULL, nlambda = 100, lmin_ratio = 1e-5, alpha = 1.0, ...) {
lambda_val = sort(as.numeric(lambda), decreasing = TRUE)
if(length(lambda_val) > 0 & any(lambda_val < 0)) {
stop("lambda must be non-negative value.")
}
if(nlambda <= 0) {
stop("nlambda must be a positive integer.")
}
if(alpha < 0 | alpha > 1) {
stop("alpha must be located in (0, 1].")
}
.self$lambda = lambda_val
.self$nlambda = nlambda
.self$lmin_ratio = lmin_ratio
.self$alpha = alpha
invisible(.self)
}
)
## options configuration
PALM_STLS$methods(
options = function(maxit = 1000L, eps_abs = 1e-6, eps_rel = 1e-6, center = TRUE, warm_start = FALSE, ...) {
if(maxit <= 0) {
stop("maxit should be positive integer.")
}
if(eps_abs < 0 | eps_rel < 0) {
stop("eps_abs/rel should be positive.")
}
.self$maxit = maxit
.self$eps_abs = eps_abs
.self$eps_rel = eps_rel
.self$center = center
.self$warm_start = warm_start
invisible(.self)
}
)
## fit
PALM_STLS$methods(
fit = function(...) {
res = .Call("PALM_STLS", .self$X, .self$y, .self$lambda, .self$nlambda, .self$lmin_ratio, .self$alpha,
list(maxit = .self$maxit,
eps_abs = .self$eps_abs,
eps_rel = .self$eps_rel,
center = .self$center,
warm_start = .self$warm_start),
PACKAGE = "PALM")
do.call(STLS_fit, res)
}
)
### STLS_fit object
STLS_fit$methods(
show_common = function() {
cat("Lambda \n")
print(.self$lambda)
cat("\n")
cat("Intercept \n")
cat(.self$intercept)
cat("\n")
cat("Coef \n")
cat(sprintf("%d x %d Sparse Matrix: \n", nrow(.self$beta), ncol(.self$beta)))
cat("\n")
cat("Niter \n")
print(.self$niter)
},
show = function() {
cat("PALM based Total Least Square fitting results \n")
show_common()
}
)
#' Fitting A sparse Total Least Square model Using PALM algorithm
#'
#' @title fitSTLS
#' @param X The observation data matrix
#' @param y The observation response vector
#' @param center Y = Y - mean(Y) and X = X - mean(X)
#' @param lambda The regularized lambda vector
#' @param nlambda The Number of lambda
#' @param lmin_ratio lambda_max / lambda_min
#' @param alpha alpha parameter for elastic net penalty (0,1](Ridge Regression --> Lasso)
#' @param maxit The maximal iteration of PALM algorithm
#' @param eps_abs Absolutely epsilon for generated varaince vector
#' @param eps_rel Relative epsilon for objective function value of generated vector
#' @param warm_start Warm start for regularizer parameter tuning
#' @return \code{STLS_fit} object
#' @import methods
#' @import Rcpp
#' @import ggplot2
#' @importClassesFrom Matrix dgCMatrix
#' @importMethodsFrom Matrix dim as.matrix coerce
#' @export
fitSTLS = function(X, y,
center = TRUE,
lambda = NULL,
nlambda = 100,
lmin_ratio = 1e-4,
alpha = 1,
eps_abs = 1e-4,
eps_rel = 1e-4,
maxit = 1000L,
warm_start = TRUE) {
fit = PALM_STLS(X, y)
fit$penalty(lambda, nlambda, lmin_ratio, alpha)
fit$options(maxit, eps_abs, eps_rel, center, warm_start)
res = fit$fit()
res
}
### cross-validation configuration
cv.PALM_STLS$methods(
initialize = function(X, y, folds, ...) {
if (nrow(X) != length(y)) {
stop("nrow(X) should be equal to size of y.")
}
.self$foldix = llply(seq(max(folds)), function(x) {
which(folds != x) - 1
})
.self$X = as.matrix(X)
.self$y = as.numeric(y)
.self$lambda = numeric(0)
.self$nlambda = 100L
.self$lmin_ratio = 1e-5
.self$alpha = 1.0
.self$maxit = 10000L
.self$eps_abs = 1e-6
.self$eps_rel = 1e-6
.self$center = TRUE
.self$warm_start = FALSE
}
)
## print cv.PALM_STLS
cv.PALM_STLS$methods(
show_common = function() {
cat(sprintf("X := %d x %d Matrix", nrow(.self$X), ncol(.self$X)))
cat(sprintf("y := %d Vector"), length(.self$y))
fields = setdiff(names(.refClassDef@fieldClasses), c("X", "y"))
for(field in fields) {
cat(field, " : ", paste(.self$field(field), collapse = " "), "\n", sep = "")
}
},
show = function()
{
cat("CV TLS model in PALM \n \n")
show_common()
}
)
## penalty parameter set-up
cv.PALM_STLS$methods(
penalty = function(lambda = NULL, nlambda = 100, lmin_ratio = 1e-5, alpha = 1.0, ...) {
lambda_val = sort(as.numeric(lambda), decreasing = TRUE)
if(length(lambda_val) > 0 & any(lambda_val < 0)) {
stop("lambda must be non-negative value.")
}
if(nlambda <= 0) {
stop("nlambda must be a positive integer.")
}
if(alpha < 0 | alpha > 1) {
stop("alpha must be located in (0, 1].")
}
.self$lambda = lambda_val
.self$nlambda = nlambda
.self$lmin_ratio = lmin_ratio
.self$alpha = alpha
invisible(.self)
}
)
## options configuration
cv.PALM_STLS$methods(
options = function(maxit = 1000L, eps_abs = 1e-6, eps_rel = 1e-6, center = TRUE, warm_start = FALSE, ...) {
if(maxit <= 0) {
stop("maxit should be positive integer.")
}
if(eps_abs < 0 | eps_rel < 0) {
stop("eps_abs/rel should be positive.")
}
.self$maxit = maxit
.self$eps_abs = eps_abs
.self$eps_rel = eps_rel
.self$center = center
.self$warm_start = warm_start
invisible(.self)
}
)
cv.PALM_STLS$methods(
fit = function(...) {
res = .Call("CVPALM_STLS", .self$X, .self$y, .self$foldix, .self$lambda, .self$nlambda, .self$lmin_ratio, .self$alpha,
list(maxit = .self$maxit,
eps_abs = .self$eps_abs,
eps_rel = .self$eps_rel,
center = .self$center,
warm_start = .self$warm_start),
PACKAGE = "PALM")
cvtest = list()
n = nrow(.self$X)
for(i in 1:length(res$cv.result)) {
ix = setdiff(seq(n), res$cv.train[[i]] + 1)
xtest = .self$X[ix, , drop =]
ytest = .self$y[ix]
tmp = xtest %*% res$cv.result[[i]]$beta
obj = vector("numeric", length(res$cv.result[[i]]$intercept))
for(j in 1:length(res$cv.result[[i]]$intercept)) {
obj[j] = crossprod(ytest - (tmp[,j] + res$cv.result[[i]]$intercept[j]))
obj[j] = 1/2 * obj[j] / (1 + sum((res$cv.result[[i]]$beta[,j])^2))
}
cvtest[[i]] = obj
}
cvtest = do.call(rbind, cvtest)
res[["cvm"]] = colMeans(cvtest)
res[["cvsd"]] = apply(cvtest, 2, sd)
res
}
)
#' Cross-validation for Sparse Total Least Square model Using PALM algorithm
#'
#' @title cv.fitSTLS
#' @param X The observation data matrix
#' @param y The observation response vector
#' @param nfolds Number of folds in cross-validation
#' @param foldid Folds id in N-fold cross-validation
#' @param center Y = Y - mean(Y) and X = X - mean(X)
#' @param lambda The regularized lambda vector
#' @param nlambda The Number of lambda
#' @param lmin_ratio lambda_max / lambda_min
#' @param alpha alpha parameter for elastic net penalty (0,1](Ridge Regression --> Lasso)
#' @param maxit The maximal iteration of PALM algorithm
#' @param eps_abs Absolutely epsilon for generated varaince vector
#' @param eps_rel Relative epsilon for objective function value of generated vector
#' @param warm_start Warm start for regularizer parameter tuning
#' @return \code{STLS_fit} object
#' @import methods
#' @import Rcpp
#' @import ggplot2
#' @importClassesFrom Matrix dgCMatrix
#' @importMethodsFrom Matrix dim as.matrix coerce
#' @export
cv.fitSTLS = function(X, y,
nfolds = 5,
foldid = NULL,
center = TRUE,
lambda = NULL,
nlambda = 100,
lmin_ratio = 1e-4,
alpha = 1,
eps_abs = 1e-4,
eps_rel = 1e-4,
maxit = 1000L,
warm_start = TRUE) {
if(nfolds >= 3 & is.null(foldid)) {
foldid = sample(seq(nfolds), size = nrow(X), replace = T)
} else if(!is.null(foldid)) {
nfolds = max(foldid)
}
cvfit = cv.PALM_STLS(X, y, foldid)
cvfit$penalty(lambda, nlambda, lmin_ratio, alpha)
cvfit$options(maxit, eps_abs, eps_rel, center, warm_start)
cvfit = cvfit$fit()
# get lambda.1se
opt.lambda = minLambda(cvfit$lambda, cvfit$cvm, cvfit$cvsd)
res = list(lambda = cvfit$lambda, cvres = cvfit$cv.result, cvtrain = cvfit$cv.train, cvm = cvfit$cvm, cvsd = cvfit$cvsd,
lambda.min = cvfit$lambda[which.min(cvfit$cvm)], lambda.1se = opt.lambda$lambda.1se, lambda.opt = opt.lambda$lambda.opt)
do.call(cv.STLS_fit, res)
}
##' plot cross-validation curive and search for lambda.min
##'
##' @rdname plot.cv.STLS_fit
##' @param x fitted cv.STLS_fit
##' @method plot cv.STLS_fit
##' @export
plot.cv.STLS_fit = function(x, ...) {
cvSTLS = data.frame(lambda = log10(x$lambda), cvm = x$cvm, cerr = x$cvsd)
ggplot(cvSTLS, aes(x = lambda, y = cvm)) + geom_errorbar(aes(ymin = cvm - cerr, ymax = cvm + cerr)) +
geom_line() + geom_point() + scale_y_log10() +
xlab(expression(log(lambda))) + ylab("objective value") +
ggtitle("cv.obj ~ lambda") +
geom_vline(xintercept = log10(x$lambda.min), linetype = "dashed") +
scale_colour_PALM() + theme_PALM()
}
### stability-selection configuration
ss.PALM_STLS$methods(
initialize = function(X, y, nperm = 100, ...) {
if (nrow(X) != length(y)) {
stop("nrow(X) should be equal to size of y.")
}
.self$X = as.matrix(X)
.self$y = as.numeric(y)
n = nrow(X)
k = floor(n / 2)
index = rep(c(0, 1), c(n-k, k))
permB = replicate(nperm, sample(index))[sample(1:n), , drop = FALSE]
.self$permix = llply(1:nperm, function(x) {
which(permB[,x] != 0)
})
.self$lambda = numeric(0)
.self$nlambda = 100L
.self$lmin_ratio = 1e-5
.self$alpha = 1.0
.self$maxit = 10000L
.self$eps_abs = 1e-6
.self$eps_rel = 1e-6
.self$center = TRUE
.self$warm_start = FALSE
}
)
## print ss.PALM_STLS
ss.PALM_STLS$methods(
show_common = function() {
cat(sprintf("X := %d x %d Matrix", nrow(.self$X), ncol(.self$X)))
cat(sprintf("y := %d Vector"), length(.self$y))
fields = setdiff(names(.refClassDef@fieldClasses), c("X", "y"))
for(field in fields) {
cat(field, " : ", paste(.self$field(field), collapse = " "), "\n", sep = "")
}
},
show = function()
{
cat("SS TLS model in PALM \n \n")
show_common()
}
)
## penalty parameter set-up
ss.PALM_STLS$methods(
penalty = function(lambda = NULL, nlambda = 100, lmin_ratio = 1e-5, alpha = 1.0, ...) {
lambda_val = sort(as.numeric(lambda), decreasing = TRUE)
if(length(lambda_val) > 0 & any(lambda_val < 0)) {
stop("lambda must be non-negative value.")
}
if(nlambda <= 0) {
stop("nlambda must be a positive integer.")
}
if(alpha < 0 | alpha > 1) {
stop("alpha must be located in (0, 1].")
}
.self$lambda = lambda_val
.self$nlambda = nlambda
.self$lmin_ratio = lmin_ratio
.self$alpha = alpha
invisible(.self)
}
)
## options configuration
ss.PALM_STLS$methods(
options = function(maxit = 1000L, eps_abs = 1e-6, eps_rel = 1e-6, center = TRUE, warm_start = FALSE, ...) {
if(maxit <= 0) {
stop("maxit should be positive integer.")
}
if(eps_abs < 0 | eps_rel < 0) {
stop("eps_abs/rel should be positive.")
}
.self$maxit = maxit
.self$eps_abs = eps_abs
.self$eps_rel = eps_rel
.self$center = center
.self$warm_start = warm_start
invisible(.self)
}
)
ss.PALM_STLS$methods(
fit = function(...) {
res = .Call("SSPALM_STLS", .self$X, .self$y, .self$permix, .self$lambda, .self$nlambda, .self$lmin_ratio, .self$alpha,
list(maxit = .self$maxit,
eps_abs = .self$eps_abs,
eps_rel = .self$eps_rel,
center = .self$center,
warm_start = .self$warm_start),
PACKAGE = "PALM")
res
}
)
#' Stability-Selection for Sparse Total Least Square model Using PALM algorithm
#'
#' @title ss.fitSTLS
#' @param X The observation data matrix
#' @param y The observation response vector
#' @param nperm Number of permutaiton of stabilioty-selection
#' @param center Y = Y - mean(Y) and X = X - mean(X)
#' @param lambda The regularized lambda vector
#' @param nlambda The Number of lambda
#' @param lmin_ratio lambda_max / lambda_min
#' @param alpha alpha parameter for elastic net penalty (0,1](Ridge Regression --> Lasso)
#' @param maxit The maximal iteration of PALM algorithm
#' @param eps_abs Absolutely epsilon for generated varaince vector
#' @param eps_rel Relative epsilon for objective function value of generated vector
#' @param warm_start Warm start for regularizer parameter tuning
#' @return \code{STLS_fit} object
#' @import methods
#' @import Rcpp
#' @import ggplot2
#' @importClassesFrom Matrix dgCMatrix
#' @importMethodsFrom Matrix dim as.matrix coerce
#' @export
ss.fitSTLS = function(X, y,
nperm = 100,
center = TRUE,
lambda = NULL,
nlambda = 100,
lmin_ratio = 1e-4,
alpha = 1,
eps_abs = 1e-4,
eps_rel = 1e-4,
maxit = 1000L,
warm_start = TRUE) {
ssfit = ss.PALM_STLS(X, y, nperm)
ssfit$penalty(lambda, nlambda, lmin_ratio, alpha)
ssfit$options(maxit, eps_abs, eps_rel, center, warm_start)
ssfit = ssfit$fit()
ssfit[["ss.prob"]] = llply(ssfit$ss.beta, `[[`, "beta")
ssfit[["ss.prob"]] = Reduce(`+`, llply(ssfit[["ss.prob"]], function(x){as.matrix(x != 0)}))
ssfit[["ss.prob"]] = rowSums(ssfit[["ss.prob"]]) / (nperm * nlambda)
ssfit
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.