R/working_correlations.R

In rdmatheus/mbcsec: Multivariate Box-Cox Symmetric Distributions Generated From Elliptical Copula

#' @name cormat.mbcsec
#'
#' @title Working Correlation Matrices for Multivariate Box-Cox Symmetric
#'     Generated by Elliptical Copula Distributions
#'
#' @param d Dimension of the correlation matrix.
#' @param p,q Order of the autoregressive and the moving average component,
#'     respectively, passed as arguments to \code{arma()}.
#' @param id 	Subject id for longitudinal/clustered data. This is a vector of
#'     the same lenght of the number of observations. Please note that data
#'     must be sorted in way that observations from the same cluster are
#'     contiguous.
#' @param type  A character string specifying the correlation structure among
#'     groups for longitudinal/clustered data. At the
#'     moment, the following are implemented:
#'     \tabular{ll}{
#'     \code{independence}  \tab Working independence. \cr
#'     \code{ar1}  \tab Autoregressive of order 1. \cr
#'     \code{ma1}  \tab Moving average of order 1. \cr
#'     \code{exchangeable1}  \tab Exchangeable. \cr
#'     \code{unstructured}  \tab Unstructured \cr
#'     }
#' @param D Matrix with values of the distances between pairs of data
#'     locations for spatial data.
#' @param smt Value of the shape parameter of the Matern correlation class.
#'     The default \code{smt = 0.5} corresponds to an exponential correlation
#'     model.
#'
#' @description Current available correlation matrices in the \code{mbcsec}
#'     package.
#'
#' @details The functions related to the association structures are a direct
#'     adaptation of the functions of the \code{\link[gcmr]{gcmr}} package.
#'     The documentation of the original functions of the package can be seen
#'     at \code{\link[gcmr]{cormat.gcmr}} documentation. The available
#'     correlation matrices are:
#'
#'  \tabular{ll}{
#'  \bold{Function}  \tab \bold{Correlation}\cr
#'  \code{ind}  \tab Working independence. \cr
#'  \code{un}  \tab Unstructured \cr
#'  \code{arma}  \tab ARMA(p, q) \cr
#'  \code{cluster}  \tab Longitudinal/clustered data \cr
#'  \code{matern}  \tab Matern spatial correlation \cr
#'  }
#'
#'  Each type of structure will require different arguments for the association
#'  matrix P to be completed. However, all functions return a list of the
#'  following components:
#'
#'  \itemize{
#'  \item{\code{npar}:}{ Number of parameters associated with the correlation
#'      structure.}
#'  \item{\code{start}:}{ A function \code{function(y)} that returns the
#'      initial value for use in optimization algorithms. Its argument is the
#'      vector/matrix of observations \code{y}.}
#'  \item{\code{P}:}{ A function \code{function(alpha, d)} which returns the
#'      association matrix itself. Its arguments are the parameter vector
#'      associated with the correlation matrix structure and the dimension
#'      of the correlation matrix.}
#'  \item{\code{name}:}{ Name of the specified correlation structure.}
#'  }
#'
#' @references Guido Masarotto, & Cristiano Varin (2017). Gaussian Copula
#'     Regression in R. Journal of Statistical Software, 77(8), 1--26.
#'
#' @author Rodrigo M. R. Medeiros <\email{rodrigo.matheus@live.com}>
#'
#' @export


### Working independence correlation -------------------------------------------
#' @rdname cormat.mbcsec
#' @export
ind <- function(){
  ans <- list()
  ans$npar <- 0
  ans$start <- function(y) NULL
  ans$P <- function(alpha = NULL, d) diag(d)
  ans$name <- "Working independence"
  ans
}

### Unstructured working correlation -------------------------------------------
#' @rdname cormat.mbcsec
#' @export
un <- function(d){
  ans <- list()
  ans$npar <- d * (d - 1) / 2
  ans$start <- function(y) unique(as.numeric(stats::cor(y)))[-1]
  ans$P <- function(alpha, d = d){
    class(alpha) <- 'dist'
    attr(alpha,'Size') <- d
    alpha <- as.matrix(alpha)
    diag(alpha) <- 1
    alpha
  }
  ans$name <- "Unstructured"
  ans
}

