R/basepcor.R

In graphpcor: Models for Correlation Matrices Based on Graphs

#' A base precision's Cholesky model for a correlation matrix
#' @description
#' The `basepcor` class contain a correlation matrix `base`,
#' the parameter vector `theta`, that generates
#' or is generated by `base`, the dimension `p`,
#' the index `iLtheta` for `theta` in the (lower) Cholesky,
#' and the Hessian around it `I0`, see details.
#' @describeIn basepcor
#' Build a `basepcor` object.
#' @param base a correlation matrix, or numeric vector as
#' the parameter(s) to define correlation matrix.
#' @param p integer (needed if `base` is vector): the dimension.
#' @param iLtheta integer vector or 'graphpcor' to specify the (vectorized)
#' position where 'theta' is placed in the initial (before the fill-in)
#' Cholesky (lower triangle) factor. If missing, default, assumes
#' the dense case as `iLtheta = which(lower.tri(...))`, giving
#' `length(theta)=p(p-1)/2`.
#' @param d0 numeric vector to specify the diagonal of the
#' Cholesky factor for the initial precision matrix `Q0`.
#' Default, if not provided, is `d0 = p:1`.
#' @param iparams integer ordered vector with length equal
#' the number of parameters used to specify common parameter values.
#' Default is `1:m`, `m=length(theta)`. Example: By setting
#' `iparams = c(1,1,2,3)`, `m=3`, the first and second parameters
#' are considered to be the same.
#' NOTE: `c(1,2,1)` is allowed, but `c(2,1,2)` is not.
#' @details
#' The Inverse Transform Parametrization - ITP,
#' is applied by starting with a
#' \deqn{\mathbf{L}^{(0)} = \left[ \begin{array}{ccccc}
#' \textrm{p} & 0 & 0 & \ldots & 0 \\
#' \theta_1 & \textrm{p-}1 & 0 & \ldots & 0 \\
#' \theta_2 & \theta_p & \textrm{p-}2 & \ddots & \vdots \\
#' \vdots & \vdots & \ddots & \ddots & 0 \\
#' \theta_{p-1} & \theta_{2p-3} & \ldots & \theta_m & 1
#' \end{array}\right] .}
#'
#' Then compute \eqn{\mathbf{Q}^{(0)}} =
#' \eqn{\mathbf{L}\mathbf{L}^{T}},
#' \eqn{\mathbf{V}^{(0)}} = \eqn{(\mathbf{Q}^{(0)})^{-1}} and
#' \eqn{s_{i}^{(0)}} = \eqn{\sqrt{\mathbf{V}_{i,i}^{(0)}}}, and
#' define \eqn{\mathbf{S}^{(0)}} = diag\eqn{(s_1^{(0)},\ldots,s_p^{(0)})}
#'  in order to have \eqn{\mathbf{C}}=
#'  \eqn{\mathbf{S}^{-1}\mathbf{V}^{(0)}\mathbf{S}^{-1}}.
#'
#' 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)}.
#' @returns a basepcor object
#' @export
basepcor <- function(
    base,
    p,
    iLtheta,
    d0,
    iparams) {
  UseMethod("basepcor")
}
#' @describeIn basepcor
#' Build a `basepcor` from the parameter vector.
#' @returns a `basepcor` object
#' @export
#' @example demo/basepcor.R
basepcor.numeric <- function(
    base,
    p,
    iLtheta,
    d0,
    iparams) {

  theta <- base

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

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

  if(missing(d0)) {
    d0 <- p:1
  } else {
    stopifnot(length(d0)==p)
    stopifnot(all(d0>0))
  }

  ## compute L0
  L0 <- Lprec0(
    theta = theta[iparams],
    p = p,
    iLtheta = iLtheta,
    d0 = d0
  )
  lq02cor <- function(l) {
    cov2cor(chol2inv(t(l)))
  }

  ## compute the base correlation matrix
  base <- lq02cor(L0)
  U0correl <- chol(base)

  ## output
  out <- list(
    base = base,
    theta = theta,
    p = p,
    d0 = d0,
    iLtheta = iLtheta,
    iparams = iparams,
    L0 = L0, ## initial precision (lower) Cholesky
    L = t(U0correl)) ## the correlation's (lower) Cholesky
  class(out) <- "basepcor"

  return(out)
}

#' @describeIn basepcor
#' Build a `basepcor` from a correlation matrix.
#' @export
basepcor.matrix <- function(
    base,
    p,
    iLtheta,
    d0,
    iparams) {
  stopifnot(all.equal(base, t(base)))
  p <- as.integer(nrow(base))
  if(missing(d0)) {
    d0 <- p:1
  } else {
    stopifnot(length(d0)==p)
    stopifnot(all(d0>0))
  }
  stopifnot((length(d0)==p) && (all(d0>0)))
  U0correl <- chol(base)
  Q <- chol2inv(U0correl)
  ilQ <-  which(
    lower.tri(matrix(1, p, p)) &
      (!is.zero(Q, tol = 0.001)))
  iLtheta <- p_iLtheta_fncheck(p, iLtheta)
  stopifnot(all(ilQ %in% iLtheta))

  ## compute theta
  LQ0 <- t(chol(Q))
  for(i in 1:p) {
      LQ0[i, ] <- (d0[i]/LQ0[i, i]) * LQ0[i, ]
  }
  theta <- LQ0[iLtheta]

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

  if(iparams[m]<m) {
    ## Check if the parameters assumed to be common actually are
    stheta <- split(theta, iparams)
    theta.diff <- all(sapply(stheta, function(x)
      all(abs(x-mean(x)))>0.001))
    if(any(theta.diff)) {
      cat("Differences found among supposed equal parameters:\n")
      print(stheta[which(theta.diff)])
      stop("Please review `iparams` definition!")
    }
    theta0 <- sapply(stheta, mean)
  } else {
    theta0 <- tapply(theta, iparams, mean)
  }

  ## output
  out <- list(
    base = base,
    theta = new("numeric", theta0),
    p = p,
    d0 = d0,
    iLtheta = iLtheta,
    iparams = iparams,
    L0 = LQ0, ## initial precision (lower) Cholesky
    L = t(U0correl)) ## the correlation's (lower) Cholesky

  class(out) <- "basepcor"

  return(out)
}
#' @describeIn basepcor
#' Print method for 'basepcor'
#' @param x a basepcor object.
#' @param ... further arguments passed on.
#' @export
print.basepcor <- function(x, ...) {
  cat("Parameters:\n", sep = "")
  print(x$theta, ...)
  cat("Base correlation matrix:\n")
  print(x$base, ...)
}