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

Achim Zeileis Achim.Zeileis at uibk.ac.at
Sat May 11 08:34:59 CEST 2013

```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)
>      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