### ARMA(p,q) working correlation for time-series ------------------------------
#' @rdname cormat.mbcsec
#' @export
arma <- function(p = 0, q = 0) {

  if (p == 0 && q == 0)
    return(ind ())

  ar <- if (p) 1:p else NULL
  ma <- if (q) (p + 1):(p + q) else NULL

  ans <- list()
  ans$npar <- p + q
  ans$start <- function(y = NULL) {
    alpha <- rep(0, p + q)
    names(alpha) <- c(if ( p ) paste("ar", 1:p, sep="") else NULL ,
                      if ( q ) paste("ma", 1:q, sep="") else NULL )
    alpha
  }
  ans$P <- function(alpha, d) {
    #if ( ( p && any(Mod(polyroot(c(1, -alpha[ar]))) < 0.99) ) ||
    #     ( q && any(Mod(polyroot(c(1, alpha[ma]))) < 0.99) ) )
    #  return( NULL )
    rho <- stats::ARMAacf(alpha[ar], alpha[ma], d - 1)
    r <- seq(1, d)
    outer(r, r, function(i, j) rho[1 + abs(i - j)] )
  }
  ans$name <- paste0("ARMA(", p, ", ", q,")")
  ans
}

### Longitudinal/Clustered data working correlation ----------------------------
# It assumes that it is not possible that all the observations inside a cluster
# can be missing
#' @rdname cormat.mbcsec
#' @export
cluster <- function(id, type = c("independence", "ar1", "ma1",
                                 "exchangeable", "unstructured")) {

  type <- match.arg(type, c("independence", "ar1", "ma1",
                            "exchangeable", "unstructured"))

  if(!length(rle(id)$values) == length(unique(id)))
    stop("data must be sorted in way that observations from the same cluster are contiguous")

  ng <- 1:length(unique(id))
  if (!(length(ng) > 1)) stop("only one strata")

  if (type == "independence") {
    ans <- ind()
    ans$id <- id
    return(ans)
  }

  ans <- list(type = type, id = id)
  ans$npar <- if (type != "unstructured") 1 else choose(max(table(id)), 2)
  data <- data.frame(id = id)
  fn <- switch(type,
               "ar1" = function(g) nlme::corAR1(g, form = ~ 1 | id),
               "ma1" = function(g) nlme::corARMA(g, form = ~ 1 |id, p = 0, q = 1),
               "exchangeable" = function(g) nlme::corCompSymm(g, form = ~ 1 | id),
               "unstructured" = function(g) nlme::corSymm(g, form = ~ 1 | id))
  ans$start <- function(y = NULL) {
    np <-  if(type != "unstructured") 1 else choose(max(table(id)), 2)
    alpha <- rep(0, np)
    names(alpha) <- switch(type, "ar1" = "ar1", "ma1" = "ma1",
                           "exchangeable" = "alpha",
                           "unstructured" = paste("alpha", 1:ans$npar, sep = ""))
    eps <- sqrt(.Machine$double.eps)
    attr(alpha, "lower") <- rep(-1 + eps, np)
    attr(alpha, "upper") <- rep(1 - eps, np)
    alpha
  }
  ans$P <- function(alpha, d) {
    q <- try(nlme::corMatrix(nlme::Initialize(fn(alpha), data = data)),
             silent = TRUE)
    if (inherits(q, "try-error")) return(NULL)
    g <- split(rep(TRUE, d), id)
    q <- try(lapply(ng, function(i) q[[i]][g[[i]],g[[i]]]), silent = TRUE)
    if (inherits(q, "try-error") ) NULL else as.matrix(Matrix::bdiag(q))
  }
  ans$name <- "Clustered"
  ans
}

### Matern working correlation for spatial data --------------------------------
## D is a distance matrix, smt is the smoothing parameter
#' @rdname cormat.mbcsec
#' @export
matern <- function(D, smt = 0.5) {
  ans <- list()
  ans$npar <- 1
  ans$start <- function(y = NULL) {
    alpha <- stats::median(D)
    names(alpha) <- c("alpha")
    attr(alpha,"lower") <- sqrt(.Machine$double.eps)
    alpha
  }
  ans$P <- function(alpha, d){
    #S <- .matern(D, alpha, smt)
    #q <- try(S[d, d],silent=TRUE)
    #if( inherits(q,"try-error") ) NULL else q
    .matern(D, alpha, smt)
  }
  ans$name <- "Matern spatial"
  ans
}

.matern <- function (u, phi, kappa)
{
  if (is.vector(u))
    names(u) <- NULL
  if (is.matrix(u))
    dimnames(u) <- list(NULL, NULL)
  uphi <- u/phi
  uphi <- ifelse(u > 0, (((2^(-(kappa - 1)))/ifelse(0, Inf,
                 gamma(kappa))) * (uphi^kappa) * besselK(x = uphi, nu = kappa)), 1)
  uphi[u > 600 * phi] <- 0
  return(uphi)
}