R/estimatePcev.R

In GreenwoodLab/pcev: Principal Component of Explained Variance

#' Estimation of PCEV
#' 
#' \code{estimatePcev} estimates the PCEV.
#' 
#' @seealso \code{\link{computePCEV}}
#' @param pcevObj A pcev object of class \code{PcevClassical}, \code{PcevBlock} or
#'   \code{PcevSingular}
#' @param shrink Should we use a shrinkage estimate of the residual variance? 
#' @param index If \code{pcevObj} is of class \code{PcevBlock}, \code{index} is a vector
#'   describing the block to which individual response variables correspond.
#' @param ... Extra parameters.
#' @return A list containing the variance components, the first PCEV, the 
#'   eigenvalues of \eqn{V_R^{-1}V_M} and the estimate of the shrinkage 
#'   parameter \eqn{\rho}
#' @export 
estimatePcev <- function(pcevObj, ...) UseMethod("estimatePcev")

#' @rdname  estimatePcev
estimatePcev.default <- function(pcevObj, ...) {
  stop(strwrap("This function should be used with a Pcev object of class 
               PcevClassical, PcevBlock or PcevSingular"))
}

#' @rdname estimatePcev
#' @importFrom stats cov na.omit
estimatePcev.PcevClassical <- function(pcevObj, shrink, index, ...) {
  #initializing parameters
  rho <- NULL
  Y <- pcevObj$Y
  N <- nrow(Y)
  p <- ncol(Y)

  # Variance decomposition
  # fit <- lm.fit(cbind(pcevObj$X, pcevObj$Z), Y)
  Yfit <- getFitted(cbind(pcevObj$X, pcevObj$Z), Y,
                    overall = pcevObj$overall)
  res <- Y - Yfit
  # fit_confounder <- lm.fit(cbind(rep_len(1, N), pcevObj$Z), Y)
  Yfit_confounder <- getFitted(cbind(rep_len(1, N), pcevObj$Z), Y,
                               overall = pcevObj$overall)
  
  Vr <- cov(res, use = "pairwise.complete.obs")
  Vm <- cov(Yfit - Yfit_confounder, Y,
            use = "pairwise.complete.obs")
  
  # Shrinkage estimate of Vr
  if (shrink) {
    out <- shrink_est(Vr, res)
    Vrs <- out$cov
    rho <- out$rho
    
    # Computing PCEV
    temp <- eigen(Vr, symmetric = TRUE)
    Ur <- temp$vectors
    diagD <- eigen(Vrs, symmetric = TRUE, only.values = TRUE)$values
  } else {
    # Computing PCEV
    temp <- eigen(Vr, symmetric = TRUE)
    Ur <- temp$vectors
    diagD <- temp$values
  }
  
  value <- 1/sqrt(diagD)
  if (any(is.nan(value))) {
    stop("The residual variance matrix is numerically singular.", call. = FALSE)
  } 
  root_Vr <- Ur %*% diag(value, nrow = length(value)) %*% t(Ur)
  mainMatrix <- root_Vr %*% Vm %*% root_Vr
  temp1 <- eigen(mainMatrix, symmetric = TRUE)
  weights <- root_Vr %*% temp1$vectors
  d <- temp1$values
  
  out <- list("residual" = Vr,
              "model" = Vm,
              "weights" = weights[,1, drop = FALSE],
              "rootVr" = root_Vr,
              "largestRoot" = d[1],
              "rho" = rho)
  if (ncol(pcevObj$X) > 2) out$otherWeights <- weights[,2:(ncol(pcevObj$X) - 1), drop = FALSE]
  
  return(out)
}

#' @rdname estimatePcev
estimatePcev.PcevBlock <- function(pcevObj, shrink, index, ...) {
  p <- ncol(pcevObj$Y)
  N <- nrow(pcevObj$Y)
  
  if (is.null(index) || p != length(index)) {
    stop("index should have length equal to number of response variables")
  }
  
  d <- length(unique(index))
  if (d > N && ncol(pcevObj$X) != 2) {
    warning("It is recommended to have a number of blocks smaller than the number of observations")
  }
  Ypcev <- matrix(NA, nrow = N, ncol = d)
  weights <- rep_len(0, p)
  rootVr <- list("first" = vector("list", d), 
                 "second" = NA)
  
  counter <- 0
  for (i in unique(index)) {
    counter <- counter + 1
    pcevObj_red <- pcevObj 
    pcevObj_red$Y <- pcevObj$Y[, index == i, drop = FALSE]
    class(pcevObj_red) <- "PcevClassical"
    result <- estimatePcev(pcevObj_red, shrink)
    weights[index == i] <- result$weights
    Ypcev[,counter] <- pcevObj_red$Y %*% weights[index == i]
    rootVr$first[[counter]] <- result$rootVr
  }
  
  pcevObj_total <- pcevObj
  pcevObj_total$Y <- Ypcev
  class(pcevObj_total) <- "PcevClassical"
  
  if (ncol(pcevObj_total$X) == 2) {
    fit_total <- lm.fit(pcevObj_total$X, pcevObj_total$Y)
    beta_total <- coefficients(fit_total)[2,]
    weight_step2 <- beta_total/c(crossprod(beta_total))
    
  } else {
    result <- estimatePcev(pcevObj_total, shrink)
    weight_step2 <- result$weights
    rootVr$second <- result$rootVr
  }
  
  counter <- 0
  for (i in unique(index)) {
    counter <- counter + 1
    weights[index == i] <- weights[index == i] * weight_step2[counter]
  }
  
  return(list("weights" = weights,
              "rootVr" = rootVr,
              "index" = index))
}

#' @rdname estimatePcev
estimatePcev.PcevSingular <- function(pcevObj, shrink, index, ...) {
  #initializing parameters
  rho <- NULL
  Y <- pcevObj$Y
  N <- nrow(Y)
  p <- ncol(Y)
  
  
  # Variance decomposition
  # fit <- lm.fit(cbind(pcevObj$X, pcevObj$Z), Y)
  Yfit <- getFitted(cbind(pcevObj$X, pcevObj$Z), Y)
  # Intercepts <- fit$coefficients[1,]
  res <- Y - Yfit
  # fit_confounder <- lm.fit(cbind(rep_len(1, N), pcevObj$Z), Y)
  Yfit_confounder <- getFitted(cbind(rep_len(1, N), pcevObj$Z), Y)
  
  Vr <- crossprod(res)
  Vm <- crossprod(Yfit - Yfit_confounder, Y)
  
  svdRes <- corpcor::fast.svd(res)
  rankVr <- corpcor::rank.condition(res)$rank
  eigVecVr <- svdRes$v[, 1:rankVr]
  eigValVrInv <- 1/svdRes$d[1:rankVr]
  Xp <- eigVecVr %*% diag(eigValVrInv)
  C <- crossprod(Xp, Vm %*% Xp)
  svdC <- corpcor::fast.svd(C)
  Xpp <- svdC$u
  singWeights <- Xp %*% Xpp
  largestRoot <- max(crossprod(singWeights,  Vm %*% singWeights))
  
  out <- list("residual" = Vr,
              "model" = Vm,
              "weights" = singWeights[,1, drop = FALSE],
              "rootVr" = NULL,
              "largestRoot" = largestRoot,
              "rho" = rho,
              "rank" = rankVr)
  if (ncol(pcevObj$X) > 2) out$otherWeights <- singWeights[, -1]
  
  return(out)
}