Model constructors

The LinearMixedModel type represents a linear mixed-effects model. Typically it is constructed from a Formula and an appropriate data type, usually a DataFrame.

MixedModels.LinearMixedModel — Type.

LinearMixedModel

Linear mixed-effects model representation

Fields

formula: the formula for the model
allterms: a vector of random-effects terms, the fixed-effects terms and the response
sqrtwts: vector of square roots of the case weights. Can be empty.
A: an nt × nt symmetric BlockMatrix of matrices representing hcat(Z,X,y)'hcat(Z,X,y)
L: a nt × nt BlockMatrix - the lower Cholesky factor of Λ'AΛ+I
optsum: an OptSummary object

Properties

θ or theta: the covariance parameter vector used to form λ
β or beta: the fixed-effects coefficient vector
λ or lambda: a vector of lower triangular matrices repeated on the diagonal blocks of Λ
σ or sigma: current value of the standard deviation of the per-observation noise
b: random effects on the original scale, as a vector of matrices
reterms: a Vector{ReMat{T}} of random-effects terms.
feterms: a Vector{FeMat{T}} of the fixed-effects model matrix and the response
u: random effects on the orthogonal scale, as a vector of matrices
lowerbd: lower bounds on the elements of θ
X: the fixed-effects model matrix
y: the response vector

source

Examples of linear mixed-effects model fits

For illustration, several data sets from the lme4 package for R are made available in .rda format in this package. These include the Dyestuff and Dyestuff2 data sets.

julia> using DataFrames, MixedModels, RData, StatsBase

julia> datf = joinpath(dirname(pathof(MixedModels)),"..","test","dat.rda")
"/home/travis/build/JuliaStats/MixedModels.jl/src/../test/dat.rda"

julia> const dat = Dict(Symbol(k)=>v for (k,v) in load(datf))
Dict{Symbol,DataFrames.DataFrame} with 62 entries:
  :bs10          => 1104×6 DataFrame…
  :Genetics      => 60×5 DataFrame…
  :Contraception => 1934×6 DataFrame…
  :Mmmec         => 354×6 DataFrame…
  :kb07          => 1790×10 DataFrame. Omitted printing of 3 columns…
  :Rail          => 18×2 DataFrame…
  :KKL           => 53765×24 DataFrame. Omitted printing of 16 columns…
  :Bond          => 21×3 DataFrame…
  :VerbAgg       => 7584×9 DataFrame. Omitted printing of 1 columns…
  :ml1m          => 1000209×3 DataFrame…
  :ergoStool     => 36×3 DataFrame…
  :s3bbx         => 2449×6 DataFrame…
  :cake          => 270×5 DataFrame…
  :Cultivation   => 24×4 DataFrame…
  :Pastes        => 60×4 DataFrame…
  :Exam          => 4059×5 DataFrame…
  :Socatt        => 1056×9 DataFrame. Omitted printing of 2 columns…
  :WWheat        => 60×3 DataFrame…
  :Pixel         => 102×5 DataFrame…
  ⋮              => ⋮

julia> describe(dat[:Dyestuff])
2×8 DataFrame. Omitted printing of 1 columns
│ Row │ variable │ mean   │ min    │ median │ max    │ nunique │ nmissing │
│     │ Symbol   │ Union… │ Any    │ Union… │ Any    │ Union…  │ Nothing  │
├─────┼──────────┼────────┼────────┼────────┼────────┼─────────┼──────────┤
│ 1   │ G        │        │ A      │        │ F      │ 6       │          │
│ 2   │ Y        │ 1527.5 │ 1440.0 │ 1530.0 │ 1635.0 │         │          │

The columns in these data sets have been renamed for convenience. The response is always named Y. Potential grouping factors for random-effects terms are named G, H, etc. Numeric covariates are named starting with U. Categorical covariates not suitable as grouping factors are named starting with A.

Models with simple, scalar random effects

The formula language in Julia is similar to that in R except that the formula must be enclosed in a call to the @formula macro. A basic model with simple, scalar random effects for the levels of G (the batch of an intermediate product, in this case) is declared and fit as

julia> fm1 = fit!(LinearMixedModel(@formula(Y ~ 1 + (1|G)),
           dat[:Dyestuff]))
Linear mixed model fit by maximum likelihood
 Y ~ 1 + (1 | G)
   logLik   -2 logLik     AIC        BIC    
 -163.66353  327.32706  333.32706  337.53065

Variance components:
            Column    Variance  Std.Dev. 
G        (Intercept)  1388.3333 37.260345
Residual              2451.2500 49.510100
 Number of obs: 30; levels of grouping factors: 6

  Fixed-effects parameters:
────────────────────────────────────────────────
              Coef.  Std. Error      z  Pr(>|z|)
────────────────────────────────────────────────
(Intercept)  1527.5     17.6946  86.33    <1e-99
────────────────────────────────────────────────

An alternative expression is

julia> fm1 = fit(MixedModel, @formula(Y ~ 1 + (1|G)), dat[:Dyestuff])
Linear mixed model fit by maximum likelihood
 Y ~ 1 + (1 | G)
   logLik   -2 logLik     AIC        BIC    
 -163.66353  327.32706  333.32706  337.53065

Variance components:
            Column    Variance  Std.Dev. 
G        (Intercept)  1388.3333 37.260345
Residual              2451.2500 49.510100
 Number of obs: 30; levels of grouping factors: 6

  Fixed-effects parameters:
────────────────────────────────────────────────
              Coef.  Std. Error      z  Pr(>|z|)
