R/getCoefVC.onePhase.R

Defines functions getCoefVC.onePhase

Documented in getCoefVC.onePhase

##' Get Variance Components' Coefficients and Mean Squares for Single-Phase or
##' Two-Phase Experiments
##' 
##' Compute the variance components' coefficients and corresponding to random
##' effects in the expected mean squares of ANOVA table in single-phase or
##' two-phase experiments. These coefficients are then inserted to a matrix
##' where the rows correspond to each source of variation and column correspond
##' to DF and every variance component. The mean squares is calculated if the
##' \code{response} argument is used.
##' 
##' The main purpose of this function is to combine the matrices presenting
##' every source of variation of the ANOVA table and the variance matrix to
##' compute the coefficients of the variance components.
##' 
##' The complication arise in giving the row names of the matrix for the source
##' of variation in the ANOVA table.
##' 
##' @param Pb a list of matrices generated by \code{\link{infoDecompMat}}
##' function.
##' @param design.df a data frame containing the experimental design. Requires
##' every column be a \code{\link{factor}}.
##' @param v.mat a list of matrix generated by \code{\link{getVMat.onePhase}}
##' or \code{\link{getVMat.twoPhase}}.
##' @param response a numeric vector contains the responses from the
##' experiment.
##' @param table.legend a logical allows users to generate a legend for the
##' variance components of the ANOVA table for large designs. Default is
##' \code{FALSE}, resulting in the use of original block factor names.
##' @param decimal a logical allows users to display the coefficients as the
##' decimals. Default is \code{FALSE}, resulting in the use of
##' \code{fractions}.
##' @param digits a integer indicating the number of decimal places. Default is
##' 2, resulting in 2 decimal places.
##' @return A matrix containing the characters.
##' @importFrom MASS fractions ginv
##' @author Kevin Chang
##' @export
##' @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)
##' 
##' trt.str <- "Trt"              
##' fT <- terms(as.formula(paste("~", trt.str, sep = "")), keep.order = TRUE)  #fixed terms
##' 
##' trtTerm <- attr(fT, "term.labels")
##' effectsMatrix <- attr(fT, "factor")        
##' 
##' T <- makeContrMat(design1, trtTerm, effectsMatrix, contr.vec = NA)
##' 
##' N =  makeOverDesMat(design1, trtTerm)
##' 
##' Replist = getTrtRep(design1, trtTerm)   
##'  
##' Rep <- Replist$Rep
##' trt.Sca <- Replist$Sca
##'     
##' effFactors = lapply(makeOrthProjectors(Z), function(z) 
##'       getEffFactor(z, T, N, Rep, trt.Sca))
##' 
##' #Now construct variance matrices
##' Pb <- effFactors[sort(1:length(effFactors), decreasing=TRUE)]
##' 
##' v.mat <- getVMat.onePhase(Z.Phase1 = Z, design.df = design.df, var.comp = NA)
##'     
##' getCoefVC.onePhase(Pb = Pb, design.df = design1, v.mat = v.mat, response = NA, 
##'     table.legend = FALSE, decimal = FALSE, digit = 2)
##' 
getCoefVC.onePhase <- function(Pb, design.df, v.mat, response, table.legend, decimal, 
    digits) {
    
    if (all(is.na(response))) {
        response <- rep(NA, nrow(design.df))
    }
    MS <- lapply(Pb, function(q) lapply(q[[1]], function(w) t(response) %*% w %*% response))
    #MS <- lapply(Pb, function(q) apply(q, 3, function(w) t(response) %*% w %*% response))
 
    
    V <- v.mat
    
    VC <- rep("1", length(V) + 2)
    names(VC) <- c("DF", names(V), "MS")
    VC <- t(data.frame(VC, stringsAsFactors = TRUE ))
    ############################################################################################################## 
    
    for (i in 1:length(Pb)) {
        tmp <- matrix(0, nrow = length(names(Pb[[i]][[1]])), ncol = length(V) + 1, dimnames = list(names(Pb[[i]][[1]]), 
            c(names(V), "MS")))
        for (j in 1:(length(names(Pb[[i]][[1]])))) {
            for (z in 1:(length(V))) {
                tmp[j, z] <- tr(Pb[[i]][[1]][[j]] %*% V[[z]])
            }
            tmp[j, (length(V) + 1)] <- as.numeric(MS[[i]][[j]])
        }
        
        
        if (nrow(tmp) == 1 && rownames(tmp) == "Residual") {
            tmp <- c(tmp[1], tmp/tmp[1])
            if (decimal) {
                VC <- rbind(VC, c(round(tmp[-length(tmp)], digits = digits), round(tmp[length(tmp)], 
                  digits = 3)))
            } else {
                VC <- rbind(VC, c(attr(fractions(tmp[-length(tmp)]), "fracs"), round(tmp[length(tmp)], 
                  digits = 3)))
            }
            
            if (grepl("Within", names(Pb[i]))) {
                rownames(VC)[nrow(VC)] <- paste(names(Pb[i]), sep = " ")
            } else {
                rownames(VC)[nrow(VC)] <- paste("Between", names(Pb[i]), sep = " ")
            }
            
        } else {
            VC <- rbind(VC, character(length = length(V) + 2))
            
            if (grepl("Within", names(Pb[i]))) {
                rownames(VC)[nrow(VC)] <- paste(names(Pb[i]), sep = " ")
            } else {
                rownames(VC)[nrow(VC)] <- paste("Between", names(Pb[i]), sep = " ")
            }
            
            rownames(tmp) <- paste("  ", rownames(tmp), sep = " ")
            
            if (decimal) {
                tmp <- t(apply(tmp, 1, function(x) c(round(c(x[1], x[-length(x)]/x[1]), 
                  digits = digits), as.character(round(x[length(x)]/x[1], digits = 3)))))
            } else {
                tmp <- t(apply(tmp, 1, function(x) c(attr(fractions(c(round(x[1]), x[-length(x)]/x[1])), 
                  "fracs"), as.character(round(x[length(x)]/x[1], digits = 3)))))
            }
            
            VC <- rbind(VC, tmp)
        }
        
    }
    
    if (all(is.na(response))) {
        VC <- VC[, -ncol(VC)]
    }
    
    if (length((which(apply(apply(VC[, -1], 2, function(x) x == VC[, "e"]), 2, all)))) > 
        1) {
        VC <- noquote(VC[-1, -2])
    } else {
        VC <- noquote(VC[-1, ])
    }
    
    if (table.legend) {
        Legend <- paste(paste(letters[1:(length(colnames(VC)) - 1)], colnames(VC)[-1], 
            sep = " = "))
        colnames(VC)[-1] <- letters[1:(length(colnames(VC)) - 1)]
        VC <- list(VC = VC, Legend = Legend)
    }
    
    return(VC)
} 

Try the infoDecompuTE package in your browser

Any scripts or data that you put into this service are public.

infoDecompuTE documentation built on April 14, 2020, 7:08 p.m.