R/linest_matrix.R

In hojsgaard/doBy: Groupwise Statistics, LSmeans, Linear Estimates, Utilities

## #############################################################################
#'
#' @title Linear estimates matrix
#' @description Generate matrix specifying linear estimate.
#' @name linest-matrix
#'
## #############################################################################
#'
#' @param object Model object
#' 
#' @param effect A vector of variables. For each configuration of
#'     these the estimate will be calculated.
#'
#' @param at Either NULL, a list or a dataframe. 1) If a list, then the list
#'     must consist of covariates (including levels of some factors)
#'     to be used in the calculations. 2) If a dataframe, the
#'     dataframe is split rowwise and the function is invoked on each
#'     row.
#'
#' @param linest_list Linear estimate list (as generated by \code{get_linest_list}).
#' 
#' @details Check this
#'
#' @seealso \code{\link{LSmeans}}, \code{\link{linest}}
#' @keywords utilities
#' 
#' @examples
#' 
#' ## Two way anova:
#' 
#' data(warpbreaks)
#'
#' ## An additive model
#' m0 <- lm(breaks ~ wool + tension, data=warpbreaks)
#'
#' ## Estimate mean for each wool type, for tension="M":
#' K <- LE_matrix(m0, at=list(wool=c("A", "B"), tension="M"))
#' K
#' 
#' ## Vanilla computation:
#' K %*% coef(m0)
#'
#' ## Alternative; also providing standard errors etc:
#' linest(m0, K)
#' esticon(m0, K)
#'
#' ## Estimate mean for each wool type when averaging over tension;
#' # two ways of doing this
#' K <- LE_matrix(m0, at=list(wool=c("A", "B")))
#' K
#' K <- LE_matrix(m0, effect="wool")
#' K
#' linest(m0, K)
#'
#' ## The linear estimate is sometimes called to "least squares mean"
#' # (LSmeans) or popupulation means.
#' # Same as
#' LSmeans(m0, effect="wool")
#'
#' ## Without mentioning 'effect' or 'at' an average across all
#' #predictors are calculated:
#' K <- LE_matrix(m0)
#' K
#' linest(m0, K)
#' 
#' ## Because the design is balanced (9 observations per combination
#' #of wool and tension) this is the same as computing the average. If
#' #the design is not balanced, the two quantities are in general not
#' #the same.
#' mean(warpbreaks$breaks)
#'
#' ## Same as 
#' LSmeans(m0)
#'
#' ## An interaction model 
#' m1 <- lm(breaks ~ wool * tension, data=warpbreaks)
#' 
#' K <- LE_matrix(m1, at=list(wool=c("A", "B"), tension="M"))
#' K
#' linest(m1, K)
#' K <- LE_matrix(m1, at=list(wool=c("A", "B")))
#' K
#' linest(m1, K)
#' K <- LE_matrix(m1, effect="wool")
#' K
#' linest(m1, K)
#' LSmeans(m1, effect="wool")
#' 
#' K <- LE_matrix(m1)
#' K
#' linest(m1, K)
#' LSmeans(m1)

#' @export
#' @rdname linest-matrix
LE_matrix <- function(object, effect=NULL, at=NULL){
  UseMethod("LE_matrix")
}

## FIXME: LE_matrix.default: Should be a check of what 'object' is
#' @export
#' @rdname linest-matrix
LE_matrix.default <- function(object, effect=NULL, at=NULL){

    effect <- .parse.effect(effect)

    ## print("HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH")
    ## str(list(effect=effect, at=at))
    
    if (!is.null(at))
        if (!inherits(at, c("list", "data.frame"))) stop("'at' must be list or data.frame\n")

    if (!inherits(at, "data.frame")){
        out <- get_linest_list(object, effect, at)
        out <- aggregate_linest_list (out)
    } else {
        ## call function for each row in at
        at.list <- split(at, 1:nrow(at))

        
        at.list <- lapply(at.list, as.list)
        out  <- lapply(at.list, function(at.arg)
            LE_matrix(object, effect, at.arg))
        out <- do.call(rbind, out)
    }
    class(out) <- c("linest_matrix_class", "matrix")
    out    
}



## OLD
## LE_matrix.default <- function(object, effect=NULL, at=NULL){
##     out <- get_linest_list(object, effect, at)
##     out <- aggregate_linest_list (out)
##     class(out) <- c("linest_matrix_class", "matrix")
##     out
## }

#' @export
#' @rdname linest-matrix
aggregate_linest_list <- function (linest_list){
    out               <- lapply( linest_list, function( mm ) apply( mm, 2, mean ) )
    out               <- do.call(rbind, out)
    attr(out, "at")   <- attr(linest_list, "at")
    attr(out, "grid") <- attr(linest_list, "grid")
    attr(out, "offset") <- attr(linest_list, "offset")
    out
}

#' @export
print.linest_matrix_class <- function(x, ...){
  prmatrix(x)
  invisible(x)
}