────────────────────────────────────────────────
(Intercept)  1527.5     17.6946  86.33    <1e-99
────────────────────────────────────────────────

(If you are new to Julia you may find that this first fit takes an unexpectedly long time, due to Just-In-Time (JIT) compilation of the code. The second and subsequent calls to such functions are much faster.)

julia> @time fit(MixedModel, @formula(Y ~ 1 + (1|G)), dat[:Dyestuff2]);
  0.648289 seconds (1.01 M allocations: 51.987 MiB)

By default, the model fit is by maximum likelihood. To use the REML criterion instead, add the optional named argument REML = true to the call to fit!

julia> fm1R = fit(MixedModel, @formula(Y ~ 1 + (1|G)),
           dat[:Dyestuff], REML=true)
Linear mixed model fit by REML
 Y ~ 1 + (1 | G)
 REML criterion at convergence: 319.6542768422538

Variance components:
            Column    Variance  Std.Dev. 
G        (Intercept)  1764.0506 42.000602
Residual              2451.2499 49.510099
 Number of obs: 30; levels of grouping factors: 6

  Fixed-effects parameters:
────────────────────────────────────────────────
              Coef.  Std. Error      z  Pr(>|z|)
────────────────────────────────────────────────
(Intercept)  1527.5     19.3834  78.80    <1e-99
────────────────────────────────────────────────

Simple, scalar random effects

A simple, scalar random effects term in a mixed-effects model formula is of the form (1|G). All random effects terms end with |G where G is the grouping factor for the random effect. The name or, more generally, the expression G should evaluate to a categorical array that has a distinct set of levels. The random effects are associated with the levels of the grouping factor.

A scalar random effect is, as the name implies, one scalar value for each level of the grouping factor. A simple, scalar random effects term is of the form, (1|G). It corresponds to a shift in the intercept for each level of the grouping factor.

Models with vector-valued random effects

The sleepstudy data are observations of reaction time, Y, on several subjects, G, after 0 to 9 days of sleep deprivation, U. A model with random intercepts and random slopes for each subject, allowing for within-subject correlation of the slope and intercept, is fit as

julia> fm2 = fit(MixedModel, @formula(Y ~ 1+U + (1+U|G)), dat[:sleepstudy])
Linear mixed model fit by maximum likelihood
 Y ~ 1 + U + (1 + U | G)
   logLik   -2 logLik     AIC        BIC    
 -875.96967 1751.93934 1763.93934 1783.09709

Variance components:
            Column    Variance   Std.Dev.    Corr.
G        (Intercept)  565.510678 23.7804684
         U             32.682124  5.7168282  0.08
Residual              654.941447 25.5918238
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  251.405      6.63226  37.91    <1e-99
U             10.4673     1.50224   6.97    <1e-11
──────────────────────────────────────────────────

Models with multiple, scalar random-effects terms

A model for the Penicillin data incorporates random effects for the plate, G, and for the sample, H. As every sample is used on every plate these two factors are crossed.

julia> fm4 = fit(MixedModel,@formula(Y ~ 1+(1|G)+(1|H)),dat[:Penicillin])
Linear mixed model fit by maximum likelihood
 Y ~ 1 + (1 | G) + (1 | H)
   logLik   -2 logLik     AIC        BIC    
 -166.09417  332.18835  340.18835  352.06760

Variance components:
            Column    Variance   Std.Dev. 
G        (Intercept)  0.71497950 0.8455646
H        (Intercept)  3.13519287 1.7706476
Residual              0.30242640 0.5499331
 Number of obs: 144; levels of grouping factors: 24, 6

  Fixed-effects parameters:
─────────────────────────────────────────────────
               Coef.  Std. Error      z  Pr(>|z|)
─────────────────────────────────────────────────
(Intercept)  22.9722    0.744596  30.85    <1e-99
─────────────────────────────────────────────────

In contrast the sample, G, grouping factor is nested within the batch, H, grouping factor in the Pastes data. That is, each level of G occurs in conjunction with only one level of H.

julia> fm5 = fit(MixedModel, @formula(Y ~ 1 + (1|G) + (1|H)), dat[:Pastes])
Linear mixed model fit by maximum likelihood
 Y ~ 1 + (1 | G) + (1 | H)
   logLik   -2 logLik     AIC        BIC    
 -123.99723  247.99447  255.99447  264.37184

Variance components:
            Column    Variance   Std.Dev.  
G        (Intercept)  8.43361634 2.90406893
H        (Intercept)  1.19918042 1.09507097
Residual              0.67800208 0.82340882
 Number of obs: 60; levels of grouping factors: 30, 10

  Fixed-effects parameters:
─────────────────────────────────────────────────
               Coef.  Std. Error      z  Pr(>|z|)
─────────────────────────────────────────────────
(Intercept)  60.0533    0.642136  93.52    <1e-99
─────────────────────────────────────────────────

In observational studies it is common to encounter partially crossed grouping factors. For example, the InstEval data are course evaluations by students, G, of instructors, H. Additional covariates include the academic department, I, in which the course was given and A, whether or not it was a service course.

julia> fm6 = fit(MixedModel,@formula(Y ~ 1+ A*I +(1|G)+(1|H)),dat[:InstEval])
Linear mixed model fit by maximum likelihood
 Y ~ 1 + A + I + A & I + (1 | G) + (1 | H)
     logLik        -2 logLik          AIC             BIC       
 -1.18792777×10⁵  2.37585553×10⁵  2.37647553×10⁵  2.37932876×10⁵

Variance components:
            Column    Variance   Std.Dev.  
