R/Lfx.R

###
### Copyright 2012  Georges Monette <georges@yorku.ca>
###
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 3 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  http://www.r-project.org/Licenses/
#

###
### Lfx.R
### part of the spida package
#
# 2012 12 05:
#  Modelled and improved from Lform: Allows estimation of effects of arbitrary
#  order evaluated over a data frame.  Correctly handles interaction terms involving
#  matrices. 
#  Lfx generates a matrix that can be used as a linear hypothesis matrix in 'wald'
#' @export
getData <- function(x,...) UseMethod("getData")



#' Hypothesis matrix generated by expressions for each column or term %%
#' ~~function to do ... ~~
#' 
#' Creates an L matrix using expressions evaluated in \code{data} for each
#' column or multi-column term of the L matrix %% ~~ A concise (1-5 lines)
#' description of what the function does. ~~
#' 
#' %% ~~ If necessary, more details than the description above ~~
#' 
#' @aliases Lfx M M.default M.factor M.formula M.M <.factor >.factor <=.factor
#' =>.factor
#' @param fit a fitted model with a 'getFix' method. %% ~~Describe \code{fit}
#' here~~
#' @param expr.list a list of expressions with one component for each column
#' (or groups of columns) of the hypothesis matrix corresponding to each term
#' of the model. A term with multiple degrees of freedom can either be
#' generated as separate single terms or with an expression that evaluates to a
#' suitable matrix.  %% ~~Describe \code{expr.list} here~~
#' @param data the data frame in which expressions are evaluated. %% ~~Describe
#' \code{data} here~~
#' @param wrap
#' 
#' %% ~~Describe \code{wrap} here~~
#' @param debug %% ~~Describe \code{debug} here~~
#' @param formula as an argument of \code{M}, a one-sided formula defining a
#' main effect or an interaction term involving factors.
#' @param expr as an argument of \code{M}, an expression that can be evaluated
#' in each row of \code{data} to form element(s) of the corresponding row of
#' \code{L}.  formula defining a main effect or an interaction term involving
#' factors.
#' @note %% ~~further notes~~
#' @author %% ~~who you are~~
#' @examples
#' 
#' ##
#' ## The estimated effect of an additional year of education
#' ##
#' 
#' data(Prestige)   # from car
#' fit <- lm( income ~ (education+I(education^2)) * type, Prestige)
#' summary(fit)
#' Lfx(fit)  # geneates expression to cut and paste and differentiate
#' 
#' pred <- expand.grid( education = 3:16, type = levels(Prestige$type))
#' 
#' Lf.ed <- Lfx( fit, 
#'               list( 0,
#'                     1 + 0* education,
#'                     2 * M(I(education)),
#'                     0 * M(type),
#'                     1 * 1 * M(type),
#'                     2 * M(I(education)) * M(type) 
#'               ),pred)
#' zw <- as.data.frame(wald(fit,Lf.ed))
#' head(zw)
#' xyplot( coef ~ education , zw, groups = type, type = 'l', auto.key = T)
#' 
#' fit2 <- lm( income ~ (education+I(education^2)) * type*women, Prestige)
#' pred <- expand.grid( education = 3:16, type = levels(Prestige$type), women=c(0,50,100))
#' Lfx(fit2)
#' zw <- Lfx(fit2,
#'     list( 0,   # differentiated wrt education
#'           1 * 1,
#'           1 * M(I(education)),
#'           0 * M(type),
#'           0 * women,
#'           1 * 1 * M(type),
#'           2 * M(I(education)) * M(type),
#'           1 * 1 * women,
#'           2 * M(I(education)) * women,
#'           0 * M(type) * women,
#'           1 * 1 * M(type) * women,
#'           2 * M(I(education)) * M(type) * women 
#'     ),pred)
#' zw    
#' zw <- as.data.frame( wald(fit2, zw))
#' summary(fit2)
#' xyplot( coef ~ education|women , zw, groups = type, type = 'l', auto.key = T,
#'         ylim = c(-30000,30000))
#' 
#' ###
#' ###  Differencing of factor
#' ###
#' 
#' pred <- expand.grid( education = 3:18, women = c(0,50,100), 
#'                      type = levels(Prestige$type),
#'                      type0 = levels(Prestige$type))
#' 
#' Lfx(fit2)
#' 
#' Lcomp <- Lfx(fit2,
#'   list( 0,
#'       0 * education,
#'       0 * M(I(education^2)),
#'       1 * M(type) - M(type0),
#'       0 * women,
#'       1 * education * (M(type)-M(type0)),
#'       1 * M(I(education^2)) * (M(type) - M(type0)),
#'       0 * education * women,
#'       0 * M(I(education^2)) * women,
#'       1 * (M(type)-M(type0)) * women,
#'       1 * education * (M(type) - M(type0)) * women,
#'       1 * M(I(education^2)) * (M(type) - M(type0)) * women 
#'   ), pred)
#' 
#' zw <- as.data.frame( wald( fit2, Lcomp))
#' zw
#' 
#' 
#' 
#' 
#' 
#' 
#' @export
Lfx <-     
  function (fit, expr.list, data = getData(fit), prefix = "1 * ", wrap = FALSE, debug = FALSE) 
  {
    # 2012-12-03:
    # 1. modified to work with M with factor instead of formula
    # 2. Lform(fit) includes comments from Lform0
    # 3. More informative error messages 
    
    Lfx0 <- function(fit, factors = getFactorNames(fit), wrap = FALSE,
                       debug = debug) {
      ts <- colnames(attr(terms(fit),'factors'))
      ts <- strsplit(ts,":")
      if(debug) disp(ts)
      # add functions as well because some might generate matrices:
      
      facl <- lapply( ts, function(x) x %in% factors | grepl(")",x))
      if(debug) disp(facl)
      ret <- lapply( seq_along(facl),  function(i){
        lapply(seq_along(facl[[i]]),
               function(j) if(facl[[i]][j]) 
                 paste("M(",ts[[i]][j],")", sep = "") else ts[[i]][j] 
        )})
      
      ret <- lapply( ret, paste, collapse = " * ")    
      ret <- lapply( ret, function(s) paste(prefix,s,sep=""))  # add prefix 
      ret <- do.call( paste, c("list( 1",ret, sep = ",\n"))
      ret <- paste(ret,"\n)")
      if (wrap) ret <- paste("with(data, do.call( cbind,\n",
                             ret, "\n ))\n")
      cat(ret)
      invisible(ret)
    }
    
    
    if (missing(expr.list)) return( Lfx0(fit, wrap = wrap, debug = debug))
    gg <- getFix(fit)
    Lsub <- do.call(cbind, eval(substitute(expr.list), data))
    L <- Lsub
    rownames(L) <- rownames(data)
    colnames(L) <- names(gg$fixed)
    
    attr(L, "data") <- data
    L
  }

