{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Complexity in fitting Linear Mixed Models"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"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.\n",
"In the language of statistical experimental design the latter variables are called _experimental factors_ whereas factors like `Subject` and `Item` are _blocking factors_.\n",
"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."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## An example data set"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The data from experiment 2 in [_Kronmueller and Barr (2007)_](https://doi.org/10.1016/j.jml.2006.05.002) are available in `.rds` (R Data Set) format in the file `kb07_exp2_rt.rds` in the [github repository](https://github.com/dalejbarr/kronmueller-barr-2007) provided by Dale Barr.\n",
"Files in this format can be loaded using the [_RData_ package](https://github.com/JuliaData/RData.jl) for [__Julia__](https://julialang.org)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Loading DataFrames support into Gadfly.jl\n",
"└ @ Gadfly /home/bates/.julia/packages/Gadfly/09PWZ/src/mapping.jl:228\n"
]
}
],
"source": [
"using BenchmarkTools, DataFrames, Distributions, FreqTables\n",
"using MixedModels, RData, Statistics, StatsBase, StatsModels\n",
"using Gadfly"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/latex": [
"\\begin{tabular}{r|cccccc}\n",
"\t& spkr & prec & load & meanRT & sdRT & n\\\\\n",
"\t\\hline\n",
"\t& Categorical… & Categorical… & Categorical… & Float64 & Float64 & Int64\\\\\n",
"\t\\hline\n",
"\t1 & new & break & no & 2345.17 & 972.693 & 223 \\\\\n",
"\t2 & new & break & yes & 2505.34 & 1012.04 & 224 \\\\\n",
"\t3 & new & maintain & no & 1760.66 & 587.445 & 224 \\\\\n",
"\t4 & new & maintain & yes & 1843.28 & 609.672 & 224 \\\\\n",
"\t5 & old & break & no & 2536.29 & 1058.78 & 224 \\\\\n",
"\t6 & old & break & yes & 2674.32 & 1119.67 & 224 \\\\\n",
"\t7 & old & maintain & no & 1770.76 & 669.901 & 222 \\\\\n",
"\t8 & old & maintain & yes & 2018.83 & 725.359 & 224 \\\\\n",
"\\end{tabular}\n"
],
"text/plain": [
"8×6 DataFrame\n",
"│ Row │ spkr │ prec │ load │ meanRT │ sdRT │ n │\n",
"│ │ \u001b[90mCategorical…\u001b[39m │ \u001b[90mCategorical…\u001b[39m │ \u001b[90mCategorical…\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mInt64\u001b[39m │\n",
"├─────┼──────────────┼──────────────┼──────────────┼─────────┼─────────┼───────┤\n",
"│ 1 │ new │ break │ no │ 2345.17 │ 972.693 │ 223 │\n",
"│ 2 │ new │ break │ yes │ 2505.34 │ 1012.04 │ 224 │\n",
"│ 3 │ new │ maintain │ no │ 1760.66 │ 587.445 │ 224 │\n",
"│ 4 │ new │ maintain │ yes │ 1843.28 │ 609.672 │ 224 │\n",
"│ 5 │ old │ break │ no │ 2536.29 │ 1058.78 │ 224 │\n",
"│ 6 │ old │ break │ yes │ 2674.32 │ 1119.67 │ 224 │\n",
"│ 7 │ old │ maintain │ no │ 1770.76 │ 669.901 │ 222 │\n",
"│ 8 │ old │ maintain │ yes │ 2018.83 │ 725.359 │ 224 │"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cellmeans = by(kb07, [:spkr, :prec, :load], meanRT = :rt_trunc => mean, sdRT = :rt_trunc => std, n = :rt_trunc => length)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An interaction plot of these means"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(cellmeans, x=:load, y=:meanRT, color=:prec, xgroup=:spkr, Geom.subplot_grid(Geom.point))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"shows that the dominant factor is the precedent with the response times under _break_ being considerably larger than under _maintain_."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Formulating a simple model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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)`.\n",
"In the _MixedModels_ package, which uses the formula specifications from the [_StatsModels_ package](https://github.com/JuliaStats/StatsModels.jl), a formula must be wrapped in a call to the `@formula` macro.\n",
"The model is created as an instance of a `LinearMixedModel` type then fit with a call to the `fit!` generic.\n",
"\n",
"(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.)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Linear mixed model fit by maximum likelihood\n",
" rt_trunc ~ 1 + spkr + prec + load + (1 | subj) + (1 | item)\n",
" logLik -2 logLik AIC BIC \n",
" -1.44115694×10⁴ 2.88231387×10⁴ 2.88371387×10⁴ 2.88755646×10⁴\n",
"\n",
"Variance components:\n",
" Variance Std.Dev. \n",
" subj 101529.32 318.63666\n",
" item 132314.19 363.75018\n",
" Residual 521145.70 721.90422\n",
" Number of obs: 1789; levels of grouping factors: 56, 32\n",
"\n",
" Fixed-effects parameters:\n",
" Estimate Std.Error z value P(>|z|)\n",
"(Intercept) 2369.24 84.3435 28.0903 <1e-99\n",
"spkr: old 136.167 34.1367 3.98886 <1e-4\n",
"prec: maintain -667.173 34.1367 -19.5441 <1e-84\n",
"load: yes 156.709 34.1366 4.59065 <1e-5\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m1 = fit!(LinearMixedModel(@formula(rt_trunc ~ 1 + spkr + prec + load + (1|subj) + (1|item)), kb07))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.)\n",
"Subsequent fits of this or similar models are much faster.\n",
"\n",
"When timing fits that take a few seconds or less the `@btime` macro will be used, as it performs the evaluation of the expression several times and takes the minimum of these evaluation times. \n",
"This is to eliminate the noise introduced by garbage collection and compilation.\n",
"\n",
"Furthermore, in what follows the formulas and models will be declared `const`, which in Julia means that the type of the value cannot be changed, although the values in, say, a vector can be.\n",
"\n",
"Also the expression is written using `$` before an argument name. This is an interpolation in a macro call. The effect is to time only the expression evaluation and not the name lookup for the arguments."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.324 ms (3082 allocations: 841.83 KiB)\n"
]
}
],
"source": [
"const f1 = @formula(rt_trunc ~ 1 + spkr + prec + load + (1|subj) + (1|item));\n",
"@btime fit!(LinearMixedModel($f1, $kb07));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The coding for the experimental factors in the model matrix for this model is the default or \"dummy\" coding, where the first level is the reference level and the estimated coefficients represent differences of the response at the indicated level relative to the first.\n",
"\n",
"In this case the model matrix for the fixed effects is"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1789×4 Array{Float64,2}:\n",
" 1.0 0.0 0.0 1.0\n",
" 1.0 1.0 1.0 0.0\n",
" 1.0 1.0 0.0 0.0\n",
" 1.0 0.0 1.0 0.0\n",
" 1.0 0.0 0.0 0.0\n",
" 1.0 1.0 1.0 1.0\n",
" 1.0 1.0 0.0 1.0\n",
" 1.0 0.0 1.0 1.0\n",
" 1.0 0.0 0.0 1.0\n",
" 1.0 1.0 1.0 0.0\n",
" 1.0 1.0 0.0 0.0\n",
" 1.0 0.0 1.0 0.0\n",
" 1.0 0.0 0.0 0.0\n",
" ⋮ \n",
" 1.0 0.0 1.0 0.0\n",
" 1.0 0.0 0.0 0.0\n",
" 1.0 1.0 1.0 1.0\n",
" 1.0 1.0 0.0 1.0\n",
" 1.0 0.0 1.0 1.0\n",
" 1.0 0.0 0.0 1.0\n",
" 1.0 1.0 1.0 0.0\n",
" 1.0 1.0 0.0 0.0\n",
" 1.0 0.0 1.0 0.0\n",
" 1.0 0.0 0.0 0.0\n",
" 1.0 1.0 1.0 1.0\n",
" 1.0 1.0 0.0 1.0"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m1.X"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In two-level factorial experiments like this it is an advantage to use a $\\pm1$ coding, which can be derived from `EffectsCoding` or `HelmertCoding` (these codings are the same for a two-level factor)."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"const contrasts = Dict(:spkr=>HelmertCoding(), :prec=>HelmertCoding(), :load=>HelmertCoding());"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.397 ms (3300 allocations: 849.63 KiB)\n"
]
},
{
"data": {
"text/plain": [
"Linear mixed model fit by maximum likelihood\n",
" rt_trunc ~ 1 + spkr + prec + load + (1 | subj) + (1 | item)\n",
" logLik -2 logLik AIC BIC \n",
" -1.44115694×10⁴ 2.88231387×10⁴ 2.88371387×10⁴ 2.88755646×10⁴\n",
"\n",
"Variance components:\n",
" Variance Std.Dev. \n",
" subj 101529.32 318.63666\n",
" item 132314.19 363.75018\n",
" Residual 521145.70 721.90422\n",
" Number of obs: 1789; levels of grouping factors: 56, 32\n",
"\n",
" Fixed-effects parameters:\n",
" Estimate Std.Error z value P(>|z|)\n",
"(Intercept) 2182.09 78.9884 27.6254 <1e-99\n",
"spkr: old 68.0833 17.0684 3.98886 <1e-4\n",
"prec: maintain -333.586 17.0684 -19.5441 <1e-84\n",
"load: yes 78.3546 17.0683 4.59065 <1e-5\n"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m1 = @btime fit!(LinearMixedModel($f1, $kb07, $contrasts))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The estimate of the `(Intercept)` coefficient is now approximately the sample mean. (It may be slightly different from the sample mean due to the \n",
"slight imbalance in the design because of two unusually long response times being dropped.)\n",
"\n",
"The coefficients for the experimental factors are half the previously estimated effects, because of the $\\pm1$ coding,\n",
"Their standard errors are also half those from the previous fit.\n",
"The estimates of the variance components, the log-likelihood and other goodness of fit criteria, and the z statistics and p-values for the effects remain the same.\n",
"\n",
"The model matrix for the fixed effects is now"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1789×4 Array{Float64,2}:\n",
" 1.0 -1.0 -1.0 1.0\n",
" 1.0 1.0 1.0 -1.0\n",
" 1.0 1.0 -1.0 -1.0\n",
" 1.0 -1.0 1.0 -1.0\n",
" 1.0 -1.0 -1.0 -1.0\n",
" 1.0 1.0 1.0 1.0\n",
" 1.0 1.0 -1.0 1.0\n",
" 1.0 -1.0 1.0 1.0\n",
" 1.0 -1.0 -1.0 1.0\n",
" 1.0 1.0 1.0 -1.0\n",
" 1.0 1.0 -1.0 -1.0\n",
" 1.0 -1.0 1.0 -1.0\n",
" 1.0 -1.0 -1.0 -1.0\n",
" ⋮ \n",
" 1.0 -1.0 1.0 -1.0\n",
" 1.0 -1.0 -1.0 -1.0\n",
" 1.0 1.0 1.0 1.0\n",
" 1.0 1.0 -1.0 1.0\n",
" 1.0 -1.0 1.0 1.0\n",
" 1.0 -1.0 -1.0 1.0\n",
" 1.0 1.0 1.0 -1.0\n",
" 1.0 1.0 -1.0 -1.0\n",
" 1.0 -1.0 1.0 -1.0\n",
" 1.0 -1.0 -1.0 -1.0\n",
" 1.0 1.0 1.0 1.0\n",
" 1.0 1.0 -1.0 1.0"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m1.X"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Model construction versus model optimization"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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!`.\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 923.172 μs (1583 allocations: 805.44 KiB)\n"
]
}
],
"source": [
"@btime LinearMixedModel($f1, $kb07, $contrasts);"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.370 ms (1716 allocations: 44.09 KiB)\n"
]
}
],
"source": [
"@btime fit!($m1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Factors affecting the time to optimize the parameters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The optimization process is summarized in the `optsum` property of the model."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Initial parameter vector: [1.0, 1.0]\n",
"Initial objective value: 28881.249730604657\n",
"\n",
"Optimizer (from NLopt): LN_BOBYQA\n",
"Lower bounds: [0.0, 0.0]\n",
"ftol_rel: 1.0e-12\n",
"ftol_abs: 1.0e-8\n",
"xtol_rel: 0.0\n",
"xtol_abs: [1.0e-10, 1.0e-10]\n",
"initial_step: [0.75, 0.75]\n",
"maxfeval: -1\n",
"\n",
"Function evaluations: 28\n",
"Final parameter vector: [0.441384, 0.503876]\n",
"Final objective value: 28823.13873586134\n",
"Return code: FTOL_REACHED\n"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m1.optsum"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"(See section 3 of [_Bates et al. 2015_](https://www.jstatsoft.org/article/view/v067i01) for details.)\n",
"Both these parameters must be non-negative (i.e. both have a lower bound of zero) and both have an initial value of one.\n",
"After 21 function evaluations an optimum is declared according to the function value tolerance, either $10^{-8}$ in absolute terms or $10^{-12}$ relative to the current value.\n",
"\n",
"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.\n",
"\n",
"Each function evaluation is of the form"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"const θ1 = m1.θ;"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 42.720 μs (57 allocations: 1.42 KiB)\n"
]
}
],
"source": [
"@btime objective(updateL!(setθ!($m1, $θ1)));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Remember that this is a \"best case\" time. During the course of the optimization some of the function evaluations could take longer due to garbage collection or other factors.\n",
"\n",
"On this machine 21 function evaluations, each taking 42 microseconds or more, gives the total function evaluation time of at least 0.85 ms., which is practically all of the time to fit the model.\n",
"\n",
"The majority of the time for the function evaluation for this model is in the call to `updateL!`"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 32.016 μs (17 allocations: 320 bytes)\n"
]
}
],
"source": [
"@btime updateL!($m1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is an operation that updates the lower Cholesky factor (often written as `L`) of a blocked sparse matrix.\n",
"\n",
"There are 4 rows and columns of blocks.\n",
"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.\n",
"Their sizes and types are"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1,1: LinearAlgebra.Diagonal{Float64,Array{Float64,1}} (56, 56) LinearAlgebra.Diagonal{Float64,Array{Float64,1}}\n",
"2,1: Array{Float64,2} (32, 56) Array{Float64,2}\n",
"2,2: LinearAlgebra.Diagonal{Float64,Array{Float64,1}} (32, 32) Array{Float64,2}\n",
"3,1: Array{Float64,2} (4, 56) Array{Float64,2}\n",
"3,2: Array{Float64,2} (4, 32) Array{Float64,2}\n",
"3,3: Array{Float64,2} (4, 4) Array{Float64,2}\n",
"4,1: Array{Float64,2} (1, 56) Array{Float64,2}\n",
"4,2: Array{Float64,2} (1, 32) Array{Float64,2}\n",
"4,3: Array{Float64,2} (1, 4) Array{Float64,2}\n",
"4,4: Array{Float64,2} (1, 1) Array{Float64,2}\n"
]
}
],
"source": [
"describeblocks(m1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are two lower-triangular blocked matrices: `A` with fixed entries determined by the model and data, and `L` which is updated for each evaluation of the objective function.\n",
"The type of the `A` block is given before the size and the type of the `L` block is after the size.\n",
"For scalar random effects, generated by a random-effects term like `(1|G)`, the (1,1) block is always diagonal for both `A` and `L`.\n",
"Its size is the number of levels of the grouping factor, `G`.\n",
"\n",
"Because subject and item are crossed, the (2,1) block of `A` is dense, as is the (2,1) block of `L`.\n",
"The (2,2) block of `A` is diagonal because, like the (1,1) block, it is generated from a scalar random effects term.\n",
"However, the (2,2) block of `L` ends up being dense as a result of \"fill-in\" in the sparse Cholesky factorization.\n",
"All the blocks associated with the fixed-effects or the response are stored as dense matrices but their dimensions are (relatively) small."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Increasing the complexity"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In general, adding more terms to a model will increase the time required to fit the model.\n",
"However, there is a big difference between adding fixed-effects terms and adding complexity to the random effects.\n",
"\n",
"Adding the two- and three-factor interactions to the fixed-effects terms increases the time required to fit the model."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"const f2 = @formula(rt_trunc ~ 1 + spkr*prec*load + (1|subj) + (1|item));"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 3.902 ms (5685 allocations: 1.65 MiB)\n"
]
},
{
"data": {
"text/plain": [
"Linear mixed model fit by maximum likelihood\n",
" rt_trunc ~ 1 + spkr + prec + load + spkr & prec + spkr & load + prec & load + spkr & prec & load + (1 | subj) + (1 | item)\n",
" logLik -2 logLik AIC BIC \n",
" -1.4409243×10⁴ 2.88184861×10⁴ 2.88404861×10⁴ 2.89008696×10⁴\n",
"\n",
"Variance components:\n",
" Variance Std.Dev. \n",
" subj 101581.912 318.71917\n",
" item 132326.015 363.76643\n",
" Residual 519722.754 720.91799\n",
" Number of obs: 1789; levels of grouping factors: 56, 32\n",
"\n",
" Fixed-effects parameters:\n",
" Estimate Std.Error z value P(>|z|)\n",
"(Intercept) 2182.0 78.9917 27.6231 <1e-99\n",
"spkr: old 68.0177 17.0451 3.99045 <1e-4\n",
"prec: maintain -333.652 17.0451 -19.5747 <1e-84\n",
"load: yes 78.4464 17.0451 4.6023 <1e-5\n",
"spkr: old & prec: maintain -21.6417 17.0451 -1.26968 0.2042\n",
"spkr: old & load: yes 18.1151 17.0451 1.06278 0.2879\n",
"prec: maintain & load: yes 4.2647 17.0451 0.250201 0.8024\n",
"spkr: old & prec: maintain & load: yes 23.2834 17.0451 1.36599 0.1719\n"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"const m2 = @btime fit!(LinearMixedModel($f2, $kb07, $contrasts))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(Notice that none of the interactions are statistically significant.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"32"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m2.optsum.feval"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"const θ2 = m2.θ;"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 48.086 μs (57 allocations: 1.42 KiB)\n"
]
}
],
"source": [
"@btime objective(updateL!(setθ!($m2, $θ2)));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Increasing complexity of the random effects"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another way in which the model can be extended is to switch to vector-valued random effects.\n",
"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`.\n",
"Instead of just one, scalar, change associated with each subject there is an entire vector of changes in the coefficients.\n",
"\n",
"A model with a random slopes for each of the experimental factors for both subject and item is specified as"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"const f3 = @formula(rt_trunc ~ 1 + spkr*prec*load + (1+spkr+prec+load|subj) + (1+spkr+prec+load|item));"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 517.561 ms (715239 allocations: 31.42 MiB)\n"
]
},
{
"data": {
"text/plain": [
"Linear mixed model fit by maximum likelihood\n",
" rt_trunc ~ 1 + spkr + prec + load + spkr & prec + spkr & load + prec & load + spkr & prec & load + (1 + spkr + prec + load | subj) + (1 + spkr + prec + load | item)\n",
" logLik -2 logLik AIC BIC \n",
" -1.43189855×10⁴ 2.8637971×10⁴ 2.8695971×10⁴ 2.8855164×10⁴\n",
"\n",
"Variance components:\n",
" Variance Std.Dev. Corr.\n",
" subj 91130.9567 301.879043\n",
" 1092.7366 33.056567 1.00\n",
" 3463.3273 58.850041 -0.62 -0.62\n",
" 4481.7460 66.945844 0.36 0.36 0.51\n",
" item 131377.4183 362.460230\n",
" 1690.3829 41.114266 0.42\n",
" 61168.0277 247.321709 -0.69 0.37\n",
" 1882.8841 43.392212 0.29 0.14 -0.13\n",
" Residual 447630.8098 669.052173\n",
" Number of obs: 1789; levels of grouping factors: 56, 32\n",
"\n",
" Fixed-effects parameters:\n",
" Estimate Std.Error z value P(>|z|)\n",
"(Intercept) 2181.64 77.3507 28.2046 <1e-99\n",
"spkr: old 67.7496 17.9606 3.77213 0.0002\n",
"prec: maintain -333.92 47.155 -7.08133 <1e-11\n",
"load: yes 78.8007 19.726 3.99476 <1e-4\n",
"spkr: old & prec: maintain -21.996 15.8191 -1.39047 0.1644\n",
"spkr: old & load: yes 18.3832 15.8191 1.16209 0.2452\n",
"prec: maintain & load: yes 4.53273 15.8191 0.286535 0.7745\n",
"spkr: old & prec: maintain & load: yes 23.6377 15.8191 1.49425 0.1351\n"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"const m3 = @btime fit!(LinearMixedModel($f3, $kb07, $contrasts))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are several interesting aspects of this model fit.\n",
"\n",
"First, the number of parameters optimized directly has increased substantially.\n",
"What was previously a 2-dimensional optimization has now become 20 dimensional."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Initial parameter vector: [1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0]\n",
"Initial objective value: 29303.38411137385\n",
"\n",
"Optimizer (from NLopt): LN_BOBYQA\n",
"Lower bounds: [0.0, -Inf, -Inf, -Inf, 0.0, -Inf, -Inf, 0.0, -Inf, 0.0, 0.0, -Inf, -Inf, -Inf, 0.0, -Inf, -Inf, 0.0, -Inf, 0.0]\n",
"ftol_rel: 1.0e-12\n",
"ftol_abs: 1.0e-8\n",
"xtol_rel: 0.0\n",
"xtol_abs: [1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10, 1.0e-10]\n",
"initial_step: [0.75, 1.0, 1.0, 1.0, 0.75, 1.0, 1.0, 0.75, 1.0, 0.75, 0.75, 1.0, 1.0, 1.0, 0.75, 1.0, 1.0, 0.75, 1.0, 0.75]\n",
"maxfeval: -1\n",
"\n",
"Function evaluations: 781\n",
"Final parameter vector: [0.451204, 0.0494081, -0.0547485, 0.0358954, 0.0, -0.0657285, -0.0891811, 0.0204788, 0.0277562, 0.0, 0.541752, 0.0259224, -0.253773, 0.0189352, 0.0557164, 0.267834, 0.000872977, 0.0226436, 0.0620244, 0.0]\n",
"Final objective value: 28637.9710101938\n",
"Return code: FTOL_REACHED\n"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m3.optsum"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and the number of function evaluations to convergence has gone from under 40 to over 650.\n",
"\n",
"The time required for each function evaluation has also increased considerably,"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"const θ3 = m3.θ;"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 592.043 μs (559 allocations: 29.98 KiB)\n"
]
}
],
"source": [
"@btime objective(updateL!(setθ!($m3, $θ3)));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"resulting in much longer times for model fitting - about half a second for the \"best case\" on this machine."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that the estimates of the fixed-effects coefficients and their standard errors have not changed substantially except for the standard error of `P` (_Precedent_), which is also the largest effect."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The parameters in the optimization for this model can be arranged as two lower-triangular 4 by 4 matrices."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 LinearAlgebra.LowerTriangular{Float64,Array{Float64,2}}:\n",
" 0.451204 ⋅ ⋅ ⋅ \n",
" 0.0494081 0.0 ⋅ ⋅ \n",
" -0.0547485 -0.0657285 0.0204788 ⋅ \n",
" 0.0358954 -0.0891811 0.0277562 0.0"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m3.λ[1]"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 LinearAlgebra.LowerTriangular{Float64,Array{Float64,2}}:\n",
" 0.541752 ⋅ ⋅ ⋅ \n",
" 0.0259224 0.0557164 ⋅ ⋅ \n",
" -0.253773 0.267834 0.0226436 ⋅ \n",
" 0.0189352 0.000872977 0.0620244 0.0"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m3.λ[2]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which generate the covariance matrices for the random effects.\n",
"The cumulative proportion of the variance in the principal components of these covariance matrices, available as"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Array{Float64,1},1}:\n",
" [0.93984, 1.0, 1.0, 1.0] \n",
" [0.853217, 0.991263, 1.0, 1.0]"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m3.rePCA"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"The random effects for item also have 99% of the variation in the first two principal directions.\n",
"\n",
"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 `P` 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 `RTtrunc ~ 1 + spkr+prec+load + (1|subj) + (1+prec|item)`."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"const f4 = @formula(rt_trunc ~ 1 + spkr+prec+load + (1|subj) + (1+prec|item));"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 17.886 ms (75285 allocations: 3.31 MiB)\n"
]
},
{
"data": {
"text/plain": [
"Linear mixed model fit by maximum likelihood\n",
" rt_trunc ~ 1 + spkr + prec + load + (1 | subj) + (1 + prec | item)\n",
" logLik -2 logLik AIC BIC \n",
" -1.43319251×10⁴ 2.86638501×10⁴ 2.86818501×10⁴ 2.87312548×10⁴\n",
"\n",
"Variance components:\n",
" Variance Std.Dev. Corr.\n",
" item 133015.242 364.71255\n",
" 63766.936 252.52116 -0.70\n",
" subj 88819.437 298.02590\n",
" Residual 462443.388 680.03190\n",
" Number of obs: 1789; levels of grouping factors: 32, 56\n",
"\n",
" Fixed-effects parameters:\n",
" Estimate Std.Error z value P(>|z|)\n",
"(Intercept) 2181.85 77.4681 28.1645 <1e-99\n",
"spkr: old 67.879 16.0785 4.22172 <1e-4\n",
"prec: maintain -333.791 47.4472 -7.03499 <1e-11\n",
"load: yes 78.5904 16.0785 4.88792 <1e-5\n"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"const m4 = @btime fit!(LinearMixedModel($f4, $kb07, $contrasts))"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"93"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m4.optsum.feval"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"const θ4 = m4.θ;"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 136.383 μs (289 allocations: 14.30 KiB)\n"
]
}
],
"source": [
"@btime objective(updateL!(setθ!($m4, $θ4)));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These two model fits can be compared with one of the information criteria, `AIC` or `BIC`, for which \"smaller is better\".\n",
"They both indicate a preference for the smaller model, `m4`.\n",
"\n",
"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.\n",
"\n",
"Because model `m4` is a special case of model `m3`, a likelihood ratio test could also be used.\n",
"The alternative model, `m3`, will always produce an objective that is less than or equal to that from the null model, `m4`.\n",
"The difference in these value is similar to the change in the residual sum of squares in a linear model fit.\n",
"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.\n",
"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"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1-element Array{Float64,1}:\n",
" 25.879117466858588"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff(objective.([m3, m4]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This difference is compared to a $\\chi^2$ distribution with degrees of freedom corresponding to the difference in the number of parameters"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1-element Array{Int64,1}:\n",
" 20"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff(dof.([m4, m3]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"producing a p-value of"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.16984148520447773"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ccdf(Chisq(20), first(diff(objective.([m3,m4]))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The \"maximal\" model as proposed by [_Barr et al., 2013_](https://www.sciencedirect.com/science/article/pii/S0749596X12001180) would include all possible interactions of experimental and grouping factors."
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [],
"source": [
"const f5 = @formula(rt_trunc ~ 1 + spkr*prec*load + (1+spkr*prec*load|subj) + (1+spkr*prec*load|item));"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 10.297 s (2568234 allocations: 107.73 MiB)\n"
]
},
{
"data": {
"text/plain": [
"Linear mixed model fit by maximum likelihood\n",
" rt_trunc ~ 1 + spkr + prec + load + spkr & prec + spkr & load + prec & load + spkr & prec & load + (1 + spkr + prec + spkr & prec + load + spkr & load + prec & load + spkr & prec & load | subj) + (1 + spkr + prec + spkr & prec + load + spkr & load + prec & load + spkr & prec & load | item)\n",
" logLik -2 logLik AIC BIC \n",
" -1.42891729×10⁴ 2.85783458×10⁴ 2.87403458×10⁴ 2.91849882×10⁴\n",
"\n",
"Variance components:\n",
" Variance Std.Dev. Corr.\n",
" subj 90608.5644 301.012565\n",
" 5143.6952 71.719559 0.43\n",
" 5644.8543 75.132245 -0.48 -0.09\n",
" 9013.7657 94.940854 -0.20 -0.76 0.55\n",
" 7685.1682 87.665091 0.22 0.22 0.39 0.20\n",
" 1817.5847 42.633141 -0.45 -0.51 0.10 0.29 0.45\n",
" 7501.2782 86.609919 -0.10 -0.14 -0.04 0.04 -0.86 -0.72\n",
" 3885.1757 62.331178 0.48 0.41 0.37 0.16 0.09 -0.78 0.39\n",
" item 130494.7324 361.240546\n",
" 1749.7106 41.829542 0.36\n",
" 62347.1395 249.694092 -0.68 0.43\n",
" 1036.0759 32.188132 -0.57 -0.79 -0.03\n",
" 2871.3208 53.584706 0.20 0.06 -0.17 -0.03\n",
" 1472.8731 38.378029 -0.30 0.03 0.31 -0.20 -0.53\n",
" 4738.0323 68.833366 0.07 0.26 0.22 0.25 -0.15 -0.03\n",
" 4990.5418 70.643767 -0.05 -0.48 -0.31 0.65 0.69 -0.72 0.11\n",
" Residual 401618.2118 633.733550\n",
" Number of obs: 1789; levels of grouping factors: 56, 32\n",
"\n",
" Fixed-effects parameters:\n",
" Estimate Std.Error z value P(>|z|)\n",
"(Intercept) 2181.86 76.9448 28.3562 <1e-99\n",
"spkr: old 67.7443 19.263 3.5168 0.0004\n",
"prec: maintain -333.925 47.6832 -7.003 <1e-11\n",
"load: yes 78.581 21.2485 3.69819 0.0002\n",
"spkr: old & prec: maintain -21.7763 20.4419 -1.06528 0.2868\n",
"spkr: old & load: yes 18.3885 17.4074 1.05636 0.2908\n",
"prec: maintain & load: yes 4.53807 22.5067 0.201632 0.8402\n",
"spkr: old & prec: maintain & load: yes 23.4181 21.2101 1.1041 0.2695\n"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"const m5 = @btime fit!(LinearMixedModel($f5, $kb07, $contrasts))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Array{Float64,1},1}:\n",
" [0.724878, 0.846927, 0.94087, 0.999998, 1.0, 1.0, 1.0, 1.0] \n",
" [0.796978, 0.94188, 0.975067, 0.996119, 0.999996, 1.0, 1.0, 1.0]"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m5.rePCA"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For both the subjects and the items practically all the variation of these 8-dimensional random effects is in the first 4 principal components.\n",
"\n",
"The dimension of $\\theta$, the parameters in the optimization, increases considerably"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"72"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"const θ5 = m5.θ;\n",
"length(θ5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of these 72 parameters, 36 are estimated from variation between items, yet there are only 32 items.\n",
"\n",
"Because the dimension of the optimization problem has gotten much larger the number of function evaluations to convergence increases correspondingly."
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2810"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m5.optsum.feval"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Also, each function evaluation requires more time"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.284 ms (561 allocations: 30.02 KiB)\n"
]
}
],
"source": [
"@btime objective(updateL!(setθ!($m5, $θ5)));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"almost all of which is for the call to `updateL!`."
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.228 ms (519 allocations: 28.88 KiB)\n"
]
}
],
"source": [
"@btime updateL!($m5);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [],
"source": [
"const f6 = @formula(rt_trunc ~ 1 + spkr*prec*load + \n",
" (1+spkr+prec+load+spkr&prec+spkr&load+prec&load|subj) + \n",
" (1+spkr+prec+load+spkr&prec+spkr&load+prec&load|item));"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.950 s (1526668 allocations: 59.76 MiB)\n"
]
},
{
"data": {
"text/plain": [
"Linear mixed model fit by maximum likelihood\n",
" rt_trunc ~ 1 + spkr + prec + load + spkr & prec + spkr & load + prec & load + 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)\n",
" logLik -2 logLik AIC BIC \n",
" -1.43019526×10⁴ 2.86039053×10⁴ 2.87339053×10⁴ 2.9090717×10⁴\n",
"\n",
"Variance components:\n",
" Variance Std.Dev. Corr.\n",
" subj 92941.31347 304.862778\n",
" 5017.82659 70.836619 0.49\n",
" 5026.50804 70.897871 -0.53 -0.18\n",
" 7197.78240 84.839745 0.24 0.24 0.33\n",
" 8159.97984 90.332607 -0.19 -0.76 0.55 0.17\n",
" 1310.88763 36.206182 -0.58 -0.52 0.43 0.52 0.43\n",
" 7451.93456 86.324588 -0.17 -0.23 -0.09 -0.93 0.03 -0.60\n",
" item 130595.36541 361.379808\n",
" 229.13521 15.137213 1.00\n",
" 62451.29181 249.902565 -0.68 -0.68\n",
" 2103.67100 45.865793 0.21 0.21 -0.10\n",
" 708.89408 26.625065 -0.59 -0.59 -0.04 -0.35\n",
" 335.37064 18.313128 -0.78 -0.78 0.69 -0.41 0.54\n",
" 4296.57645 65.548276 0.11 0.11 0.27 -0.32 0.10 0.42\n",
" Residual 417488.17855 646.133251\n",
" Number of obs: 1789; levels of grouping factors: 56, 32\n",
"\n",
" Fixed-effects parameters:\n",
" Estimate Std.Error z value P(>|z|)\n",
"(Intercept) 2181.91 77.2928 28.2292 <1e-99\n",
"spkr: old 67.7697 18.1705 3.72966 0.0002\n",
"prec: maintain -333.9 47.6945 -7.00081 <1e-11\n",
"load: yes 78.5329 20.6802 3.79749 0.0001\n",
"spkr: old & prec: maintain -21.7283 20.0317 -1.08469 0.2781\n",
"spkr: old & load: yes 18.3631 16.349 1.12319 0.2614\n",
"prec: maintain & load: yes 4.51266 22.3772 0.201663 0.8402\n",
"spkr: old & prec: maintain & load: yes 23.37 15.2775 1.5297 0.1261\n"
]
},
"execution_count": 51,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"const m6 = @btime fit!(LinearMixedModel($f6, $kb07, $contrasts))"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"56"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"const θ6 = m6.θ;\n",
"length(θ6)"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.568 ms (561 allocations: 30.02 KiB)\n"
]
}
],
"source": [
"@btime objective(updateL!(setθ!($m6, $θ6)));"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.538 ms (519 allocations: 28.88 KiB)\n"
]
}
],
"source": [
"@btime updateL!($m6);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Summary of goodness of fit"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Apply the goodness of fit measures to `m1` to `m6` creating a data frame"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"