G        (Intercept)  0.10541797 0.32468133
H        (Intercept)  0.25841639 0.50834672
Residual              1.38472777 1.17674457
 Number of obs: 73421; levels of grouping factors: 2972, 1128

  Fixed-effects parameters:
──────────────────────────────────────────────────────
                    Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────────
(Intercept)    3.22961      0.064053   50.42    <1e-99
A: 1           0.252025     0.0686507   3.67    0.0002
I: 5           0.129536     0.101294    1.28    0.2010
I: 10         -0.176751     0.0881352  -2.01    0.0449
I: 12          0.0517102    0.0817524   0.63    0.5270
I: 6           0.0347319    0.085621    0.41    0.6850
I: 7           0.14594      0.0997984   1.46    0.1436
I: 4           0.151689     0.0816897   1.86    0.0633
I: 8           0.104206     0.118751    0.88    0.3802
I: 9           0.0440401    0.0962985   0.46    0.6474
I: 14          0.0517546    0.0986029   0.52    0.5997
I: 1           0.0466719    0.101942    0.46    0.6471
I: 3           0.0563461    0.0977925   0.58    0.5645
I: 11          0.0596536    0.100233    0.60    0.5517
I: 2           0.00556281   0.110867    0.05    0.9600
A: 1 & I: 5   -0.180757     0.123179   -1.47    0.1423
A: 1 & I: 10   0.0186492    0.110017    0.17    0.8654
A: 1 & I: 12  -0.282269     0.0792937  -3.56    0.0004
A: 1 & I: 6   -0.494464     0.0790278  -6.26    <1e-9
A: 1 & I: 7   -0.392054     0.110313   -3.55    0.0004
A: 1 & I: 4   -0.278547     0.0823727  -3.38    0.0007
A: 1 & I: 8   -0.189526     0.111449   -1.70    0.0890
A: 1 & I: 9   -0.499868     0.0885423  -5.65    <1e-7
A: 1 & I: 14  -0.497162     0.0917162  -5.42    <1e-7
A: 1 & I: 1   -0.24042      0.0982071  -2.45    0.0144
A: 1 & I: 3   -0.223013     0.0890548  -2.50    0.0123
A: 1 & I: 11  -0.516997     0.0809077  -6.39    <1e-9
A: 1 & I: 2   -0.384773     0.091843   -4.19    <1e-4
──────────────────────────────────────────────────────

Simplifying the random effect correlation structure

MixedEffects.jl estimates not only the variance of the effects for each random effect level, but also the correlation between the random effects for different predictors. So, for the model of the sleepstudy data above, one of the parameters that is estimated is the correlation between each subject's random intercept (i.e., their baseline reaction time) and slope (i.e., their particular change in reaction time over days of sleep deprivation). In some cases, you may wish to simplify the random effects structure by removing these correlation parameters. This often arises when there are many random effects you want to estimate (as is common in psychological experiments with many conditions and covariates), since the number of random effects parameters increases as the square of the number of predictors, making these models difficult to estimate from limited data.

A model with uncorrelated random effects for the intercept and slope by subject is fit as

julia> fm3 = fit!(zerocorr!(LinearMixedModel(@formula(Y ~ 1+U+(1+U|G)),
           dat[:sleepstudy])))
Linear mixed model fit by maximum likelihood
 Y ~ 1 + U + (1 + U | G)
   logLik   -2 logLik     AIC        BIC    
 -876.00163 1752.00326 1762.00326 1777.96804

Variance components:
            Column    Variance  Std.Dev.   Corr.
G        (Intercept)  584.258971 24.17145
         U             33.632805  5.79938   .  
Residual              653.115782 25.55613
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  251.405      6.70771  37.48    <1e-99
U             10.4673     1.51931   6.89    <1e-11
──────────────────────────────────────────────────

Note that the use of zerocorr! requires the model to be constructed, then altered to eliminate the correlation of the random effects, then fit with a call to the mutating function, fit!.

MixedModels.zerocorr! — Function.

zerocorr!(m::LinearMixedModel[, trmnms::Vector{Symbol}])

Rewrite the random effects specification for the grouping factors in trmnms to zero correlation parameter.

The default for trmnms is all the names of random-effects terms.

A random effects term is in the zero correlation parameter configuration when the off-diagonal elements of λ are all zero - hence there are no correlation parameters in that term being estimated.

source

The special syntax zerocorr can be applied to individual random effects terms inside the @formula:

MixedModels.zerocorr — Function.

zerocorr(term::RandomEffectsTerm)

Remove correlations between random effects in term.

source

julia> fit(MixedModel, @formula(Y ~ 1 + U + zerocorr(1+U|G)), dat[:sleepstudy])
Linear mixed model fit by maximum likelihood
 Y ~ 1 + U + MixedModels.ZeroCorr((1 + U | G))
   logLik   -2 logLik     AIC        BIC    
 -876.00163 1752.00326 1762.00326 1777.96804

Variance components:
            Column    Variance  Std.Dev.   Corr.
G        (Intercept)  584.258971 24.17145
         U             33.632805  5.79938   .  
Residual              653.115782 25.55613
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  251.405      6.70771  37.48    <1e-99
U             10.4673     1.51931   6.89    <1e-11
──────────────────────────────────────────────────

Alternatively, correlations between parameters can be removed by including them as separate random effects terms:

julia> fit(MixedModel, @formula(Y ~ 1+U+(1|G)+(0+U|G)), dat[:sleepstudy])
Linear mixed model fit by maximum likelihood
 Y ~ 1 + U + (1 | G) + (0 + U | G)
   logLik   -2 logLik     AIC        BIC    
 -876.00163 1752.00326 1762.00326 1777.96804

