R/spls.cox.R

In bigPLScox: Partial Least Squares for Cox Models with Big Matrices

spls.cox=function (x, y, K, eta, kappa = 0.5, select = "pls2", fit = "regression", 
                   scale.x = TRUE, scale.y = FALSE, eps = 1e-04, maxstep = 100, 
                   trace = FALSE) 
{
  force(K)
  x <- as.matrix(x)
  n <- nrow(x)
  p <- ncol(x)
  ip <- c(1:p)
  y <- as.matrix(y)
  q <- ncol(y)
  one <- matrix(1, 1, n)
  mu <- one %*% y/n
  y <- scale(y, drop(mu), FALSE)
  meanx <- drop(one %*% x)/n
  x <- scale(x, meanx, FALSE)
  if (scale.x) {
    normx <- sqrt(drop(one %*% (x^2))/(n - 1))
    if (any(normx < .Machine$double.eps)) {
      stop("Some of the columns of the predictor matrix have zero variance.")
    }
    x <- scale(x, FALSE, normx)
  }
  else {
    normx <- rep(1, p)
  }
  if (scale.y) {
    normy <- sqrt(drop(one %*% (y^2))/(n - 1))
    if (any(normy < .Machine$double.eps)) {
      stop("Some of the columns of the response matrix have zero variance.")
    }
    y <- scale(y, FALSE, normy)
  }
  else {
    normy <- rep(1, q)
  }
  betahat <- matrix(0, p, q)
  betamat <- list()
  x1 <- x
  y1 <- y
  type <- correctp.cox(x, y, eta, K, kappa, select, fit)
  eta <- type$eta
  K <- type$K
  kappa <- type$kappa
  select <- type$select
  fit <- type$fit
  if (is.null(colnames(x))) {
    xnames <- c(1:p)
  }
  else {
    xnames <- colnames(x)
  }
  new2As <- list()
  if (trace) {
    cat("The variables that join the set of selected variables at each step:\n")
  }
  for (k in 1:K) {
    Z <- t(x1) %*% y1
    what <- spls.dv(Z, eta, kappa, eps, maxstep)
    A <- unique(ip[what != 0 | betahat[, 1] != 0])
    new2A <- ip[what != 0 & betahat[, 1] == 0]
    xA <- x[, A, drop = FALSE]
    plsfit <- pls.cox(X=xA, Y=y, ncomp = min(k, length(A)), mode = fit, 
                      scale.X = FALSE, scale.Y=FALSE)
    predplsfit <- predict.pls.cox(plsfit,newdata=xA,scale.X = FALSE, scale.Y=FALSE)
    betahat <- matrix(0, p, q)
    betahat[A, ] <- matrix(predplsfit$B.hat[,,plsfit$ncomp], length(A), q)
    betamat[[k]] <- betahat
    #        pj <- plsfit$projection
    if (select == "pls2") {
      y1 <- y - predplsfit$predict[,,plsfit$ncomp]
    }
    #        if (select == "simpls") {
    #            pw <- pj %*% solve(t(pj) %*% pj) %*% t(pj)
    #            x1 <- x
    #            x1[, A] <- x[, A, drop = FALSE] - x[, A, drop = FALSE] %*% 
    #                pw
    #        }
    new2As[[k]] <- new2A
    if (trace) {
      if (length(new2A) <= 10) {
        cat(paste("- ", k, "th step (K=", k, "):\n", 
                  sep = ""))
        cat(xnames[new2A])
        cat("\n")
      }
      else {
        cat(paste("- ", k, "th step (K=", k, "):\n", 
                  sep = ""))
        nlines <- ceiling(length(new2A)/10)
        for (i in 0:(nlines - 2)) {
          cat(xnames[new2A[(10 * i + 1):(10 * (i + 1))]])
          cat("\n")
        }
        cat(xnames[new2A[(10 * (nlines - 1) + 1):length(new2A)]])
        cat("\n")
      }
    }
  }
  if (!is.null(colnames(x))) {
    rownames(betahat) <- colnames(x)
  }
  if (q > 1 & !is.null(colnames(y))) {
    colnames(betahat) <- colnames(y)
  }
  object <- list(x = x, y = y, betahat = betahat, A = A, betamat = betamat, 
                 new2As = new2As, mu = mu, meanx = meanx, normx = normx, 
                 normy = normy, eta = eta, K = K, kappa = kappa, select = select, 
                 fit = fit, projection = NA, plsmod=plsfit)
  class(object) <- "spls"
  object
}

ust=function (b, eta) 
{
  b.ust <- matrix(0, length(b), 1)
  if (eta < 1) {
    valb <- abs(b) - eta * max(abs(b))
    b.ust[valb >= 0] <- valb[valb >= 0] * (sign(b))[valb >= 
                                                      0]
  }
  return(b.ust)
}