# 2012 12 03: Replaced with following functions
# M <- function(form, ..., keep.intercept = FALSE) {
#   ret <- model.matrix(form,...)
#   if( !keep.intercept ){
#     pos <- grep( '^\\(Intercept\\)$', colnames(ret))
#     ret <- if( length(pos) ) ret[, -pos] 
#   }
#   ret
# }

# To do: return object of type "M" and
# define '*' method to do column by column multiplication
# -- DONE

#' @export
M <- function(x,...) UseMethod("M")

#' @export
M.formula <- function (form, ..., keep.intercept = FALSE) { 
  ret <- model.matrix(form, ...)
  if (!keep.intercept) {
    pos <- grep("^\\(Intercept\\)$", colnames(ret))
    ret <- if (length(pos)) 
      ret[, -pos]
  }
  class(ret) <- 'M'
  ret
}

#' @export
M.factor <- function( x, base = NULL) {
  ret <- contrasts(x)[x,]
  if( ! is.null(base)) {
    if ( !is.factor(base)) stop("base must be a factor")
    if (!identical(levels(x), levels(base))) warning("'x' and 'base' should have identical levels")
    ret <- ret - contrasts(base)[base,]
  }
  class(ret) <- 'M'
  ret
}

#' @export
M.M <- function(x, ...) x

#' @export
M.default <- function(x,...) {
  x <- cbind(x)
  class(x) <- c('M',attr(x,"class"))
  x
}

#' @export
"*.M" <- function (e1, e2) {
  # column by column multiplication
  # to expand interaction terms
  # LHS runs fastest
  if( !inherits(e1,"M") && is.numeric(e1) && length(e1) == 1) {
    return( M( e1 * unclass(e2)))    
  }
  if( !inherits(e2,"M") && is.numeric(e2) && length(e2) == 1) {
    return( M( unclass(e1) * e2))    
  }
  e1 <- cbind(e1)
  e2 <- cbind(e2)
  n1 <- ncol(e1)
  n2 <- ncol(e2)
  ret <- c(e1) * c(apply(e2, 2, rep, n1))
  ret <- matrix(ret, nrow = nrow(e1))
  #  rownames(ret) <- rownames(e1)
  if (!is.null(c1 <- colnames(e1)) && !is.null(c2 <- colnames(e2))) {
    cn <- paste(rep(c1, length(c2)), rep(c2, each = length(c1)), 
                sep = ":")
    colnames(ret) <- cn
  }
  class(ret) <- "M"
  ret
}

# To make is easier to generate pairwise comparisons without
# redundancies
# e.g.
#  > pred <- expand.grid( type = levels(Prestige$type),
#                        type0 = levels(Prestige$type),
#                        education = seq(6,18,3))
#  > pred <- subset( pred, type0 < type)


#' @export
"<.factor" <- function(x,y) {
  if ( ! identical( levels(x), levels(y))) warning( "Different 'levels' in x and y")
  as.numeric(x) < as.numeric(y)
}


#' @export
">.factor" <- function(x,y) {
  if ( ! identical( levels(x), levels(y))) warning( "Different 'levels' in x and y")
  as.numeric(x) > as.numeric(y)
}

#' @export
"<=.factor" <- function(x,y) {
  if ( ! identical( levels(x), levels(y))) warning( "Different 'levels' in x and y")
  as.numeric(x) <= as.numeric(y)
}

#' @export
">=.factor" <- function(x,y) {
  if ( ! identical( levels(x), levels(y))) warning( "Different 'levels' in x and y")
  as.numeric(x) >= as.numeric(y)
}
gmonette/spida15 documentation built on May 17, 2019, 7:26 a.m.