Retired stats. prof., R and Julia enthusiast

Complexity in fitting Linear Mixed Models

Linear mixed-effects models are increasingly used for the analysis of data from experiments in fields like psychology where several subjects are each exposed to each of several different items. In addition to a response, which here will be assumed to be on a continuous scale, such as a response time, a number of experimental conditions are systematically varied during the experiment. In the language of statistical experimental design the latter variables are called experimental factors whereas factors like Subject and Item are blocking factors. That is, these are known sources of variation that usually are not of interest by themselves but still should be accounted for when looking for systematic variation in the response.

An example data set

The data from experiment 2 in Kronmueller and Barr (2007) are available in .rds (R Data Set) format in the file kb07_exp2_rt.rds in the github repository provided by Dale Barr. Files in this format can be loaded using the RData package for Julia.

Loading the data in Julia

Install some packages not part of the standard nextjournal environment. The version of MixedModels should be at least v"2.1.0".

using Pkg
Pkg.add(["GLM", "MixedModels", "RData", "StatsModels"])

Attach the packages to be used.

using BenchmarkTools, CSV, DataFrames, GLM, MixedModels
using RData, Statistics, StatsBase, StatsModels, StatsPlots
gr();  # plot backend
kb07 = load("/kronmueller-barr-2007/kb07_exp2_rt.rds")
describe(kb07)

The blocking factors are subj and item with 56 and 32 levels respectively. There are three experimental factors each with two levels: spkr (speaker), prec (precedence), and load (cognitive load). The response time, rt_raw, is measured in milliseconds. A few very large values, e.g. the maximum which is nearly 16 seconds, which could skew the results, are truncated in the rt_trunc column. In addition, three erroneously recorded responses (values of 300 ms.) have been dropped, resulting in a slight imbalance in the data.

A table of mean responses and standard deviations for combinations of the experimental factors, as shown in Table 3 of the paper and on the data repository can be reproduced as

cellmeans = by(kb07, [:spkr, :prec, :load], 
  meanRT = :rt_trunc => mean, sdRT = :rt_trunc => std, n = :rt_trunc => length,
  semean = :rt_trunc => x -> std(x)/sqrt(length(x))
)

The data are slightly imbalanced because 3 unrealistically low response times were removed.

An interaction plot of the cell means shows that the main effect of prec is the dominant effect. (Need to fix this plot to show levels of prec in different colors/symbols.)

@df cellmeans scatter(:load, :meanRT, group=:spkr,
  layout=(1,2), link=:both, xlabel="Cognitive load",
  ylabel="Mean response time (ms)", label="")

Loading the data in R

kb07 <- readRDS("/kronmueller-barr-2007/kb07_exp2_rt.rds")
str(kb07)

The positions of the missing observations can be determined from

(subjitemtbl <- xtabs(~ subj + item, kb07))
table(subjitemtbl)

All of the experimental factors vary within subject and within item, as can be verified by examining the frequency tables for the experimental and grouping factors. For example

xtabs(~ spkr + subj, kb07)

Formulating a simple model

Installing the required R package

For R lme4 and its dependencies must be installed

if (is.na(match("lme4", rownames(installed.packages())))) {
   install.packages("lme4", repos="https://cloud.R-project.org")
}
require(lme4, quietly=TRUE)

Formula and model for simple, scalar random effects

A simple model with main-effects for each of the experimental factors and with random effects for subject and for item is described by the formula rt_trunc ~ 1 + spkr + prec + load + (1|subj) + (1|item). In the MixedModels package, which uses the formula specifications from the StatsModels package, a formula must be wrapped in a call to the @formula macro.

f1 = @formula(rt_trunc ~ 1 + spkr + prec + load + (1|subj) + (1|item));

In R a formula can stand alone.

f1 <- rt_trunc ~ 1 + spkr + prec + load + (1|subj) + (1|item)

For the MixedModels package the model is constructed as a LinearMixedModel then fit with a call to fit!

(By convention, the names in Julia of mutating functions, which modify the value of one or more of their arguments, end in ! as a warning to the user that arguments, usually just the first argument, can be overwritten with new values.)

