For more years than I care to admit, I have been developing algorithms and implementations to fit what statisticians call mixed-effects models, most recently in Julia. These models are most useful when applied to very large data sets, which makes creating an effective, flexible implementation difficult. Evaluation of the log-likelihood to be maximized by the parameter estimates requires updating large structures, preferably in place, and some intricate linear algebra. Efficiency is best achieved by taking into account the special structure of the matrices representing the model. The Julia implementation can now out-perform earlier implementations that I and others have created, in some cases by orders of magnitude. The Julia type system, multiple dispatch, and, most of all, the fact that one can achieve both flexibility and performance in a single language,have been fundamental to achieving this performance.
Mixed-effects models are similar to linear models (LMs) or generalized linear models (GLMs) in that fitting the model requires estimation of coefficients in a linear predictor expression. In a linear mixed model (LMM) or a generalized linear mixed model (GLMM) there are two types of coefficients, or effects; fixed-effects, which are associated with factors that have reproducible levels, such as male/female, and random-effects, which are associated with factors whose levels are a sample from a population, such as subject or item.
Random-effects can be viewed as nuisance parameters in that they account for sources of variability in the observations but are often not by themselves of interest. Frequently the number of random effects grows with the number of observations (you observe a few observations on each of many subjects). Because their estimates will be ill-defined they are regularized. The model-fitting involves selecting the regularization parameters to maximize a log-likelihood.
In an LMM the $n$-dimensional response, $\bf y$, is considered to be the realization of random variable $\mathcal{Y}$ with conditional distribution $$ {\mathcal Y}|({\mathcal B} = \bf b) \sim {\mathcal N}({\bf X\beta + Z b},\sigma^2\bf I_n) . $$ The linear predictor, $\bf X\beta+Z b$, is based on two, known model matrices; the dense matrix $\bf X$ of size $n\times p$, and sparse matrix $\bf Z$ of size $n\times q$.
For example, we define a model for the InstEval data set from the lme4 package
In [1]:
using HDF5,ReTerms,StatsBase
inst = h5open("/var/tmp/dat.h5","r") do io g2dict(io,"inst") end;
dump(inst)
where Y, the response, is the evaluation by student, S, of instructor, D, in a course from department, P. The service factor, R, indicates if the course was a service course or not. We create random effects terms for student, instructor and department.
In [2]:
revec = [reterm(inst[s]) for s in [:S,:D,:P]];
for i in eachindex(revec) println(size(revec[i])) end
The three sparse matrices derived from these terms are, in this case, the indicator columns for the 2972 students, for the 1128 instructors, and for the 14 departments. The matrix $\bf Z$, of size $73421\times 4114$ is the horizontal concatenation of these three. It is (comparatively) large but very sparse -- each row has only three nonzeros.
The fixed effects in the model will be the intercept term and an indicator of a service course. The intercept corresponds to a typical evaluation in a non-service course, after adjusting for the student, instructor and department. The coefficient for R is the change in typical evaluation for service compared to non-service courses.
In [3]:
X = hcat(ones(length(inst[:R])),convert(Vector{Float64},inst[:R] .== 2));
size(X)
Out[3]:
The unconditional distribution of the random effects is also multivariate normal $$ {\mathcal B}\sim\mathcal{N}(\bf 0,\Sigma_{\theta}) $$ The positive-definite covariance matrix, $\Sigma$, depends on a small parameter vector $\theta$. For scalar random -effects such as these (a scalar effect for each student, each instuctor and each department), the elements of $\mathcal B$ are independent within and between groups. That is, $\Sigma$ is diagonal. Furthermore the variances of the effects are constant within group. For technical reasons we write $$ \Sigma_{\theta} = \sigma^2\Lambda_{\theta}\Lambda_{\theta}' $$ where $\sigma^2$ is the same scale parameter as in the conditional distribution of $\mathcal Y$. $\Lambda_{\theta}$ is the lower Cholesky factor of the relative covariance matrix, $\Sigma/\sigma^2$.
Since the 1970's it has been recognized that, given a value of the covariance parameters, $\theta$, the conditional estimate of the fixed-effects, $\widehat{\bf\beta}$, and the "BLUPS" of the random effects, $\tilde{\bf b}$ are the solution to "Henderson's mixed-model equations" $$ \begin{bmatrix} \bf X'X & \bf X'Z \\ \bf Z'X & \bf Z'Z+\Sigma_{\theta}^{-1} \end{bmatrix} \begin{bmatrix} \widehat{\beta}\ \tilde{\bf b}
= \begin{bmatrix} \bf \Lambda_\theta'Z'y\\ \bf X'y \end{bmatrix} $$ corresponding to the _pseudo data_ representation of the penalized least squares problem $$ \rho^2(\theta)= \min{\beta,\bf u} \left| \begin{bmatrix} \bf y\\ \bf 0 \end{bmatrix} - \begin{bmatrix} \bf Z\Lambda\theta & \bf X \ \bf I_q & \bf 0 \end{bmatrix} \begin{bmatrix} \bf u\\ \beta \end{bmatrix} \right|^2 . $$
In this formulation it is natural to form the Cholesky factor of the system matrix for which the upper left block is $$ \bf L_\theta L_\theta' = \bf \Lambda_{\theta}'Z'Z\Lambda_{\theta} + I $$ providing straightforward evaluation of $$ \log\left(\left|\bf \Lambda_{\theta}'Z'Z\Lambda_{\theta} + I\right|\right) = 2\log\left(\left|\bf L_\theta\right|\right). $$
The profiled log-likelihood can be written as a function of $\theta$ only $$ -2\ell(\theta|{\bf y}) = -2\log(|\bf L_\theta|)+n\left(1+\frac{2\pi\rho^2(\theta)}n\right) . $$
In [4]:
m1 = LMM(X,revec,inst[:Y]);
fieldnames(m1)
Out[4]:
In [5]:
size(m1.A)
Out[5]:
In [6]:
A = m1.A;
for j in 1:4, i in j:4 println(i,", ",j," ",size(A[i,j])," ",typeof(A[i,j])) end
In [7]:
L = m1.L
for j in 1:4, i in j:4 println(i,", ",j," ",size(L[i,j])," ",typeof(L[i,j])) end
In [8]:
fit(m1,true);
In [9]:
gc(); @time fit(LMM(X,revec,inst[:Y]));