Prediction and simulation in Mixed-Effects Models

We recommend the MixedModelsSim.jl package and associated documentation for useful tools in constructing designs to simulate. For now, we'll use the sleep study data as a starting point.

using DataFrames
using MixedModels
using MixedModelsDatasets
using StatsBase
# use a DataFrame to make it easier to change things later
slp = DataFrame(MixedModelsDatasets.dataset(:sleepstudy))
slpm = fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj)), slp)

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + (1 | subj)
   logLik   -2 logLik     AIC       AICc        BIC    
  -897.0393  1794.0786  1802.0786  1802.3072  1814.8505

Variance components:
            Column    Variance Std.Dev.
subj     (Intercept)  1296.8699 36.0121
Residual               954.5278 30.8954
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  251.405     9.50618   26.45    <1e-99
days          10.4673    0.801735  13.06    <1e-38
──────────────────────────────────────────────────

Prediction

The simplest form of prediction are the fitted values from the model: they are indeed the model's predictions for the observed data.

predict(slpm) ≈ fitted(slpm)

true

When generalizing to new data, we need to consider what happens if there are new, previously unobserved levels of the grouping variable(s). MixedModels.jl provides three options:

1. :error: error on encountering unobserved levels
2. :population: use population values (i.e. only the fixed effects) for observations with unobserved levels
3. :missing: return missing for observations with unobserved levels.

Providing either no prediction (:error, :missing) or providing the population-level values seem to be the most reasonable ways for predicting new values. For simulating new values based on previous estimates of the variance components, use simulate.

In the case where there are no new levels of the grouping variable, all three of these methods provide the same results:

predict(slpm, slp; new_re_levels=:population) ≈ fitted(slpm)

true

predict(slpm, slp; new_re_levels=:missing) ≈ fitted(slpm)

true

predict(slpm, slp; new_re_levels=:error) ≈ fitted(slpm)

true

In the case where there are new levels of the grouping variable, these methods differ.

# create a new level
slp2 = transform(slp, :subj => ByRow(x -> (x == "S308" ? "NEW" : x)) => :subj)

180×3 DataFrame
 Row │ subj    days  reaction 
     │ String  Int8  Float32  
─────┼────────────────────────
   1 │ NEW        0   249.56
   2 │ NEW        1   258.705
   3 │ NEW        2   250.801
   4 │ NEW        3   321.44
   5 │ NEW        4   356.852
   6 │ NEW        5   414.69
   7 │ NEW        6   382.204
   8 │ NEW        7   290.149
  ⋮  │   ⋮      ⋮       ⋮
 174 │ S372       3   310.632
 175 │ S372       4   287.173
 176 │ S372       5   329.608
 177 │ S372       6   334.482
 178 │ S372       7   343.22
 179 │ S372       8   369.142
 180 │ S372       9   364.124
              165 rows omitted

try
  predict(slpm, slp2; new_re_levels=:error)
catch e
  show(e)
end

ArgumentError("New level encountered in subj")

predict(slpm, slp2; new_re_levels=:missing)

180-element Vector{Union{Missing, Float64}}:
    missing
    missing
    missing
    missing
    missing
    missing
    missing
    missing
    missing
    missing
   ⋮
 279.9221289700728
 290.3894144756822
 300.85669998129174
 311.32398548690117
 321.7912709925107
 332.25855649812013
 342.72584200372967
 353.1931275093391
 363.66041301494863

predict(slpm, slp2; new_re_levels=:population)

180-element Vector{Float64}:
 251.40510605320847
 261.8723915588179
 272.3396770644274
 282.8069625700369
 293.27424807564637
 303.74153358125585
 314.20881908686533
 324.6761045924748
 335.1433900980843
 345.6106756036938
   ⋮
 279.9221289700728
 290.3894144756822
 300.85669998129174
 311.32398548690117
 321.7912709925107
 332.25855649812013
 342.72584200372967
 353.1931275093391
 363.66041301494863

For generalized linear mixed models, there is an additional keyword argument to predict: type specifies whether the predictions are returned on the scale of the linear predictor (:linpred) or on the level of the response (:response) (i.e. the level at which the values were originally observed).

cbpp = DataFrame(MixedModelsDatasets.dataset(:cbpp))
cbpp.rate = cbpp.incid ./ cbpp.hsz
gm = fit(MixedModel, @formula(rate ~ 1 + period + (1|herd)), cbpp, Binomial(), wts=float(cbpp.hsz))
predict(gm, cbpp; type=:response) ≈ fitted(gm)

false

logit(x) = log(x / (1 - x))
predict(gm, cbpp; type=:linpred) ≈ logit.(fitted(gm))

