On Fri, 1 Jul 2016, Faradj Koliev wrote:
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
This probably means that not all rows of "mydata" have been used for the
estimation of the model, e.g., due to missing values or something like
that. Hence, the mismatch seems to occur.
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]>
The variable "tets" is presumably one of the regressors in your "model"
and apparently it cannot be found when extracting the model scores. Maybe
you haven't stored the model frame and deleted the data.
Hard to say without a simple reproducible example.
Z
What can I do to fix it? What am I doing wrong now?
1 jul 2016 kl. 14:57 skrev Achim Zeileis
<achim.zeil...@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@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.
______________________________________________
R-help@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.
