Dear Achim Zeileis,
Many thanks for your quick and informative answer.
I?m sure that the vcovCL should work, however, I experience some problems.
> coeftest(model, vcov=vcovCL(model, cluster=mydata$ID))
First I got this error:
Error in vcovCL(model, cluster = mydata$ID) :
length of 'cluster' does not match number of observations
After checking the observations I got this error:
Error in vcovCL(model, cluster = mydata$ID) : object 'tets' not found
Called from: vcovCL(model, cluster = mydata$ID)
Browse[1]>
What can I do to fix it? What am I doing wrong now?
> 1 jul 2016 kl. 14:57 skrev Achim Zeileis <Achim.Zeileis at
uibk.ac.at>:
>
> On Fri, 1 Jul 2016, Faradj Koliev wrote:
>
>> Dear all,
>>
>> I use ?polr? command (library: MASS) to estimate an ordered logistic
regression.
>>
>> My model: summary( model<- polr(y ~ x1+x2+x3+x4+x1*x2
,data=mydata, Hess = TRUE))
>>
>> But how do I get robust clustered standard errors?
>> I??ve tried coeftest(resA, vcov=vcovHC(resA, cluster=lipton$ID))
>
> The vcovHC() function currently does not (yet) have a "cluster"
argument. We are working on it but it's not finished yet.
>
> As an alternative I include the vcovCL() function below that computes the
usual simple clustered sandwich estimator. This can be applied to
"polr" objects and plugged into coeftest(). So
>
> coeftest(resA, vcov=vcovCL(resA, cluster=lipton$ID))
>
> should work.
>
>> and summary(a <- robcov(model,mydata$ID)).
>
> The robcov() function does in principle what you want by I'm not sure
whether it works with polr(). But for sure it works with lrm() from the
"rms" package.
>
> Hope that helps,
> Z
>
> vcovCL <- function(object, cluster = NULL, adjust = NULL)
> {
> stopifnot(require("sandwich"))
>
> ## cluster specification
> if(is.null(cluster)) cluster <- attr(object, "cluster")
> if(is.null(cluster)) stop("no 'cluster' specification
found")
> cluster <- factor(cluster)
>
> ## estimating functions and dimensions
> ef <- estfun(object)
> n <- NROW(ef)
> k <- NCOL(ef)
> if(n != length(cluster))
> stop("length of 'cluster' does not match number of
observations")
> m <- length(levels(cluster))
>
> ## aggregate estimating functions by cluster and compute meat
> ef <- sapply(levels(cluster), function(i) colSums(ef[cluster == i, ,
> drop = FALSE]))
> ef <- if(NCOL(ef) > 1L) t(ef) else matrix(ef, ncol = 1L)
> mt <- crossprod(ef)/n
>
> ## bread
> br <- try(bread(object), silent = TRUE)
> if(inherits(br, "try-error")) br <- vcov(object) * n
>
> ## put together sandwich
> vc <- 1/n * (br %*% mt %*% br)
>
> ## adjustment
> if(is.null(adjust)) adjust <- class(object)[1L] == "lm"
> adj <- if(adjust) m/(m - 1L) * (n - 1L)/(n - k) else m/(m - 1L)
>
> ## return
> return(adj * vc)
> }
>
>
>> Neither works for me. So I wonder what am I doing wrong here?
>>
>>
>> All suggestions are welcome ? thank you!
>> [[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.
[[alternative HTML version deleted]]