Model constructors

The LinearMixedModel type represents a linear mixed-effects model. Typically it is constructed from a Formula and an appropriate Table 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
reterms: a Vector{AbstractReMat{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.
parmap : Vector{NTuple{3,Int}} of (block, row, column) mapping of θ to λ
dims : NamedTuple{(:n, :p, :nretrms),NTuple{3,Int}} of dimensions. p is the rank of X, which may be smaller than size(X, 2).
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 .arrow format in this package. Often, for convenience, we will convert these to DataFrames. These data sets include the dyestuff and dyestuff2 data sets.

MixedModels.dataset — Function

dataset(nm)

Return, as an Arrow.Table, the test data set named nm, which can be a String or Symbol

source

using DataFrames, MixedModels
using StatsBase: describe
dyestuff = DataFrame(MixedModels.dataset(:dyestuff))
describe(dyestuff)

2 rows × 8 columns

	variable	mean	min	median	max	nunique	nmissing	eltype
	Symbol	Union…	Any	Union…	Any	Union…	Nothing	DataType
1	batch		A		F	6		String
2	yield	1527.5	1440	1530.0	1635			Int16

Models with simple, scalar random effects

The formula language in Julia is similar to that in R which is based on (Wilkinson and Rogers 1973). In Julia a formula must be enclosed in a call to the @formula macro.

StatsModels.@formula — Macro

@formula(ex)

Capture and parse a formula expression as a Formula struct.

A formula is an abstract specification of a dependence between left-hand and right-hand side variables as in, e.g., a regression model. Each side specifies at a high level how tabular data is to be converted to a numerical matrix suitable for modeling. This specification looks something like Julia code, is represented as a Julia Expr, but uses special syntax. The @formula macro takes an expression like y ~ 1 + a*b, transforms it according to the formula syntax rules into a lowered form (like y ~ 1 + a + b + a&b), and constructs a Formula struct which captures the original expression, the lowered expression, and the left- and right-hand-side.

Operators that have special interpretations in this syntax are

~ is the formula separator, where it is a binary operator (the first argument is the left-hand side, and the second is the right-hand side.
+ concatenates variables as columns when generating a model matrix.
& representes an interaction between two or more variables, which corresponds to a row-wise kronecker product of the individual terms (or element-wise product if all terms involved are continuous/scalar).
* expands to all main effects and interactions: a*b is equivalent to a+b+a&b, a*b*c to a+b+c+a&b+a&c+b&c+a&b&c, etc.
1, 0, and -1 indicate the presence (for 1) or absence (for 0 and -1) of an intercept column.

The rules that are applied are

The associative rule (un-nests nested calls to +, &, and *).
The distributive rule (interactions & distribute over concatenation +).
The * rule expands a*b to a+b+a&b (recursively).
Subtraction is converted to addition and negation, so x-1 becomes x + -1 (applies only to subtraction of literal 1).
Single-argument & calls are stripped, so &(x) becomes the main effect x.

A basic model with simple, scalar random effects for the levels of batch (the batch of an intermediate product, in this case) is declared and fit as

fm = @formula(yield ~ 1 + (1|batch))
fm1 = fit(MixedModel, fm, dyestuff)

Linear mixed model fit by maximum likelihood
 yield ~ 1 + (1 | batch)
   logLik   -2 logLik     AIC       AICc        BIC    
  -163.6635   327.3271   333.3271   334.2501   337.5307

Variance components:
            Column    Variance Std.Dev.
batch    (Intercept)  1388.3333 37.2603
Residual              2451.2500 49.5101
 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 subsequent calls to such functions are much faster.)

using BenchmarkTools
dyestuff2 = MixedModels.dataset(:dyestuff2)
@benchmark fit(MixedModel, $fm, $dyestuff2)

BenchmarkTools.Trial: 
  memory estimate:  57.91 KiB
  allocs estimate:  1129
  --------------
  minimum time:     295.894 μs (0.00% GC)
  median time:      350.094 μs (0.00% GC)
  mean time:        388.377 μs (5.24% GC)
  maximum time:     47.754 ms (95.23% GC)
  --------------
  samples:          10000
  evals/sample:     1

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

fm1reml = fit(MixedModel, fm, dyestuff, REML=true)

Linear mixed model fit by REML
 yield ~ 1 + (1 | batch)
 REML criterion at convergence: 319.6542768422538

Variance components:
            Column    Variance Std.Dev.
batch    (Intercept)  1764.0506 42.0006
Residual              2451.2499 49.5101
 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
────────────────────────────────────────────────

Float-point type in the model

The type of fm1

typeof(fm1)

LinearMixedModel{Float64}

includes the floating point type used internally for the various matrices, vectors, etc. that represent the model. At present, this will always be Float64 because the parameter estimates are optimized using the NLopt package which calls compiled C code that only allows for optimization with respect to a Float64 parameter vector.

So in theory other floating point types, such as BigFloat or Float32, can be used to define a model but in practice only Float64 works at present.

In theory, theory and practice are the same. In practice, they aren't. – Anon

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, reaction, on several subjects, subj, after 0 to 9 days of sleep deprivation, days. A model with random intercepts and random slopes for each subject, allowing for within-subject correlation of the slope and intercept, is fit as

sleepstudy = MixedModels.dataset(:sleepstudy)
fm2 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy)

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + (1 + days | subj)
   logLik   -2 logLik     AIC       AICc        BIC    
  -875.9697  1751.9393  1763.9393  1764.4249  1783.0971

Variance components:
            Column    Variance Std.Dev.   Corr.
subj     (Intercept)  565.51069 23.78047
         days          32.68212  5.71683 +0.08
Residual              654.94145 25.59182
 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
days          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, and for the sample. As every sample is used on every plate these two factors are crossed.

penicillin = MixedModels.dataset(:penicillin)
fm3 = fit(MixedModel, @formula(diameter ~ 1 + (1|plate) + (1|sample)), penicillin)

Linear mixed model fit by maximum likelihood
 diameter ~ 1 + (1 | plate) + (1 | sample)
   logLik   -2 logLik     AIC       AICc        BIC    
  -166.0942   332.1883   340.1883   340.4761   352.0676

Variance components:
            Column   Variance Std.Dev. 
plate    (Intercept)  0.714979 0.845565
sample   (Intercept)  3.135194 1.770648
Residual              0.302426 0.549933
 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 cask grouping factor is nested within the batch grouping factor in the Pastes data.

pastes = DataFrame(MixedModels.dataset(:pastes))
describe(pastes)

3 rows × 8 columns

	variable	mean	min	median	max	nunique	nmissing	eltype
	Symbol	Union…	Any	Union…	Any	Union…	Nothing	DataType
1	batch		A		J	10		String
2	cask		a		c	3		String
3	strength	60.0533	54.2	59.3	66.0			Float64

This can be expressed using the solidus (the "/" character) to separate grouping factors, read "cask nested within batch":

fm4a = fit(MixedModel, @formula(strength ~ 1 + (1|batch/cask)), pastes)

Linear mixed model fit by maximum likelihood
 strength ~ 1 + (1 | batch) + (1 | batch & cask)
   logLik   -2 logLik     AIC       AICc        BIC    
  -123.9972   247.9945   255.9945   256.7217   264.3718

Variance components:
                Column   Variance Std.Dev. 
batch & cask (Intercept)  8.433616 2.904069
batch        (Intercept)  1.199180 1.095071
Residual                  0.678002 0.823409
 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
─────────────────────────────────────────────────

If the levels of the inner grouping factor are unique across the levels of the outer grouping factor, then this nesting does not need to expressed explicitly in the model syntax. For example, defining sample to be the combination of batch and cask, yields a naming scheme where the nesting is apparent from the data even if not expressed in the formula. (That is, each level of sample occurs in conjunction with only one level of batch.) As such, this model is equivalent to the previous one.

pastes.sample = (string.(pastes.cask, "&",  pastes.batch))
fm4b = fit(MixedModel, @formula(strength ~ 1 + (1|sample) + (1|batch)), pastes)

Linear mixed model fit by maximum likelihood
 strength ~ 1 + (1 | sample) + (1 | batch)
   logLik   -2 logLik     AIC       AICc        BIC    
  -123.9972   247.9945   255.9945   256.7217   264.3718

Variance components:
            Column   Variance Std.Dev. 
sample   (Intercept)  8.433617 2.904069
batch    (Intercept)  1.199179 1.095070
Residual              0.678002 0.823409
 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, s, of instructors, d. Additional covariates include the academic department, dept, in which the course was given and service, whether or not it was a service course.

insteval = MixedModels.dataset(:insteval)
fm5 = fit(MixedModel, @formula(y ~ 1 + service * dept + (1|s) + (1|d)), insteval)

Linear mixed model fit by maximum likelihood
 y ~ 1 + service + dept + service & dept + (1 | s) + (1 | d)
    logLik     -2 logLik       AIC         AICc          BIC     
 -118792.7767  237585.5534  237647.5534  237647.5804  237932.8763

Variance components:
            Column   Variance Std.Dev. 
s        (Intercept)  0.105418 0.324681
d        (Intercept)  0.258416 0.508347
Residual              1.384728 1.176745
 Number of obs: 73421; levels of grouping factors: 2972, 1128

  Fixed-effects parameters:
────────────────────────────────────────────────────────────────
                              Coef.  Std. Error      z  Pr(>|z|)
────────────────────────────────────────────────────────────────
(Intercept)              3.27628      0.0793647  41.28    <1e-99
service: Y               0.0116044    0.0699321   0.17    0.8682
dept: D02               -0.0411091    0.120331   -0.34    0.7326
dept: D03                0.00967413   0.108411    0.09    0.9289
dept: D04                0.105017     0.0944964   1.11    0.2664
dept: D05                0.0828643    0.11148     0.74    0.4573
dept: D06               -0.01194      0.0978342  -0.12    0.9029
dept: D07                0.0992679    0.110598    0.90    0.3694
dept: D08                0.0575337    0.127935    0.45    0.6529
dept: D09               -0.00263181   0.107085   -0.02    0.9804
dept: D10               -0.223423     0.099838   -2.24    0.0252
dept: D11                0.0129816    0.110639    0.12    0.9066
dept: D12                0.00503825   0.0944243   0.05    0.9574
dept: D14                0.0050827    0.109041    0.05    0.9628
dept: D15               -0.0466719    0.101942   -0.46    0.6471
service: Y & dept: D02  -0.144352     0.0929527  -1.55    0.1204
service: Y & dept: D03   0.0174078    0.0927237   0.19    0.8511
service: Y & dept: D04  -0.0381262    0.0810901  -0.47    0.6382
service: Y & dept: D05   0.0596632    0.123952    0.48    0.6303
service: Y & dept: D06  -0.254044     0.080781   -3.14    0.0017
service: Y & dept: D07  -0.151634     0.11157    -1.36    0.1741
service: Y & dept: D08   0.0508942    0.112189    0.45    0.6501
service: Y & dept: D09  -0.259448     0.0899448  -2.88    0.0039
service: Y & dept: D10   0.25907      0.111369    2.33    0.0200
service: Y & dept: D11  -0.276577     0.0819621  -3.37    0.0007
service: Y & dept: D12  -0.0418489    0.0792928  -0.53    0.5977
service: Y & dept: D14  -0.256742     0.0931016  -2.76    0.0058
service: Y & dept: D15   0.24042      0.0982071   2.45    0.0144
────────────────────────────────────────────────────────────────

Simplifying the random effect correlation structure

MixedModels.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 per day 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.

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

fm2zerocorr_fm = fit(MixedModel, @formula(reaction ~ 1 + days + zerocorr(1 + days|subj)), sleepstudy)

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + MixedModels.ZeroCorr((1 + days | subj))
   logLik   -2 logLik     AIC       AICc        BIC    
  -876.0016  1752.0033  1762.0033  1762.3481  1777.9680

Variance components:
            Column    Variance Std.Dev.   Corr.
subj     (Intercept)  584.25897 24.17145
         days          33.63281  5.79938   .  
Residual              653.11578 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
days          10.4673     1.51931   6.89    <1e-11
──────────────────────────────────────────────────

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

fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj) + (days|subj)), sleepstudy)

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + (1 | subj) + (days | subj)
   logLik   -2 logLik     AIC       AICc        BIC    
  -876.0016  1752.0033  1762.0033  1762.3481  1777.9680

