Correlated Random Effects in MixedModels.jl

Specific Domains Statistics
1

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:

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.

2

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.

function prepcre!(data, fixedvars, randomvars)
    for f in fixedvars
        for r in randomvars
            transform!(groupby(data, r),
                f => mean => "$(f)_$(r)")                                                                                                                          
        end
    end
end