corLevSymm | R Documentation |
Compute the correlation matrix for a general symmetric correlation structure.
corLevSymm(par, nlevels, levels, lowerSQRT = FALSE, compGrad = TRUE,
cov = 0, impl = c("C", "R"))
par |
A numeric vector with length |
nlevels |
Number of levels. |
levels |
Character representing the levels. |
lowerSQRT |
Logical. When |
compGrad |
Logical. Should the gradient be computed? This is only possible for the C implementation. |
cov |
Integer |
impl |
A character telling which of the C and R implementations should be chosen. |
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.
A correlation matrix (or its Cholesky root) with the optional
gradient
attribute.
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.
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.
The corSymm
correlation structure in the nlme
package.
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")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.