R/cancor.R

In friendly/candisc: Visualizing Generalized Canonical Discriminant and Canonical Correlation Analysis

# canonical correlation analysis, as a more general method, with method functions
# Updated to implement weights
#
# compare with :
#      stats::cancor (very basic),
#      yacca::cca (fairly complete, but very messy return structure)
#      CCA::cc (fairly complete, but very messy return structure, no longer maintained)

cancor <- function(x, ...) {
  UseMethod("cancor", x)
}

cancor.formula <- function(formula, data, subset, weights,
                           na.rm = TRUE,           # DONE: now just use na.rm
                           method = "gensvd",      # "qr" not implemented
                           # 		contrasts = NULL,  # would it make any sense to allow factors? should we test for factors in X?
                           ...) {

  # remove the intercept [solution by John Fox]
  formula <- update(formula, . ~ . - 1)
  cl <- match.call()
  cl$formula <- formula
  mf <- match.call(expand.dots = FALSE)
  mf$formula <- formula
  m <- match(c("formula", "data", "subset", "weights"), names(mf), 0L)
  mf <- mf[c(1L, m)]


  mf[[1L]] <- as.name("model.frame")
  mf <- eval(mf, parent.frame())
  mt <- attr(mf, "terms")

  y <- model.response(mf, "numeric")
  w <- as.vector(model.weights(mf))
  if (!is.null(w) && !is.numeric(w)) {
    stop("'weights' must be a numeric vector")
  }

  x <- model.matrix(mt, mf, contrasts)
  z <- cancor.default(x, y, weights = w, na.rm = na.rm, ...)

  z$call <- cl
  z$terms <- mt
  z
}


# var-cov matrix allowing weights
# Without weights, honors na.rm and use
# With weights, na.rm --> use='complete'

Var <- function(x, na.rm = TRUE, use, weights) {
  if (missing(weights) || is.null(weights)) {
    res <- var(x, na.rm = na.rm, use = use)
  } else {
    # cov.wt doesn't allow for missing data, and doesn't allow x, y=NULL, ...
    if (na.rm) {
      if (!pmatch(use, "complete")) warning("Use of weights only supports use='complete'")
      OK <- complete.cases(x)
      x <- x[OK, ]
    }
    res <- cov.wt(x, wt = weights)$cov
  }
  res
}

# DONE:  add weights component to result
# DONE:  add a method=argument
# TODO:  should replace X, Y by x, y throughout to avoid another copy
# DONE:  allow weights, by use of cov.wt() as in dataEllipse()
# FIXED: row.names should refer to original x, because as.matrix() strips them

cancor.default <- function(x, y,
                           weights,
                           X.names = colnames(x),
                           Y.names = colnames(y),
                           row.names = rownames(x),
                           xcenter = TRUE, ycenter = TRUE, # not yet implemented (used implicitly)
                           xscale = FALSE, yscale = FALSE, # not yet implemented
                           ndim = min(p, q),
                           set.names = c("X", "Y"),
                           prefix = c("Xcan", "Ycan"), # s/b: paste0(set.names, "can")
                           na.rm = TRUE,
                           use = if (na.rm) "complete" else "pairwise",
                           method = "gensvd",
                           ...) {
  X <- as.matrix(x)
  Y <- as.matrix(y)
  p <- ncol(X)
  q <- ncol(Y)
  n <- length(complete.cases(X, Y)) # DONE: take account of 0 weights
  if (!missing(weights)) n <- n - sum(weights == 0)

  C <- Var(cbind(X, Y), na.rm = TRUE, use = use, weights = weights)
  Cxx <- C[1:p, 1:p]
  Cyy <- C[-(1:p), -(1:p)]
  Cxy <- C[1:p, -(1:p)]

  res <- gensvd(Cxy, Cxx, Cyy, nu = ndim, nv = ndim)
  names(res) <- c("cor", "xcoef", "ycoef")
  colnames(res$xcoef) <- paste(prefix[1], 1:ndim, sep = "")
  colnames(res$ycoef) <- paste(prefix[2], 1:ndim, sep = "")

  scores <- can.scores(X, Y, res$xcoef, res$ycoef)
  colnames(scores$xscores) <- paste(prefix[1], 1:ndim, sep = "")
  colnames(scores$yscores) <- paste(prefix[2], 1:ndim, sep = "")

  structure <- can.structure(X, Y, scores, use = use)
  result <- list(
    cancor = res$cor,
    names = list(
      X = X.names, Y = Y.names,
      row.names = row.names, set.names = set.names
    ),
    ndim = ndim,
    dim = list(p = p, q = q, n = n),
    coef = list(X = res$xcoef, Y = res$ycoef),
    scores = list(X = scores$xscores, Y = scores$yscores),
    X = X, Y = Y,
    weights = if (missing(weights)) NULL else weights,
    structure = structure
  )
  class(result) <- "cancor"
  return(result)
}

# scores on canonical variates
can.scores <- function(X, Y, xcoef, ycoef) {
  X.aux <- scale(X, center = TRUE, scale = FALSE) # TODO: incorporate xscale, yscale, etc here
  Y.aux <- scale(Y, center = TRUE, scale = FALSE)
  X.aux[is.na(X.aux)] <- 0
  Y.aux[is.na(Y.aux)] <- 0

  xscores <- X.aux %*% xcoef
  yscores <- Y.aux %*% ycoef
  return(list(xscores = xscores, 
              yscores = yscores))
}

