title: R: Clustering standard errors in MASS::polr()
body: <p>I am trying to estimate an ordinal logistic regression with clustered standard errors using the <code>MASS</code> package's <code>polr()</code> function. There is no built-in clustering feature, so I am looking for (a) packages or (b) manual methods for calculating clustered standard errors <strong>using the model output</strong>. I plan to use <code>margins</code> package to estimate marginal effects from the model.</p>
  
  <p>Here is an example:</p>
  
  <pre><code>library(MASS)
  set.seed(1)
  obs &lt;- 500
  
  # Create data frame
  dat &lt;- 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 &lt;- MASS::polr(y ~x, data = dat, Hess = TRUE)
  
  summary(m1)
  </code></pre>
  
  <p>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 <a href="https://stackoverflow.com/questions/16498849/logistic-regression-with-robust-clustered-standard-errors-in-r">logistic regression</a>), it's unclear how to cluster errors in <strong>ordered</strong> 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 <code>MASS::polr()</code> which is compatible with <code>predict()</code>.</p>