Nothing
### calc_d2Omega.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: sep 16 2021 (13:18)
## Version:
## Last-Updated: jul 26 2023 (14:27)
## By: Brice Ozenne
## Update #: 217
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
## * calc_d2Omega
##' @title Second Derivative of the Residual Variance-Covariance Matrix
##' @description Second derivative of the residual variance-covariance matrix for given parameter values.
##' @noRd
##'
##' @param structure [structure]
##' @param param [named numeric vector] values of the parameters.
##' @param Omega [list of matrices] Residual Variance-Covariance Matrix for each pattern.
##' @param dOmega [list of matrices] First derivative of the residual Variance-Covariance Matrix for each pattern.
##' @param Jacobian [matrix] Jacobian of the reparametrisation.
##' @param dJacobian [array] First derivative of the Jacobian of the reparametrisation.
##' @param transform.sigma,transform.k,transform.rho [character] Transformation used on the variance/correlation coefficients.
##' Only active if \code{"log"}, \code{"log"}, \code{"atanh"}: then the derivative is directly computed on the transformation scale instead of using the Jacobian.
##'
##' @keywords internal
##'
##' @examples
##' data(gastricbypassL, package = "LMMstar")
##' gastricbypassL$gender <- c("M","F")[as.numeric(gastricbypassL$id) %% 2+1]
##' dd <- gastricbypassL[!duplicated(gastricbypassL[,c("time","gender")]),]
##'
##' ## independence
##' Sid1 <- skeleton(IND(~1, var.time = "time"), data = dd)
##' Sid4 <- skeleton(IND(~1|id, var.time = "time"), data = dd)
##' Sdiag1 <- skeleton(IND(~visit), data = dd)
##' Sdiag4 <- skeleton(IND(~visit|id), data = dd)
##' Sdiag24 <- skeleton(IND(~visit+gender|id, var.time = "time"), data = dd)
##' param24 <- setNames(c(1,2,2,3,3,3,5,3),Sdiag24$param$name)
##'
##' .calc_d2Omega(Sid4, param = c(sigma = 2))
##' .calc_d2Omega(Sdiag1, param = c(sigma = 1, k.visit2 = 2, k.visit3 = 3, k.visit4 = 4))
##' .calc_d2Omega(Sdiag4, param = c(sigma = 1, k.visit2 = 2, k.visit3 = 3, k.visit4 = 4))
##' .calc_d2Omega(Sdiag24, param = param24)
##'
##' ## compound symmetry
##' Scs4 <- skeleton(CS(~1|id, var.time = "time"), data = gastricbypassL)
##' Scs24 <- skeleton(CS(gender~time|id), data = gastricbypassL)
##'
##' .calc_d2Omega(Scs4, param = c(sigma = 1,rho=0.5))
##' .calc_d2Omega(Scs4, param = c(sigma = 2,rho=0.5))
##' .calc_d2Omega(Scs24, param = c("sigma:F" = 2, "sigma:M" = 1,
##' "rho:F"=0.5, "rho:M"=0.25))
##'
##' ## unstructured
##' Sun4 <- skeleton(UN(~visit|id), data = gastricbypassL)
##' param4 <- setNames(c(1,1.1,1.2,1.3,0.5,0.45,0.55,0.7,0.1,0.2),Sun4$param$name)
##' Sun24 <- skeleton(UN(gender~visit|id), data = gastricbypassL)
##' param24 <- setNames(c(param4,param4*1.1),Sun24$param$name)
##'
##' .calc_d2Omega(Sun4, param = param4)
##' .calc_d2Omega(Sun24, param = param24)
`.calc_d2Omega` <-
function(object, param, Omega, dOmega, Jacobian, dJacobian,
transform.sigma, transform.k, transform.rho) UseMethod(".calc_d2Omega")
## * calc_d2Omega.ID
.calc_d2Omega.ID <- function(object, param, Omega, dOmega, Jacobian = NULL, dJacobian = NULL,
transform.sigma = NULL, transform.k = NULL, transform.rho = NULL){
## ** prepare
type <- stats::setNames(object$param$type,object$param$name)
name.sigma <- object$param$name[type=="sigma"]
name.k <- object$param$name[type=="k"]
name.rho <- object$param$name[type=="rho"]
name.paramVar <- c(name.sigma,name.k,name.rho)
param <- param[name.paramVar]
name.param <- names(param)
if(missing(Omega)){
Omega <- .calc_Omega(object, param = param, keep.interim = TRUE)
}
if(missing(dOmega)){
dOmega <- .calc_dOmega(object, param = param, Omega = Omega, Jacobian = Jacobian,
transform.sigma = transform.sigma, transform.k = transform.k, transform.rho = transform.rho)
}
Upattern <- object$Upattern
n.Upattern <- NROW(Upattern)
X.var <- object$var$Xpattern
X.cor <- object$cor$Xpattern
if(identical(transform.sigma,"log") && identical(transform.k,"log") && identical(transform.rho,"atanh")){
Jacobian <- NULL
dJacobian <- NULL
}else{
transform.sigma <- "none"
transform.k <- "none"
transform.rho <- "none"
}
if(!is.null(Jacobian)){
test.nooffdiag <- all(abs(c(Jacobian[lower.tri(Jacobian,diag=FALSE)],Jacobian[upper.tri(Jacobian,diag=FALSE)]))<1e-10)
if(test.nooffdiag){
JacobianM1 <- Jacobian
diag(JacobianM1) <- 1/diag(Jacobian)
## range(JacobianM1 - solve(Jacobian))
}else{
JacobianM1 <- solve(Jacobian)
}
}
## ** loop over covariance patterns
out <- lapply(1:n.Upattern, function(iPattern){ ## iPattern <- 1
iPattern.var <- Upattern[iPattern,"var"]
iPattern.cor <- Upattern[iPattern,"cor"]
iNtime <- Upattern[iPattern,"n.time"]
iName.param <- Upattern[iPattern,"param"][[1]]
if(is.null(iName.param)){return(NULL)}
iOmega.sd <- attr(Omega[[iPattern]],"sd")
iOmega.var <- tcrossprod(iOmega.sd)
iOmega.cor <- attr(Omega[[iPattern]],"cor")
iOmega <- Omega[[iPattern]] ; attr(iOmega,"sd") <- NULL; attr(iOmega,"cor") <- NULL; attr(iOmega,"time") <- NULL;
iScore <- stats::setNames(vector(mode = "list", length = length(iName.param)), iName.param)
iPair <- object$pair.vcov[[Upattern[iPattern,"name"]]]
n.iPair <- NCOL(iPair)
iHess <- lapply(1:n.iPair, function(iPair){matrix(0, nrow = iNtime, ncol = iNtime)})
for(iiPair in 1:n.iPair){ ## iiPair <- 2
## name of parameters
iCoef1 <- iPair[1,iiPair]
iCoef2 <- iPair[2,iiPair]
## type of parameters
iType1 <- type[iCoef1]
iType2 <- type[iCoef2]
## indicators
iMindicator.var1 <- attr(X.var[[iPattern.var]],"Mindicator.param")[[iCoef1]]
iMindicator.var2 <- attr(X.var[[iPattern.var]],"Mindicator.param")[[iCoef2]]
iIndicator.cor2 <- attr(X.cor[[iPattern.cor]],"indicator.param")[[iCoef2]]
if(iType1 == "sigma"){
if(iType2 == "sigma"){
if(iCoef1!=iCoef2){
stop("Cannot compute the Hessian with interacting sigma coefficients. \n")
}
if(transform.sigma == "log"){
iHess[[iiPair]] <- 4 * iOmega
}else{ ## no transformation (other transformations are made through jacobian)
iHess[[iiPair]] <- 2 * iOmega / param[iCoef1]^2
}
}else if(iType2 == "k"){
if(transform.k == "log"){
iHess[[iiPair]] <- iMindicator.var1 * iMindicator.var2 * iOmega
}else{ ## no transformation (other transformations are made through jacobian)
iHess[[iiPair]] <- iMindicator.var1 * iMindicator.var2 * iOmega / (param[iCoef1]*param[iCoef2])
}
}else if(iType2 == "rho"){
if(transform.rho == "atanh"){
iHess[[iiPair]][iIndicator.cor2] <- 2 * iOmega.var[iIndicator.cor2] * (1-param[iCoef2]^2)
}else{ ## no transformation (other transformations are made through jacobian)
iHess[[iiPair]][iIndicator.cor2] <- 2 * iOmega.var[iIndicator.cor2] / param[iCoef1]
}
}
}else if(iType1 == "k"){
if(iType2 == "k"){
if(transform.k == "log"){
iHess[[iiPair]] <- iMindicator.var1 * iMindicator.var2 * iOmega
}else{ ## no transformation (other transformations are made through jacobian)
if(iCoef1==iCoef2){
iHess[[iiPair]] <- iMindicator.var1*(iMindicator.var1-1) * iOmega / param[iCoef1]^2
}else{
iHess[[iiPair]] <- iMindicator.var1 * iMindicator.var2 * iOmega / (param[iCoef1]*param[iCoef2])
}
}
}else if(iType2 == "rho"){
if(transform.k == "log"){
iHess[[iiPair]][iIndicator.cor2] <- iMindicator.var1[iIndicator.cor2] * iOmega.var[iIndicator.cor2] * (1-param[iCoef2]^2)
}else{ ## no transformation (other transformations are made through jacobian)
iHess[[iiPair]][iIndicator.cor2] <- iMindicator.var1[iIndicator.cor2] * iOmega.var[iIndicator.cor2] / param[iCoef1]
}
}
}else if(transform.rho == "atanh" && iType1 == "rho" && iType2 == "rho" && iCoef1 == iCoef2){
iHess[[iiPair]][iIndicator.cor2] <- - 2 * iOmega.var[iIndicator.cor2] * param[iCoef2] * (1 - param[iCoef2]^2)
}
}
## apply transformation
if(!is.null(Jacobian) || !is.null(dJacobian)){
iParamVar <- Upattern[iPattern,"param"][[1]]
n.iParamVar <- length(iParamVar)
if(any(abs(JacobianM1[iParamVar,setdiff(name.paramVar,iParamVar),drop=FALSE])>1e-10) || any(abs(dJacobian[iParamVar,setdiff(name.paramVar,iParamVar),,drop=FALSE])>1e-10)){
stop("Something went wrong when computing the derivative of the residual variance covariance matrix. \n",
"Contact the package manager with a reproducible example generating this error message. \n")
}
M.iScore <- do.call(cbind,lapply(dOmega[[iPattern]],as.double)) %*% JacobianM1[iParamVar,iParamVar,drop=FALSE]
iHess2 <- vector(mode = "list", length = n.iPair)
for(iP in 1:n.iParamVar){ ## iP <- 2
## d/d theta_1 =
## [dOmega_[11]/d2 theta_1] ... [dOmega_[11]/d theta_1 d theta_p] %*% Jacobian + [dOmega_[11]/d theta_1] ... [dOmega_[11]/d theta_p] %*% dJacobian/d theta_1
## [dOmega_[ij]/d2 theta_1] ... [dOmega_[ij]/d theta_1 d theta_p] %*% Jacobian + [dOmega_[ij]/d theta_1] ... [dOmega_[ij]/d theta_p] %*% dJacobian/d theta_1
## [dOmega_[mm]/d2 theta_1] ... [dOmega_[mm]/d theta_1 d theta_p] %*% Jacobian + [dOmega_[mm]/d theta_1] ... [dOmega_[mm]/d theta_p] %*% dJacobian/d theta_1
iCoef1 <- iParamVar[iP]
## find all pairs involving the current parameter plus another parameter
## e.g. if the current parameter is sigma we want (sigma,sigma) (sigma,k2) (sigma,k3) (sigma k4) (sigma rho)
## e.g. if the current parameter is k4 we want (sigma,k4) (k2,k4) (k3,k4) (k4 k4) (k4 rho)
iIndex.pair <- which(colSums(iPair == iCoef1)>0)
## name of the coefficient of the other pair
iCoef.pairCoef <- apply(iPair[,iIndex.pair,drop=FALSE], 2, function(iCol){
iOut <- setdiff(iCol,iCoef1)
if(length(iOut)==0){return(iCoef1)}else{return(iOut)}
})
## reorder the pairs according to the order of the parameters
iIndex.pair <- iIndex.pair[match(iCoef.pairCoef, iParamVar)]
## apply transformation
M.iHess2 <- (do.call(cbind,lapply(iHess[iIndex.pair],as.double)) %*% Jacobian[iParamVar,iParamVar,drop=FALSE] + M.iScore %*% dJacobian[iParamVar,iParamVar,iCoef1]) * Jacobian[iCoef1,iCoef1]
## convert back to time format
iHess2[iIndex.pair] <- lapply(1:NCOL(M.iScore), function(iCol){matrix(M.iHess2[,iCol], nrow = iNtime, ncol = iNtime, byrow = FALSE)})
}
iHess <- iHess2
}
return(iHess)
})
## ** export
out <- stats::setNames(out,Upattern$name)
return(out)
}
## * calc_d2Omega.IND
.calc_d2Omega.IND <- .calc_d2Omega.ID
## * calc_d2Omega.CS
.calc_d2Omega.CS <- .calc_d2Omega.ID
## * calc_d2Omega.RE
.calc_d2Omega.RE <- .calc_d2Omega.ID
## * calc_d2Omega.TOEPLITZ
.calc_d2Omega.TOEPLITZ <- .calc_d2Omega.ID
## * calc_d2Omega.UN
.calc_d2Omega.UN <- .calc_d2Omega.ID
## * calc_d2Omega.CUSTOM
.calc_d2Omega.CUSTOM <- function(object, param, Omega, dOmega, Jacobian = NULL, dJacobian = NULL,
transform.sigma = NULL, transform.k = NULL, transform.rho = NULL){
Upattern <- object$Upattern
n.Upattern <- NROW(Upattern)
X.var <- object$var$Xpattern
X.cor <- object$cor$Xpattern
FCT.sigma <- object$FCT.sigma
FCT.rho <- object$FCT.rho
dFCT.sigma <- object$dFCT.sigma
dFCT.rho <- object$dFCT.rho
d2FCT.sigma <- object$d2FCT.sigma
d2FCT.rho <- object$d2FCT.rho
name.sigma <- object$param[object$param$type=="sigma","name"]
name.rho <- object$param[object$param$type=="rho","name"]
pair.varcoef <- object$pair.vcov
if(!is.null(FCT.sigma) && is.null(d2FCT.sigma) || !is.null(FCT.rho) && is.null(d2FCT.rho) ){
## second derivative
## unlist(.calc_dOmega.CUSTOM(object, param = param, Omega = Omega))
vec.dOmega <- numDeriv::jacobian(func = function(x){
unlist(.calc_dOmega.CUSTOM(object, param = x, Omega = Omega))
}, x = param[c(name.sigma,name.rho)])
## indicator of pattern
vec.pattern <- unlist(lapply(names(dOmega), function(iName){ ## iName <- names(dOmega)[1]
iParam <- names(dOmega[[iName]])
iNtime <- Upattern[Upattern$name==iName,"n.time"]
iOut <- lapply(iParam, function(iP){matrix(iName, nrow = iNtime, ncol = iNtime)})
return(iOut)
}))
## indicator of param
vec.param <- unlist(lapply(names(dOmega), function(iName){ ## iName <- names(dOmega)[1]
iParam <- names(dOmega[[iName]])
iNtime <- Upattern[Upattern$name==iName,"n.time"]
iOut <- lapply(iParam, function(iP){matrix(iP, nrow = iNtime, ncol = iNtime)})
return(iOut)
}))
## matrix to list of matrices according to param and pattern
vec.patternXparam <- paste0(vec.pattern,"_",vec.param)
list.d2Omega <- by(data = data.frame(pattern = vec.pattern, param = vec.param, vec.dOmega), INDICES = vec.patternXparam, FUN = function(idOmega){ ## idOmega <- vec.dOmega[1:16,]
iOut <- apply(idOmega[,-(1:2),drop=FALSE], MARGIN = 2, simplify = FALSE, function(iVec){
iNtime <- sqrt(length(iVec))
matrix(iVec, nrow = iNtime, ncol = iNtime)
})
names(iOut) <- c(name.sigma,name.rho)
return(iOut)
}, simplify = FALSE)[unique(vec.patternXparam)]
## normalize to expected output
out <- stats::setNames(vector(mode = "list", length = n.Upattern), Upattern$name)
for(iPattern in 1:n.Upattern){
iIndex.pattern <- vec.pattern[!duplicated(vec.patternXparam)]==Upattern$name[iPattern]
iParam.pattern <- vec.param[!duplicated(vec.patternXparam)][iIndex.pattern]
iList.d2Omega <- list.d2Omega[iIndex.pattern]
out[[iPattern]] <- apply(pair.varcoef[[Upattern$name[iPattern]]], MARGIN = 2, simplify = FALSE, function(iCol){ ## iCol <- pair.varcoef[[Upattern$name[iPattern]]][,1]
iList.d2Omega[[which(iParam.pattern==iCol[1])]][[iCol[2]]]
})
}
}else{
out <- stats::setNames(vector(mode = "list", length = n.Upattern), Upattern$name)
for(iPattern in 1:n.Upattern){ ## iPattern <- 1
iPattern.var <- Upattern$var[iPattern]
iNtime <- Upattern$n.time[iPattern]
iX.var <- X.var[[iPattern.var]]
iOmega.sd <- attr(Omega[[iPattern]], "sd")
idOmega.sd <- dFCT.sigma(p = param[name.sigma], n.time = iNtime, X = iX.var)
id2Omega.sd <- d2FCT.sigma(p = param[name.sigma], n.time = iNtime, X = iX.var)
if(iNtime > 1 && !is.na(Upattern$cor[iPattern])){
iPattern.cor <- Upattern$cor[iPattern]
iX.cor <- X.cor[[iPattern.cor]]
iOmega.cor <- attr(Omega[[iPattern]], "cor")
idOmega.cor <- dFCT.rho(p = param[name.rho], n.time = iNtime, X = iX.cor)
id2Omega.cor <- d2FCT.rho(p = param[name.rho], n.time = iNtime, X = iX.cor)
}
out[[iPattern]] <- apply(pair.varcoef[[Upattern$name[iPattern]]], MARGIN = 2, simplify = FALSE, function(iCol){ ## iCol <- pair.varcoef[[Upattern$name[iPattern]]][,1]
iDeriv <- matrix(0, iNtime, iNtime)
if(iCol[1] %in% name.sigma && iCol[2] %in% name.sigma){
if(iCol[1]==iCol[2]){
iDeriv1 <- idOmega.sd[[iCol[1]]]
iDeriv2 <- id2Omega.sd[[iCol[1]]]
## diagonal sigma terms: f(a)^2 --> 2f(a)f'(a) --> 2[f'(a)f'(a)+f(a)f''(a)]
iDeriv <- iDeriv + 2*diag(iDeriv1^2 + iOmega.sd*iDeriv2, nrow = iNtime, ncol = iNtime)
## off-diagonal sigma terms: f(a)f(b) --> f''(a)f(b)
if(iNtime > 1 && !is.null(X.cor)){
iDeriv <- iDeriv + iOmega.cor * (iDeriv2 %*% t(iOmega.sd) + iOmega.sd %*% t(iDeriv2))
}
}else if(iNtime > 1 && !is.null(X.cor)){
iDeriv1 <- idOmega.sd[[iCol[1]]]
iDeriv2 <- idOmega.sd[[iCol[2]]]
## off-diagonal sigma terms: f(a)f(b) --> f'(a)f'(b)
iDeriv <- iDeriv + iOmega.cor * (iDeriv1 %*% t(iDeriv2) + iDeriv2 %*% t(iDeriv1))
}
}else if(iCol[1] %in% name.rho && iCol[2] %in% name.rho){
## diagonal sigma terms: 0
## off-diagonal sigma terms: f(a)f(b)f(c) --> f(a)f(b)f''(c)
if(iCol[1]==iCol[2]){
iDeriv <- iDeriv + id2Omega.cor[[iCol[1]]] * tcrossprod(iOmega.sd)
}
}else{
iDeriv1 <- idOmega.sd[[iCol[iCol %in% name.sigma]]]
iDeriv2 <- idOmega.cor[[iCol[iCol %in% name.rho]]]
## diagonal sigma terms: 0
## off-diagonal sigma terms: f(a)f(b)f(c) --> f'(a)f(b)f'(c)
iDeriv <- iDeriv + iDeriv2 * (iDeriv1 %*% t(iOmega.sd) + iOmega.sd %*% t(iDeriv1))
}
return(iDeriv)
})
}
}
return(out)
}
##----------------------------------------------------------------------
### calc_d2Omega.R ends here
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.