Variance components:
            Column    Variance  Std.Dev.   Corr.
G        (Intercept)  584.258971 24.17145
         U             33.632805  5.79938   .  
Residual              653.115782 25.55613
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  251.405      6.70771  37.48    <1e-99
U             10.4673     1.51931   6.89    <1e-11
──────────────────────────────────────────────────

Note that it is necessary to explicitly block the inclusion of an intercept term by adding 0 in the random-effects term (0+U|G).

Finally, for predictors that are categorical, MixedModels.jl will estimate correlations between each level. Notice the large number of correlation parameters if we treat U as a categorical variable by giving it contrasts:

julia> fit(MixedModel, @formula(Y ~ 1+U+(1+U|G)), dat[:sleepstudy],
           contrasts=Dict(:U=>DummyCoding()))
Linear mixed model fit by maximum likelihood
 Y ~ 1 + U + (1 + U | G)
   logLik   -2 logLik     AIC        BIC    
  -805.3996  1610.7992  1742.7992  1953.5344

Variance components:
            Column     Variance   Std.Dev.   Corr.
G        (Intercept)   956.231148 30.922987
         U: 1.0        496.640794 22.285439 -0.30
         U: 2.0        914.755039 30.244918 -0.57  0.75
         U: 3.0       1264.283956 35.556771 -0.37  0.72  0.87
         U: 4.0       1480.122126 38.472355 -0.32  0.58  0.67  0.91
         U: 5.0       2288.746783 47.840848 -0.25  0.46  0.45  0.70  0.85
         U: 6.0       3842.541380 61.988236 -0.27  0.30  0.48  0.70  0.77  0.75
         U: 7.0       1805.494419 42.491110 -0.16  0.22  0.47  0.50  0.62  0.64  0.71
         U: 8.0       3147.273280 56.100564 -0.20  0.28  0.36  0.56  0.73  0.90  0.73  0.74
         U: 9.0       3068.800196 55.396753  0.05  0.25  0.16  0.38  0.58  0.78  0.38  0.53  0.85
Residual                18.991546  4.357929
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  256.652       7.36064  34.87    <1e-99
U: 1.0         7.84395     5.44989   1.44    0.1501
U: 2.0         8.71009     7.27529   1.20    0.2312
U: 3.0        26.3402      8.50577   3.10    0.0020
U: 4.0        31.9976      9.18364   3.48    0.0005
U: 5.0        51.8667     11.3694    4.56    <1e-5
U: 6.0        55.5264     14.6828    3.78    0.0002
U: 7.0        62.0988     10.1201    6.14    <1e-9
U: 8.0        79.9777     13.3026    6.01    <1e-8
U: 9.0        94.1994     13.1377    7.17    <1e-12
───────────────────────────────────────────────────

Separating the 1 and U random effects into separate terms removes the correlations between the intercept and the levels of U, but not between the levels themselves:

julia> fit(MixedModel, @formula(Y ~ 1+U+(1|G)+(0+U|G)), dat[:sleepstudy],
           contrasts=Dict(:U=>DummyCoding()))
Linear mixed model fit by maximum likelihood
 Y ~ 1 + U + (1 | G) + (0 + U | G)
   logLik   -2 logLik     AIC        BIC    
  -805.3993  1610.7986  1744.7986  1958.7267

Variance components:
            Column     Variance    Std.Dev.    Corr.
G        (Intercept)   254.8329728 15.9634887
         U: 0.0        715.0961607 26.7412820   .  
         U: 1.0        796.0939722 28.2151373   .    0.65
         U: 2.0        561.6178331 23.6984774   .    0.26  0.69
         U: 3.0       1167.4873418 34.1685139   .    0.32  0.68  0.85
         U: 4.0       1449.4232150 38.0712912   .    0.32  0.58  0.63  0.90
         U: 5.0       2270.1876345 47.6464861   .    0.26  0.45  0.40  0.69  0.84
         U: 6.0       3508.7575522 59.2347664   .    0.11  0.22  0.39  0.65  0.73  0.72
         U: 7.0       2106.4678189 45.8962724   .    0.40  0.38  0.52  0.54  0.65  0.65  0.68
         U: 8.0       3160.4372900 56.2177667   .    0.23  0.31  0.32  0.55  0.72  0.89  0.71  0.74
         U: 9.0       3973.1028107 63.0325536   .    0.47  0.51  0.36  0.53  0.70  0.83  0.42  0.63  0.87
Residual                 4.5694002  2.1376155
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  256.652       7.35791  34.88    <1e-99
U: 1.0         7.84395     5.45552   1.44    0.1505
U: 2.0         8.71009     7.28243   1.20    0.2317
U: 3.0        26.3402      8.52339   3.09    0.0020
U: 4.0        31.9976      9.20582   3.48    0.0005
U: 5.0        51.8667     11.3981    4.55    <1e-5
U: 6.0        55.5265     14.7068    3.78    0.0002
U: 7.0        62.0988     10.1208    6.14    <1e-9
U: 8.0        79.9777     13.3242    6.00    <1e-8
U: 9.0        94.1994     13.1547    7.16    <1e-12
───────────────────────────────────────────────────

But using zerocorr on the individual terms (or zerocorr! on the constructed model object as above) does remove the correlations between the levels:

julia> fit(MixedModel, @formula(Y ~ 1+U+zerocorr(1+U|G)), dat[:sleepstudy],
           contrasts=Dict(:U=>DummyCoding()))
