corLevSymm {kergp}

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"))

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.

Value

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

Note

This function is essentially for internal use and the corresponding correlation or covariance kernels are created as covQual objects by using the q1Symm creator.

The parameters \theta_{ij} are used in row order rather than in the column order as in the reference or in the nlme package. This order simplifies the computation of the gradients.

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")
}