Model constructors

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
  • reterms: a Vector{ReMat{T}} of random-effects terms.
  • feterms: a Vector{FeMat{T}} of the fixed-effects model matrix 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
  • 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 DataFrames.DataFrame…
  :Genetics      => 60×5 DataFrames.DataFrame…
  :Contraception => 1934×6 DataFrames.DataFrame. Omitted printing of 1 columns…
  :Mmmec         => 354×6 DataFrames.DataFrame…
  :kb07          => 1790×10 DataFrames.DataFrame. Omitted printing of 4 columns…
  :Rail          => 18×2 DataFrames.DataFrame…
  :KKL           => 53765×24 DataFrames.DataFrame. Omitted printing of 18 colum…
  :Bond          => 21×3 DataFrames.DataFrame…
  :VerbAgg       => 7584×9 DataFrames.DataFrame. Omitted printing of 4 columns…
  :ml1m          => 1000209×3 DataFrames.DataFrame…
  :ergoStool     => 36×3 DataFrames.DataFrame…
  :s3bbx         => 2449×6 DataFrames.DataFrame. Omitted printing of 1 columns…
  :cake          => 270×5 DataFrames.DataFrame…
  :Cultivation   => 24×4 DataFrames.DataFrame…
  :Pastes        => 60×4 DataFrames.DataFrame…
  :Exam          => 4059×5 DataFrames.DataFrame…
  :Socatt        => 1056×9 DataFrames.DataFrame. Omitted printing of 5 columns…
  :WWheat        => 60×3 DataFrames.DataFrame…
  :Pixel         => 102×5 DataFrames.DataFrame…
  ⋮              => ⋮

julia> describe(dat[:Dyestuff])
2×8 DataFrames.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:
──────────────────────────────────────────────────
             Estimate  Std.Error  z value  P(>|z|)
──────────────────────────────────────────────────
(Intercept)    1527.5    17.6946   86.326   <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:
──────────────────────────────────────────────────
             Estimate  Std.Error  z value  P(>|z|)
──────────────────────────────────────────────────
(Intercept)    1527.5    17.6946   86.326   <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.715776 seconds (1.05 M allocations: 55.982 MiB, 3.47% gc time)

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:
──────────────────────────────────────────────────
             Estimate  Std.Error  z value  P(>|z|)
──────────────────────────────────────────────────
(Intercept)    1527.5    19.3834  78.8045   <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.51069 23.780469
         U             32.68212  5.716828  0.08
Residual              654.94145 25.591824
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
             Estimate  Std.Error   z value  P(>|z|)
───────────────────────────────────────────────────
(Intercept)  251.405     6.63226  37.9064    <1e-99
U             10.4673    1.50224   6.96781   <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.71497949 0.8455646
H        (Intercept)  3.13519360 1.7706478
Residual              0.30242640 0.5499331
 Number of obs: 144; levels of grouping factors: 24, 6

  Fixed-effects parameters:
──────────────────────────────────────────────────
             Estimate  Std.Error  z value  P(>|z|)
──────────────────────────────────────────────────
(Intercept)   22.9722   0.744596  30.8519   <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:
──────────────────────────────────────────────────
             Estimate  Std.Error  z value  P(>|z|)
──────────────────────────────────────────────────
(Intercept)   60.0533   0.642136  93.5212   <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.10541798 0.32468136
H        (Intercept)  0.25841636 0.50834669
Residual              1.38472777 1.17674456
 Number of obs: 73421; levels of grouping factors: 2972, 1128

  Fixed-effects parameters:
────────────────────────────────────────────────────────
                Estimate  Std.Error     z value  P(>|z|)