Linear mixed model fit by maximum likelihood
 Y ~ 1 + U + MixedModels.ZeroCorr((1 + U | G))
   logLik   -2 logLik     AIC        BIC    
  -882.9138  1765.8276  1807.8276  1874.8797

Variance components:
            Column    Variance   Std.Dev.   Corr.
G        (Intercept)   958.51186 30.959843
         U: 1.0          0.00000  0.000000   .  
         U: 2.0          0.00000  0.000000   .     .  
         U: 3.0          0.00000  0.000000   .     .     .  
         U: 4.0          0.00000  0.000000   .     .     .     .  
         U: 5.0        519.82399 22.799649   .     .     .     .     .  
         U: 6.0       1703.51659 41.273679   .     .     .     .     .     .  
         U: 7.0        609.09627 24.679876   .     .     .     .     .     .     .  
         U: 8.0       1273.52942 35.686544   .     .     .     .     .     .     .     .  
         U: 9.0       1753.95967 41.880302   .     .     .     .     .     .     .     .     .  
Residual               434.82413 20.852437
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  256.652       8.79816  29.17    <1e-99
U: 1.0         7.84395     6.95081   1.13    0.2591
U: 2.0         8.71009     6.95081   1.25    0.2102
U: 3.0        26.3402      6.95081   3.79    0.0002
U: 4.0        31.9976      6.95081   4.60    <1e-5
U: 5.0        51.8666      8.78595   5.90    <1e-8
U: 6.0        55.5264     11.9563    4.64    <1e-5
U: 7.0        62.0988      9.0638    6.85    <1e-11
U: 8.0        79.9777     10.9117    7.33    <1e-12
U: 9.0        94.1994     12.0729    7.80    <1e-14
───────────────────────────────────────────────────

julia> fit(MixedModel, @formula(Y ~ 1+U+(1|G)+zerocorr(0+U|G)),dat[:sleepstudy],
           contrasts=Dict(:U=>DummyCoding()))
Linear mixed model fit by maximum likelihood
 Y ~ 1 + U + (1 | G) + MixedModels.ZeroCorr((0 + U | G))
   logLik   -2 logLik     AIC        BIC    
  -878.9843  1757.9686  1801.9686  1872.2137

Variance components:
            Column      Variance    Std.Dev.    Corr.
G        (Intercept)  1135.34423225 33.6948695
         U: 0.0        776.23600556 27.8610123   .  
         U: 1.0        358.01572563 18.9213035   .     .  
         U: 2.0        221.49645908 14.8827571   .     .     .  
         U: 3.0          0.38236102  0.6183535   .     .     .     .  
         U: 4.0         44.71468158  6.6869037   .     .     .     .     .  
         U: 5.0        670.50512881 25.8941138   .     .     .     .     .     .  
         U: 6.0       1740.07339494 41.7141870   .     .     .     .     .     .     .  
         U: 7.0        908.93184419 30.1484965   .     .     .     .     .     .     .     .  
         U: 8.0       1458.06668413 38.1846394   .     .     .     .     .     .     .     .     .  
         U: 9.0       2028.19601076 45.0354972   .     .     .     .     .     .     .     .     .     .  
Residual               180.63063839 13.4398898
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  256.652      10.7812   23.81    <1e-99
U: 1.0         7.84395     9.11505   0.86    0.3895
U: 2.0         8.71009     8.68905   1.00    0.3161
U: 3.0        26.3402      7.95082   3.31    0.0009
U: 4.0        31.9976      8.10422   3.95    <1e-4
U: 5.0        51.8666     10.0222    5.18    <1e-6
U: 6.0        55.5264     12.6438    4.39    <1e-4
U: 7.0        62.0988     10.6626    5.82    <1e-8
U: 8.0        79.9777     12.0082    6.66    <1e-10
U: 9.0        94.1994     13.2617    7.10    <1e-11
───────────────────────────────────────────────────

Fitting generalized linear mixed models

To create a GLMM representation

MixedModels.GeneralizedLinearMixedModel — Type.

GeneralizedLinearMixedModel

Generalized linear mixed-effects model representation

Fields

LMM: a LinearMixedModel - the local approximation to the GLMM.
β: the pivoted and possibly truncated fixed-effects vector
β₀: similar to β. Used in the PIRLS algorithm if step-halving is needed.
θ: covariance parameter vector
b: similar to u, equivalent to broadcast!(*, b, LMM.Λ, u)
u: a vector of matrices of random effects
u₀: similar to u. Used in the PIRLS algorithm if step-halving is needed.
resp: a GlmResp object
η: the linear predictor
wt: vector of prior case weights, a value of T[] indicates equal weights.

The following fields are used in adaptive Gauss-Hermite quadrature, which applies only to models with a single random-effects term, in which case their lengths are the number of levels in the grouping factor for that term. Otherwise they are zero-length vectors.

devc: vector of deviance components
devc0: vector of deviance components at offset of zero
sd: approximate standard deviation of the conditional density
mult: multiplier

Properties

In addition to the fieldnames, the following names are also accessible through the . extractor

theta: synonym for θ
beta: synonym for β
σ or sigma: common scale parameter (value is NaN for distributions without a scale parameter)
lowerbd: vector of lower bounds on the combined elements of β and θ
formula, trms, A, L, and optsum: fields of the LMM field
X: fixed-effects model matrix
y: response vector

source

the distribution family for the response, and possibly the link function, must be specified.

julia> verbaggform = @formula(r2 ~ 1+a+g+b+s+m+(1|id)+(1|item));