Variance components:
            Column    Variance Std.Dev.   Corr.
subj     (Intercept)  584.25897 24.17145
         days          33.63281  5.79938   .  
Residual              653.11578 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
days          10.4673     1.51931   6.89    <1e-11
──────────────────────────────────────────────────

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

fit(MixedModel, @formula(reaction ~ 1 + days + (1 + days|subj)), sleepstudy,
    contrasts = Dict(:days => DummyCoding()))

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + (1 + days | subj)
   logLik   -2 logLik     AIC       AICc        BIC    
  -805.3993  1610.7985  1742.7985  1821.0640  1953.5337

Variance components:
            Column    Variance  Std.Dev.   Corr.
subj     (Intercept)   956.16843 30.92197
         days: 1       497.04350 22.29447 -0.30
         days: 2       915.85978 30.26318 -0.57 +0.75
         days: 3      1267.14813 35.59702 -0.37 +0.72 +0.87
         days: 4      1484.70937 38.53193 -0.32 +0.58 +0.67 +0.91
         days: 5      2296.94899 47.92650 -0.25 +0.46 +0.45 +0.70 +0.85
         days: 6      3849.41364 62.04364 -0.28 +0.30 +0.48 +0.70 +0.77 +0.75
         days: 7      1805.20141 42.48766 -0.16 +0.22 +0.47 +0.50 +0.63 +0.64 +0.71
         days: 8      3153.67996 56.15763 -0.20 +0.29 +0.36 +0.56 +0.73 +0.90 +0.73 +0.74
         days: 9      3075.28984 55.45530 +0.05 +0.25 +0.16 +0.38 +0.59 +0.78 +0.38 +0.53 +0.85
