R/basecor.R

In graphpcor: Models for Correlation Matrices Based on Graphs

#' A base Cholesky model of a correlation matrix
#' @description
#' The `basecor` class contain a correlation matrix `base`,
#' the parameter vector `theta`, that generates
#' or is generated by `base`, the dimention `p`,
#' the index `iLtheta` for `theta` in the (lower) Cholesky,
#' and the Hessian around it `I0`, see details.

#' @describeIn basecor
#' Build a `basecor` object.
#' @param base numeric vector/matrix used to define the base
#' correlation matrix. If numeric vector with length 'm',
#' 'm' should be 'p(p-1)/2' in the dense model case and
#' 'length(iLtheta)' in the sparse model case.
#' @param p integer with the dimension,
#' the number of rows/columns of the correlation matrix.
#' @param parametrization character to specify the
#' parametrization used. The available ones are
#' "cpc" (or "CPC") or "sap" (or "SAP").
#' See Details. The default is "cpc".
#' @param iparams integer ordered vector with length equal
#' the number of parameters used to specify common parameter values.
#' If missing, assumed to be `1:length(theta)`. Example: By setting
#' `iparams = c(1,1,2,3)` the first and second parameters are
#' considered to be the same.
#' @param iLtheta integer vector to specify the (vectorized) position
#' where 'theta' will be placed in the (lower triangle)  Cholesky
#' factorization of the correlation matrix. Default
#' (missing, or if NULL) and assumes `iLtheta = which(lower.tri(...))`.
#' @details
#' For 'parametrization' = "CPC" or 'parametrization' = "cpc":
#' The Canonical Partial Correlation - CPC parametrization,
#'  Lewandowski, Kurowicka, and Joe (2009), compute
#' \eqn{r[k]} = tanh(\eqn{\theta[k]}), for \eqn{k=1,...,m},
#' and the two \eqn{p\times p} matrices
#' \deqn{A = \left[
#' \begin{array}{ccccc}
#'   1 & & & & \\
#'   r_1 & 1 & & & \\
#'   r_2 & r_p & 1 & & \\
#'   \vdots & \vdots & \ddots & \ddots & \\
#'   r_{p-1} & r_{2p-3} & \ldots & r_m & 1
#' \end{array} \right]
#' \textrm{ and } B = \left[
#' \begin{array}{ccccc}
#'   1 & & & & \\
#'   \sqrt{1-r_1^2} & 1 & & & \\
#'   \sqrt{1-r_2^2} & \sqrt{1-r_p^2} & 1 & & \\
#'   \vdots & \vdots & \ddots & \ddots & \\
#'   \sqrt{1-r_{p-1}^2} & \sqrt{r_{2p-3}^2} & \ldots & \sqrt{1-r_m^2} & 1
#' \end{array} \right] }
#'
#' The matrices \eqn{A} and \eqn{B} are then used
#' to build the Cholesky factor of the correlation matrix,
#' given as
#' \deqn{L = \left[
#' \begin{array}{ccccc}
#'   1 & 0 & 0 & \ldots & 0\\
#'   A_{2,1} & B_{2,1} & 0 & \ldots & 0\\
#'   A_{3,1} & A_{3,2}B_{3,1} & B_{3,1}B_{3,2} & & \vdots \\
#'   \vdots & \vdots & \ddots & \ddots & 0\\
#'   A_{p,1} & A_{p,2}B_{p,1} & \ldots &
#'   A_{p,p-1}\prod_{k=1}^{p-1}B_{p,k} & \prod_{k=1}^{p-1}B_{p,k}
#' \end{array} \right]}
#' Note: The determinant of the correlation matriz is
#' \deqn{\prod_{i=2}^p\prod_{j=1}^{i-1}B_{i,j} = \prod_{i=2}^pL_{i,i}}

