R's basic lm() and summary.lm functions
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.