Residual                19.24598  4.38702
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  256.652       7.36136  34.86    <1e-99
days: 1        7.84395     5.45454   1.44    0.1504
days: 2        8.71009     7.28145   1.20    0.2316
days: 3       26.3402      8.51678   3.09    0.0020
days: 4       31.9976      9.19904   3.48    0.0005
days: 5       51.8666     11.3906    4.55    <1e-05
days: 6       55.5264     14.6968    3.78    0.0002
days: 7       62.0988     10.1206    6.14    <1e-09
days: 8       79.9777     13.317     6.01    <1e-08
days: 9       94.1994     13.1525    7.16    <1e-12
───────────────────────────────────────────────────

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

fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj) + (days|subj)), sleepstudy,
    contrasts = Dict(:days => DummyCoding()))

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + (1 | subj) + (days | subj)
   logLik   -2 logLik     AIC       AICc        BIC    
  -827.7748  1655.5496  1769.5496  1823.7463  1951.5482

Variance components:
            Column    Variance Std.Dev.  Corr.
subj     (Intercept)   789.8157 28.1037
         days: 1         0.0000  0.0000   .  
         days: 2       357.6745 18.9123   .     NaN
         days: 3       684.5752 26.1644   .     NaN +0.81
         days: 4       949.5146 30.8142   .     NaN +0.57 +0.91
         days: 5      1752.0021 41.8569   .     NaN +0.26 +0.66 +0.87
         days: 6      3355.6664 57.9281   .     NaN +0.45 +0.72 +0.80 +0.76
         days: 7      1538.9135 39.2290   .     NaN +0.39 +0.42 +0.59 +0.62 +0.71
         days: 8      2736.2930 52.3096   .     NaN +0.22 +0.52 +0.75 +0.93 +0.72 +0.75
         days: 9      2768.2331 52.6140   .     NaN -0.05 +0.28 +0.57 +0.80 +0.34 +0.52 +0.87