m1 = fit(MixedModel, f1, kb07)
Linear mixed model fit by maximum likelihood rt_trunc ~ 1 + spkr + prec + load + (1 | subj) + (1 | item) logLik -2 logLik AIC BIC -1.44115694×10⁴ 2.88231387×10⁴ 2.88371387×10⁴ 2.88755646×10⁴ Variance components: Column Variance Std.Dev. subj (Intercept) 101529.32 318.63666 item (Intercept) 132314.19 363.75017 Residual 521145.70 721.90422 Number of obs: 1789; levels of grouping factors: 56, 32 Fixed-effects parameters: ─────────────────────────────────────────────────────── Estimate Std.Error z value P(>|z|) ─────────────────────────────────────────────────────── (Intercept) 2369.24 84.3435 28.0903 <1e-99 spkr: old 136.167 34.1367 3.98886 <1e-4 prec: maintain -667.173 34.1367 -19.5441 <1e-84 load: yes 156.709 34.1366 4.59065 <1e-5 ───────────────────────────────────────────────────────

The first fit of such a model can take several seconds because the Just-In-Time (JIT) compiler must analyze and compile a considerable amount of code. (All of the code in the MixedModels package is Julia code.) Subsequent fits of this or similar models are much faster.

The comparable model fit with lme4 in R is

m1 <- lmer(f1, kb07, REML=FALSE)
summary(m1, corr=FALSE)

The estimated coefficients for the fixed-effects are different from those in the Julia fit because the reference level for load is different. Conversion to a factor (CategoricalArray) is done slightly differently.

Assigning contrasts

For two-level experimental factors, such as prec, spkr and load, in a (nearly) balanced design such as this it is an advantage to use a encoding of the levels in the model matrix. This is known as "assigning contrasts" in both R and Julia.

In R one assigns a contrasts specification to the factor itself. The "Helmert" contrast type produces a encoding with two levels.

contrasts(kb07$spkr) <- contr.helmert(2)
contrasts(kb07$prec) <- contr.helmert(2)
contrasts(kb07$load) <- contr.helmert(2)
system.time(m1 <- lmer(f1, kb07, REML=FALSE))
summary(m1, corr=FALSE)

The change in coding results in estimated coefficients (and standard errors) for the experimental factors being half the previous estimate, because of the coding. Furthermore, the (Intercept) estimate now is now a "typical" response over all the conditions. It would be the sample mean if the experiment were completely balanced.

mean(kb07$rt_trunc)

In Julia the contrasts are specified in a dictionary that is passed as an argument to the model constructor.

const HC = HelmertCoding();
const contrasts = Dict(:spkr => HC, :prec => HC, :load=> HC);
m1 = fit(MixedModel, f1, kb07, contrasts=contrasts)
Linear mixed model fit by maximum likelihood rt_trunc ~ 1 + spkr + prec + load + (1 | subj) + (1 | item) logLik -2 logLik AIC BIC -1.44115694×10⁴ 2.88231387×10⁴ 2.88371387×10⁴ 2.88755646×10⁴ Variance components: Column Variance Std.Dev. subj (Intercept) 101529.32 318.63666 item (Intercept) 132314.19 363.75018 Residual 521145.70 721.90422 Number of obs: 1789; levels of grouping factors: 56, 32 Fixed-effects parameters: ──────────────────────────────────────────────────────── Estimate Std.Error z value P(>|z|) ──────────────────────────────────────────────────────── (Intercept) 2182.09 78.9884 27.6254 <1e-99 spkr: old 68.0833 17.0684 3.98886 <1e-4 prec: maintain -333.586 17.0684 -19.5441 <1e-84 load: yes 78.3546 17.0683 4.59065 <1e-5 ────────────────────────────────────────────────────────

The model matrix for the fixed effects is

Int.(m1.X)    # convert to integers for cleaner printing
1789×4 Array{Int64,2}: 1 -1 -1 1 1 1 1 -1 1 1 -1 -1 1 -1 1 -1 1 -1 -1 -1 1 1 1 1 1 1 -1 1 1 -1 1 1 1 -1 -1 1 1 1 1 -1 1 1 -1 -1 1 -1 1 -1 1 -1 -1 -1 ⋮ 1 -1 1 -1 1 -1 -1 -1 1 1 1 1 1 1 -1 1 1 -1 1 1 1 -1 -1 1 1 1 1 -1 1 1 -1 -1 1 -1 1 -1 1 -1 -1 -1 1 1 1 1 1 1 -1 1

An advantage of the encoding is thatis nearly a multiple of the identity matrix. Furthermore, any interactions of these two-level factors will also have an encoding.