────────────────────────────────────────────────────────
(Intercept)    3.22961    0.064053   50.4209      <1e-99
A: 1           0.252025   0.0686507   3.67112     0.0002
I: 5           0.129536   0.101294    1.27882     0.2010
I: 10         -0.176751   0.0881352  -2.00545     0.0449
I: 12          0.0517102  0.0817523   0.632522    0.5270
I: 6           0.0347319  0.085621    0.405647    0.6850
I: 7           0.14594    0.0997984   1.46235     0.1436
I: 4           0.151689   0.0816897   1.85689     0.0633
I: 8           0.104206   0.118751    0.877517    0.3802
I: 9           0.0440401  0.0962985   0.457329    0.6474
I: 14          0.0517546  0.0986029   0.524879    0.5997
I: 1           0.0466719  0.101942    0.457828    0.6471
I: 3           0.0563461  0.0977925   0.57618     0.5645
I: 11          0.0596536  0.100233    0.59515     0.5517
I: 2           0.0055628  0.110867    0.0501756   0.9600
A: 1 & I: 5   -0.180757   0.123179   -1.46744     0.1423
A: 1 & I: 10   0.0186492  0.110017    0.169513    0.8654
A: 1 & I: 12  -0.282269   0.0792937  -3.55979     0.0004
A: 1 & I: 6   -0.494464   0.0790278  -6.25683     <1e-9 
A: 1 & I: 7   -0.392054   0.110313   -3.55403     0.0004
A: 1 & I: 4   -0.278547   0.0823727  -3.38154     0.0007
A: 1 & I: 8   -0.189526   0.111449   -1.70056     0.0890
A: 1 & I: 9   -0.499868   0.0885423  -5.64553     <1e-7 
A: 1 & I: 14  -0.497162   0.0917162  -5.42065     <1e-7 
A: 1 & I: 1   -0.24042    0.0982071  -2.4481      0.0144
A: 1 & I: 3   -0.223013   0.0890548  -2.50422     0.0123
A: 1 & I: 11  -0.516997   0.0809077  -6.38997     <1e-9 
A: 1 & I: 2   -0.384773   0.091843   -4.18946     <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.258970 24.17145
         U             33.632805  5.79938   .  
Residual              653.115782 25.55613
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
             Estimate  Std.Error   z value  P(>|z|)
───────────────────────────────────────────────────
(Intercept)  251.405     6.70771  37.48      <1e-99
U             10.4673    1.51931   6.88951   <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.258970 24.17145
         U             33.632805  5.79938   .  
Residual              653.115782 25.55613
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
             Estimate  Std.Error   z value  P(>|z|)
───────────────────────────────────────────────────
(Intercept)  251.405     6.70771  37.48      <1e-99
U             10.4673    1.51931   6.88951   <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.258970 24.17145
         U             33.632805  5.79938   .  
Residual              653.115782 25.55613
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
             Estimate  Std.Error   z value  P(>|z|)
───────────────────────────────────────────────────
(Intercept)  251.405     6.70771  37.48      <1e-99
U             10.4673    1.51931   6.88951   <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.39922 1610.79843 1742.79843 1953.53359

Variance components:
            Column     Variance   Std.Dev.    Corr.
G        (Intercept)   955.578308 30.9124297
         U: 1.0        497.183571 22.2976136 -0.30
         U: 2.0        915.645818 30.2596401 -0.57  0.75
         U: 3.0       1267.226692 35.5981276 -0.37  0.72  0.87
         U: 4.0       1484.909119 38.5345185 -0.31  0.59  0.67  0.91
         U: 5.0       2296.038533 47.9169963 -0.25  0.46  0.45  0.70  0.85
         U: 6.0       3850.672299 62.0537855 -0.28  0.30  0.48  0.70  0.77  0.75
         U: 7.0       1806.262575 42.5001479 -0.16  0.22  0.47  0.50  0.63  0.64  0.71
         U: 8.0       3151.792889 56.1408309 -0.20  0.29  0.36  0.56  0.73  0.90  0.73  0.74
         U: 9.0       3075.767360 55.4596011  0.06  0.25  0.16  0.38  0.59  0.78  0.38  0.53  0.85
Residual                19.261211  4.3887597
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
────────────────────────────────────────────────────
              Estimate  Std.Error   z value  P(>|z|)
