-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Frank
E Harrell Jr
Sent: Sunday, December 24, 2006 10:29 AM
To: John Fox
Cc: r-help at stat.math.ethz.ch; 'ivo welch'; 'Dirk Eddelbuettel'
Subject: Re: [R] extend summary.lm for hccm?
John Fox wrote:
Dear Dirk and Ivo,
It's true that the sandwich package has more extensive
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
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",
"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),
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