R/intercorr_cont.R

In SimCorrMix: Simulation of Correlated Data with Multiple Variable Types Including Continuous and Count Mixture Distributions

#' @title Calculate Intermediate MVN Correlation for Continuous Variables Generated by Polynomial Transformation Method
#'
#' @description This function finds the intermediate correlation for standard normal random variables
#'     which are used in Fleishman's third-order (\doi{10.1007/BF02293811}) or Headrick's fifth-order
#'     (\doi{10.1016/S0167-9473(02)00072-5}) polynomial transformation method (PMT) using \code{\link[nleqslv]{nleqslv}}.  It is used in
#'     \code{\link[SimCorrMix]{intercorr}} and \code{\link[SimCorrMix]{intercorr2}} and would not ordinarily be called by the user.  The
#'     correlations are found pairwise so that eigen-value or principal components decomposition should be done on the resulting \code{Sigma}
#'     matrix. The \bold{Comparison of Correlation Methods 1 and 2} vignette contains the equations which were derived by Headrick and Sawilowsky
#'     (\doi{10.1007/BF02294317}) or Headrick (\doi{10.1016/S0167-9473(02)00072-5}).
#' @param method the method used to generate the continuous variables.  "Fleishman" uses Fleishman's third-order polynomial transformation
#'     and "Polynomial" uses Headrick's fifth-order transformation.
#' @param constants a matrix with each row a vector of constants c0, c1, c2, c3 (if \code{method} = "Fleishman") or
#'     c0, c1, c2, c3, c4, c5 (if \code{method} = "Polynomial"), like that returned by
#'     \code{\link[SimMultiCorrData]{find_constants}}
#' @param rho_cont a matrix of target correlations among continuous variables, does not have to be symmetric
#' @import nleqslv
#' @export
#' @keywords correlation, continuous, Fleishman, Headrick
#' @seealso \code{\link[SimCorrMix]{intercorr}}, \code{\link[SimCorrMix]{intercorr2}}, \code{\link[nleqslv]{nleqslv}}
#' @return the intermediate matrix of correlations with the same dimensions as \code{rho_cont}
#' @references Please see additional references for \code{\link[SimCorrMix]{SimCorrMix}}.
#'
#' Fialkowski AC (2018). SimMultiCorrData: Simulation of Correlated Data with Multiple Variable Types. R package version 0.2.2.
#'     \url{https://CRAN.R-project.org/package=SimMultiCorrData}.
#'
#' Headrick TC (2002). Fast Fifth-order Polynomial Transforms for Generating Univariate and Multivariate
#'     Non-normal Distributions. Computational Statistics & Data Analysis, 40(4):685-711. \doi{10.1016/S0167-9473(02)00072-5}.
#'     (\href{http://www.sciencedirect.com/science/article/pii/S0167947302000725}{ScienceDirect})
#'
#' Headrick TC, Kowalchuk RK (2007). The Power Method Transformation: Its Probability Density Function, Distribution
#'     Function, and Its Further Use for Fitting Data. Journal of Statistical Computation and Simulation, 77:229-249. \doi{10.1080/10629360600605065}.
#'
#' Headrick TC, Sawilowsky SS (1999). Simulating Correlated Non-normal Distributions: Extending the Fleishman Power
#'     Method. Psychometrika, 64:25-35. \doi{10.1007/BF02294317}.
#'
intercorr_cont <- function(method = c("Fleishman", "Polynomial"),
                           constants = NULL, rho_cont = NULL) {
  Sigma <- matrix(1, nrow(rho_cont), ncol(rho_cont))
  if (method == "Fleishman") {
    intercorrfleish <- function(r, c, rho) {
      f <-  numeric(1)
      a1 <- c[1, 1]
      b1 <- c[1, 2]
      d1 <- c[1, 4]
      a2 <- c[2, 1]
      b2 <- c[2, 2]
      d2 <- c[2, 4]
      f[1] <- r * (b1 * b2 + 3 * b2 * d1 + 3 * b1 * d2 + 9 * d1 * d2 +
                     2 * a1 * a2 * r + 6 * d1 * d2 * r^2) - rho
      return(f)
    }
    if (isSymmetric(rho_cont) & all(diag(rho_cont) == 1)) {
      for (i in 1:nrow(rho_cont)) {
        for (j in 1:ncol(rho_cont)) {
          if (j > i) {
            Sigma[i, j] <- nleqslv(x = rho_cont[i, j], intercorrfleish,
              c = constants[c(i, j), ], rho = rho_cont[i, j],
              method = "Broyden", control = list(ftol = 1e-5))$x
            Sigma[j, i] <- Sigma[i, j]
          }
        }
      }
    } else {
      for (i in 1:nrow(rho_cont)) {
        for (j in 1:ncol(rho_cont)) {
          Sigma[i, j] <- nleqslv(x = rho_cont[i, j], intercorrfleish,
            c = constants[c(i, j), ], rho = rho_cont[i, j],
            method = "Broyden", control = list(ftol = 1e-5))$x
        }
      }
    }
  }
  if (method == "Polynomial") {
    intercorrpoly <- function(r, c, rho) {
      f <- numeric(1)
      c0a <- c[1, 1]
      c1a <- c[1, 2]
      c2a <- c[1, 3]
      c3a <- c[1, 4]
      c4a <- c[1, 5]
      c5a <- c[1, 6]
      c0b <- c[2, 1]
      c1b <- c[2, 2]
      c2b <- c[2, 3]
      c3b <- c[2, 4]
      c4b <- c[2, 5]
      c5b <- c[2, 6]
      f[1] <- (3 * c4a * c0b + 3 * c4a * c2b + 9 * c4a * c4b + c0a *
                 (c0b + c2b + 3 * c4b) + c1a * c1b * r + 3 * c3a * c1b * r +
                 15 * c5a * c1b * r + 3 * c1a * c3b * r + 9 * c3a * c3b * r +
                 45 * c5a * c3b * r + 15 * c1a * c5b * r + 45 * c3a * c5b * r +
                 225 * c5a * c5b * r + 12 * c4a * c2b * r^2 +
                 72 * c4a * c4b * r^2 + 6 * c3a * c3b * r^3 +
                 60 * c5a * c3b * r^3 + 60 * c3a * c5b * r^3 +
                 600 * c5a * c5b * r^3 + 24 * c4a * c4b * r^4 +
                 120 * c5a * c5b * r^5 + c2a * (c0b + c2b + 3 * c4b +
                 2 * c2b * r^2 + 12 * c4b * r^2)) - rho
      return(f)
    }
    if (isSymmetric(rho_cont) & all(diag(rho_cont) == 1)) {
      for (i in 1:nrow(rho_cont)) {
        for (j in 1:ncol(rho_cont)) {
          if (j > i) {
            Sigma[i, j] <- nleqslv(x = rho_cont[i, j], intercorrpoly,
              c = constants[c(i, j), ], rho = rho_cont[i, j],
              method = "Broyden", control = list(ftol = 1e-5))$x
            Sigma[j, i] <- Sigma[i, j]
          }
        }
      }
    } else {
      for (i in 1:nrow(rho_cont)) {
        for (j in 1:ncol(rho_cont)) {
          Sigma[i, j] <- nleqslv(x = rho_cont[i, j], intercorrpoly,
            c = constants[c(i, j), ], rho = rho_cont[i, j],
            method = "Broyden", control = list(ftol = 1e-5))$x
        }
      }
    }
  }
  return(Sigma)
}