extend summary.lm for hccm?

Dear Frank,

If I remember Freedman's recent paper correctly, he argues that sandwich
variance estimator, though problematic in general, is not problematic in the
case that White described -- an otherwise correctly specified linear model
with heteroscedasticity estimated by least-squares. 

Regards,
 John

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

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