Multivariate partly reduced-rank regression

rrr_comb.R

##' Multivariate partly reduced-rank regression
##'
##' @usage
##' rrr_comb(Y, X, Z, ...)
##'
##' @param Y a matrix of response (n by q)
##' @param X a matrix of covariates whose association with the responses has reduced rank (n by p1)
##' @param Z a matrix of covariates whose association with the responses has full rank (n by p2)
##' All other parameters will be passed to \code{rrr}.
##'
##' @return
##'
##' S3 \code{rrr_comb} object, a list consisting of
##' \item{rrr}{the object returned by \code{rrr}}
##' \item{X, Y, Z}{original input matrices}
##' \item{D}{least-squares regression coefficient matrix for the full-rank portion (p2 by q))}
##'
##' @references
##'
##' Reinsel, G. C. & Velu, R. P. (1998)
##' \emph{Multivariate Reduced-Rank Regression : Theory and Applications}, Springer, Chapter 3.
##'
##' @examples
##' library(rrpack)
##' p <- 50; q <- 50; n <- 100; nrank <- 3
##' mydata <- rrr.sim1(n, p, q, nrank, s2n = 1, sigma = NULL,
##'                    rho_X = 0.5, rho_E = 0.3)
##' mydata$Z <- matrix(1, n, 1) #intercept term
##' rcombfit <- with(mydata, rrr_comb(Y, X, Z, maxrank = 10))
##' sse <- sum((rcombfit$Y - rcombfit$Z %*% rcombfit$D - rcombfit$rrr$X %*% rcombfit$rrr$coef) ^ 2)
##' summary(rcombfit)
##' @export
##'
##
## Copyright 2018, Nir Krakauer, mail@nirkrakauer.net
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.

rrr_comb <- function(Y, X, Z, ...)
{
  ## record function call
  Call <- match.call()

  library(pracma) #for mldivide
  rfit <- rrr(Y -  Z %*% mldivide(Z, Y), X - Z %*% mldivide(Z, X), ...)

  if (rfit$rank != 0) {
    D <- mldivide(Z, Y - rfit$X %*% rfit$coef)
  } else {
    D <- mldivide(Z, Y)
  }

  rval <- list(
    call = Call,
    Y = Y,
    X = X,
    Z = Z,
    rrr = rfit,
    D = D
  )
  class(rval) <- "rrr_comb"
}