Int.(m1.X'm1.X)
4×4 Array{Int64,2}: 1789 -1 -1 3 -1 1789 -3 1 -1 -3 1789 1 3 1 1 1789

A benchmark of this fit shows that it is quite fast - on the order of a few milliseconds.

m1bmk = @benchmark fit(MixedModel, $f1, $kb07, contrasts = $contrasts)
BenchmarkTools.Trial: memory estimate: 846.20 KiB allocs estimate: 3131 -------------- minimum time: 3.319 ms (0.00% GC) median time: 3.531 ms (0.00% GC) mean time: 4.067 ms (8.99% GC) maximum time: 20.826 ms (81.80% GC) -------------- samples: 1227 evals/sample: 1

Model construction versus model optimization

The m1 object is created in the call to the constructor function, LinearMixedModel, then the parameters are optimized or fit in the call to fit!. Usually the process of fitting a model will take longer than creating the numerical representation but, for simple models like this, the creation time can be a significant portion of the overall running time.

bm1construct = @benchmark LinearMixedModel($f1, $kb07, contrasts=$contrasts)
BenchmarkTools.Trial: memory estimate: 802.22 KiB allocs estimate: 1418 -------------- minimum time: 1.530 ms (0.00% GC) median time: 1.654 ms (0.00% GC) mean time: 2.065 ms (15.21% GC) maximum time: 16.490 ms (80.65% GC) -------------- samples: 2410 evals/sample: 1
bm1fit = @benchmark fit!($m1)
BenchmarkTools.Trial: memory estimate: 43.89 KiB allocs estimate: 1712 -------------- minimum time: 1.670 ms (0.00% GC) median time: 1.767 ms (0.00% GC) mean time: 1.903 ms (1.42% GC) maximum time: 19.660 ms (85.89% GC) -------------- samples: 2618 evals/sample: 1

Factors affecting the time to optimize the parameters

The optimization process is summarized in the optsum property of the model.

m1.optsum
Initial parameter vector: [1.0, 1.0] Initial objective value: 28881.249730604657 Optimizer (from NLopt): LN_BOBYQA Lower bounds: [0.0, 0.0] ftol_rel: 1.0e-12 ftol_abs: 1.0e-8 xtol_rel: 0.0 xtol_abs: [1.0e-10, 1.0e-10] initial_step: [0.75, 0.75] maxfeval: -1 Function evaluations: 28 Final parameter vector: [0.4413835700252912, 0.5038759520997933] Final objective value: 28823.13873586132 Return code: FTOL_REACHED

For this model there are two parameters to be optimized because the objective function, negative twice the log-likelihood, can be profiled with respect to all the other parameters. (See section 3 of Bates et al. 2015 for details.) Both these parameters must be non-negative (i.e. both have a lower bound of zero) and both have an initial value of one. After 28 function evaluations an optimum is declared according to the function value tolerance, either in absolute terms or relative to the current value.

The optimization itself has a certain amount of setup and summary time but the majority of the time is spent in the evaluation of the objective - the profiled log-likelihood.

Each function evaluation is of the form

θ1 = m1.θ;
m1objective = @benchmark objective(updateL!(setθ!($m1, $θ1)))
BenchmarkTools.Trial: memory estimate: 1.42 KiB allocs estimate: 57 -------------- minimum time: 52.548 μs (0.00% GC) median time: 54.541 μs (0.00% GC) mean time: 58.507 μs (0.00% GC) maximum time: 2.534 ms (0.00% GC) -------------- samples: 10000 evals/sample: 1

On this machine 28 function evaluations, each taking around 60 microseconds, gives the total function evaluation time of at least 1.7 ms., which is practically all of the time to fit the model.

The majority of the time for the function evaluation for this model is in the call to updateL!

m1update = @benchmark updateL!($m1)
BenchmarkTools.Trial: memory estimate: 320 bytes allocs estimate: 17 -------------- minimum time: 43.119 μs (0.00% GC) median time: 44.110 μs (0.00% GC) mean time: 46.043 μs (0.00% GC) maximum time: 853.151 μs (0.00% GC) -------------- samples: 10000 evals/sample: 1

This is an operation that updates the lower Cholesky factor (often written as L) of a blocked sparse matrix.

There are 4 rows and columns of blocks. The first row and column correspond to the random effects for subject, the second to the random effects for item, the third to the fixed-effects parameters and the fourth to the response. Their sizes and types are

describeblocks(m1)

There are two lower-triangular blocked matrices in the model representation: A with fixed entries determined by the model and data, and L which is updated for each evaluation of the objective function. The type of the A block is given before the size and the type of the L block is after the size. For scalar random effects, generated by a random-effects term like (1|G), the (1,1) block is always diagonal in both A and L. Its size is the number of levels of the grouping factor, G.

Because subject and item are crossed, the (2,1) block of A is dense, as is the (2,1) block of L. The (2,2) block of A is diagonal because, like the (1,1) block, it is generated from a scalar random effects term. However, the (2,2) block of L ends up being dense as a result of "fill-in" in the sparse Cholesky factorization. All the blocks associated with the fixed-effects or the response are stored as dense matrices but their dimensions are (relatively) small.

Increasing the model complexity

In general, adding more terms to a model will increase the time required to fit the model. However, there is a big difference between adding fixed-effects terms and adding complexity to the random effects.

Increasing the complexity of the fixed effects

Adding the two- and three-factor interactions to the fixed-effects terms increases the time required to fit the model.

f2 = @formula(rt_trunc ~ 1 + spkr*prec*load + (1|subj) + (1|item));
m2 = fit(MixedModel, f2, kb07, contrasts=contrasts)

(Notice that none of the interactions are statistically significant.)

m2bmk = @benchmark fit(MixedModel, $f2, $kb07,contrasts=$contrasts)
BenchmarkTools.Trial: memory estimate: 1.65 MiB allocs estimate: 5354 -------------- minimum time: 4.997 ms (0.00% GC) median time: 5.276 ms (0.00% GC) mean time: 6.250 ms (11.08% GC) maximum time: 33.497 ms (70.26% GC) -------------- samples: 799 evals/sample: 1

In this case, the increase in fitting time is more because the number of function evaluations to determine the optimum increases than because of increased evaluation time for the objective function.

m2.optsum.feval
32
θ2 = m2.θ;
m2objective = @benchmark objective(updateL!(setθ!($m2, $θ2)))
BenchmarkTools.Trial: memory estimate: 1.42 KiB allocs estimate: 57 -------------- minimum time: 57.826 μs (0.00% GC) median time: 61.344 μs (0.00% GC) mean time: 65.811 μs (0.00% GC) maximum time: 2.933 ms (0.00% GC) -------------- samples: 10000 evals/sample: 1
f2 <- rt_trunc ~ 1 + spkr*prec*load + (1|subj) + (1|item)
system.time(m2 <- lmer(f2, kb07, REML=FALSE))
summary(m2, corr=FALSE)

Increasing complexity of the random effects

Another way in which the model can be extended is to switch to vector-valued random effects. Sometimes this is described as having random slopes, so that a subject not only brings their own shift in the typical response but also their own shift in the change due to, say, Load versus No Load. Instead of just one, scalar, change associated with each subject there is an entire vector of changes in the coefficients.

A model with a random slopes for each of the experimental factors for both subject and item is specified as

f3 = @formula(
  rt_trunc ~ 1 + spkr*prec*load + (1+spkr+prec+load|subj) +
    (1+spkr+prec+load|item)
  );
m3 = fit(MixedModel, f3, kb07, contrasts=contrasts)
m3bmk = @benchmark fit(MixedModel, $f3, $kb07, contrasts=$contrasts)
BenchmarkTools.Trial: memory estimate: 22.32 MiB allocs estimate: 367842 -------------- minimum time: 670.530 ms (1.28% GC) median time: 691.396 ms (1.31% GC) mean time: 699.907 ms (1.29% GC) maximum time: 753.922 ms (1.14% GC) -------------- samples: 8 evals/sample: 1
f3 <- 
  rt_trunc ~ 1 + spkr*prec*load + (1+spkr+prec+load|subj) + 
		(1+spkr+prec+load|item)
system.time(m3 <- lmer(f3, kb07, REML=FALSE,
                       control=lmerControl(calc.derivs=FALSE)))
summary(m3, corr=FALSE)

There are several interesting aspects of this model fit.

First, the number of parameters optimized directly has increased substantially. What was previously a 2-dimensional optimization has now become 20 dimensional.

m3.optsum

and the number of function evaluations to convergence has gone from under 40 to over 600.

The time required for each function evaluation has also increased considerably,

θ3 = m3.θ;
m3objective = @benchmark objective(updateL!(setθ!($m3, $θ3)))
BenchmarkTools.Trial: memory estimate: 29.98 KiB allocs estimate: 559 -------------- minimum time: 950.854 μs (0.00% GC) median time: 990.975 μs (0.00% GC) mean time: 1.064 ms (1.13% GC) maximum time: 12.447 ms (90.22% GC) -------------- samples: 4672 evals/sample: 1

resulting in much longer times for model fitting - about three-quarters of a second in Julia and over 6 seconds in R.

Notice that the estimates of the fixed-effects coefficients and their standard errors have not changed substantially except for the standard error of prec, which is also the largest effect.

The parameters in the optimization for this model can be arranged as two lower-triangular 4 by 4 matrices.

m3.λ[1]
4×4 LinearAlgebra.LowerTriangular{Float64,Array{Float64,2}}: 0.451492 ⋅ ⋅ ⋅ 0.0502586 0.040191 ⋅ ⋅ -0.0552264 0.0726916 0.017086 ⋅ 0.0351312 0.0850068 0.0320935 0.0
m3.λ[2]
4×4 LinearAlgebra.LowerTriangular{Float64,Array{Float64,2}}: 0.541793 ⋅ ⋅ ⋅ 0.026789 0.0543733 ⋅ ⋅ -0.253091 0.268958 0.0 ⋅ 0.019991 0.00613279 -0.0456205 0.0381363

which generate the covariance matrices for the random effects. The cumulative proportion of the variance in the principal components of these covariance matrices, available as

m3.rePCA
(subj = [0.9317024826383811, 0.9991564337515368, 1.0, 1.0], item = [0.8536519268476426, 0.991942454066671, 1.0, 1.0])

show that 93% of the variation in the random effects for subject is in the first principal direction and 99% in the first two principal directions. The random effects for item also have 99% of the variation in the first two principal directions.

Furthermore the estimates of the standard deviations of the "slope" random effects are much smaller than the those of the intercept random effects except for the prec coefficient random effect for item, which suggests that the model could be reduced to rt_trunc ~ 1 + spkr*prec*load + (1|subj) + (1+prec|item) or even rt_trunc ~ 1 + spkr+prec+load + (1|subj) + (1+prec|item).

f4 = @formula(rt_trunc ~ 1 + spkr+prec+load + (1|subj) + (1+prec|item));
m4 = fit(MixedModel, f4, kb07, contrasts=contrasts)
m4bmk = @benchmark fit(MixedModel, $f4, $kb07, contrasts=$contrasts)
BenchmarkTools.Trial: memory estimate: 2.54 MiB allocs estimate: 35776 -------------- minimum time: 19.224 ms (0.00% GC) median time: 20.152 ms (0.00% GC) mean time: 21.900 ms (5.27% GC) maximum time: 45.768 ms (0.00% GC) -------------- samples: 229 evals/sample: 1
m4.optsum.feval
93
f4 <- rt_trunc ~ 1 + spkr+prec+load + (1+prec|item) + (1|subj)
system.time(m4 <- lmer(f4, kb07, REML=FALSE,  control=lmerControl(calc.derivs=FALSE)))
summary(m4, corr=FALSE)
θ4 = m4.θ;
m4objective = @benchmark objective(updateL!(setθ!($m4, $θ4)))
BenchmarkTools.Trial: memory estimate: 14.30 KiB allocs estimate: 289 -------------- minimum time: 157.573 μs (0.00% GC) median time: 171.155 μs (0.00% GC) mean time: 186.610 μs (3.18% GC) maximum time: 10.499 ms (96.83% GC) -------------- samples: 10000 evals/sample: 1

These two model fits can be compared with one of the information criteria, AIC or BIC, for which "smaller is better". They both indicate a preference for the smaller model, m4.

These criteria are values of the objective, negative twice the log-likelihood at convergence, plus a penalty that depends on the number of parameters being estimated.

Because model m4 is a special case of model m3, a likelihood ratio test could also be used. The alternative model, m3, will always produce an objective that is less than or equal to that from the null model, m4. The difference in these value is similar to the change in the residual sum of squares in a linear model fit. This objective would be called the deviance if there was a way of defining a saturated model but it is not clear what this should be. However, if there was a way to define a deviance then the difference in the deviances would be the same as the differences in these objectives, which is

diff(objective.([m3, m4]))
1-element Array{Float64,1}: 26.7107862699595

This difference is compared to a distribution with degrees of freedom corresponding to the difference in the number of parameters

diff(dof.([m4, m3]))
1-element Array{Int64,1}: 20

producing a p-value of about 14%.

pchisq(26.7108, 20, lower.tail=FALSE)

Going maximal

The "maximal" model as proposed by Barr et al., 2013 would include all possible interactions of experimental and grouping factors.

f5 = @formula(rt_trunc ~ 1 + spkr*prec*load + (1+spkr*prec*load|subj) +      (1+spkr*prec*load|item));
m5 = fit(MixedModel, f5, kb07, contrasts=contrasts)
m5bmk = @benchmark fit(MixedModel, $f5, $kb07, contrasts=$contrasts)
BenchmarkTools.Trial: memory estimate: 93.77 MiB allocs estimate: 1638022 -------------- minimum time: 12.313 s (0.30% GC) median time: 12.313 s (0.30% GC) mean time: 12.313 s (0.30% GC) maximum time: 12.313 s (0.30% GC) -------------- samples: 1 evals/sample: 1

As is common in models with high-dimensional vector-valued random effects, the dominant portion of the variation is in the first few principal components

m5.rePCA
(subj = [0.7252516043443922, 0.8471888880119028, 0.9411190095298798, 0.9999933288201127, 0.9999999801038078, 0.9999999989133693, 1.0, 1.0], item = [0.7968023853414774, 0.9418350206829007, 0.9750907771940907, 0.9961087873841261, 0.9999942056296, 0.999999881579346, 1.0, 1.0])

For both the subjects and the items practically all the variation of these 8-dimensional random effects is in the first 4 principal components.

The dimension of , the parameters in the optimization, increases considerably

θ5 = m5.θ;
length(θ5)
72

Of these 72 parameters, 36 are estimated from variation between items, yet there are only 32 items.

Because the dimension of the optimization problem has gotten much larger the number of function evaluations to convergence increases correspondingly.

m5.optsum.feval
2862

Also, each function evaluation requires more time

m5objective = @benchmark objective(updateL!(setθ!($m5, $θ5)))
BenchmarkTools.Trial: memory estimate: 30.02 KiB allocs estimate: 561 -------------- minimum time: 3.418 ms (0.00% GC) median time: 3.500 ms (0.00% GC) mean time: 3.671 ms (0.19% GC) maximum time: 14.367 ms (0.00% GC) -------------- samples: 1359 evals/sample: 1

almost all of which is for the call to updateL!.

@benchmark updateL!($m5)
BenchmarkTools.Trial: memory estimate: 28.88 KiB allocs estimate: 519 -------------- minimum time: 3.374 ms (0.00% GC) median time: 3.478 ms (0.00% GC) mean time: 3.672 ms (0.17% GC) maximum time: 17.212 ms (0.00% GC) -------------- samples: 1359 evals/sample: 1

This model takes a long time to fit using lme4.

f5 <- rt_trunc ~ 1 + spkr*prec*load + (1+spkr*prec*load|subj) +
    (1+spkr*prec*load|item)
system.time(m5 <- lmer(f5, kb07, REML=FALSE,        control=lmerControl(calc.derivs=FALSE)))
summary(m5, corr=FALSE)

A model of intermediate complexity

To provide more granularity in the plots of execution time shown below, fit one more model without random effects for the third-order interaction of the experimental factors.

f6 = @formula(rt_trunc ~ 1 + spkr*prec*load + 
    (1+spkr+prec+load+spkr&prec+spkr&load+prec&load|subj) + 
    (1+spkr+prec+load+spkr&prec+spkr&load+prec&load|item));
m6 = fit(MixedModel, f6, kb07, contrasts=contrasts)
m6bmk = @benchmark fit(MixedModel, $f6, $kb07, contrasts=$contrasts)
BenchmarkTools.Trial: memory estimate: 61.91 MiB allocs estimate: 1057606 -------------- minimum time: 6.796 s (0.26% GC) median time: 6.796 s (0.26% GC) mean time: 6.796 s (0.26% GC) maximum time: 6.796 s (0.26% GC) -------------- samples: 1 evals/sample: 1
θ6 = m6.θ;
length(θ6)
56
m6objective = @benchmark objective(updateL!(setθ!($m6, $θ6)))
BenchmarkTools.Trial: memory estimate: 30.02 KiB allocs estimate: 561 -------------- minimum time: 3.155 ms (0.00% GC) median time: 3.268 ms (0.00% GC) mean time: 3.528 ms (0.38% GC) maximum time: 25.832 ms (74.16% GC) -------------- samples: 1414 evals/sample: 1
@benchmark updateL!($m6)
BenchmarkTools.Trial: memory estimate: 28.88 KiB allocs estimate: 519 -------------- minimum time: 3.142 ms (0.00% GC) median time: 3.290 ms (0.00% GC) mean time: 3.491 ms (0.18% GC) maximum time: 12.883 ms (0.00% GC) -------------- samples: 1428 evals/sample: 1

Summary of goodness of fit

Apply the goodness of fit measures to m1 to m6 creating a data frame

const mods = [m1, m2, m3, m4, m5, m6];
gofsumry = DataFrame(dof=dof.(mods), deviance=deviance.(mods),
    AIC = aic.(mods), AICc = aicc.(mods), BIC = bic.(mods))

6 rows × 5 columns

dofdevianceAICAICcBIC
Int64Float64Float64Float64Float64
1728823.128837.128837.228875.6
21128818.528840.528840.628900.9
32928637.128695.128696.128854.3
4928663.928681.928682.028731.3
58128578.328740.328748.129185.0
66528603.928733.928738.929090.7

Here dof or degrees of freedom is the total number of parameters estimated in the model and deviance is simply negative twice the log-likelihood at convergence, without a correction for a saturated model. All the information criteria are on a scale of "smaller is better" and all would select model 4 as "best".

map(nm -> argmin(getproperty(gofsumry,nm)), (AIC=:AIC, AICc=:AICc, BIC=:BIC))
(AIC = 4, AICc = 4, BIC = 4)

The benchmark times are recorded in nanoseconds. The median time in seconds to fit each model is evaluated as

fits = [median(b.times) for b in [m1bmk, m2bmk, m3bmk, m4bmk, m5bmk, m6bmk]] ./ 10^9;
nfe(m) = length(coef(m));
nre1(m) = MixedModels.vsize(first(m.reterms));
nlv1(m) = MixedModels.nlevs(first(m.reterms));
nre2(m) = MixedModels.vsize(last(m.reterms));
nlv2(m) = MixedModels.nlevs(last(m.reterms));
(m) = sum(MixedModels., m.reterms);
nev = [m.optsum.feval for m in mods];
dimsumry = DataFrame(p = nfe.(mods), q1 = nre1.(mods), n1 = nlv1.(mods),
    q2 = nre2.(mods), n2 = nlv2.(mods),  = .(mods),
    npar = dof.(mods), nev = nev, fitsec = fits)

In this table, p is the dimension of the fixed-effects vector, q1 is the dimension of the random-effects for each of the n1 levels of the first grouping factor while q2 and n2 are similar for the second grouping factor. is the dimension of the parameter vector and npar is the total number of parameters being estimated, and is equal to nθ + p + 1.

nev is the number of function evaluations to convergence. Because this number will depend on the number of parameters in the model and several other factors such as the setting of the convergence criteria, only its magnitude should be considered reproducible. As described above, fitsec is the median time, in seconds, to fit the model.

Relationships between model complexity and fitting time

As would be expected, the time required to fit a model increases with the complexity of the model. In particular, the number of function evaluations to convergence increases with the number of parameters being optimized.

@df dimsumry scatter(:nθ, :nev, xlabel="# of parameters in the optimization",
    ylabel="Function evaluations to convergence", label="")

The relationship is more-or-less linear, based on this small sample.

Finally, consider the time to fit the model versus the number of parameters in the optimization. On a log-log scale these form a reasonably straight line.

@df dimsumry scatter(:nθ, :fitsec, xlabel="# of parameters in the optimization",
  ylabel = "Median benchmark time (sec) to fit model", label="", yscale=:log10,
  xscale=:log10)

with a slope of about 2.2

coef(lm(@formula(log(fitsec) ~ 1 + log()), dimsumry))
2-element Array{Float64,1}: -6.980738425978648 2.2136833389993336

or close to a quadratic relationship between time to fit and the number of parameters in the optimization. However, the number of parameters to optimize is itself quadratic in the dimension of the random effects vector in each random-effects term. Thus, increasing the dimension of the vector-valued random effects can cause a considerable increase in the amount of time required to fit the model.

Furthermore, the estimated covariance matrix for high-dimensional vector-valued random effects is singular in m3, m5 and m6, as is often the case.