#' For 'parametrization' = "SAP" or 'parametrization' = "sap":
#' The Standard Angles Parametrization - SAP, as described in
#' Rapisarda, Brigo and Mercurio (2007), compute
#' \eqn{x[k] = \pi/(1+\exp(-\theta[k]))}, for \eqn{k=1,...,m},
#' and the two \eqn{p\times p} matrices
#' \deqn{A = \left[
#' \begin{array}{ccccc}
#'   1 & & & & \\
#'   \cos(x_1) & 1 & & & \\
#'   \cos(x_2) & \cos(x_p) & 1 & & \\
#'   \vdots & \vdots & \ddots & \ddots & \\
#'   \cos(x_{p-1}) & \cos(x_{2p-3}) & \ldots & \cos(x_m) & 1
#' \end{array} \right] \textrm{ and } B = \left[
#' \begin{array}{ccccc}
#'   1 & & & & \\
#'   \sin(x_1) & 1 & & & \\
#'   \sin(x_2) & \sin(x_p) & 1 & & \\
#'   \vdots & \vdots & \ddots & \ddots & \\
#'   \sin(x_{p-1}) & \sin(x_{2p-3}) & \ldots & \sin(x_m) & 1
#' \end{array} \right]}
#'
#' The decomposition of the Hessian matrix around the base model,
#' `I0`, formally \eqn{\mathbf{I}(\theta_0)}, is numerically computed.
#' This element has the following attributes:
#' 'h.5' as \eqn{\mathbf{I}^{1/2}(\theta_0)}, and
#' 'hneg.5' as \eqn{\mathbf{I}^{-1/2}(\theta_0)}.
#' @references
#' Rapisarda, Brigo and Mercurio (2007).
#'   Parameterizing correlations: a geometric interpretation.
#'   IMA Journal of Management Mathematics (2007) 18, 55-73.
#'   <doi 10.1093/imaman/dpl010>
#'
#' Lewandowski, Kurowicka and Joe (2009)
#' Generating Random Correlation Matrices Based
#' on Vines and Extended Onion Method.
#' Journal of Multivariate Analysis 100: 1989–2001.
#' <doi: 10.1016/j.jmva.2009.04.008>
#' @export
basecor <- function(
    base,
    p,
    parametrization = "cpc",
    iparams,
    iLtheta) {
  UseMethod("basecor")
}
#' @describeIn basecor
#' Build a `basecor` from the parameter vector.
#' @returns a `basecor` object
#' @export
#' @example demo/basecor.R
basecor.numeric <- function(
    base,
    p,
    parametrization = "cpc",
    iparams,
    iLtheta) {

  parametrization <- match.arg(
    arg = tolower(parametrization),
    choices = c("cpc", "sap")
  )
  theta <- base

  ## check p and iLtheta
  iLtheta <- p_iLtheta_fncheck(p, iLtheta)
  if(missing(p)) {
    p <- attr(iLtheta, "p")
  }
  stopifnot(p>1)

  m <- length(iLtheta)
  ## check iparams
  iparams <- m_iparams_fncheck(m, iparams)

  L <- cholcor(
    theta = theta[iparams],
    p = p,
    parametrization = parametrization,
    iLtheta = iLtheta)
  base <- tcrossprod(L)

  out <- list(
    base = base,
    theta = theta,
    p = p,
    parametrization = parametrization,
    iLtheta = iLtheta,
    iparams = iparams,
    L = L)

  class(out) <- "basecor"
  return(out)
}
#' @describeIn basecor
#' Build a `basecor` from a correlation matrix.
#' @export
basecor.matrix <- function(
    base,
    p,
    parametrization = "cpc",
    iparams,
    iLtheta) {

  parametrization <- match.arg(
    arg = tolower(parametrization),
    choices = c("cpc", "sap")
  )

  stopifnot(all.equal(base, t(base)))
  L0 <- t(chol(base))

  p <- as.integer(nrow(base))
  if(missing(iLtheta)) {
      iLtheta <- which(lower.tri(diag(p)))
  }
  m <- length(iLtheta)

  ## check iparams
  iparams <- m_iparams_fncheck(m, iparams)
  m0 <- iparams[m]

  theta <- stats::optim(
    rep(0.0, m0), function(x)
      mean((tcrossprod(cholcor(
        theta = x[iparams],
        p = p,
        iLtheta = iLtheta,
        parametrization = parametrization))-base)^2),
      method = 'BFGS')$par

  out <- list(
    base = base,
    theta = theta,
    p = p,
    parametrization = parametrization,
    iparams = iparams,
    iLtheta = iLtheta,
    L = L0)
  class(out) <- "basecor"
  return(out)
}
#' @describeIn basecor
#' Print method for 'basecor'
#' @param x a basecor object.
#' @param ... further arguments passed on.
#' @export
print.basecor <- function(x, ...) {
  cat("Parameters (", toupper(x$parametrization),
      " parametrization):\n", sep = "")
  print(x$theta, ...)
  cat("Base correlation matrix:\n")
  print(x$base, ...)
}