Residual               135.0177 11.6197
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  256.652       7.16796  35.81    <1e-99
days: 1        7.84395     3.87324   2.03    0.0429
days: 2        8.71009     5.90532   1.47    0.1402
days: 3       26.3402      7.28244   3.62    0.0003
days: 4       31.9976      8.23121   3.89    0.0001
days: 5       51.8667     10.5988    4.89    <1e-06
days: 6       55.5265     14.1925    3.91    <1e-04
days: 7       62.0988     10.0248    6.19    <1e-09
days: 8       79.9777     12.9236    6.19    <1e-09
days: 9       94.1994     12.992     7.25    <1e-12
───────────────────────────────────────────────────

(Notice that the variance component for days: 1 is estimated as zero, so the correlations for this component are undefined and expressed as NaN, not a number.)

An alternative is to force all the levels of days as indicators using fulldummy encoding.

MixedModels.fulldummy — Function

fulldummy(term::CategoricalTerm)

Assign "contrasts" that include all indicator columns (dummy variables) and an intercept column.

This will result in an under-determined set of contrasts, which is not a problem in the random effects because of the regularization, or "shrinkage", of the conditional modes.

The interaction of fulldummy with complex random effects is subtle and complex with numerous potential edge cases. As we discover these edge cases, we will document and determine their behavior. Until such time, please check the model summary to verify that the expansion is working as you expected. If it is not, please report a use case on GitHub.

source

fit(MixedModel, @formula(reaction ~ 1 + days + (1 + fulldummy(days)|subj)), sleepstudy,
    contrasts = Dict(:days => DummyCoding()))

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + (1 + days | subj)
   logLik   -2 logLik     AIC       AICc        BIC    
  -805.3991  1610.7982  1764.7982  1882.5630  2010.6559

Variance components:
            Column    Variance  Std.Dev.   Corr.
subj     (Intercept)   691.51895 26.29675
         days: 0       472.08523 21.72752 -0.18
         days: 1       427.62891 20.67919 -0.08 +0.45
         days: 2       426.31157 20.64731 -0.29 -0.02 +0.53
         days: 3       358.23780 18.92717 +0.36 -0.53 +0.21 +0.60
         days: 4       351.41637 18.74610 +0.66 -0.81 -0.26 -0.09 +0.66
         days: 5      1155.74302 33.99622 +0.37 -0.45 -0.17 -0.22 +0.26 +0.67
         days: 6      2377.88790 48.76359 +0.27 -0.48 -0.36 -0.07 +0.38 +0.60 +0.56
         days: 7      1182.36885 34.38559 +0.26 -0.10 -0.19 +0.07 -0.01 +0.22 +0.37 +0.50
         days: 8      1983.74678 44.53927 +0.31 -0.36 -0.29 -0.24 +0.10 +0.50 +0.83 +0.56 +0.57
         days: 9      2395.74371 48.94633 +0.44 -0.10 -0.06 -0.31 -0.05 +0.38 +0.69 +0.10 +0.35 +0.80
Residual                19.13253  4.37407
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  256.652       7.35933  34.87    <1e-99
days: 1        7.84395     5.4559    1.44    0.1505
days: 2        8.71009     7.28391   1.20    0.2318
days: 3       26.3402      8.52324   3.09    0.0020
days: 4       31.9976      9.20578   3.48    0.0005
days: 5       51.8666     11.3968    4.55    <1e-05
days: 6       55.5264     14.7108    3.77    0.0002
days: 7       62.0988     10.1325    6.13    <1e-09
days: 8       79.9777     13.3263    6.00    <1e-08
days: 9       94.1994     13.1581    7.16    <1e-12
───────────────────────────────────────────────────

This fit produces a better fit as measured by the objective (negative twice the log-likelihood is 1610.8) but at the expense of adding many more parameters to the model. As a result, model comparison criteria such, as AIC and BIC, are inflated.

But using zerocorr on the individual terms does remove the correlations between the levels:

fit(MixedModel, @formula(reaction ~ 1 + days + zerocorr(1 + days|subj)), sleepstudy,
    contrasts = Dict(:days => DummyCoding()))

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + MixedModels.ZeroCorr((1 + days | subj))
   logLik   -2 logLik     AIC       AICc        BIC    
  -882.9138  1765.8276  1807.8276  1813.6757  1874.8797