false

Simulation

In contrast to predict, simulate and simulate! introduce randomness. This randomness occurs both at the level of the observation-level (residual) variance and at the level of the random effects, where new conditional modes are sampled based on the specified covariance parameter (θ; see Details of the parameter estimation), which defaults to the estimated value of the model. For reproducibility, we specify a pseudorandom generator here; if none is provided, the global PRNG is taken as the default.

The simplest example of simulate takes a fitted model and generates a new response vector based on the existing model matrices combined with noise.

using Random
ynew = simulate(MersenneTwister(42), slpm)

180-element Vector{Float64}:
 283.0086324277372
 296.92519287408305
 321.96087087300515
 322.3364895747736
 396.98425562676067
 317.1515220750256
 348.5132549850892
 378.2824424540446
 302.30014992830576
 425.39058801095723
   ⋮
 314.68558962508183
 350.0319301979975
 293.6498136375726
 375.4019457774969
 324.58236629946265
 311.10277355039045
 338.7622249672039
 345.3555959218867
 301.287184629893

The simulated response can also be placed in a pre-allocated vector:

ynew2 = zeros(nrow(slp))
simulate!(MersenneTwister(42), ynew2, slpm)
ynew2 ≈ ynew

true

Or even directly replace the previous response vector in a model, at which point the model must be refit to the new values:

slpm2 = deepcopy(slpm)
refit!(simulate!(MersenneTwister(42), slpm2))

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + (1 | subj)
   logLik   -2 logLik     AIC       AICc        BIC    
  -903.1987  1806.3974  1814.3974  1814.6259  1827.1692

Variance components:
            Column    Variance Std.Dev.
subj     (Intercept)  2420.4298 49.1979
Residual               964.3708 31.0543
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  263.59      12.3684    21.31    <1e-99
days           9.35638    0.805859  11.61    <1e-30
───────────────────────────────────────────────────

This inplace simulation actually forms the basis of parametricbootstrap.

Finally, we can also simulate the response from entirely new data.

df = DataFrame(days = repeat(1:10, outer=20), subj=repeat(1:20, inner=10))
df[!, :subj] = string.("S", lpad.(df.subj, 2, "0"))
df[!, :reaction] .= 0
df

200×3 DataFrame
 Row │ days   subj    reaction 
     │ Int64  String  Int64    
─────┼─────────────────────────
   1 │     1  S01            0
   2 │     2  S01            0
   3 │     3  S01            0
   4 │     4  S01            0
   5 │     5  S01            0
   6 │     6  S01            0
   7 │     7  S01            0
   8 │     8  S01            0
  ⋮  │   ⋮      ⋮        ⋮
 194 │     4  S20            0
 195 │     5  S20            0
 196 │     6  S20            0
 197 │     7  S20            0
 198 │     8  S20            0
 199 │     9  S20            0
 200 │    10  S20            0
               185 rows omitted

ysim = simulate(MersenneTwister(42), slpm, df)

200-element Vector{Float64}:
 262.9357776809343
 276.8523381272802
 301.8880161262023
 302.26363482797075
 376.9114008799578
 297.0786673282227
 328.4404002382863
 358.2095877072418
 282.2272951815029
 405.3177332641544
   ⋮
 336.4096599426081
 274.0198071455625
 232.79273085362502
 281.77385811683087
 296.8405259840217
 379.16601424258704
 305.03738166252464
 358.3844522027118
 381.7236757584809

Note that this is a convenience method for creating a new model and then using the parameters from the old model to call simulate on that model. In other words, this method incurs the cost of constructing a new model and then discarding it. If you could re-use that model (e.g., fitting that model as part of a simulation study), it often makes sense to do these steps to perform these steps explicitly and avoid the unnecessary construction and discarding of an intermediate model:

msim = LinearMixedModel(@formula(reaction ~ 1 + days + (1|subj)), df)
simulate!(MersenneTwister(42), msim; θ=slpm.θ, β=slpm.β, σ=slpm.σ)
response(msim) ≈ ysim

true

fit!(msim)

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + (1 | subj)
   logLik   -2 logLik     AIC       AICc        BIC    
  -996.0696  1992.1392  2000.1392  2000.3443  2013.3325

Variance components:
            Column    Variance Std.Dev.
subj     (Intercept)   663.5395 25.7593
Residual              1012.9884 31.8275
 Number of obs: 200; levels of grouping factors: 20

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  259.607      7.53746   34.44    <1e-99
days           9.46755    0.783538  12.08    <1e-32
───────────────────────────────────────────────────

For simulating from generalized linear mixed models, there is no type option because the observation-level always occurs at the level of the response and not of the linear predictor.