#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.