Variance components:
            Column    Variance Std.Dev.  Corr.
subj     (Intercept)   958.5617 30.9606
         days: 1         0.0000  0.0000   .  
         days: 2         0.0000  0.0000   .     .  
         days: 3         0.0000  0.0000   .     .     .  
         days: 4         0.0000  0.0000   .     .     .     .  
         days: 5       519.6940 22.7968   .     .     .     .     .  
         days: 6      1703.8588 41.2778   .     .     .     .     .     .  
         days: 7       608.8685 24.6753   .     .     .     .     .     .     .  
         days: 8      1273.1085 35.6806   .     .     .     .     .     .     .     .  
         days: 9      1753.9456 41.8801   .     .     .     .     .     .     .     .     .  
Residual               434.8523 20.8531
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  256.652       8.7984   29.17    <1e-99
days: 1        7.84395     6.95104   1.13    0.2591
days: 2        8.71009     6.95104   1.25    0.2102
days: 3       26.3402      6.95104   3.79    0.0002
days: 4       31.9976      6.95104   4.60    <1e-05
days: 5       51.8666      8.78572   5.90    <1e-08
days: 6       55.5264     11.9572    4.64    <1e-05
days: 7       62.0988      9.06328   6.85    <1e-11
days: 8       79.9777     10.9108    7.33    <1e-12
days: 9       94.1994     12.073     7.80    <1e-14
───────────────────────────────────────────────────

fit(MixedModel, @formula(reaction ~ 1 + days + (1|subj) + zerocorr(days|subj)), sleepstudy,
    contrasts = Dict(:days => DummyCoding()))

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + (1 | subj) + MixedModels.ZeroCorr((days | subj))
   logLik   -2 logLik     AIC       AICc        BIC    
  -882.9138  1765.8276  1807.8276  1813.6757  1874.8797

Variance components:
            Column    Variance Std.Dev.  Corr.
subj     (Intercept)   958.5617 30.9606
         days: 1         0.0000  0.0000   .  
         days: 2         0.0000  0.0000   .     .  
         days: 3         0.0000  0.0000   .     .     .  
         days: 4         0.0000  0.0000   .     .     .     .  
         days: 5       519.6940 22.7968   .     .     .     .     .  
         days: 6      1703.8588 41.2778   .     .     .     .     .     .  
         days: 7       608.8685 24.6753   .     .     .     .     .     .     .  
         days: 8      1273.1085 35.6806   .     .     .     .     .     .     .     .  
         days: 9      1753.9456 41.8801   .     .     .     .     .     .     .     .     .  
Residual               434.8523 20.8531
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  256.652       8.7984   29.17    <1e-99
days: 1        7.84395     6.95104   1.13    0.2591
days: 2        8.71009     6.95104   1.25    0.2102
days: 3       26.3402      6.95104   3.79    0.0002
days: 4       31.9976      6.95104   4.60    <1e-05
days: 5       51.8666      8.78572   5.90    <1e-08
days: 6       55.5264     11.9572    4.64    <1e-05
days: 7       62.0988      9.06328   6.85    <1e-11
days: 8       79.9777     10.9108    7.33    <1e-12
days: 9       94.1994     12.073     7.80    <1e-14
───────────────────────────────────────────────────

fit(MixedModel, @formula(reaction ~ 1 + days + zerocorr(1 + fulldummy(days)|subj)), sleepstudy,
    contrasts = Dict(:days => DummyCoding()))

Linear mixed model fit by maximum likelihood
 reaction ~ 1 + days + MixedModels.ZeroCorr((1 + days | subj))
   logLik   -2 logLik     AIC       AICc        BIC    
  -878.9843  1757.9686  1801.9686  1808.4145  1872.2137

Variance components:
            Column     Variance    Std.Dev.    Corr.
subj     (Intercept)  1135.3845922 33.6954684
         days: 0       776.0604926 27.8578623   .  
         days: 1       357.7267988 18.9136670   .     .  
         days: 2       221.1165324 14.8699876   .     .     .  
         days: 3         0.0024532  0.0495301   .     .     .     .  
         days: 4        44.5173907  6.6721354   .     .     .     .     .  
         days: 5       670.5498913 25.8949781   .     .     .     .     .     .  
         days: 6      1740.1206146 41.7147530   .     .     .     .     .     .     .  
         days: 7       908.9721651 30.1491652   .     .     .     .     .     .     .     .  
         days: 8      1458.0940000 38.1849971   .     .     .     .     .     .     .     .     .  
         days: 9      2028.4320380 45.0381176   .     .     .     .     .     .     .     .     .     .  
Residual               180.8469951 13.4479365
 Number of obs: 180; levels of grouping factors: 18

  Fixed-effects parameters:
───────────────────────────────────────────────────
                 Coef.  Std. Error      z  Pr(>|z|)
