extend summary.lm for hccm?
John Fox wrote:
Dear Dirk and Ivo, It's true that the sandwich package has more extensive facilities for sandwich variance estimation than the hccm function in the car package, but I think that the thrust of Ivo's question is to get a model summary that automatically uses the adjusted standard errors rather than having to compute them and use them "manually." Here's one approach, which could be modified slightly if Ivo wants "hc0" as the default; it could also be adapted to use the sandwich package. I hope this helps, John
Another approach:
library(Design) # also requires library(Hmisc)
f <- ols(y ~ x1 + x2, x=TRUE, y=TRUE)
f <- robcov(f) # sandwich; also allows clustering. Also see bootcov
anova(f) # all later steps use sandwich variance matrix
summmary(f)
contrast(f, list(x1=.5), list(x1=.2))
BUT note that sandwich covariance matrix estimators can have poor mean
squared error (a paper by Bill Gould in Stata Technical Bulletin years
ago related to logistic regression showed an example with a 100-fold
increase in the variance of a variance estimate) and can give you the
right estimate of the wrong quantity (reference below).
Frank Harrell
@Article{free06so,
author = {Freedman, David A.},
title = {On the so-called ``{Huber} sandwich
estimator'' and
``robust standard errors''},
journal = The American Statistician,
year = 2006,
volume = 60,
pages = {299-302},
annote = {nice summary of derivation of sandwich
estimators;questions why we should be interested in getting the right
variance of the wrong parameters when the model doesn't fit}
}
----------- snip --------------
summaryHCCM <- function(model, ...) UseMethod("summaryHCCM")
summaryHCCM.lm <- function(model, type=c("hc3", "hc0", "hc1", "hc2", "hc4"),
...){
if (!require(car)) stop("Required car package is missing.")
type <- match.arg(type)
V <- hccm(model, type=type)
sumry <- summary(model)
table <- coef(sumry)
table[,2] <- sqrt(diag(V))
table[,3] <- table[,1]/table[,2]
table[,4] <- 2*pt(abs(table[,3]), df.residual(model), lower.tail=FALSE)
sumry$coefficients <- table
print(sumry)
cat("Note: Heteroscedasticity-consistant standard errors using
adjustment",
type, "\n")
}
--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
--------------------------------
-----Original Message----- From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Dirk Eddelbuettel Sent: Saturday, December 23, 2006 11:16 PM To: ivo welch Cc: r-help at stat.math.ethz.ch Subject: Re: [R] extend summary.lm for hccm? On 23 December 2006 at 20:46, ivo welch wrote: | I wonder whether it is possible to extend the summary method for the | lm function, so that it uses an option "hccm" (well, model "hc0"). In | my line of work, it is pretty much required in reporting of almost all | linear regressions these days, which means that it would be very nice | not to have to manually library car, then sqrt the diagonal, and | recompute T-stats; instead, I would love to get everything in the same | format as the current output---except errors heteroskedasticity | adjusted. | | easy or hard? Did you consider the 'sandwich' package? A simple
> install.packages("sandwich")
> library(sandwich)
> ?vcovHC
> ?vcovHAC
should get you there. Dirk --
Frank E Harrell Jr Professor and Chair School of Medicine
Department of Biostatistics Vanderbilt University