julia> gm1 = fit(MixedModel, verbaggform, dat[:VerbAgg], Bernoulli())
Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
  r2 ~ 1 + a + g + b + s + m + (1 | id) + (1 | item)
  Distribution: Distributions.Bernoulli{Float64}
  Link: GLM.LogitLink()

  Deviance: 8135.8329

Variance components:
        Column     Variance   Std.Dev.  
id   (Intercept)  1.793533889 1.33922884
item (Intercept)  0.117159757 0.34228607

 Number of obs: 7584; levels of grouping factors: 316, 24

Fixed-effects parameters:
─────────────────────────────────────────────────────
                  Coef.  Std. Error       z  Pr(>|z|)
─────────────────────────────────────────────────────
(Intercept)   0.553166     0.38537     1.44    0.1512
a             0.0574327    0.016753    3.43    0.0006
g: M          0.320676     0.191209    1.68    0.0935
b: scold     -1.05976      0.184166   -5.75    <1e-8
b: shout     -2.10365      0.186524  -11.28    <1e-28
s: self      -1.05425      0.1512     -6.97    <1e-11
m: do        -0.707208     0.151013   -4.68    <1e-5
─────────────────────────────────────────────────────

The canonical link, which is GLM.LogitLink for the Bernoulli distribution, is used if no explicit link is specified.

Note that, in keeping with convention in the GLM package, the distribution family for a binary (i.e. 0/1) response is the Bernoulli distribution. The Binomial distribution is only used when the response is the fraction of trials returning a positive, in which case the number of trials must be specified as the case weights.

Optional arguments to fit!

An alternative approach is to create the GeneralizedLinearMixedModel object then call fit! on it. In this form optional arguments fast and/or nAGQ can be passed to the optimization process.

As the name implies, fast=true, provides a faster but somewhat less accurate fit. These fits may suffice for model comparisons.

julia> gm1a = fit!(GeneralizedLinearMixedModel(verbaggform, dat[:VerbAgg],
           Bernoulli()), fast=true)
Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
  r2 ~ 1 + a + g + b + s + m + (1 | id) + (1 | item)
  Distribution: Distributions.Bernoulli{Float64}
  Link: GLM.LogitLink()

  Deviance: 8136.1709

Variance components:
        Column    Variance   Std.Dev.  
id   (Intercept)  1.79270002 1.33891748
item (Intercept)  0.11875573 0.34460953

 Number of obs: 7584; levels of grouping factors: 316, 24

Fixed-effects parameters:
─────────────────────────────────────────────────────
                  Coef.  Std. Error       z  Pr(>|z|)
─────────────────────────────────────────────────────
(Intercept)   0.548543    0.385673     1.42    0.1549
a             0.0543802   0.0167462    3.25    0.0012
g: M          0.304244    0.191141     1.59    0.1114
b: scold     -1.01749     0.185216    -5.49    <1e-7
b: shout     -2.02067     0.187522   -10.78    <1e-26
s: self      -1.01255     0.15204     -6.66    <1e-10
m: do        -0.679102    0.151857    -4.47    <1e-5
─────────────────────────────────────────────────────

julia> deviance(gm1a) - deviance(gm1)
0.33801218210646766

julia> @time fit!(GeneralizedLinearMixedModel(verbaggform, dat[:VerbAgg],
           Bernoulli()));
  3.650444 seconds (278.72 k allocations: 20.271 MiB)

julia> @time fit!(GeneralizedLinearMixedModel(verbaggform, dat[:VerbAgg],
           Bernoulli()), fast=true);
  0.620628 seconds (50.41 k allocations: 9.654 MiB)

The optional argument nAGQ=k causes evaluation of the deviance function to use a k point adaptive Gauss-Hermite quadrature rule. This method only applies to models with a single, simple, scalar random-effects term, such as

julia> contraform = @formula(use ~ 1+a+abs2(a)+l+urb+(1|d));

julia> @time gm2 = fit!(GeneralizedLinearMixedModel(contraform,
           dat[:Contraception], Bernoulli()), nAGQ=9)
  3.717716 seconds (5.00 M allocations: 261.568 MiB, 1.77% gc time)
Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 9)
  use ~ 1 + a + :(abs2(a)) + l + urb + (1 | d)
  Distribution: Distributions.Bernoulli{Float64}
  Link: GLM.LogitLink()

  Deviance: 2372.4589

Variance components:
     Column    Variance   Std.Dev. 
d (Intercept)  0.22909357 0.4786372

 Number of obs: 1934; levels of grouping factors: 60

Fixed-effects parameters:
──────────────────────────────────────────────────────
                   Coef.   Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────────
(Intercept)  -1.03543     0.174556     -5.93    <1e-8
a             0.00353269  0.00923289    0.38    0.7020
abs2(a)      -0.00456312  0.000725329  -6.29    <1e-9
l: 1          0.815134    0.162222      5.02    <1e-6
l: 2          0.916551    0.18514       4.95    <1e-6
l: 3+         0.915354    0.185812      4.93    <1e-6
urb: Y        0.696706    0.119944      5.81    <1e-8
──────────────────────────────────────────────────────

julia> @time deviance(fit!(GeneralizedLinearMixedModel(contraform,
           dat[:Contraception], Bernoulli()), nAGQ=9, fast=true))
  0.068797 seconds (22.07 k allocations: 2.454 MiB)
2372.513592622964

julia> @time deviance(fit!(GeneralizedLinearMixedModel(contraform,
           dat[:Contraception], Bernoulli())))
  0.243641 seconds (68.63 k allocations: 5.027 MiB)