# canonical structure coefficients: correlations
can.structure <- function(X, Y, scores, use = "complete.obs") {
  xscores <- scores$xscores
  yscores <- scores$yscores
  X.xscores <- cor(X, xscores, use = use)
  Y.xscores <- cor(Y, xscores, use = use)
  X.yscores <- cor(X, yscores, use = use)
  Y.yscores <- cor(Y, yscores, use = use)

  return(list(
    X.xscores = X.xscores,
    Y.xscores = Y.xscores,
    X.yscores = X.yscores,
    Y.yscores = Y.yscores
  ))
}

# TODO: replace with equivalent qr stuff
gensvd <- function(Rxy, Rxx, Ryy, nu = p, nv = q) {
  p <- dim(Rxy)[1]
  q <- dim(Rxy)[2]

  if (missing(Rxx)) Rxx <- diag(p)
  if (missing(Ryy)) Ryy <- diag(q)

  if (dim(Rxx)[1] != dim(Rxx)[2]) stop("Rxx must be square")
  if (dim(Ryy)[1] != dim(Ryy)[2]) stop("Ryy must be square")

  s <- min(p, q)
  if (max(abs(Rxx - t(Rxx))) / max(abs(Rxx)) > 1e-10) {
    warning("Rxx not symmetric.")
    Rxx <- (Rxx + t(Rxx)) / 2
  }
  if (max(abs(Ryy - t(Ryy))) / max(abs(Ryy)) > 1e-10) {
    warning("Ryy not symmetric.")
    Ryy <- (Ryy + t(Ryy)) / 2
  }

  Rxxinv <- solve(chol(Rxx))
  Ryyinv <- solve(chol(Ryy))
  Dform <- t(Rxxinv) %*% Rxy %*% Ryyinv
  if (p >= q) {
    result <- svd(Dform, nu = nu, nv = nv)
    values <- result$d
    Xmat <- Rxxinv %*% result$u
    Ymat <- Ryyinv %*% result$v
  } else {
    result <- svd(t(Dform), nu = nv, nv = nu)
    values <- result$d
    Xmat <- Rxxinv %*% result$v
    Ymat <- Ryyinv %*% result$u
  }
  gsvdlist <- list(values = values, Xmat = Xmat, Ymat = Ymat)
  return(gsvdlist)
}


# vector of stars, of max. width, corresponding to proportions
stars <- function(p, width = 30) {
  p <- p / sum(p)
  reps <- round(p * width / max(p))
  res1 <- sapply(reps, function(x) paste(rep("*", x), sep = "", collapse = ""))
  res2 <- sapply(reps, function(x) paste(rep(" ", width - x), sep = "", collapse = ""))
  res <- paste(res1, res2, sep = "")
  res
}

# DONE: move the printout of coefficients to a summary method
print.cancor <- function(x, digits = max(getOption("digits") - 2, 3), ...) {
  names <- x$names
  cat("\nCanonical correlation analysis of:\n")
  cat("\t", x$dim$p, " ", names$set.names[1], " variables: ", paste(names$X, collapse = ", "), "\n")
  cat("  with\t", x$dim$q, " ", names$set.names[2], " variables: ", paste(names$Y, collapse = ", "), "\n")
  cat("\n")
  canr <- x$cancor
  lambda <- canr^2 / (1 - canr^2)
  pct <- 100 * lambda / sum(lambda)
  cum <- cumsum(pct)
  # DONE: add stars column, showing pct
  stwidth <- getOption("width") - 40
  scree <- stars(pct, width = min(30, stwidth))
  canrdf <- data.frame("CanR" = canr, 
                       "CanRSQ" = canr^2, 
                       "Eigen" = lambda, 
                       "percent" = pct, 
                       "cum" = cum, 
                       "scree" = scree)
  print(canrdf, digits = 4)

  tests <- Wilks.cancor(x)
  print(tests, digits = digits)

  invisible(x)
}

# moved printout of coefficients here
summary.cancor <- function(object, digits = max(getOption("digits") - 2, 3), ...) {
  names <- object$names
  print(object, digits = digits, ...)

  cat("\nRaw canonical coefficients\n")
  cat("\n  ", names$set.names[1], " variables: \n")
  print(object$coef$X, digits = digits)
  cat("\n  ", names$set.names[2], " variables: \n")
  print(object$coef$Y, digits = digits)
}



# extractor functions
scores <- function(x, ...) {
  UseMethod("scores")
}

# TODO: check rownames
scores.cancor <- function(x, type = c("x", "y", "both", "list", "data.frame"), ...) {
  type <- match.arg(type)
  switch(type,
    x = x$scores$X,
    y = x$scores$Y,
    both = x$scores,
    list = x$scores,
    data.frame = data.frame(x$scores$X, x$scores$Y)
  )
}

coef.cancor <- function(object, type = c("x", "y", "both", "list"), standardize = FALSE, ...) {
  coef <- object$coef
  if (standardize) {
    coef$X <- diag(sqrt(diag(cov(object$X)))) %*% coef$X
    coef$Y <- diag(sqrt(diag(cov(object$Y)))) %*% coef$Y
    rownames(coef$X) <- rownames(object$coef$X)
    rownames(coef$Y) <- rownames(object$coef$Y)
  }
  type <- match.arg(type)
  switch(type,
    x = coef$X,
    y = coef$Y,
    both = list(coef$X, coef$Y),
    list = list(coef$X, coef$Y)
  )
}