Run_simulation/compCL_test/compCL2.R

In jiji6454/Rpac_compReg: Compositional Data Regression

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