2372.7285823652946

julia> @time deviance(fit!(GeneralizedLinearMixedModel(contraform,
           dat[:Contraception], Bernoulli()), fast=true))
  0.038119 seconds (11.59 k allocations: 1.823 MiB)
2372.7844291358947

Extractor functions

LinearMixedModel and GeneralizedLinearMixedModel are subtypes of StatsBase.RegressionModel which, in turn, is a subtype of StatsBase.StatisticalModel. Many of the generic extractors defined in the StatsBase package have methods for these models.

Model-fit statistics

The statistics describing the quality of the model fit include

StatsBase.loglikelihood — Method.

loglikelihood(obj::StatisticalModel)

Return the log-likelihood of the model.

StatsBase.aic — Function.

aic(obj::StatisticalModel)

Akaike's Information Criterion, defined as $-2 \log L + 2k$, with $L$ the likelihood of the model, and k its number of consumed degrees of freedom (as returned by dof).

StatsBase.bic — Function.

bic(obj::StatisticalModel)

Bayesian Information Criterion, defined as $-2 \log L + k \log n$, with $L$ the likelihood of the model, $k$ its number of consumed degrees of freedom (as returned by dof), and $n$ the number of observations (as returned by nobs).

StatsBase.dof — Method.

dof(obj::StatisticalModel)

Return the number of degrees of freedom consumed in the model, including when applicable the intercept and the distribution's dispersion parameter.

StatsBase.nobs — Method.

nobs(obj::StatisticalModel)

Return the number of independent observations on which the model was fitted. Be careful when using this information, as the definition of an independent observation may vary depending on the model, on the format used to pass the data, on the sampling plan (if specified), etc.

julia> loglikelihood(fm1)
-163.66352994057004

julia> aic(fm1)
333.3270598811401

julia> bic(fm1)
337.5306520261266

julia> dof(fm1)   # 1 fixed effect, 2 variances
3

julia> nobs(fm1)  # 30 observations
30

julia> loglikelihood(gm1)
-4067.916429761946

In general the deviance of a statistical model fit is negative twice the log-likelihood adjusting for the saturated model.

StatsBase.deviance — Method.

deviance(obj::StatisticalModel)

Return the deviance of the model relative to a reference, which is usually when applicable the saturated model. It is equal, up to a constant, to $-2 \log L$, with $L$ the likelihood of the model.

Because it is not clear what the saturated model corresponding to a particular LinearMixedModel should be, negative twice the log-likelihood is called the objective.

MixedModels.objective — Function.

objective(m::LinearMixedModel)

Return negative twice the log-likelihood of model m

source

This value is also accessible as the deviance but the user should bear in mind that this doesn't have all the properties of a deviance which is corrected for the saturated model. For example, it is not necessarily non-negative.

julia> objective(fm1)
327.3270598811401

julia> deviance(fm1)
327.3270598811401

The value optimized when fitting a GeneralizedLinearMixedModel is the Laplace approximation to the deviance or an adaptive Gauss-Hermite evaluation.

MixedModels.deviance! — Function.

deviance!(m::GeneralizedLinearMixedModel, nAGQ=1)

Update m.η, m.μ, etc., install the working response and working weights in m.LMM, update m.LMM.A and m.LMM.R, then evaluate the deviance.

source

julia> MixedModels.deviance!(gm1)
8135.83285952388

Fixed-effects parameter estimates

The coef and fixef extractors both return the maximum likelihood estimates of the fixed-effects coefficients.

StatsBase.coef — Function.

coef(obj::StatisticalModel)

Return the coefficients of the model.

MixedModels.fixef — Function.

fixef(m::MixedModel, permuted=true)

Return the fixed-effects parameter vector estimate of m.

If permuted is true the vector elements are permuted according to m.trms[end - 1].piv and truncated to the rank of that term.

source

julia> show(coef(fm1))
[1527.4999999999993]
julia> show(fixef(fm1))
[1527.4999999999993]
julia> show(fixef(gm1))
[0.5531660445886255, 0.057432668911874415, 0.3206764775020489, -1.0597583684026035, -2.1036493549432858, -1.0542481028081037, -0.7072083631129692]

An alternative extractor for the fixed-effects coefficient is the β property. Properties whose names are Greek letters usually have an alternative spelling, which is the name of the Greek letter.

julia> show(fm1.β)
[1527.4999999999993]
julia> show(fm1.beta)
[1527.4999999999993]
julia> show(gm1.β)
[0.5531660445886255, 0.057432668911874415, 0.3206764775020489, -1.0597583684026035, -2.1036493549432858, -1.0542481028081037, -0.7072083631129692]

A full list of property names is returned by propertynames

julia> propertynames(fm1)
(:formula, :sqrtwts, :A, :L, :optsum, :θ, :theta, :β, :beta, :λ, :lambda, :stderror, :σ, :sigma, :σs, :sigmas, :b, :u, :lowerbd, :X, :y, :rePCA, :reterms, :feterms, :objective, :pvalues)

julia> propertynames(gm1)
(:A, :L, :theta, :beta, :coef, :λ, :lambda, :σ, :sigma, :X, :y, :lowerbd, :σρs, :σs, :LMM, :β, :β₀, :θ, :b, :u, :u₀, :resp, :η, :wt, :devc, :devc0, :sd, :mult)

The variance-covariance matrix of the fixed-effects coefficients is returned by

StatsBase.vcov — Function.

vcov(obj::StatisticalModel)

Return the variance-covariance matrix for the coefficients of the model.