────────────────────────────────────────────────────
(Intercept)  256.652      7.3592   34.875     <1e-99
U: 1.0         7.84395    5.45541   1.43783   0.1505
U: 2.0         8.71009    7.28075   1.19632   0.2316
U: 3.0        26.3402     8.51714   3.09261   0.0020
U: 4.0        31.9976     9.19973   3.4781    0.0005
U: 5.0        51.8667    11.3885    4.5543    <1e-5 
U: 6.0        55.5265    14.6992    3.77752   0.0002
U: 7.0        62.0988    10.1236    6.13404   <1e-9 
U: 8.0        79.9777    13.3131    6.00743   <1e-8 
U: 9.0        94.1994    13.1536    7.16152   <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.39935 1610.79870 1744.79870 1958.72681

Variance components:
            Column     Variance   Std.Dev.    Corr.
G        (Intercept)   181.213462 13.4615550
         U: 0.0        782.422065 27.9718084   .  
         U: 1.0        863.277684 29.3815875   .    0.69
         U: 2.0        627.634941 25.0526434   .    0.34  0.73
         U: 3.0       1231.455404 35.0920989   .    0.37  0.71  0.87
         U: 4.0       1513.458193 38.9031900   .    0.36  0.60  0.65  0.91
         U: 5.0       2335.558526 48.3276166   .    0.30  0.48  0.43  0.70  0.85
         U: 6.0       3570.086607 59.7502017   .    0.14  0.25  0.41  0.66  0.74  0.73
         U: 7.0       2171.047045 46.5944959   .    0.43  0.42  0.55  0.56  0.67  0.67  0.69
         U: 8.0       3225.633405 56.7946600   .    0.26  0.34  0.35  0.57  0.73  0.90  0.72  0.75
         U: 9.0       4044.254162 63.5944507   .    0.49  0.52  0.38  0.54  0.70  0.83  0.43  0.64  0.88
Residual                10.760307  3.2802908
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
────────────────────────────────────────────────────
              Estimate  Std.Error   z value  P(>|z|)
────────────────────────────────────────────────────
(Intercept)  256.652      7.35752  34.8829    <1e-99
U: 1.0         7.84395    5.45364   1.4383    0.1504
U: 2.0         8.71009    7.28356   1.19586   0.2318
U: 3.0        26.3402     8.51764   3.09243   0.0020
U: 4.0        31.9976     9.19943   3.47822   0.0005
U: 5.0        51.8666    11.3926    4.55267   <1e-5 
U: 6.0        55.5264    14.7034    3.77643   0.0002
U: 7.0        62.0988    10.1198    6.13637   <1e-9 
U: 8.0        79.9777    13.32      6.00433   <1e-8 
U: 9.0        94.1994    13.1575    7.1594    <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.56168 30.960647
         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.69396 22.796797   .     .     .     .     .  
         U: 6.0       1703.85878 41.277824   .     .     .     .     .     .  
         U: 7.0        608.86850 24.675261   .     .     .     .     .     .     .  
         U: 8.0       1273.10854 35.680647   .     .     .     .     .     .     .     .  
         U: 9.0       1753.94556 41.880133   .     .     .     .     .     .     .     .     .  
Residual               434.85234 20.853113
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
────────────────────────────────────────────────────
              Estimate  Std.Error   z value  P(>|z|)
────────────────────────────────────────────────────
(Intercept)  256.652      8.7984   29.1703    <1e-99
U: 1.0         7.84395    6.95104   1.12846   0.2591
U: 2.0         8.71009    6.95104   1.25306   0.2102
U: 3.0        26.3402     6.95104   3.78939   0.0002
U: 4.0        31.9976     6.95104   4.60329   <1e-5 
U: 5.0        51.8666     8.78572   5.90352   <1e-8 
U: 6.0        55.5264    11.9572    4.64375   <1e-5 
U: 7.0        62.0988     9.06328   6.85169   <1e-11
U: 8.0        79.9777    10.9108    7.33015   <1e-12
U: 9.0        94.1994    12.073     7.80246   <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.9844  1757.9688  1801.9688  1872.2139

