Nothing
################################################################################
#' (Co)variance parameter transformations.
#'
#' Converts lists of (co)variance parameters either between \code{list} and
#' \code{vector} format or between the theta and nu scales.
#'
#' \itemize{
#' \item{\code{stTrans} }{Transform start parameters into lower triangle
#' matrices of class \code{dsCMatrix}.}
#' \item{\code{conTrans} }{Transformation of starting constraints to correct
#' format.}
#' \item{\code{start2theta} }{Converts lists of starting values for (co)variance
#' parameters to a theta object used to structure the (co)variance components
#' within gremlin.}
#' \item{\code{matlist2vech} }{Converts a \code{list} of (co)variance parameter
#' matrices to a vector with a \dQuote{skel} attribute.}
#' \item{\code{vech2matlist} }{Converts a vector of (co)variance parameters to
#' a list of covariance matrices.}
#' \item{\code{theta2nu_trans} }{Transforms theta to nu scale by taking the
#' Cholesky factor of each covariance matrix and then replacing the diagonals
#' with their (natural) logarithms. Done to ensure matrices are positive
#' definite.}
#' \item{\code{nu2theta_trans} }{Back transformation from
#' \code{theta2nu_trans}: exponentiates the diagonal elements of each matrix
#' then calculates the cross-product.}
#' \item{\code{theta2nu_lambda} }{Transformation that factors out a residual
#' variance so that \code{nu} contains the \sQuote{lambda} parameterization:
#' ratios of variance parameters with the residual variance.}
#' \item{\code{nu2theta_lambda} }{Back transformation from
#' \code{theta2nu_lambda}.}
#' \item{\code{nuVar2thetaVar_lambda} }{Transformation of Sampling Variances
#' from \code{lambda} Scale for \code{theta}.}
#' \item{\code{nuAI2thetaAIinv_lambda} }{Transform AI matrix from \code{lambda}
#' Scale to AI-inverse of \code{theta}.}
#' \item{\code{nu2theta_noTrans} }{Structures \code{theta} when not
#' transformed.}
#' }
#'
#' @aliases stTrans conTrans vech2matlist start2theta matlist2vech theta2nu_trans
#' nu2theta_trans theta2nu_lambda nu2theta_lambda nuVar2thetaVar_lambda
#' nuAI2thetaAIinv_lambda nu2theta_noTrans
#' @param x,theta,nu A \code{list} of matrices containing the (co)variance
#' parameters of the model.
#' @param object An object of \code{class} \sQuote{gremlin}.
#' @param vech A \code{vector} of (co)variance parameters.
#' @param skeleton An example structure to map \code{vech} onto.
#' @param Gstart,Rstart A \code{list} of starting (co)variance values for the
#' G-structure (random effects terms) or R-structure (residual).
#' @param Gcon,Rcon A \code{list} of starting (co)variance constraints for the
#' G-structure (random effects terms) or R-structure (residual).
#' @param name An (optional) character \code{vector} containing the (co)variance
#' component names.
#' @param thetaG,thetaR A \code{vector} indexing the G-structure or R-structure
#' components, respectively.
#' @param sigma2e A \code{numeric} estimate of the factored out residual
#' variance from the mixed model equations (i.e., the \sQuote{lambda} scale)
#' \eqn{\sigma^{2}_{e}}.
#' @return Functions are specified to mostly return either a \code{list} of
#' matrices (structure as defined by the \dQuote{skel} attribute or in
#' the \code{skeleton} object) or a \code{vector} containing the (co)variance
#' parameters of the model. Additional list elements returned can be:
#' \describe{
#' \item{thetaG }{A \code{vector} indexing the G-structure components.}
#' \item{thetaR }{A \code{vector} indexing the R-structure components.}
#' }
#' Alternatively, \code{nuVar2thetaVar_lambda} and \code{nuAI2thetaAIinv_lambda}
#' return a \code{vector} and \code{matrix}, respectively, holding the sampling
#' (co)variances of the model (co)variance parameters both on the \code{theta}
#' scale. These are elements of the inverse Average Information matrix.
#' @author \email{matthewwolak@@gmail.com}
#' @import Matrix
#' @examples
#' # User-specified starting parameters
#' thetaOut <- start2theta(Gstart = list(matrix(1), matrix(2)),
#' Rstart = matrix(3))
#' ## convert to a vector and then back into a matrix list
#' thetav <- matlist2vech(thetaOut$theta)
#' theta <- vech2matlist(thetav, attr(thetav, "skel"))
#' identical(thetaOut$theta, theta) #<-- should be TRUE
#' # lambda parameterization transformation
#' nu <- theta2nu_lambda(theta, thetaOut$thetaG, thetaOut$thetaR)
#' # back-transform from (lambda scale) nu to theta
#' ## For example, when the sigma2e estimate=0.5
#' theta2 <- nu2theta_lambda(nu, sigma2e = 0.5, thetaOut$thetaG, thetaOut$thetaR)
#' @name covFun
NULL
# Transformation of starting parameters to correct matrix format.
#' @rdname covFun
#' @export
stTrans <- function(x){
if(is.numeric(x) && !is.matrix(x)) x <- as.matrix(x)
if(!isSymmetric(x)) stop(cat("Element", x, "must be a symmetric matrix or a number\n"))
x <- as(x, "symmetricMatrix")
x@uplo <- "L"
x <- as(x, "dsCMatrix")
x
}
# Transformation of starting constraints to correct format.
#' @rdname covFun
#' @export
conTrans <- function(Gcon, Rcon){
l2symat <- function(x){
if(is.character(x) && !is.matrix(x)) x <- as.matrix(x)
if(!isSymmetric(x)){
stop(cat(x, "must be a symmetric matrix or a single character\n"))
}
x[lower.tri(x, diag = TRUE)]
}
c(G = sapply(Gcon, FUN = l2symat), R. = l2symat(Rcon))
}
# Starting Parameters to theta List.
#' @rdname covFun
#' @export
start2theta <- function(Gstart, Rstart, name = NULL){
theta <- c(G = sapply(Gstart, FUN = stTrans), R. = stTrans(Rstart))
thetaGorR <- sapply(strsplit(names(theta), ".", fixed = TRUE), FUN = "[[", i = 1)
#FIXME ensure grep() is best way and won't mess up when multiple Gs and **Rs**
thetaG <- grep("G", thetaGorR, ignore.case = FALSE)
thetaR <- grep("R", thetaGorR, ignore.case = FALSE)
thNames <- if(!is.null(name)) name else seq_len(length(thetaG))
names(theta) <- c(paste0(rep("G.", length(thetaG)), thNames),
paste0("ResVar", seq(length(thetaR))))
return(list(theta = theta, thetaG = thetaG, thetaR = thetaR))
}
# theta List to Vector
#' @rdname covFun
#' @export
matlist2vech <- function(theta){
thetav <- sapply(theta, FUN = slot, name = "x")
skel <- lapply(seq(length(theta)),
FUN = function(i){mapply(slot, theta[i], c("i", "p", "Dim"))})
names(skel) <- names(theta)
attr(thetav, "skel") <- skel
thetav
}
# Vector to List of Matrices.
#' @rdname covFun
#' @export
vech2matlist <- function(vech, skeleton){
newmatlist <- vector("list", length = length(skeleton))
si <- 1
for(s in 1:length(skeleton)){
lss <- length(skeleton[[s]][[1]])
newmatlist[[s]] <- sparseMatrix(i = skeleton[[s]][[1]],
p = skeleton[[s]][[2]],
x = vech[seq(from = si, by = 1, length.out = lss)],
dims = skeleton[[s]][[3]], symmetric = TRUE, index1 = FALSE)
si <- si + lss
}
if(!is.null(names(skeleton))) names(newmatlist) <- names(skeleton)
newmatlist
}
############################
# Log-diagonals of Cholesky
############################
#FIXME: change G to cholesky of G with log(diagonals)
### when to do cholesky factorization of G?
#### is it easier to do direct product of cholesky's than direct product then cholesky?
#### If the latter then save the symbolic factorization and just do updates?
#XXX Don't need for EM algorithm
#TODO FIXME need to take log of just diagonals
#TODO also convert to lambda scale
# Transform theta to nu With log-Cholesky Parameterization
#' @rdname covFun
#' @export
#' @import Matrix
theta2nu_trans <- function(theta){
list(G = lapply(theta$G, FUN = function(x){log(chol(x))}),
R = log(chol(theta$R)))
}
# Back-transform nu to theta, Reversing log-Cholesky Parameterization
#' @rdname covFun
#' @export
#' @import Matrix
#TODO FIXME when converting back, do exp() of just diagonals
#TODO also back-convert from lambda scale
nu2theta_trans <- function(nu){
list(G = lapply(nu[[1L]], FUN = function(x){exp(crossprod(x))}),
R = exp(crossprod(nu[[2L]])))
}
######################
# Lambda scale
######################
# Transform theta to nu on the lambda Scale
#' @rdname covFun
#' @export
#' @import Matrix
theta2nu_lambda <- function(theta, thetaG, thetaR){
nu <- theta
cRinv <- solve(chol(theta[[thetaR]]))
# Meyer 1991, p.77 (<-- ?p70/eqn4?)
nu[thetaG] <- lapply(thetaG,
FUN = function(x){as(as(crossprod(cRinv, nu[[x]]) %*% cRinv, "symmetricMatrix"),
"dsCMatrix")})
#TODO what to do if R is a matrix
nu[[thetaR]] <- as(as(crossprod(cRinv, nu[[thetaR]]) %*% cRinv, "symmetricMatrix"),
"dsCMatrix")
nu
}
# Back-transform lambda Scale nu to theta
#' @rdname covFun
#' @export
#' @import Matrix
nu2theta_lambda <- function(nu, sigma2e, thetaG, thetaR){ #TODO FIXME
cR <- chol(matrix(sigma2e))
theta <- nu
theta[thetaG] <- lapply(thetaG,
FUN = function(x){as(as(crossprod(cR, nu[[x]]) %*% cR, "symmetricMatrix"),
"dsCMatrix")})
theta[[thetaR]] <- as(as(crossprod(cR, nu[[thetaR]]) %*% cR, "symmetricMatrix"),
"dsCMatrix")
theta
}
#############################
#XXX Works on AI matrix XXX
# Back-transform Sampling Variances from lambda Scale for theta
#' @rdname covFun
#' @export
nuVar2thetaVar_lambda <- function(object){
if(!object$grMod$lambda){
stop("object must be of type lambda = TRUE")
}
s2e <- object$grMod$sigma2e
se <- sqrt(s2e)
nuv <- c(matlist2vech(object$grMod$nu))
s2e_ind <- length(nuv) #FIXME assumes always at the end (find which ==1.0)
nuv[s2e_ind] <- s2e
invAI <- solve(object$grMod$AI) #<-- TODO some check about invertibility?
s2e_var <- diag(invAI)[s2e_ind]
c(nuv[-s2e_ind]^2 * s2e_var +
2 * s2e * nuv[-s2e_ind] * invAI[s2e_ind, -s2e_ind] +
s2e^2 * diag(invAI)[-s2e_ind],
s2e_var)
}
#XXX Works on AI matrix XXX
# Transform AI matrix from lambda Scale to AI-inverse of theta
#' @rdname covFun
#' @export
#XXX See Lynch & Walsh (1998) p.818 for variance of a product (eqn. A1.18b,c)
#XXX for psi: cov(psi, theta_i) = nu_i * cov(s2e, psi) + s2e * cov(psi, nu_i)
nuAI2thetaAIinv_lambda <- function(object){
if(!object$grMod$lambda){
stop("object must be of type lambda = TRUE")
}
s2e <- object$grMod$sigma2e
thetav <- c(object$grMod$thetav)
nuv <- c(matlist2vech(object$grMod$nu))
s2e_ind <- length(nuv) #FIXME assumes always at the end (find which ==1.0)
invAIl <- solve(object$grMod$AI) #<-- TODO some check about invertibility?
invAI <- matrix(NA, nrow = length(thetav), ncol = length(thetav))
# cov(theta_i, s2e)
invAI[s2e_ind, -s2e_ind] <- invAI[-s2e_ind, s2e_ind] <- nuv[-s2e_ind] *
invAIl[s2e_ind, s2e_ind] + s2e * invAIl[s2e_ind, -s2e_ind]
# var(theta)
diag(invAI) <- nuVar2thetaVar_lambda(object)
# cov(theta_r, theta_c)
for(c in seq(length(nuv))[-s2e_ind]){
for(r in seq(c+1, length(nuv), 1)){
if(r == s2e_ind) next
invAI[r, c] <- invAI[c, r] <- nuv[r] * nuv[c] * invAIl[s2e_ind, s2e_ind] +
thetav[r] * invAIl[s2e_ind, c] +
thetav[c] * invAIl[s2e_ind, r] +
s2e^2 * invAIl[r, c]# - invAIl[s2e_ind, r] * invAIl[s2e_ind, c]
} #<-- end for r
} #<-- end for c
invAI
}
######################
# No transformations
######################
## Mostly ensures correct formatting of theta
## TODO combine with above function (reconcile with `vech2matlist()`)
#' @rdname covFun
#' @export
#' @import Matrix
# theta2nu_noTrans <- function(theta) theta
# Back-transform nu to theta When No Transformation Occurred
nu2theta_noTrans <- function(nu, thetaG, thetaR){
theta <- nu
theta[thetaG] <- lapply(thetaG,
FUN = function(x){as(as(matrix(nu[[x]],1,1), "symmetricMatrix"),
"dsCMatrix")})
theta[[thetaR]] <- as(as(matrix(nu[[thetaR]], 1, 1), "symmetricMatrix"),
"dsCMatrix")
theta
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.