kergp (version 0.5.7)

corLevSymm: Correlation Matrix for a General Symmetric Correlation Structure

Description

Compute the correlation matrix for a general symmetric correlation structure.

Usage

corLevSymm(par, nlevels, levels, lowerSQRT = FALSE, compGrad = TRUE,
           cov = 0, impl = c("C", "R"))

Value

A correlation matrix (or its Cholesky root) with the optional

gradient attribute.

Arguments

par

A numeric vector with length npCor + npVar where npCor = nlevels * (nlevels - 1) / 2 is the number of correlation parameters, and npVar is the number of variance parameters, which depends on the value of cov. The value of npVar is 0, 1 or nlevels corresponding to the values of cov: 0, 1 and 2. The correlation parameters are assumed to be located at the head of par i.e. at indices 1 to npCor. The variance parameter(s) are assumed to be at the tail, i.e. at indices npCor + 1 to npCor + npVar.

nlevels

Number of levels.

levels

Character representing the levels.

lowerSQRT

Logical. When TRUE the (lower) Cholesky root \(\mathbf{L}\) of the correlation or covariance matrix \(\mathbf{C}\) is returned instead of the correlation matrix.

compGrad

Logical. Should the gradient be computed? This is only possible for the C implementation.

cov

Integer 0, 1 or 2. If cov is 0, the matrix is a correlation matrix (or its Cholesky root). If cov is 1 or 2, the matrix is a covariance (or its Cholesky root) with constant variance vector for code = 1 and with arbitrary variance for code = 2. The variance parameters par are located at the tail of the par vector, so at locations npCor + 1 to npCor + nlevels when code = 2 where npCor is the number of correlation parameters, i.e. nlevels * (nlevels - 1) / 2.

impl

A character telling which of the C and R implementations should be chosen.

Details

The correlation matrix with dimension \(n\) is the general symmetric correlation matrix as described by Pinheiro and Bates and implemented in the nlme package. It depends on \(n \times (n - 1) / 2\) parameters \(\theta_{ij}\) where the indices \(i\) and \(j\) are such that \(1 \leq j < i \leq n\). The parameters \(\theta_{ij}\) are angles and are to be taken to be in \([0, \pi)\) for a one-to-one parameterisation.

References

Jose C. Pinheiro and Douglas M. Bates (1996). "Unconstrained Parameterizations for Variance-Covariance matrices". Statistics and Computing, 6(3) pp. 289-296.

Jose C. Pinheiro and Douglas M. Bates (2000) Mixed-Effects Models in S and S-PLUS, Springer.

See Also

The corSymm correlation structure in the nlme package.

Examples

checkGrad <- TRUE
nlevels <- 12
npar <- nlevels * (nlevels - 1) / 2
par <- runif(npar, min = 0, max = pi)
##============================================================================
## Compare R and C implementations for 'lowerSQRT = TRUE'
##============================================================================
tR <- system.time(TR <- corLevSymm(nlevels = nlevels,
                                   par = par, lowerSQRT = TRUE, impl = "R"))
tC <- system.time(T <- corLevSymm(nlevels = nlevels, par = par,
                                  lowerSQRT = TRUE))
tC2 <- system.time(T2 <- corLevSymm(nlevels = nlevels, par = par,
                                    lowerSQRT = TRUE, compGrad = FALSE))
## time
rbind(R = tR, C = tC, C2 = tC2)

## results
max(abs(T - TR))
max(abs(T2 - TR))

##============================================================================
## Compare R and C implementations for 'lowerSQRT = FALSE'
##============================================================================
tR <- system.time(TRF <- corLevSymm(nlevels = nlevels, par = par,
                                    lowerSQRT = FALSE, impl = "R"))
tC <- system.time(TCF <- corLevSymm(nlevels = nlevels, par = par,
                                    compGrad = FALSE, lowerSQRT = FALSE))
tC2 <- system.time(TCF2 <- corLevSymm(nlevels = nlevels, par = par,
                                      compGrad = TRUE, lowerSQRT = FALSE))
rbind(R = tR, C = tC, C2 = tC2)
max(abs(TCF - TRF))
max(abs(TCF2 - TRF))

##===========================================================================
## Compare the gradients
##===========================================================================

if (checkGrad) {

    library(numDeriv)

    ##==================
    ## lower SQRT case
    ##==================
    JR <- jacobian(fun = corLevSymm, x = par, nlevels = nlevels,
                   lowerSQRT = TRUE, method = "complex", impl = "R")
    J <- attr(T, "gradient")

    ## redim and compare.
    dim(JR) <- dim(J)
    max(abs(J - JR))
    nG <- length(JR)
    plot(1:nG, as.vector(JR), type = "p", pch = 21, col = "SpringGreen3",
         cex = 0.8, ylim = range(J, JR),
         main = "gradient check, lowerSQRT = TRUE")
    points(x = 1:nG, y = as.vector(J), pch = 16, cex = 0.6, col = "orangered")

    ##==================
    ## Symmetric case
    ##==================
    JR <- jacobian(fun = corLevSymm, x = par, nlevels = nlevels,
                   lowerSQRT = FALSE, impl = "R", method = "complex")
    J <- attr(TCF2, "gradient")

    ## redim and compare.
    dim(JR) <- dim(J)
    max(abs(J - JR))
    nG <- length(JR)
    plot(1:nG, as.vector(JR), type = "p", pch = 21, col = "SpringGreen3",
         cex = 0.8,
         ylim = range(J, JR),
         main = "gradient check, lowerSQRT = FALSE")
    points(x = 1:nG, y = as.vector(J), pch = 16, cex = 0.6, col = "orangered")
}