Variance components:
            Column      Variance      Std.Dev.    Corr.
G        (Intercept)  1134.018890418 33.67519696
         U: 0.0        775.482297289 27.84748278   .  
         U: 1.0        356.726068226 18.88719323   .     .  
         U: 2.0        220.598457629 14.85255728   .     .     .  
         U: 3.0          0.041430896  0.20354581   .     .     .     .  
         U: 4.0         43.482811564  6.59414980   .     .     .     .     .  
         U: 5.0        670.667693798 25.89725263   .     .     .     .     .     .  
         U: 6.0       1740.087547279 41.71435661   .     .     .     .     .     .     .  
         U: 7.0        904.947269466 30.08234149   .     .     .     .     .     .     .     .  
         U: 8.0       1456.355135727 38.16222132   .     .     .     .     .     .     .     .     .  
         U: 9.0       2026.627987970 45.01808512   .     .     .     .     .     .     .     .     .     .  
Residual               181.681891819 13.47894253
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
─────────────────────────────────────────────────────
              Estimate  Std.Error    z value  P(>|z|)
─────────────────────────────────────────────────────
(Intercept)  256.652     10.7785   23.8114     <1e-99
U: 1.0         7.84395    9.11523   0.860533   0.3895
U: 2.0         8.71009    8.69049   1.00226    0.3162
U: 3.0        26.3402     7.95434   3.31143    0.0009
U: 4.0        31.9976     8.10462   3.94807    <1e-4 
U: 5.0        51.8667    10.0264    5.17301    <1e-6 
U: 6.0        55.5265    12.6468    4.39056    <1e-4 
U: 7.0        62.0988    10.6557    5.82775    <1e-8 
U: 8.0        79.9777    12.0074    6.6607     <1e-10
U: 9.0        94.1994    13.2612    7.10338    <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.64355891 1.28201362
item (Intercept)  0.10735667 0.32765327

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

Fixed-effects parameters:
─────────────────────────────────────────────────────
              Estimate  Std.Error    z value  P(>|z|)
─────────────────────────────────────────────────────
(Intercept)   0.553512  0.368905     1.50042   0.1335
a             0.05742   0.0160373    3.58041   0.0003
g: M          0.320766  0.183041     1.75243   0.0797
b: scold     -1.05987   0.176294    -6.01194   <1e-8 
b: shout     -2.10387   0.178552   -11.783     <1e-31
s: self      -1.05433   0.144737    -7.28446   <1e-12
m: do        -0.707051  0.144558    -4.89112   <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.636122430 1.27911001
item (Intercept)  0.108383391 0.32921633

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

Fixed-effects parameters:
──────────────────────────────────────────────────────
               Estimate  Std.Error    z value  P(>|z|)
──────────────────────────────────────────────────────
(Intercept)   0.548543   0.368446     1.4888    0.1365
a             0.0543802  0.0159982    3.39915   0.0007
g: M          0.304244   0.182603     1.66614   0.0957
b: scold     -1.01749    0.176943    -5.75038   <1e-8 
b: shout     -2.02067    0.179146   -11.2795    <1e-28
s: self      -1.01255    0.145248    -6.97114   <1e-11
m: do        -0.679102   0.145074    -4.68108   <1e-5 
──────────────────────────────────────────────────────

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

julia> @time fit!(GeneralizedLinearMixedModel(verbaggform, dat[:VerbAgg],
    Bernoulli()));
  3.930082 seconds (15.28 M allocations: 135.041 MiB, 1.75% gc time)

