R/makeOrthProjectors.R

In kcha193/infoDecompuTE: Information Decomposition of Two-Phase Experiments

##' Construct Orthogonal Projector Matrices
##' 
##' Construct a list of orthogonal projector matrices corresponding to all
##' strata of the experiment.
##' 
##' The strata decomposition is performed within this function. The first step
##' is to convert the list of block design matrices generated by
##' \code{\link{makeBlkDesMat}} to projection matrices using
##' \code{\link{projMat}}. The second step is to use these projection matrices
##' to project the raw data vector from one stratum to next stratum of the
##' experiment; the resulting matrix corresponds to each stratum is the
##' orthogonal projector matrix of the given stratum.
##' 
##' @param BlkDesList a list of block design matrices generated by
##' \code{\link{makeBlkDesMat}}.
##' @return A list containing matrices.
##' @author Kevin Chang
##' @examples
##' 
##' design1 <- local({ 
##'   Ani = as.factor(LETTERS[c(1,2,3,4,
##'                             5,6,7,8)])
##'   Trt = as.factor(letters[c(1,1,1,1,
##'                             2,2,2,2)])
##'   data.frame(Ani, Trt, stringsAsFactors = TRUE )
##' })
##' 
##' blk.str = "Ani"
##'     
##' rT = terms(as.formula(paste("~", blk.str, sep = "")), keep.order = TRUE) 
##' blkTerm = attr(rT,"term.labels")
##'      
##' Z = makeBlkDesMat(design1, blkTerm)
##' Pb = makeOrthProjectors(Z)
##' 
##' @export makeOrthProjectors
makeOrthProjectors <- function(BlkDesList) {
    
    initial <- diag(nrow(BlkDesList[[1]])) - K(nrow(BlkDesList[[1]]))
    n <- nrow(BlkDesList$e)
    Q <- lapply(BlkDesList, function(z) projMat(z))
    
    Q <- Q[sort(1:length(Q), decreasing = TRUE)]
    
    elementToRemove <- numeric(0)
    
    cusumQ <- P <- NULL
    P[[1]] <- Q[[1]] %*% initial
    
    if (all(P[[1]] < 1e-06)) {
        elementToRemove <- 1  #revord the elements that are all zeros
        P[[1]] <- matrix(0, nrow = n, ncol = n)  #make the projectors which has less than 1e-6 to zero, to avoid rounding error
    }
    
    cusumQ[[1]] <- initial - P[[1]]
    
    for (i in 2:(length(Q))) {
        P[[i]] <- Q[[i]] %*% cusumQ[[i - 1]]
        
        if ((all(P[[i]] < 1e-06)) || tr(P[[i]]) <= 0) {
            elementToRemove <- c(elementToRemove, i)  #revord the elements that are all zeros
            P[[i]] <- matrix(0, nrow = n, ncol = n)  #make the projectors which has less than 1e-6 to zero, to avoid rounding error
        }
        cusumQ[[i]] <- cusumQ[[i - 1]] - P[[i]]
    }
    
    P <- P[sort(1:length(P), decreasing = TRUE)]
    names(P) <- names(BlkDesList)
    
    if (length(elementToRemove) > 0) 
        P <- P[-(length(Q) - elementToRemove + 1)]
    
    return(P)
}