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

[Achim Zeileis]

On Fri, 10 May 2013, ivo welch wrote:

> 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.

If you use NeweyWest() from "sandwich", the code becomes much shorter:

ols2 <- function (..., newey.west= 0, stdcoefs = TRUE)
{
   ## model and summary
   m <- lm(...)
   rval <- summary(m)

   ## Newey-West standard errors and t-statistic
   nw <- sqrt(diag(NeweyWest(m, lag = newey.west, prewhite = FALSE)))
   rval$coefficients <- cbind(rval$coefficients,
     "NW" = nw, "t (NW)" = coef(m)/nw)

   ## standardized coefficients
   if(stdcoefs) rval$coefficients <- cbind(rval$coefficients, "Std. Coef" =
     coef(m) * apply(model.matrix(m), 2, sd) /
     sd(model.response(model.frame(m))))

   return(rval)
}

The ols2(y ~ x + z) produces output equivalent to ols(y ~ x + z).

Personally, I always just use

   m <- lm(....)
   coeftest(m, vcov = NeweyWest)

to get the coefficient tests with Newey-West standard errors. By default, 
this uses automatic lag selection and prewhitening but you can switch both 
off if you want:

   coeftest(m, vcov = NeweyWest(m, lag = 0, prewhite = FALSE))

Best,
Z

>
> 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
> }
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>