julia> @time fit!(GeneralizedLinearMixedModel(verbaggform, dat[:VerbAgg],
    Bernoulli()), fast=true);
  0.635499 seconds (2.44 M allocations: 27.878 MiB, 1.27% gc time)

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)
  4.904342 seconds (9.16 M allocations: 319.795 MiB, 5.72% 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.22222789 0.47141053

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

Fixed-effects parameters:
─────────────────────────────────────────────────────────
                Estimate    Std.Error    z value  P(>|z|)
─────────────────────────────────────────────────────────
(Intercept)  -1.0354      0.171933     -6.02213    <1e-8 
a             0.00353186  0.00909421    0.388363   0.6977
abs2(a)      -0.00456294  0.000714433  -6.38679    <1e-9 
l: 1          0.815183    0.159785      5.10173    <1e-6 
l: 2          0.916476    0.18236       5.02565    <1e-6 
l: 3+         0.915342    0.183021      5.00131    <1e-6 
urb: Y        0.696715    0.118142      5.89727    <1e-8 
─────────────────────────────────────────────────────────

julia> @time deviance(fit!(GeneralizedLinearMixedModel(contraform,
    dat[:Contraception], Bernoulli()), nAGQ=9, fast=true))
  0.083751 seconds (421.33 k allocations: 5.277 MiB)
2372.5135926229646

julia> @time deviance(fit!(GeneralizedLinearMixedModel(contraform,
    dat[:Contraception], Bernoulli())))
  0.296034 seconds (1.47 M allocations: 14.747 MiB)
2372.7285827581495

julia> @time deviance(fit!(GeneralizedLinearMixedModel(contraform,
    dat[:Contraception], Bernoulli()), fast=true))
  0.074024 seconds (280.48 k allocations: 3.856 MiB, 19.05% gc time)
2372.784429135894

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.9164280257423

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.832856051482

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.5535116086379206, 0.05742001139348312, 0.32076602931752635, -1.059867665970283, -2.103869984409917, -1.0543329105826895, -0.7070511487154889]

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.5535116086379206, 0.05742001139348312, 0.32076602931752635, -1.059867665970283, -2.103869984409917, -1.0543329105826895, -0.7070511487154889]

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.

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.136091    -0.00513612   -0.00895404   …  -0.0104975    -0.0104986  
 -0.00513612   0.000257194   6.59009e-5      -1.35602e-5   -9.38556e-6 
 -0.00895404   6.59009e-5    0.0335039       -7.37176e-5   -4.81924e-5 
 -0.0155523   -1.31704e-5   -8.4825e-5        0.000243582   0.000157714
 -0.0157104   -2.6628e-5    -0.000148816      0.000603881   0.000477019
 -0.0104975   -1.35602e-5   -7.37176e-5   …   0.0209489     0.000227073
 -0.0104986   -9.38556e-6   -4.81924e-5       0.000227073   0.0208971

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.632257825314581, 1.502235453639816]
julia> show(StatsBase.stderror(gm1))
[0.36890512061819736, 0.01603727875112492, 0.1830407884432926, 0.1762937885165568, 0.17855176330630151, 0.1447372582156264, 0.14455818519253938]

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)
───────────────────────────────────────────────────
             Estimate  Std.Error   z value  P(>|z|)
───────────────────────────────────────────────────
(Intercept)  251.405     6.63226  37.9064    <1e-99
U             10.4673    1.50224   6.96781   <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.51069 23.780469
         U             32.68212  5.716828  0.08
Residual              654.94145 25.591824


julia> VarCorr(gm1)
Variance components:
        Column    Variance   Std.Dev.  
id   (Intercept)  1.64355891 1.28201362
item (Intercept)  0.10735667 0.32765327

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.9414513956141

julia> sdest(fm2)
25.59182391693906

julia> fm2.σ
25.59182391693906

Conditional modes of the random effects

The ranef extractor

MixedModels.ranef — Function.
ranef(m::LinearMixedModel; uscale=false) #, named=true)

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.62822143006434 0.36951603177972425 … 53.57982460798641 -42.49434365460919]

julia> fm1.b
1-element Array{Array{Float64,2},1}:
 [-16.62822143006434 0.36951603177972425 … 53.57982460798641 -42.49434365460919]

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]