compCL2 <- function(y, Z, Zc = NULL, intercept = TRUE,
pf = rep(1, times = p),
lam = NULL, nlam = 100, lambda.factor = ifelse(n < p, 0.05, 0.001),
dfmax = p, pfmax = min(dfmax * 1.5, p),
u = 1,
mu_ratio = 1.01, tol = 1e-10,
outer_maxiter = 1e+08, outer_eps = 1e-8,
inner_maxiter = 1e+3, inner_eps = 1e-6
) {
#u <- 1
this.call <- match.call()
y <- drop(y)
Z <- as.matrix(Z)
# if( any(abs(rowSums(Z) - 1) > 1e-10) ) {
# message("Z is transformed into compositional data by deviding rowSums")
# Z <- Z / rowSums(Z)
# }
# if(any(Z == 0)) stop("There is zero entry in compositional data")
n <- length(y)
p <- dim(Z)[2]
inter <- as.integer(intercept)
Znames <- colnames(Z)
if (is.null(Znames)) Znames <- paste0("Z", seq(p))
Zc <- as.matrix(cbind(Zc, rep(inter, n)))
if(inter == 1) {
Zc_proj <- solve(crossprod(Zc)) %*% t(Zc)
} else {
if(dim(Zc)[2] == 1) Zc_proj <- t(Zc) else Zc_proj <- rbind(solve(crossprod(Zc[, -ncol(Zc)])) %*% t(Zc[, -ncol(Zc)]), rep(inter, n))
}
beta_c <- as.vector(Zc_proj %*% y)
beta.ini <- c(rep(0, times = p), beta_c)
m <- dim(Zc)[2]
if(m > 1) {
Zcnames <- colnames(Z)
if (is.null(Zcnames)) Zcnames <- paste0("Zc", seq(m-1))
} else {
Zcnames <- NULL
}
if(is.null(lam)) {
lam0 <- t(Z) %*% (y - Zc %*% beta_c) / n
lam0 <- max(abs(lam0) / pf)
lam <- exp(seq(from=log(lam0), to=log(lambda.factor * lam0),length=nlam))
}
pfmax <- as.integer(pfmax)
dfmax <- as.integer(dfmax)
inner_maxiter <- as.integer(inner_maxiter)
outer_maxiter <- as.integer(outer_maxiter)
inner_eps <- as.double(inner_eps)
outer_eps <- as.double(outer_eps)
tol <- as.double(tol)
u <- as.double(u)
mu_ratio <- as.double(mu_ratio)
#estar <- as.double(estar)
output <- ALM_B(y = y, Z = Z, Zc = Zc, Zc_proj = Zc_proj, beta = beta.ini,
lambda = lam, pf = pf, dfmax = dfmax, pfmax = pfmax,
inner_eps = inner_eps, outer_eps = outer_eps,
inner_maxiter = inner_maxiter, outer_maxiter = outer_maxiter,
u_ini = u, mu_ratio = mu_ratio, tol = tol
)
output$call <- this.call
class(output) <- "compCL"
output$dim <- dim(output$beta)
output$lam <- drop(output$lam)
dimnames(output$beta) <- list(c(Znames, Zcnames, "Intercept"), paste0("L", seq(along = output$lam)) )
return(output)
}
cv.compCL2 <- function(y, Z, Zc = NULL, intercept = FALSE,
lam = NULL,
nfolds = 10, foldid,
trim = 0.1, ...) {
y <- drop(y)
n <- length(y)
if (missing(foldid)) foldid <- sample(rep(seq(nfolds), length = n)) else nfolds <- max(foldid)
if (nfolds < 3) stop("nfolds must be bigger than 3; nfolds=10 recommended")
compCL.object <- compCL2(y = y, Z = Z, Zc = Zc, intercept = intercept,
lam = lam, ...)
cv.lam <- compCL.object$lam
outlist <- as.list(seq(nfolds))
###Now fit the nfold models and store them
for (i in seq(nfolds)) {
which <- foldid == i
outlist[[i]] <- compCL2(y = y[!which], Z = Z[!which, , drop = FALSE], Zc = Zc[!which, , drop = FALSE], intercept = intercept,
lam = cv.lam, ...)
}
cvstuff <- cv.test2(outlist, y, Z = Z, Zc = Zc, foldid, lam = cv.lam, trim = trim)
cvm <- drop(cvstuff$cvm)
cvsd <- drop(cvstuff$cvsd)
if(trim > 0) {
cvm.trim <- drop(cvstuff$cvmtrim)
cvsd.trim <- drop(cvstuff$cvsdtrim)
}
Ftrim = list(cvm = cvm, cvsd = cvsd, cvupper = cvm + cvsd, culo = cvm - cvsd)
lam.min <- getmin(lam = cv.lam, cvm = cvm, cvsd = cvsd)
Ftrim <- c(Ftrim, lam.min)
if(trim > 0) {
Ttrim <- list(cvm = cvm.trim, cvsd = cvsd.trim, cvupper = cvm.trim + cvsd.trim,
cvlo = cvm.trim-cvsd.trim)
lam.min <- getmin(lam = cv.lam, cvm = cvm.trim, cvsd = cvsd.trim)
Ttrim <- c(Ttrim, lam.min)
} else {
Ttrim <- NULL
}
result <- list(compCL.fit = compCL.object,
lam = cv.lam,
Ftrim = Ftrim,
Ttrim = Ttrim,
foldid = foldid)
class(result) <- "cv.compCL"
return(result)
}
GIC.compCL2 <- function(y, Z, Zc = NULL, intercept = FALSE,
lam = NULL, ...) {
digits = 5
y <- drop(y)
n <- length(y)
compCL.object <- compCL2(y = y, Z = Z, Zc = Zc, intercept = intercept, lam = lam, ...)
#print(compCL.object)
lam <- compCL.object$lam
GIC <- GIC.test2(object = compCL.object, y = y, Z = Z, Zc = Zc, intercept = intercept)
#print(GIC.test2(object = compCL.object, y = y, Z = Z, Zc = Zc, intercept = intercept))
GIC1 <- round(GIC, digits = digits)
GIC.min <- min(GIC1)
idmin <- GIC1 <= GIC.min
lam.min <- max(lam[idmin])
result <- list(compCL.fit = compCL.object,
lam = lam,
GIC = GIC,
lam.min = lam.min)
class(result) <- "GIC.compCL"
return(result)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.