A MixedModels.LinearMixedModel can be constructed without being fit.
For the purposes of examples, load the dataset stored with the MixedModels package.
In [1]:
using DataFrames, MixedModels, RData
const dat = convert(Dict{Symbol,DataFrame},
load(Pkg.dir("MixedModels", "test", "dat.rda")))
Out[1]:
A model for the sleepstudy data with random-effects for the intercept and for the days of sleep deprivation, U, by subject, G is declared as
In [2]:
m = lmm(@formula(Y ~ 1 + U + (1+U|G)), dat[:sleepstudy]);
The fields in m
In [3]:
fieldnames(m)
Out[3]:
include trms, a vector of AbstractTerms, and two blocked arrays, A and L. The last two terms are always MatrixTerms representing the fixed-effects model matrix, $\bf X$, and the response vector, $\bf y$. The terms preceding these two are always random-effects terms.
In [4]:
typeof.(m.trms)
Out[4]:
The A field is a blocked matrix
In [5]:
typeof(m.A)
Out[5]:
whose rows and columns of blocks correspond to the terms.
In [6]:
nblocks(m.A)
Out[6]:
The second last diagonal block in A is $\bf X'X$.
In [7]:
m.A[Block(2,2)]
Out[7]:
This can be verified by noting
In [8]:
dat[:sleepstudy]
Out[8]:
In [9]:
X = m.trms[end - 1].x # the model matrix X
Out[9]:
In [10]:
X'X
Out[10]:
A pivoted Cholesky decomposition of this matrix returns the computational rank, which is also the rank of $\bf X$.
In [11]:
Lpiv = cholfact(Symmetric(m.A[Block(2,2)], :L), Val{true})
Out[11]:
In [12]:
rank(Lpiv)
Out[12]:
An optional argument, tol, in the call to cholfact specifies the tolerance for determining the rank. It defaults to zero which would correspond to exact rank deficiency. A more reasonable value may be $\sqrt{\epsilon}$ where $\epsilon$ is the relative machine precision.
In [13]:
√eps()
Out[13]:
As explained in the documentation for dpstrf, this value should be negated to be used as a relative tolerance. Otherwise it is used as an absolute tolerance which usually is not what you want.
In [14]:
Lpiv = cholfact(Symmetric(m.A[Block(2,2)], :L), Val{true}, tol = -√eps())
Out[14]: