# Correlated Random Effects in MixedModels.jl
**URL:** https://discourse.julialang.org/t/correlated-random-effects-in-mixedmodels-jl/54612
**Category:** Statistics
**Created:** [February 4, 2021, 11:27am UTC](https://discourse.julialang.org/t/correlated-random-effects-in-mixedmodels-jl/54612 "2021-02-04T11:27:11Z")
**Posts on this page:** 2
**Page:** 1
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [February 4, 2021, 11:27am UTC](https://discourse.julialang.org/t/correlated-random-effects-in-mixedmodels-jl/54612/1 "2021-02-04T11:27:11Z")
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:
y\_{it} = \beta x\_{it} + \gamma \bar x\_i + r\_i + u\_{ir}
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:
```julia
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
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
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.
---
### Author:  [@danielw2904](https://discourse.julialang.org/u/danielw2904)
#### Post date: [February 10, 2021, 8:40am UTC](https://discourse.julialang.org/t/correlated-random-effects-in-mixedmodels-jl/54612/2 "2021-02-10T08:40:57Z")
I have used this small function so far and it seems to work fine. It creates the avgs. for each of the fixed effect variables wrt each of the random effect variables. It should probably better return the formula for the model. I’ll work on that when I have time.
```julia
function prepcre!(data, fixedvars, randomvars)
for f in fixedvars
for r in randomvars
transform!(groupby(data, r),
f => mean => "$(f)_$(r)")
end
end
end
```