How does one get multiway clustered standard errors in R for plm objects, where the clustering is not at the level of the panel's time/group IDs?
The package plm provides support to calculate cluster-robust standard errors using the function plm::vcovHC. Unfortunately, this function only supports clustering at the group or time IDs of the panel. In certain cases one would want to cluster standard errors at a different level than the panel's unit of observation. An example is a regression with individual fixed effects where variables at a higher level of aggregation are used as independent variables.
A similar question was asked in 2014 and a bootstrap function was recommended. This question differs because (1) I would like to cluster on two variables, not just one and (2) I would prefer not to use the bootstrap.
The package multiwaycov provides something very close to what I want, but unfortunately does not support plm objects.
