[R] sandwich package: HAC estimators

[Achim Zeileis]

On Sat, 28 May 2016, T.Riedle wrote:

> Dear R users,
>
> I am running a logistic regression using the rms package and the code 
> looks as follows:
>
> crisis_bubble4<-lrm(stock.market.crash~crash.MA+bubble.MA+MP.MA+UTS.MA+UPR.MA+PPI.MA+RV.MA,data=Data_logitregression_movingaverage)
>
> Now, I would like to calculate HAC robust standard errors using the 
> sandwich package assuming the NeweyWest estimator which looks as 
> follows:
>
> coeftest(crisis_bubble4,df=Inf,vcov=NeweyWest)
>
> Error in match.arg(type) :
>
>  'arg' should be one of "li.shepherd", "ordinary", "score", 
> "score.binary", "pearson", "deviance", "pseudo.dep", "partial", 
> "dfbeta", "dfbetas", "dffit", "dffits", "hat", "gof", "lp1"
>
> As you can see, it doesn't work.

Yes. The "sandwich" package relies on two methods being available: bread() 
and estfun(). See vignette("sandwich-OOP", package = "sandwich") for the 
background details.

For objects of class "lrm" no such methods are available. But as "lrm" 
objects inherit from "glm" the corresponding methods are called. However, 
"lrm" objects are actually too different from "glm" objects (despite the 
inheritance) resulting in the error.

It is easy to add these methods, though, because "lrm" brings all the 
necessary information:

bread.lrm <- function(x, ...) vcov(x) * nobs(x)
estfun.lrm <- function(x, ...) residuals(x, "score")

> Therefore, I did the same using the glm() instead of lrm():
>
> crisis_bubble4<-glm(stock.market.crash~crash.MA+bubble.MA+MP.MA+UTS.MA+UPR.MA+PPI.MA+RV.MA,family=binomial("logit"),data=Data_logitregression_movingaverage)
>
> If I use the coeftest() function, I get following results.
>
> coeftest(crisis_bubble4,df=Inf,vcov=NeweyWest)
>
> z test of coefficients:
>
>              Estimate Std. Error z value Pr(>|z|)
> (Intercept)   -5.26088    5.01706 -1.0486  0.29436
> crash.MA       0.49219    2.41688  0.2036  0.83863
> bubble.MA     12.12868    5.85228  2.0725  0.03822 *
> MP.MA        -20.07238  499.37589 -0.0402  0.96794
> UTS.MA       -58.18142   77.08409 -0.7548  0.45038
> UPR.MA      -337.57985  395.35639 -0.8539  0.39318
> PPI.MA       729.37693  358.60868  2.0339  0.04196 *
> RV.MA        116.00106   79.52421  1.4587  0.14465
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Some of these coefficients and standard errors are suspiciously large. It 
might make sense to check for quasi-complete separation.

> I am unsure whether the coeftest from the lmtest package is appropriate 
> in case of a logistic regression.

Yes, this is ok. (Whether or not the application of HAC standard errors is 
the best way to go is a different matter though.)

> Is there another function for logistic regressions? Furthermore, I would 
> like to present the regression coefficients, the F-statistic and the HAC 
> estimators in one single table. How can I do that?

Running first coeftest() and then lrtest() should get you close to what 
you want - even though it is not a single table.

> I thought it would be useful to incorporate the HAC consistent 
> covariance matrix into the logistic regression directly and generate an 
> output of coefficients and the corresponding standard errors. Is there 
> such a function in R?

Not with HAC standard errors, I think.

> Thanks for your support.
>
> Kind regards
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>