[R] sandwich package: HAC estimators
Achim.Zeileis at uibk.ac.at
Sat May 28 20:50:47 CEST 2016
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:
> Now, I would like to calculate HAC robust standard errors using the
> sandwich package assuming the NeweyWest estimator which looks as
> 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
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
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():
> If I use the coeftest() function, I get following results.
> 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
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
More information about the R-help