vcov(m::LinearMixedModel)

Returns the variance-covariance matrix of the fixed effects. If corr=true, then correlation of fixed effects is returned instead.

source

julia> vcov(fm2)
2×2 Array{Float64,2}:
 43.9868   -1.37039
 -1.37039   2.25671

julia> vcov(gm1)
7×7 Array{Float64,2}:
  0.14851     -0.00560477   -0.00977107   …  -0.0114559    -0.0114571
 -0.00560477   0.000280662   7.19137e-5      -1.47974e-5   -1.02422e-5
 -0.00977107   7.19137e-5    0.036561        -8.04374e-5   -5.25878e-5
 -0.0169722   -1.43718e-5   -9.25577e-5       0.000265793   0.000172104
 -0.0171447   -2.90572e-5   -0.000162381      0.000658948   0.000520539
 -0.0114559   -1.47974e-5   -8.04374e-5   …   0.0228615     0.000247791
 -0.0114571   -1.02422e-5   -5.25878e-5       0.000247791   0.0228049

The standard errors are the square roots of the diagonal elements of the estimated variance-covariance matrix of the fixed-effects coefficient estimators.

StatsBase.stderror — Function.

stderror(obj::StatisticalModel)

Return the standard errors for the coefficients of the model.

julia> show(StatsBase.stderror(fm2))
[6.632257775855945, 1.5022355309693658]
julia> show(StatsBase.stderror(gm1))
[0.3853701596029146, 0.016752979637955016, 0.19120935781772874, 0.18416579838061983, 0.18652438070920965, 0.15120011466690528, 0.15101307761446386]

Finally, the coeftable generic produces a table of coefficient estimates, their standard errors, and their ratio. The p-values quoted here should be regarded as approximations.

StatsBase.coeftable — Function.

coeftable(obj::StatisticalModel; level::Real=0.95)

Return a table of class CoefTable with coefficients and related statistics. level determines the level for confidence intervals (by default, 95%).

julia> coeftable(fm2)
──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  251.405      6.63226  37.91    <1e-99
U             10.4673     1.50224   6.97    <1e-11
──────────────────────────────────────────────────

Covariance parameter estimates

The covariance parameters estimates, in the form shown in the model summary, are a VarCorr object

MixedModels.VarCorr — Type.

VarCorr

Information from the fitted random-effects variance-covariance matrices.

Members

σρ: a NamedTuple of NamedTuples as returned from σρs
s: the estimate of the scale parameter in the distribution of the conditional dist'n of Y

The main purpose of defining this type is to isolate the logic in the show method.

source

julia> VarCorr(fm2)
Variance components:
            Column    Variance   Std.Dev.    Corr.
G        (Intercept)  565.510678 23.7804684
         U             32.682124  5.7168282  0.08
Residual              654.941447 25.5918238


julia> VarCorr(gm1)
Variance components:
        Column     Variance   Std.Dev.  
id   (Intercept)  1.793533889 1.33922884
item (Intercept)  0.117159757 0.34228607

Individual components are returned by other extractors

MixedModels.varest — Function.

varest(m::LinearMixedModel)

Returns the estimate of σ², the variance of the conditional distribution of Y given B.

source

MixedModels.sdest — Function.

sdest(m::LinearMixedModel)

Return the estimate of σ, the standard deviation of the per-observation noise.

source

julia> varest(fm2)
654.94144673283

julia> sdest(fm2)
25.591823825839963

julia> fm2.σ
25.591823825839963

Conditional modes of the random effects

The ranef extractor

MixedModels.ranef — Function.

ranef(m::MixedModel; uscale=false, named=false)

Return, as a Vector{Vector{T}} (Vector{NamedVector{T}} if named=true), the conditional modes of the random effects in model m.

If uscale is true the random effects are on the spherical (i.e. u) scale, otherwise on the original scale.

source

julia> ranef(fm1)
1-element Array{Array{Float64,2},1}:
 [-16.62822143006428 0.3695160317797815 … 53.57982460798647 -42.49434365460914]

julia> fm1.b
1-element Array{Array{Float64,2},1}:
 [-16.62822143006428 0.3695160317797815 … 53.57982460798647 -42.49434365460914]

returns the conditional modes of the random effects given the observed data. That is, these are the values that maximize the conditional density of the random effects given the observed data. For a LinearMixedModel these are also the conditional mean values.

These are sometimes called the best linear unbiased predictors or BLUPs but that name is not particularly meaningful.

At a superficial level these can be considered as the "estimates" of the random effects, with a bit of hand waving, but pursuing this analogy too far usually results in confusion.

The corresponding conditional variances are returned by

MixedModels.condVar — Function.

condVar(m::LinearMixedModel)

Return the conditional variances matrices of the random effects.

The random effects are returned by ranef as a vector of length k, where k is the number of random effects terms. The ith element is a matrix of size vᵢ × ℓᵢ where vᵢ is the size of the vector-valued random effects for each of the ℓᵢ levels of the grouping factor. Technically those values are the modes of the conditional distribution of the random effects given the observed data.

This function returns an array of k three dimensional arrays, where the ith array is of size vᵢ × vᵢ × ℓᵢ. These are the diagonal blocks from the conditional variance-covariance matrix,

s² Λ(Λ'Z'ZΛ + I)⁻¹Λ'

source

julia> condVar(fm1)
1-element Array{Array{Float64,3},1}:
 [362.3104715146578]

[362.3104715146578]

[362.3104715146578]

[362.3104715146578]

[362.3104715146578]

[362.3104715146578]