R/vcov.lmerMod.R

In merDeriv: Case-Wise and Cluster-Wise Derivatives for Mixed Effects Models

vcov.lmerMod <- function(object, ...) {
  
  if (!is(object, "lmerMod")) stop("vcov.lmerMod() only works for lmer() models.")
  
  dotdotdot <- list(...)
  if("full" %in% names(dotdotdot)){
    full <- dotdotdot$full
  } else {
    full <- FALSE
  }
  if ("information" %in% names(dotdotdot)) {
    information <- dotdotdot$information
  } else {
    information <- "expected"
  }
  if(!(full %in% c("TRUE", "FALSE"))) stop("invalid 'full' argument supplied")
  if(!(information %in% c("expected", "observed"))) stop("invalid 'information' argument supplied")

  if("ranpar" %in% names(dotdotdot)){
    ranpar <- dotdotdot$ranpar
  } else {
    ranpar <- "var"
  }  

  ## preparation for short cuts:
  ## get all elements by getME and exclude multiple random effect models.
  parts <- getME(object, "ALL")
  
  ## shortcut, y-X\beta
  yXbe <- parts$y - tcrossprod(parts$X, t(parts$beta))
  ## total variance
  uluti <- length(parts$theta)
  Zlam <- tcrossprod(parts$Z, parts$Lambdat)
  V <- (tcrossprod(Zlam, Zlam) + Diagonal(parts$n, 1)) * (parts$sigma)^2
  M <- solve(chol(V))
  invV <- tcrossprod(M, M)
  LambdaInd <- parts$Lambda
  LambdaInd@x[] <- parts$Lind 
  invVX <- crossprod(parts$X, invV)
  Pmid <- solve(crossprod(parts$X, t(invVX)))
  P <- invV - tcrossprod(crossprod(invVX, Pmid), t(invVX))
  
  ## block 1 fixhes
  fixvar <- solve(tcrossprod(crossprod(parts$X, invV), t(parts$X)))
  if (full == FALSE) {
    fixvar  
  } else {
    fixhes <- tcrossprod(crossprod(parts$X, invV), t(parts$X))
    
    ## length = (variance covariance parameter in G + residual variance.)
    uluti <- length(parts$theta)
    devV <- vector("list", (uluti + 1))
    
    ## get d(G)/d(sigma), faciliated by d(Lambda). 
    devLambda <- vector("list", uluti)
    score_varcov <- matrix(NA, nrow = length(parts$y), ncol = uluti)
    for (i in 1:uluti) {
      devLambda[[i]] <- forceSymmetric(LambdaInd==i, uplo = "L")
      devV[[i]] <- tcrossprod(tcrossprod(parts$Z, t(devLambda[[i]])), parts$Z)
    }
    devV[[(uluti+1)]] <- Diagonal(nrow(parts$X), 1)
    
    ## Block 4: sigma's second derivative.
    ranhes <- matrix(NA, nrow = (uluti + 1), ncol = (uluti + 1))
    ## combinations (allow repeat) to indicate entries
    entries <- rbind(matrix(rep(1: (uluti + 1), each = 2),
      (uluti + 1), 2, byrow = TRUE), t(combn((uluti + 1), 2)))
    entries <- entries[order(entries[,1], entries[,2]), ]
    ## ML estimates
    if (parts$devcomp$dims[['REML']] == 0) {
      if(information == "expected") {
        ranhes[lower.tri(ranhes, diag = TRUE)] <- 
          apply(entries, 1, function(x) as.numeric((1/2) *
                lav_matrix_trace(tcrossprod(tcrossprod(crossprod(invV,
                devV[[x[1]]]), invV), t(devV[[x[2]]])))))
      }
      if(information == "observed") {
        ranhes[lower.tri(ranhes, diag = TRUE)] <- unlist(apply(entries, 1, 
          function(x) as.vector(-as.numeric((1/2) *
          lav_matrix_trace(tcrossprod(tcrossprod(crossprod(invV,
          devV[[x[1]]]), invV), t(devV[[x[2]]])))) +
          tcrossprod((tcrossprod((crossprod(yXbe,
          tcrossprod(tcrossprod(crossprod(invV,
          devV[[x[1]]]), invV), t(devV[[x[2]]])))), invV)), t(yXbe)))))
      }      
    }
    ## REML estimates
    if (parts$devcomp$dims[['REML']] > 0) {
      if(information == "expected") {
        ## this is excessively slow and could be improved:
        ranhes[lower.tri(ranhes, diag = TRUE)] <- apply(entries, 1, 
          function(x) as.numeric((1/2) * lav_matrix_trace(tcrossprod(
          tcrossprod(crossprod(P, devV[[x[1]]]), P), t(devV[[x[2]]])))))
      }
      if(information == "observed") {
        ranhes[lower.tri(ranhes, diag = TRUE)] <- apply(entries, 1, 
          function(x) -as.numeric((1/2) * lav_matrix_trace(
          tcrossprod(tcrossprod(crossprod(P, devV[[x[1]]]), P), 
          t(devV[[x[2]]])))) + tcrossprod((tcrossprod((crossprod(
          yXbe, tcrossprod(tcrossprod(crossprod(invV,
          devV[[x[1]]]), P), t(devV[[x[2]]])))), invV)), t(yXbe)))             
      }      
    }
    ranhes <- forceSymmetric(ranhes, uplo = "L")
    
    ## block 2 and block 3: second derivative of sigma and beta.
    if (information == "expected"){
      varcov_beta <- matrix(0, length(devV), length(parts$beta))
    }
    if (information == "observed"){
      varcov_beta <- matrix(NA, length(devV), length(parts$beta))
      for (j in 1 : (length(devV))) {
        varcov_beta[j,] <- as.vector(tcrossprod(crossprod(parts$X, 
        (tcrossprod(crossprod(invV, devV[[j]]), invV))), t(yXbe)))
      }
    }
    ## Organize full_varcov
    ## reprameterize sd and var
    if (ranpar == "var") {
        ranhes <- ranhes
        varcov_beta <- varcov_beta
    } else if (ranpar == "sd") {
        ## varcov_beta reparameterization
        sdcormat <- as.data.frame(VarCorr(object,comp = "Std.Dev"),
          order = "lower.tri")
        sdcormat$sdcor2[which(is.na(sdcormat$var2))] <-
          sdcormat$sdcor[which(is.na(sdcormat$var2))]*2
        sdcormat$sdcor2[which(!is.na(sdcormat$var2))] <-
          sdcormat$vcov[which(!is.na(sdcormat$var2))]/
            sdcormat$sdcor[which(!is.na(sdcormat$var2))]
        varcov_beta <- sweep(varcov_beta, MARGIN = 1, sdcormat$sdcor2, `*`)
        ## ranhes reparameterization
        weight <- apply(entries, 1, function(x)
          sdcormat$sdcor2[x[1]] * sdcormat$sdcor2[x[2]])
        ranhes[lower.tri(ranhes, diag = TRUE)] <- weight *
          ranhes[lower.tri(ranhes, diag = TRUE)]
        ranhes <- forceSymmetric(ranhes, uplo = "L")
    } else {
        stop("ranpar needs to be var or sd for lmerMod object.")
    }

    full_varcov <- solve(rbind(cbind(fixhes, t(varcov_beta)),
                               cbind(varcov_beta, ranhes)))
    
    colnames(full_varcov) <- c(names(parts$fixef), paste("cov",
                             names(parts$theta), sep = "_"), "residual")
    
    callingFun <- try(deparse(sys.call(-2)), silent = TRUE)
    if(length(callingFun) > 1) callingFun <- paste(callingFun, collapse="")
    if(!inherits(callingFun, "try-error") &
       grepl("summary.merMod", callingFun)){
      return(fixvar)
    } else {
      return(full_varcov)
    }
  }
}