# Correlation Matrix for a General Symmetric Correlation Structure Compute the correlation matrix for a general symmetric correlation structure. ```r 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 $L$ of the correlation or covariance matrix $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. ## Returns A correlation matrix (or its Cholesky root) with the optional `gradient` attribute. ## 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 * (n - 1) / 2$ parameters $theta[i, j]$ where the indices $i$ and $j$ are such that $1 <= j < i <= n$. The parameters $theta[i, j]$ are angles and are to be taken to be in $[0, pi)$ for a one-to-one parameterisation. ## 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[i, j]$ 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. ## Examples ```r 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") } ``` ## 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.

corLevSymm function