#' @export
summary.linest_matrix_class <- function(object, ...){
    print(object)
    cat("at: \n"); print(attr(object, "at"))
    cat("grid: \n"); print(attr(object, "grid"))
    invisible(object)      
}


.parse.effect <- function(effect=NULL){
    ## All returns c("A", "B")
    ## parse.effect(~A+B)
    ## parse.effect(y~A+B)
    ## parse.effect(c("A", "B"))

    if (is.null(effect)) return(NULL)
    if (!inherits(effect, c("character", "formula")))
        stop("'effect' must be a character vector or a formula")
    
    if (inherits(effect, "character")){
        if (identical(effect, character(0)) | identical(effect, ""))
            NULL
        else
            effect
    } else {
        all.vars(effect[[length(effect)]])
    }
}

## This is the workhorse for generating the "contrast matrix"
#' @export
#' @rdname linest-matrix
get_linest_list <- function(object, effect=NULL, at=NULL){

    pr <- FALSE
    ##cat(".get_linest_list\n")

    effect <- .parse.effect(effect)
    
    trms     <- delete.response( terms(object) )
    fact.lev <- get_xlevels( object )            ## factor levels
    ## cat("+++ fact.lev:\n"); print(fact.lev)
    cov.ave  <- .get_covariate_ave(object, at)  ## average of covariates (except those mentioned in 'at')
    ## cat("+++ cov.ave:\n"); print(cov.ave)
    vartype  <- get_vartypes( object )           ## which are factors and which are numerics
    ## cat("+++ vartype:\n"); print(vartype)
    at.factor.name <- intersect(vartype$factor, names(at))
    cov.ave.name   <- names( cov.ave )
    effect         <- setdiff( effect, at.factor.name )

    if (is.null(effect))
    {
        new.fact.lev <- if (length(at.factor.name) > 0) at[ at.factor.name ]
                        else NULL        
    } else {
        zz  <- set_xlevels(fact.lev, at=at)
        ## print(zz)
        new.fact.lev  <- zz[c(effect, at.factor.name)]
    }

    ## str(list(effect=effect, at=at, fact.lev=fact.lev, new.fact.lev=new.fact.lev, 
    ##          vartype=vartype, cov.ave=cov.ave, at.factor.name=at.factor.name,
    ##          cov.ave.name=cov.ave.name))

    if (is.null(new.fact.lev)){
        ##cat("case: No 'effect' and no 'at'-factors; hence just a global average. \n")
        if (length(fact.lev) > 0)
        {
            ##cat("there are factors in the model\n")
            newdata <- expand.grid(fact.lev)
            if (length( cov.ave.name ) > 0){
                ##cat("there are covariates in the model\n")
                newdata[, cov.ave.name] <- cov.ave ## Augment with covariates
            }
        } else {
            ## cat("No factors in the model\n")
            if (length( cov.ave.name ) > 0){
                ##cat("yes there are covariates\n")
                newdata <- data.frame(matrix(unlist(cov.ave), nrow=1L))                
                names(newdata) <- cov.ave.name
            } else {
                ##cat("there are no factors or covariates\n")
                newdata <- data.frame(1)
            }
        }

        XXlist <- list(get_X(object, newdata=newdata, at=NULL))
        ## cat("XXlist:\n"); print(XXlist)
        attr(XXlist, "at")   <- at[intersect(vartype$numeric, names(at))]
        attr(XXlist, "grid") <- NULL
    }
    else
    {
        if(pr)cat("The general case; there are 'effect' factors or 'at' factors.\n")
        grid.data <- expand.grid(new.fact.lev)
        grid.data <- as.data.frame(lapply(grid.data, as.character), stringsAsFactors=FALSE)
        
        ## nfl <<- new.fact.lev
        ## gd <<- grid.data
        ##str(list(new.fact.lev=new.fact.lev, grid.data=grid.data))

        ## cat("HHHH grid.data\n")
        ## print(grid.data)
        
        XXlist    <- list()
        for (ii in 1:nrow(grid.data)){
            config    <- grid.data[ii, ,drop=FALSE ]
            fact.lev2 <- set_xlevels(fact.lev, at=config)
            newdata   <- expand.grid(fact.lev2)
            newdata[, cov.ave.name]  <- cov.ave
            XX             <- get_X(object, newdata=newdata, at=at)
            XXlist[[ ii ]] <- XX
        }

        grid.data[, names(cov.ave) ] <- cov.ave
        attr(XXlist, "at") <- at
        attr(XXlist, "grid") <- grid.data
        attr(XXlist, "offset") <- attr(XX, "offset")  ## FIXME: reference to XX; what is offset?
    }
    class(XXlist) <- "linest_list_class"
    XXlist
}



## --------------------------------------------------------------------

setOldClass("linest_matrix_class")

setAs("linest_matrix_class", "matrix",
      function(from){
          attr(from, "at") <- attr(from, "grid") <- NULL
          class(from) <- "matrix"
          from
      })

setAs("linest_matrix_class", "Matrix",
      function(from){
          attr(from, "at") <- attr(from, "grid") <- NULL
          class(from) <- "matrix"
          as(from, "Matrix")
      })