I run an OLS regression with multiple control variables. How is it possible to include clustered standard errors for groups using LinearModelFit?
5
-
$\begingroup$ 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, $\endgroup$– Cameron MurrayOct 28, 2016 at 13:15
-
$\begingroup$ 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. $\endgroup$– JimBOct 28, 2016 at 15:14
-
$\begingroup$ @Jim Baldwin Number 1. in the below thread shows a better explanation of the statistics. stats.stackexchange.com/questions/49050/… $\endgroup$– YvesOct 28, 2016 at 15:30
-
$\begingroup$ 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? $\endgroup$– JimBOct 28, 2016 at 15:53
-
$\begingroup$ 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) $\endgroup$– YvesOct 28, 2016 at 19:09
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
2
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`}}*)
-
$\begingroup$ 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. $\endgroup$– YvesOct 29, 2016 at 12:20
-
$\begingroup$ @Yves if you could provide a pointer to documentation of the approach you seek you might get better solutions. $\endgroup$ Oct 29, 2016 at 15:50
$\begingroup$
$\endgroup$
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
answered Oct 28, 2016 at 19:32