R/infoDecompMat.R

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

##' Construct the Matrix from Information Decomposition
##' 
##' Perform the information decomposition for either the block or treatment
##' effects within a single stratum.
##' 
##' The main purpose of this function is to construct a list of resultant
##' matrices associated with each source of variation after the information
##' decomposition.
##' 
##' This list of matrices are then used to compute the coefficient of the
##' variance components in the expected mean squares.
##' 
##' @param z a matrix containing the orthogonal projector for a single stratum
##' generated by \code{\link{makeOrthProjectors}}.
##' @param T a list of contrast matrices generated by
##' \code{\link{makeContrMat}}.
##' @param N a matrix containing the design matrix generated by
##' \code{\link{makeOverDesMat}}.
##' @return A list of 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)
##' 
##' trt.str = "Trt"              
##' fT <- terms(as.formula(paste("~", trt.str, sep = "")), keep.order = TRUE)  #fixed terms
##' 
##' trtTerm <- attr(fT, "term.labels")
##' effectsMatrix <- attr(fT, "factors")        
##' 
##' T <- makeContrMat(design1, trtTerm, effectsMatrix, contr.vec = NA)
##' 
##' N =  makeOverDesMat(design1, trtTerm)
##' 
##' 
##' infoDecompMat(Pb[[1]], T, N)
##' 
##' 
##' @export infoDecompMat
infoDecompMat <- function(z, T, N) {
    if (!is.matrix(z)) 
        return(z)
    
    nEffect <- length(T)
    PNTginvATNP <- T
    
   PNTginvATNP[[1]] <- z %*% N %*% T[[1]] %*% invInfMat(C = z, N = N, T = T[[1]]) %*% 
        T[[1]] %*% t(N) %*% t(z)
    
    # PNTginvATNP[[1]] = z %*% N %*% T[[1]] %*% ginv(t(T[[1]]) %*% t(N) %*% z %*% N %*%
    # T[[1]]) %*% t(N) %*% t(z)
    
    newZ <- (z %*% t(z)) - PNTginvATNP[[1]]
    
    if (nEffect != 1) {
        for (i in 2:nEffect) {
            PNTginvATNP[[i]] <- newZ %*% N %*% t(T[[i]]) %*% invInfMat(C = newZ, N = N, 
                T = T[[i]]) %*% T[[i]] %*% t(N) %*% t(newZ)
            
            # PNTginvATNP[[i]] = newZ %*% N %*% t(T[[i]]) %*% ginv(t(T[[i]]) %*% t(N) %*% z %*%
            # N %*% T[[i]]) %*% T[[i]] %*% t(N) %*% t(newZ)
            newZ <- (newZ %*% t(newZ)) - PNTginvATNP[[i]]
        }
    }
    
    PNTginvATNP$Residual <- newZ
    elementToRemove <- numeric(0)
    for (i in 1:length(PNTginvATNP)) {
        if (all(PNTginvATNP[[i]] < 1e-06)) 
            elementToRemove <- c(elementToRemove, i)
    }
    
    if (length(elementToRemove) > 0) 
        PNTginvATNP <- PNTginvATNP[-elementToRemove]
    
    return(PNTginvATNP)
}