───────────────────────────────────────────────────
(Intercept)  256.652      10.7814   23.81    <1e-99
days: 1        7.84395     9.11495   0.86    0.3895
days: 2        8.71009     8.68866   1.00    0.3161
days: 3       26.3402      7.95039   3.31    0.0009
days: 4       31.9976      8.10443   3.95    <1e-04
days: 5       51.8667     10.023     5.17    <1e-06
days: 6       55.5265     12.6444    4.39    <1e-04
days: 7       62.0988     10.6634    5.82    <1e-08
days: 8       79.9777     12.0089    6.66    <1e-10
days: 9       94.1994     13.2627    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.

verbagg = MixedModels.dataset(:verbagg)
verbaggform = @formula(r2 ~ 1 + anger + gender + btype + situ + mode + (1|subj) + (1|item));
gm1 = fit(MixedModel, verbaggform, verbagg, Bernoulli())

Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
  r2 ~ 1 + anger + gender + btype + situ + mode + (1 | subj) + (1 | item)
  Distribution: Bernoulli{Float64}
  Link: LogitLink()


   logLik    deviance     AIC       AICc        BIC    
 -4067.9164  8135.8329  8153.8329  8153.8566  8216.2370
Variance components:
        Column   Variance Std.Dev. 
subj (Intercept)  1.793471 1.339205
item (Intercept)  0.117137 0.342253

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

Fixed-effects parameters:
──────────────────────────────────────────────────────
                   Coef.  Std. Error       z  Pr(>|z|)
──────────────────────────────────────────────────────
(Intercept)   -0.152791    0.385217    -0.40    0.6916
anger          0.0573903   0.0167527    3.43    0.0006
gender: M      0.320786    0.191206     1.68    0.0934
btype: scold  -1.05991     0.184151    -5.76    <1e-08
btype: shout  -2.10398     0.186509   -11.28    <1e-28
situ: self    -1.05438     0.151188    -6.97    <1e-11
mode: want     0.70702     0.151001     4.68    <1e-05
──────────────────────────────────────────────────────

The canonical link, which is 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. The optional arguments fast and/or nAGQ can be passed to the optimization process via both fit and fit! (i.e these optimization settings are not used nor recognized when constructing the model).

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

gm1a = fit(MixedModel, verbaggform, verbagg, Bernoulli(), fast = true)
deviance(gm1a) - deviance(gm1)
@benchmark fit(MixedModel, $verbaggform, $verbagg, Bernoulli())
@benchmark fit(MixedModel, $verbaggform, $verbagg, Bernoulli(), fast = true)

BenchmarkTools.Trial: 
  memory estimate:  11.63 MiB
  allocs estimate:  83551
  --------------
  minimum time:     434.721 ms (0.00% GC)
  median time:      442.984 ms (0.00% GC)
  mean time:        444.607 ms (0.00% GC)
  maximum time:     468.463 ms (0.00% GC)
  --------------
  samples:          12
  evals/sample:     1

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

contraception = MixedModels.dataset(:contra)
contraform = @formula(use ~ 1 + age + abs2(age) + livch + urban + (1|dist));
bernoulli = Bernoulli()
deviances = Dict{Symbol,Float64}()
b = @benchmarkable deviances[:default] = deviance(fit(MixedModel, $contraform, $contraception, $bernoulli));
run(b)
b = @benchmarkable deviances[:fast] = deviance(fit(MixedModel, $contraform, $contraception, $bernoulli, fast = true));
run(b)
b = @benchmarkable deviances[:nAGQ] = deviance(fit(MixedModel, $contraform, $contraception, $bernoulli, nAGQ=9));
run(b)
b = @benchmarkable deviances[:nAGQ_fast] = deviance(fit(MixedModel, $contraform, $contraception, $bernoulli, nAGQ=9, fast=true));
run(b)
sort(deviances)

OrderedCollections.OrderedDict{Symbol, Float64} with 4 entries:
  :default   => 2372.73
  :fast      => 2372.78
  :nAGQ      => 2372.46
  :nAGQ_fast => 2372.51

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(model::StatisticalModel)

Return the log-likelihood of the model.

StatsBase.aic — Function

