R/mLogLik.R

In LMMsolver: Linear Mixed Model Solver

# aux function
logLikelihood_aux <- function(y,
                          X,
                          Z,
                          lGinv,
                          lRinv,
                          theta)
{
  thetaMatrix <- theta
  Ntot <- length(y)
  p <- ncol(X)
  q <- ncol(Z)
  Nres <- length(lRinv)
  Nvarcomp <- length(lGinv)
  NvarcompTot <- Nres + Nvarcomp
  dimMME <- p + q
  W <- spam::cbind.spam(X, Z)
  lWtRinvW <- lapply(X = lRinv, FUN = function(x) {
    d <- spam::diag.spam(x)
    # if x is diagonal calculate x %*% W in a more efficient way
    if (isTRUE(all.equal(spam::diag.spam(d), x))) {
      xW <- W
      xW@entries <- xW@entries * rep(d, times = diff(xW@rowpointers))
    } else {
      xW <- x %*% W
    }
    MatrixProduct(t(W), xW)
  })
  lWtRinvY <- lapply(X = lRinv, FUN = function(x) {
    spam::crossprod.spam(W, x %*% y) })
  lYtRinvY <- lapply(X = lRinv, FUN= function(x) {
    quadForm(y, x) })

  ## Extend lGinv with extra zeros for fixed effect.
  lQ <- lapply(X = lGinv, FUN = function(x) {
    zero <- spam::spam(0, ncol = p, nrow = p)
    return(spam::bdiag.spam(zero, x))
  })
  lC <- c(lQ, lWtRinvW)
  ## Remove some extra zeros.
  lWtRinvW <- lapply(X = lWtRinvW, FUN = spam::cleanup)
  lQ <- lapply(X = lQ, FUN = spam::cleanup)
  lGinv <- lapply(X = lGinv, FUN = spam::cleanup)
  lRinv <- lapply(X = lRinv, FUN = spam::cleanup)
  lC <- lapply(X = lC, FUN = spam::cleanup)

  ## Check the stucture of Rinv, don't allow for overlapping penalties
  M <- sapply(lRinv, FUN = function(x) {
    abs(spam::diag.spam(x)) > getOption("spam.eps")})
  rSums <- rowSums(M)
  #if (!isTRUE(all.equal(rSums, rep(1, nrow(M))))) {
  if (max(rSums) > 1) {
    stop("Overlapping penalties for residual part.\n", call. = FALSE)
  }
  ## Calculate number of elements per group for residuals.
  nR <- colSums(M)
  EDmax_phi <- nR
  Rinv <- Reduce("+", lRinv)
  logdetRinvConstant <- as.numeric(spam::determinant.spam(Rinv)$modulus)

  ## Make ADchol for Ginv and C:
  if (Nvarcomp > 0) {
    ADcholGinv <- ADchol(lGinv)
  } else {
    ADcholGinv <- NULL
  }
  ADcholC <- ADchol(lC)

  ## Initialize values for loop.
  nRowTheta <- nrow(thetaMatrix)
  logL <- rep(NA, nRowTheta)
  deriv <- matrix(data=NA, nrow = nRowTheta, ncol = NvarcompTot)
  for (i in seq_len(nRowTheta)) {
    theta <- thetaMatrix[i, ]
    if (Nvarcomp > 0) {
      psi <- theta[seq_len(Nvarcomp)]
      phi <- theta[-(seq_len(Nvarcomp))]
    } else {
      psi <- NULL
      phi <- theta
    }
    ## calculate logdet for Rinv.
    logdetR <- -sum(nR * log(phi)) - logdetRinvConstant

    ## calculated logdet and dlogdet for Ginv and C.
    ## Ginv, if exists
    if (!is.null(ADcholGinv)) {
      dlogdetGinv <- dlogdet(ADcholGinv, psi)
      logdetG <- -attr(dlogdetGinv, which = "logdet")
    } else {
      logdetG <- 0
    }
    ## update the expressions including Rinv.
    YtRinvY <- sum(phi * unlist(lYtRinvY))
    WtRinvY <- as.vector(linearSum(theta = phi, matrixList = lWtRinvY))

    ## matrix C.
    dlogdetC <- dlogdet(ADcholC, theta, WtRinvY)
    logdetC <- attr(dlogdetC, which = "logdet")
    a <- attr(dlogdetC, which = "x.coef")
    if (all(!is.na(dlogdetC))) {
      if (!is.null(ADcholGinv)) {
        EDmax_psi <- psi * dlogdetGinv
      } else {
        EDmax_psi <- NULL
      }
      EDmax <- c(EDmax_psi, EDmax_phi)
      ED <- EDmax - theta * dlogdetC

      ## to make sure ED is always positive.
      ED <- pmax(ED, .Machine$double.eps)

      ## calculate Sum of Squares.
      SS_all <- calcSumSquares(lYtRinvY, lWtRinvY, lWtRinvW, lQ, a, Nvarcomp)

      ## Johnson and Thompson 1995, see below eq [A5]: Py = R^{-1} r.
      yPy <- YtRinvY - sum(a * WtRinvY)

      ## calculate REMLlogL, see e.g. Smith 1995.
      logL[i] <- -0.5 * (logdetR + logdetG + logdetC + yPy)
      deriv[i,] <- (1/theta)*ED - SS_all
    }
  }
  df <- cbind(thetaMatrix, deriv, logL)
  df <- as.data.frame((df))
  names(df) <- c(paste0("theta",1:NvarcompTot),
                 paste0("dL_dtheta",1:NvarcompTot),
                 "logL")
  return(df)
}

#' Function to obtain restricted log-likelihood and the first derivatives of the log-likelihood,
#' given values for the penalty parameters
#'
#' @param object an object of class LMMsolve
#' @param theta a matrix with values of precision parameters theta.

#' @returns A data.frame with logL and the first derivatives of log-likelihood
#'
#' @export
mLogLik <- function(object, theta) {
  len <- length(object$theta)
  if (!is.matrix(theta)) {
    stop("theta should be matrix")
  }
  if (length(object$theta) != ncol(theta)) {
    stop("theta has wrong number of columns")
  }
  X <- object$X
  Z <- object$Z
  lGinv <- object$lGinv
  lRinv <- object$lRinv
  y <- object$y
  return(logLikelihood_aux(y, X, Z, lGinv, lRinv, theta))
}