In [1]:
if !isdir(Pkg.dir("ReTerms"))
Pkg.clone("git@github.com:dmbates/ReTerms.jl")
end
#Pkg.update() # to ensure the latest versions
In the ReTerms package a linear mixed-effects model is represented by a vector of terms, trms, a vector of parameterized lower triangular matrices, Λ, a blocked symmetric matrix, A, and a blocked upper triangular Cholesky factor, R.
The k elements of Λ correspond to the grouping factors for the random effects. The trms vector is of length k+1. The first k elements are implicit representations blocks in $\bf Z$. The last element is the horizontal concatenation of the fixed-effects model matrix, $\bf X$ and the response vector $\bf y$.
The symmetric matrix A, of which only the upper triangle is stored, represents
$$
\begin{bmatrix}
\bf Z'Z & \bf Z'X & \bf Z'y\\
\bf X'Z & \bf X'X & \bf X'y\\
\bf y'Z & \bf y'X & \bf y'y
\end{bmatrix}
$$
but the blocking scheme corresponds to the trms vector. The same blocking scheme is used for the upper Cholesky factor R.
A couple of familiar examples may be helpful.
Dyestuff example
In [2]:
using DataFrames,ReTerms
include(Pkg.dir("ReTerms","test","data.jl"));
dump(ds) # the Dyestuff data
In [3]:
m1 = LMM([ReMat(ds[:Batch])],ds[:Yield]);
length(m1.trms)
Out[3]:
In [4]:
m1.trms[2]'
Out[4]:
In [5]:
m1[:θ] # parameter vector
Out[5]:
In [6]:
m1.Λ
Out[6]:
In [7]:
m1.A[2,2]
Out[7]:
For small examples and when running parallel simulations it is best to allow only one BLAS thread.
In [8]:
blas_set_num_threads(1)
fit(m1,true);
In [9]:
m1.Λ
Out[9]:
In [10]:
m1.R[1,1]
Out[10]:
In [11]:
m1.R[1,2]
Out[11]:
In [12]:
m1.R[2,2]
Out[12]:
In [13]:
pwrss(m1)
Out[13]:
In [14]:
√pwrss(m1)
Out[14]:
The steps in evaluating the objective function are to split the parameter vector over the Λ vector, copy A into R (upper triangle) and update its blocks, except for the last column, to $\bf\Lambda'Z'Z\Lambda+I$, then perform the blocked Cholesky factorization. As seen above, the lower right element of the lower right block of R is the square root of the penalized residual sum-of-squares.
During the iterations there is no need to solve for the coefficients and there are no operations that involve vectors of length n or matrices with a dimension n.
sleepstudy example
In [15]:
dump(slp)
In [16]:
const X = hcat(ones(size(slp,1)),slp[:Days]);
X'
Out[16]:
In [17]:
m2 = fit(LMM([VectorReMat(slp[:Subject],X')],X,slp[:Reaction]),true);
In [18]:
gc(); @time fit(LMM([VectorReMat(slp[:Subject],X')],X,slp[:Reaction]));
In [19]:
gc(); @time fit(LMM([VectorReMat(slp[:Subject],X')],X,slp[:Reaction]));
In [20]:
UpperTriangular(m2.R[end,end])
Out[20]:
In [21]:
m2.R[1,1] # a homogeneous block diagonal matrix - all blocks are the same size
Out[21]:
Only the upper triangle of the blocks on the diagonal of m2.R[1,1] is used. The diagonal blocks are stored as a 3-dimensional array, of size $2\times2\times18$ in this case. The [1,1] block of A and of R is always diagonal or homogeneous block diagonal with small blocks. Thus we order multiple random-effects terms by decreasing number of columns in the relevant block of $\bf Z$.
penicillin example
In [22]:
dump(pen)
In [23]:
m3 = fit(LMM([ReMat(pen[s]) for s in [:Plate,:Sample]],pen[:Diameter]),true);
In [24]:
m3.A[1,2]
Out[24]:
In [25]:
m3.A[2,2]
Out[25]:
In [26]:
m3.R[1,2]
Out[26]:
In [27]:
UpperTriangular(m3.R[2,2])
Out[27]:
In [28]:
UpperTriangular(m3.R[3,3])
Out[28]:
In this example the product
In [29]:
m3.A[1,2]'m3.A[1,2]
Out[29]:
is dense, hence m3.R[2,2] becomes dense even though m3.A[2,2] is diagonal.
In [30]:
gc(); @time fit(LMM([ReMat(pen[s]) for s in [:Plate,:Sample]],pen[:Diameter]));
In [31]:
gc(); @time fit(LMM([ReMat(pen[s]) for s in [:Plate,:Sample]],pen[:Diameter]));
In [32]:
dump(psts)
In [33]:
m4 = fit(LMM([ReMat(psts[s]) for s in [:Sample,:Batch]],psts[:Strength]),true);
In [34]:
m4.R[1,2]
Out[34]:
In [35]:
m4.R[1,2]'m4.R[1,2]
Out[35]:
Because this product is diagonal, m4.R[2,2] is also diagonal.
In [36]:
m4.R[2,2]
Out[36]:
In working on some simulations of data from a model and fitting that model and a couple of others to it, I realized that I was essentially doing a parametric bootstrap. I wrote a bootstrap that uses a callback to accumulate the desired result. Suppose that I want to save the fixed-effects, their standard errors and the deviance from 10,000 parametric bootstrap fits of model m2.
In [37]:
N = 10000;
const betamat = zeros(2,N);
const stderrmat = zeros(2,N);
const devvec = zeros(N);
function saveresults(i::Integer,m::LMM) # the callback function
copy!(sub(betamat,:,i),coef(m))
copy!(sub(stderrmat,:,i),stderr(m))
devvec[i] = objective(m)
end
Out[37]:
In [38]:
gc(); @time ReTerms.bootstrap(m2,N,saveresults);
In [39]:
using Gadfly
In [40]:
plot(x=sub(betamat,1,:),Geom.density)
Out[40]:
There is still quite a bit to do but I think that 10,000 bootstrap fits in under 20 seconds is pretty good. The next stage is parallelization.