2

I am trying to estimate an ordinal logistic regression with clustered standard errors using the MASS package's polr() function. There is no built-in clustering feature, so I am looking for (a) packages or (b) manual methods for calculating clustered standard errors using the model output. I plan to use margins package to estimate marginal effects from the model.

Here is an example:

library(MASS)
set.seed(1)
obs <- 500

# Create data frame
dat <- data.frame(y = as.factor(round(rnorm(n = obs, mean = 5, sd = 1), 0)),
              x = sample(x = 1:obs, size = obs, replace = T),
              clust = rep(c(1,2), 250))

# Estimate and summarize model
m1 <- MASS::polr(y ~x, data = dat, Hess = TRUE)

summary(m1)

While many questions on Stack Overflow ask about how to cluster standard errors in R for ordinary least squares models (and in some cases for logistic regression), it's unclear how to cluster errors in ordered logistic regression (i.e. proportional odds logistic regression). Additionally, the existing SO questions focus on packages that have other severe drawbacks (e.g. the classes of model outputs are not compatible with other standard packages for analysis and presentation of results) rather than using MASS::polr() which is compatible with predict().

socialscientist
  • 3,759
  • 5
  • 23
  • 58

1 Answers1

1

This is essentially following an answer offered by Achim Zeleis on rhelp in 2016.

library(lmtest)
library("sandwich")

coeftest(m1, vcov=vcovCL(m1, factor(dat$clust) ))

t test of coefficients:

    Estimate Std. Error t value  Pr(>|t|)    
x 0.00093547 0.00023777  3.9343 9.543e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
IRTFM
  • 258,963
  • 21
  • 364
  • 487