[R] R's basic lm() and summary.lm functions

[ivo welch]

I ended up wrapping my own new "ols()" function in the end.  it is my
replacement for lm() and summary.lm.  this way, I don't need to alter
internals.  in case someone else needs it, it is included.  of course,
feel free to ignore.


docs[["ols"]] <- c(Rd= '
@TITLE ols.R
@AUTHOR ivo.welch at gmail.com
@DATE 2013
@DESCRIPTION
  adds newey-west and stdandardized coefficients to the lm function,
and adds the summary.lm information at the same time.
@USAGE ols(..., newey.west=0, stdcoefs=TRUE)
@ARGUMENTS
@DETAILS
@SEEALSO
@EXAMPLES
', test= '
  x <- rnorm(12); y <- rnorm(12); z <- rnorm(12); x[2] <-NA;
  ols( y ~ x + z )
', changes= '
')

ols <- function (..., x = FALSE, newey.west=(0), stdcoefs=TRUE) {

  ## R is painfully error-tolerant. I prefer reasonable and immediate
error warnings.
  stopifnot( (is.vector(newey.west))&(length(newey.west)==1)|(is.numeric(newey.west))
)
  stopifnot( (is.vector(stdcoefs))&(length(stdcoefs)==1)|(is.logical(stdcoefs))
)
  stopifnot( (is.vector(x))&(length(x)==1)|(is.logical(x)) )
  ## I wish I could check lm()'s argument, but I cannot.

  lmo <- lm(..., x=TRUE)
  ## note that both the x matrix and the residuals from the model have
their NA's omitted by default

  if (newey.west>=0) {
    resids <- residuals(lmo)
    diagband.matrix <- function(m, ar.terms) {
      nomit <- m - ar.terms - 1
      mm <- matrix(TRUE, nrow = m, ncol = m)
      mm[1:nomit, (ncol(mm) - nomit + 1):ncol(mm)] <-
(lower.tri(matrix(TRUE, nrow = nomit, ncol = nomit)))
      mm[(ncol(mm) - nomit + 1):ncol(mm), 1:nomit] <-
(upper.tri(matrix(TRUE, nrow = nomit, ncol = nomit)))
      mm
    }
    invx <- chol2inv(chol(crossprod(lmo$x)))
    invx.x <- invx %*% t(lmo$x)
    if (newey.west==0)
      resid.matrix <- diag(resids^2)
    else {
      full <- resids %*% t(resids)
      maskmatrix <- diagband.matrix(length(resids), newey.west)
      resid.matrix <- full * maskmatrix
    }
    vmat <- invx.x %*% resid.matrix %*% t(invx.x)

    nw <- newey.west  ## the number of AR terms
    nw.se <- sqrt(diag(vmat))  ## the standard errors
  }

  if (stdcoefs) stdcoefs.v <- lmo$coefficients*apply(lmo$x,2,sd)/sd(lmo$model$y)

  full.x.matrix <- if (x) lmo$x else NULL
  lmo <- summary(lmo)  ## the summary.lm object
  if (x) lmo$x <- full.x.matrix

  if (stdcoefs) {
    lmo$coefficients <- cbind(lmo$coefficients, stdcoefs.v )
    colnames(lmo$coefficients)[ncol(lmo$coefficients)] <- "stdcoefs"
  }

  if (newey.west>=0) {
    lmo$coefficients <- cbind(lmo$coefficients, nw.se)
    colnames(lmo$coefficients)[ncol(lmo$coefficients)] <-
paste0("se.nw(",newey.west,")")
    lmo$coefficients <- cbind(lmo$coefficients, lmo$coefficients[,1]/nw.se)
    colnames(lmo$coefficients)[ncol(lmo$coefficients)] <-
paste0("Tse.nw(",newey.west,")")
  }

  lmo
}