Clustered standard errors in LinearModelFit

I run an OLS regression with multiple control variables. How is it possible to include clustered standard errors for groups using LinearModelFit?

• I believe it is not an option. But would love it if someone came up with a neat function to do this, taking a column from the data representing the groups, – Cameron Murray Oct 28 '16 at 13:15
• Is "clustered variance" jargon in some non-statistical field? Is this a Stata term? (Stata is a statistical software package.) My point is that a giving a more complete definition so that one determine how to efficiently calculate that statistic would get you responses from a larger group of folks. – JimB Oct 28 '16 at 15:14
• @Jim Baldwin Number 1. in the below thread shows a better explanation of the statistics. stats.stackexchange.com/questions/49050/… – Yves Oct 28 '16 at 15:30
• Sorry, but I found no clarity from that link. Are you wanting to fit "mixed models" (models with 2 or more random terms) or wanting to allow for different variances for a grouping factor? – JimB Oct 28 '16 at 15:53
• I'd like to allow for different standard errors for a grouping factor to account for clustered errors (errors being independent across clusters but correlated within clusters) – Yves Oct 28 '16 at 19:09

Here is an alternative answer. @george2079 's answer finds estimates after the fit. The problem with that is that the statement model = LinearModelFit[all, x, x] assumes a constant variance which is contrary to the assumption of two different variances. To allow for two different variances currently one needs to use the LogLikelihood function. I've shamelessly borrowed from @george2079's answer:

truth[x_] = 2 x + 1;
SeedRandom[123];
cluster1 = {#, truth[#] + RandomVariate[NormalDistribution[0, 1]]} & /@ RandomReal[{0, 5}, 20]
cluster2 = {#, truth[#] + RandomVariate[NormalDistribution[0, 2]]} & /@ RandomReal[{0, 5}, 20]

(* Log of the likelihood *)
logL =
LogLikelihood[NormalDistribution[0, σ1], cluster1[[All, 2]] - (a + b cluster1[[All, 1]])] +
LogLikelihood[NormalDistribution[0, σ2], cluster2[[All, 2]] - (a + b cluster2[[All, 1]])];

(* Maximum likelihood estimates *)
sol = FindMaximum[{logL, σ1 > 0 && σ2 > 0}, {a, b, σ1, σ2}]
(* {-68.43789185435824,{a -> 0.2845791362978213,
b -> 2.1193090163306247,σ1 -> 0.9514019056833594,σ2 -> 1.8848273579627017}}*)

• thank you. Problem with this solution is that it seems to get rather inconvenient for the case of 10 explanatory variables and 20 clusters for example. Unfortunately there's no simple solution like in SAS or Stata. – Yves Oct 29 '16 at 12:20
• @Yves if you could provide a pointer to documentation of the approach you seek you might get better solutions. – george2079 Oct 29 '16 at 15:50

something like this?

truth[x_] = 2 x + 1;
cluster1 = {#, truth[#] + RandomVariate[NormalDistribution[0, 1]]} & /@
RandomReal[{0, 5}, 20];
cluster2 = {#, truth[#] + RandomVariate[NormalDistribution[0, 2]]} & /@
RandomReal[{0, 5}, 20];
all = cluster1~Join~cluster2;
model = LinearModelFit[all, x, x]
Show[{Plot[model[x], {x, 0, 5}], ListPlot[{cluster1, cluster2}]}]


variance calculation:

Sqrt[(model@#[[1]] - #[[2]] & /@ cluster1)^2 // Mean]
Sqrt[(model@#[[1]] - #[[2]] & /@ cluster2)^2 // Mean]
`

1.11905

2.34978