Reproducing a simple clustered-robust standard error OLS with Stan

[martina...@gmail.com]

Dear all,

I am trying to reproduce a simple OLS model fitted with clustered-robust standard errors with Stan.

In R, my model is the following:

mymodel <- (lm(Y ~ X1 + X2 + I(X1*X2), data = mydata));

In order to calculate the clustered-robust SE I use the following code:

clusteredSE <- function(fm, dfcw, cluster)

{

      library(sandwich)

      library(lmtest)

      library(zoo)

      M <- length(unique(cluster))

      N <- length(cluster)

      dfc <- (M / (M - 1)) * ((N - 1) / (N - fm$rank))

      u <- apply(estfun(fm), 2, function(x) tapply(x, cluster, sum))

      vcovCL <- dfc * sandwich(fm, meat = crossprod(u) / N) * dfcw

      coeftest(fm, vcovCL)

}

clusteredSE(mymodel, 1, mydata$countries)

Now, in Stan, I understood how to specify both simple OLS models and also simple Hierarchical models. However, I did not understand nor find an example of how should I proceed for estimating what in the frequentist world would be a clustered-robust standard error OLS. Could any of you point me in the right direction?

Thanks in advance for your time.

Martina Z. Andrews

[Rachael Meager]

Hi Martina,

I think folks here would suggest you use a hierarchical model to create a cluster-specific random effects structure rather than using the frequentist cluster-robust SE formulas. That's definitely what I would do in your position.

If you do want to implement the cluster SE formula, I suspect you can explicitly parameterise the variance-covariance matrix (or the sum of its elements, which is what the formula actually specifies, if I recall correctly) to take the form that you want, but it might take a bit of work to code up in stan.

I am sorry I can't be more specific -- I have never implemented this myself so I am just giving you generalised advice at this point!

Cheers

Rachael  

[martina...@gmail.com]

Hi Rachael,

Many thanks for dropping an answer. What exactly do you mean by "use a hierarchical model to create a cluster-specific random effects structure"? Do you mean just moving towards the use of a random-intercept and random slope model (for all slopes)? So, actually my current interest lays precisely in estimating both the hierarchical model and the cluster-robust SE within the Bayesian framework. The hierarchical model I know how to specify in Stan.

Best,

Martina

[CoreySparks]

I believe what Rachael is referring to is to include a cluster-specific random intercept term, in addition to any other random effect terms in your model. By including a random intercept for the clustering variables, you are effectively controlling for within cluster homogeneity.

Cheers,

Corey

[Bob Carpenter]

I don't know what any of those models are because I don't know
what all those R functions do or even what "clustered-robust standard
errors" are.

We typically do Bayesian models and only compute standard errors on
our posterior mean estimates. In general, we do full Bayes (event
estimation, probability estimation, deicison theory) rather than
hypothesis testing.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Rachael Meager]

Corey's right! Include a cluster-specific random effect for the intercept. You can also do cluster-specific REs on the slopes. If your variables only vary at the cluster level the slopes correction may be harder but even then I think you might be able to do it as long as you constrain the functional form that governs the REs. 

Cluster-Robust SEs actually do less controlling for intracluster correlation than cluster REs do. Once you have the cluster REs I don't think the cluster-robust SE makes any sense, there should be no residual correlation to pick up...

Like Bob said, Stan is really for fitting Bayesian models, and cluster-robust SEs are a specific kind of frequentist corrective to standard errors on a frequentist estimate. I don't think smushing these things together will produce anything coherent, let alone useful. It is an interest of mine to understand the relationship between these kinds of frequentist SE correctives and the bayesian modelling tools we use to combat the same problems, but there is not much work on it yet. If I were you I would first try to understand the relationship between these things before smushing them together, I suspect you'll find that it won't make sense to combine them, but I am willing to be shown a proof to the contrary. :) 

Cheers

Rachael