aic(model::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(model::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(model::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(model::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.

loglikelihood(fm1)

-163.66352994057004

aic(fm1)

333.3270598811401

bic(fm1)

337.5306520261266

dof(fm1)   # 1 fixed effect, 2 variances

nobs(fm1)  # 30 observations

loglikelihood(gm1)

-4067.916430206548

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

StatsBase.deviance — Method

deviance(model::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.

objective(fm1)

327.3270598811401

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

MixedModels.deviance!(gm1)

8135.832860413099

Fixed-effects parameter estimates

The coef and fixef extractors both return the maximum likelihood estimates of the fixed-effects coefficients. They differ in their behavior in the rank-deficient case. The associated coefnames and fixefnames return the corresponding coefficient names.

StatsBase.coef — Function

coef(model::StatisticalModel)

Return the coefficients of the model.

StatsBase.coefnames — Function

coefnames(model::StatisticalModel)

Return the names of the coefficients.

coefnames(term::AbstractTerm)

Return the name(s) of column(s) generated by a term. Return value is either a String or an iterable of Strings.

MixedModels.fixef — Function

fixef(m::MixedModel)

Return the fixed-effects parameter vector estimate of m.

In the rank-deficient case the truncated parameter vector, of length rank(m) is returned. This is unlike coef which always returns a vector whose length matches the number of columns in X.

source

MixedModels.fixefnames — Function

fixefnames(m::MixedModel)

Return a (permuted and truncated in the rank-deficient case) vector of coefficient names.

source

coef(fm1)
coefnames(fm1)

1-element Vector{String}:
 "(Intercept)"

fixef(fm1)
fixefnames(fm1)

1-element Vector{String}:
 "(Intercept)"

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.

fm1.β

1-element Vector{Float64}:
 1527.4999999999993

fm1.beta

1-element Vector{Float64}:
 1527.4999999999993

gm1.β

7-element Vector{Float64}:
 -0.15279123679334106
  0.05739028136275959
  0.3207863096008982
 -1.059914203911305
 -2.10398032428924
 -1.0543786940179338
  0.7070198407561885

A full list of property names is returned by propertynames

propertynames(fm1)

(:formula, :sqrtwts, :A, :L, :optsum, :θ, :theta, :β, :beta, :λ, :lambda, :stderror, :σ, :sigma, :σs, :sigmas, :b, :u, :lowerbd, :X, :y, :corr, :vcov, :PCA, :rePCA, :reterms, :feterms, :allterms, :objective, :pvalues)

propertynames(gm1)

(:A, :L, :theta, :beta, :coef, :fixef, :λ, :lambda, :σ, :sigma, :X, :y, :lowerbd, :objective, :σρs, :σs, :corr, :vcov, :PCA, :rePCA, :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(model::StatisticalModel)

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

vcov(m::MixedModel; corr=false)

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

source

vcov(fm2)

2×2 Matrix{Float64}:
 43.9868   -1.37039
 -1.37039   2.25671

vcov(gm1)

7×7 Matrix{Float64}:
  0.148392    -0.00561483   -0.00982334   …  -0.0112063    -0.011346
 -0.00561483   0.000280653   7.19115e-5      -1.47963e-5    1.02409e-5
 -0.00982334   7.19115e-5    0.0365599       -8.04431e-5    5.25884e-5
 -0.0167974   -1.43708e-5   -9.25637e-5       0.000265804  -0.0001721
 -0.0166214   -2.90551e-5   -0.000162392      0.000658969  -0.000520526
 -0.0112063   -1.47963e-5   -8.04431e-5   …   0.0228577    -0.000247783
 -0.011346     1.02409e-5    5.25884e-5      -0.000247783   0.0228012

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(model::StatisticalModel)

Return the standard errors for the coefficients of the model.

stderror(fm2)

2-element Vector{Float64}:
 6.63225782531458
 1.502235453639816

stderror(gm1)

7-element Vector{Float64}:
 0.3852173194197659
 0.016752703440790974
 0.19120637608556615
 0.18415061521167464
 0.1865094286245418
 0.1511877797756335
 0.151000703244838

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(model::StatisticalModel; level::Real=0.95)

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

The returned CoefTable object implements the Tables.jl interface, and can be converted e.g. to a DataFrame via using DataFrames; DataFrame(coeftable(model)).

coeftable(fm2)

──────────────────────────────────────────────────
                Coef.  Std. Error      z  Pr(>|z|)
──────────────────────────────────────────────────
(Intercept)  251.405      6.63226  37.91    <1e-99
days          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 per-observation dispersion parameter

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

source

VarCorr(fm2)

Variance components:
            Column    Variance Std.Dev.   Corr.
subj     (Intercept)  565.51069 23.78047
         days          32.68212  5.71683 +0.08
Residual              654.94145 25.59182

VarCorr(gm1)

Variance components:
        Column   Variance Std.Dev. 
subj (Intercept)  1.793471 1.339205
item (Intercept)  0.117137 0.342253

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

varest(m::GeneralizedLinearMixedModel)

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

For models with a dispersion parameter ϕ, this is simply ϕ². For models without a dispersion parameter, this value is missing. This differs from disperion, which returns 1 for models without a dispersion parameter.

For Gaussian models, this parameter is often called σ².

source

MixedModels.sdest — Function

sdest(m::LinearMixedModel)

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

source

sdest(m::GeneralizedLinearMixedModel)

Return the estimate of the dispersion, i.e. the standard deviation of the per-observation noise.

For models with a dispersion parameter ϕ, this is simply ϕ. For models without a dispersion parameter, this value is missing. This differs from disperion, which returns 1 for models without a dispersion parameter.

For Gaussian models, this parameter is often called σ.

source

varest(fm2)

654.9414513956141

sdest(fm2)

25.59182391693906

fm2.σ

25.59182391693906

Conditional modes of the random effects

The ranef extractor

MixedModels.ranef — Function

ranef(m::MixedModel; uscale=false)

Return, as a Vector{Matrix{T}}, 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.

For a named variant, see @raneftables.

source

ranef(fm1)

1-element Vector{Matrix{Float64}}:
 [-16.62822143006434 0.36951603177972425 … 53.57982460798641 -42.49434365460919]

fm1.b

1-element Vector{Matrix{Float64}}:
 [-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.

To obtain tables associating the values of the conditional modes with the levels of the grouping factor, use

MixedModels.raneftables — Function

raneftables(m::LinearMixedModel; uscale = false)

Return the conditional means of the random effects as a NamedTuple of columntables

source

as in

DataFrame(only(raneftables(fm1)))

6 rows × 2 columns

	batch	(Intercept)
	String	Float64
1	A	-16.6282
2	B	0.369516
3	C	26.9747
4	D	-21.8014
5	E	53.5798
6	F	-42.4943

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

condVar(fm1)

1-element Vector{Array{Float64, 3}}:
 [362.3104715146577]

[362.3104715146577]

[362.3104715146577]

[362.3104715146577]

[362.3104715146577]

[362.3104715146577]

Case-wise diagnostics and residual degrees of freedom

The leverage values

StatsBase.leverage — Function

leverage(model::RegressionModel)

Return the diagonal of the projection matrix of the model.

leverage(fm1)

30-element Vector{Float64}:
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 ⋮
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486
 0.15650534392640486

are used in diagnostics for linear regression models to determine cases that exert a strong influence on their own predicted response.

The documentation refers to a "projection". For a linear model without random effects the fitted values are obtained by orthogonal projection of the response onto the column span of the model matrix and the sum of the leverage values is the dimension of this column span. That is, the sum of the leverage values is the rank of the model matrix and n - sum(leverage(m)) is the degrees of freedom for residuals. The sum of the leverage values is also the trace of the so-called "hat" matrix, H. (The name "hat matrix" reflects the fact that $\hat{\mathbf{y}} = \mathbf{H} \mathbf{y}$. That is, H puts a hat on y.)

For a linear mixed model the sum of the leverage values will be between p, the rank of the fixed-effects model matrix, and p + q where q is the total number of random effects. This number does not represent a dimension (or "degrees of freedom") of a linear subspace of all possible fitted values because the projection is not an orthogonal projection. Nevertheless, it is a reasonable measure of the effective degrees of freedom of the model and n - sum(leverage(m)) can be considered the effective residual degrees of freedom.

For model fm1 the dimensions are

n, p, q, k = size(fm1)

(30, 1, 6, 1)

which implies that the sum of the leverage values should be in the range [1, 7]. The actual value is

sum(leverage(fm1))

4.695160317792145

For model fm2 the dimensions are

n, p, q, k = size(fm2)

(180, 2, 36, 1)

providing a range of [2, 38] for the effective degrees of freedom for the model. The observed value is

sum(leverage(fm2))

28.611525857134506

When a model converges to a singular covariance, such as

fm3 = fit(MixedModel, @formula(yield ~ 1+(1|batch)), MixedModels.dataset(:dyestuff2))

Linear mixed model fit by maximum likelihood
 yield ~ 1 + (1 | batch)
   logLik   -2 logLik     AIC       AICc        BIC    
   -81.4365   162.8730   168.8730   169.7961   173.0766

Variance components:
            Column   Variance Std.Dev.
batch    (Intercept)   0.00000 0.00000
Residual              13.34610 3.65323
 Number of obs: 30; levels of grouping factors: 6

  Fixed-effects parameters:
───────────────────────────────────────────────
              Coef.  Std. Error     z  Pr(>|z|)
───────────────────────────────────────────────
(Intercept)  5.6656    0.666986  8.49    <1e-16
───────────────────────────────────────────────

the effective degrees of freedom is the lower bound.

sum(leverage(fm3))

0.9999999999999998

Models for which the estimates of the variances of the random effects are large relative to the residual variance have effective degrees of freedom close to the upper bound.

fm4 = fit(MixedModel, @formula(diameter ~ 1+(1|plate)+(1|sample)),
    MixedModels.dataset(:penicillin))

Linear mixed model fit by maximum likelihood
 diameter ~ 1 + (1 | plate) + (1 | sample)
   logLik   -2 logLik     AIC       AICc        BIC    
  -166.0942   332.1883   340.1883   340.4761   352.0676

Variance components:
            Column   Variance Std.Dev. 
plate    (Intercept)  0.714979 0.845565
sample   (Intercept)  3.135194 1.770648
Residual              0.302426 0.549933
 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
─────────────────────────────────────────────────

sum(leverage(fm4))

27.46531792571996

Also, a model fit by the REML criterion generally has larger estimates of the variance components and hence a larger effective degrees of freedom.

fm4r = fit(MixedModel, @formula(diameter ~ 1+(1|plate)+(1|sample)),
    MixedModels.dataset(:penicillin), REML=true)

Linear mixed model fit by REML
 diameter ~ 1 + (1 | plate) + (1 | sample)
 REML criterion at convergence: 330.8605889909948

Variance components:
            Column   Variance Std.Dev. 
plate    (Intercept)  0.716908 0.846704
sample   (Intercept)  3.730909 1.931556
Residual              0.302415 0.549923
 Number of obs: 144; levels of grouping factors: 24, 6

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

sum(leverage(fm4r))

27.472361770234787