I recently read up on Fixed Effects vs. Random Effects and learned about the Correlated Random Effects model in Wooldridge’s Introductory Econometrics. The model looks as follows:

where \bar x_i is the time average of the explanatory variables for each i for which we have random effect. Then \beta is the same as the one obtained from fixed effects and we have a test of whether fixed or random effects should be used (through \gamma).

I was wondering if there is an easy way to implement this in MixedModels.jl? Doing this by hand does not scale well. Her is an example where the parameter of interest is `NDI`

and we want state and year effects accounted for:

```
using DataFrames, RDatasets, FixedEffectModels, MixedModels, Statistics
df = dataset("plm", "Cigar");
categorical!(df, :State);
categorical!(df, :Year);
transform!(groupby(df, :State), :NDI => mean => :NDI_stateavg);
transform!(groupby(df, :Year), :NDI => mean => :NDI_yearavg);
femod = reg(df, @formula(Sales ~ NDI + fe(State) + fe(Year)));
cremod = fit(MixedModel,
@formula(Sales ~ NDI + NDI_stateavg + NDI_yearavg + (1|State)+(1|Year)), df);
```

```
julia> cremod
Linear mixed model fit by maximum likelihood
Sales ~ 1 + NDI + NDI_stateavg + NDI_yearavg + (1 | State) + (1 | Year)
logLik -2 logLik AIC AICc BIC
-5652.7170 11305.4341 11319.4341 11319.5157 11356.0430
Variance components:
Column Variance Std.Dev.
State (Intercept) 587.15400 24.23126
Year (Intercept) 56.30887 7.50392
Residual 170.89733 13.07277
Number of obs: 1380; levels of grouping factors: 46, 30
Fixed-effects parameters:
─────────────────────────────────────────────────────────
Coef. Std. Error z Pr(>|z|)
─────────────────────────────────────────────────────────
(Intercept) 27.5711 24.9605 1.10 0.2693
NDI -0.00684385 0.000443253 -15.44 <1e-53
NDI_stateavg 0.0143558 0.00326266 4.40 <1e-4
NDI_yearavg 0.00529594 0.000541579 9.78 <1e-21
─────────────────────────────────────────────────────────
```

```
julia> femod
Fixed Effect Model
======================================================================
Number of obs: 1380 Degrees of freedom: 77
R2: 0.832 R2 Adjusted: 0.822
F Statistic: 238.031 p-value: 0.000
R2 within: 0.154 Iterations: 2
Converged: true
======================================================================
Estimate Std.Error t value Pr(>|t|) Lower 95% Upper 95%
----------------------------------------------------------------------
NDI -0.00684385 0.000443591 -15.4283 0.000 -0.00771408 -0.00597361
======================================================================
```

I can try to contribute this to the package if it doesn’t exist and it is interesting to others.