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\le j<i\le 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")
}