[R] Robust SE for lrm object
Achim Zeileis
Achim.Zeileis at uibk.ac.at
Sat Mar 6 07:22:08 CET 2010
On Sat, 6 Mar 2010, David Winsemius wrote:
> On Mar 5, 2010, at 11:54 PM, Patrick Shea wrote:
>
>>
>> I'm trying to obtain the robust standard errors for a multinomial ordered
>> logit model:
>>
>> mod6 <- lrm(wdlshea ~ initdesch + concap + capasst + qualrat +
>> terrain,data=full2)
>>
>> The model is fine but when I try to get the RSE I get an error.
>>
>> coeftest(mod6, vcov = vcovHAC(mod6))
>>
>> Error in match.arg(type) :
>> 'arg' should be one of ?ordinary?, ?score?, ?score.binary?, ?pearson?,
>> ?deviance?, ?pseudo.dep?, ?partial?, ........etc.
>>
>> I'm a novice R user and am not sure how to address this problem. I have
>> also tried to use alternatives (zelig, polr) but have had no luck. Any
>> assistance on generating RSE for a multinomial order logit model would be
>> appreciated
>
> Have you loaded the library that contains the vcovHAC function?
That is in the "sandwich" package. However, I doubt that it makes sense in
this context. Using HAC covariances would imply that you have time series
data and want to correct for heteroskedasticity and autocorrelation. I'm
not even sure whether sandwich standard errors would be terribly useful.
Both would require that you correctly specified the estimating functions
of your proportional odds logistic regression but misspecified a few other
aspects of the remaining likelihood. Not sure whether that can be obtained
for an ordinal multinomial response.
> (And do you know whether coeftest works with Design/rms objects?)
It does (unlike its own summary function in some situations):
library("rms")
library("lmtest")
data("BankWages", package = "AER")
fm <- lrm(job ~ ., data = BankWages)
summary(fm)
coeftest(fm)
The reason why vcovHAC() or sandwich() do not work is that bread() and
estfun() methods would need to be available for "lrm" objects which is
currently not the case (dito for "polr" objects). In principle they could
be written, see
vignette("sandwich-OOP", package = "sandwich")
but as I said above I'm not sure whether it would be very useful.
Z
> --
>
> David Winsemius, MD
> West Hartford, CT
>
> ______________________________________________
> R-help at r-project.org mailing list
> 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.
More information about the R-help
mailing list