
# Representation of mixed-effects models

The representation of a linear mixed-effects model
```jl
type LinearMixedModel{T <: AbstractFloat} <: MixedModel
    formula::Formula
    mf::ModelFrame
    wttrms::Vector
    trms::Vector
    sqrtwts::Diagonal{T}
    Λ::Vector
    A::Hermitian
    L::LowerTriangular
    optsum::OptSummary{T}
end
```
underlies all the other types of mixed-effects models. The members of this struct are
* `formula`: the formula for the model
* `mf`: the model frame, the `terms` component for labelling fixed effects
* `wttrms`: a length `nt` vector of weighted model matrices. The last two elements are `X` and `y`.
* `trms`: a vector of unweighted model matrices.  If `isempty(sqrtwts)` the same object as `wttrms`
* `Λ`: a length `nt - 2` vector of lower triangular matrices
* `sqrtwts`: the `Diagonal` matrix of the square roots of the case weights.  Allowed to be size 0
* `A`: a Hermitian blocked matrix representing `hcat(Z,X,y)'hcat(Z,X,y)`
* `L`: a LowerTriangular blocked matrix - the Cholesky factor of `Λ'AΛ+I`
* `opt`: an [`OptSummary`](@ref) object

If there are no case weights then the size of `sqrtwts` is $0\times 0$ and `wttrms` is the same as `trms`.  To describe the other components, it is helpful to consider a few simple examples.

## A model with a single, scalar random-effects term

```{julia;term=true}
using Feather, LinearAlgebra, MixedModels
```

```{julia;term=true}
dsfilename = Pkg.dir("MixedModels", "test", "data", "Dyestuff.feather")
dyestuff = Feather.read(dsfilename, nullable = false);
```

```{julia;term=true}
m1 = lmm(Yield ~ 1 + (1 | Batch), dyestuff);
```

In the model formula there are three terms: the _response term_, `Yield`, the _Intercept_ term, `1`, which is part of the fixed-effects, and the scalar random-effects term, `(1 | Batch)`.  Random-effects terms consist of an expression on the left-hand side of the `|`, which is evaluated as a model matrix, and an expression on the right hand side: `Batch`, in this case, which is evaluated to a categorical array called the _grouping factor_ for the term.  When the left-hand side evaluates to a model matrix with a single column, as it does here, the term is said to be a _scalar random-effects term_.

The `trms` member contains these three terms in the order `random-effects`, `fixed-effects`, `response`.

```{julia;term=true}
m1.trms
```

The `ScalarReMat` is a compact representation of the _indicator matrix_ for the grouping factor.  Its columns are indicators of the 6 levels of the grouping factor, `Batch`.

```{julia;term=true}
full(m1.trms[1])
```

The `A` member is a blocked matrix with $3$ blocks in each dimension, corresponding to the 3 terms.  The $(i, j)$ block is `m1.trms[i]'m1.trms[j]`  Because of symmetry only the blocks in the lower triangle are stored explicitly.  (Generally the term _Hermitian_ is used in Julia rather than _symmetric_ even for real-valued matrices, where they are synonyms.)

```{julia;term=true}
m1.A[1, 1]
```

```{julia;term=true}
m1.A[2, 1]
```

```{julia;term=true}
m1.A[2, 2]
```

```{julia;term=true}
m1.A[3, 1]
```

```{julia;term=true}
m1.A[3, 2]
```

```{julia;term=true}
m1.A[3, 3]
```

Blocks on the diagonal of `A` will be positive-semi-definite Hermitian matrices.  More importantly, diagonal blocks corresponding to scalar random-effects terms are diagonal matrices.

The `Λ` member is a vector of lower-triangular matrices whose values are modified during the iterations.  They correspond to the random-effects terms.  In the case of a scalar random-effects term, the corresponding element of `Λ` is a multiple of the identity.  The multiple is the parameter $\theta$ over which the log-likelihood is optimized.

```{julia;term=true}
m1.Λ
```

```{julia;term=true}
setθ!(m1, [0.75258]);
m1.Λ
```

After setting a value of $\theta$ the blocked lower Cholesky is updated.
```jl
function cholBlocked!{T}(m::LinearMixedModel{T})
    A, Λ, L = m.A.data.blocks, m.Λ, m.L.data.blocks
    n = LinAlg.checksquare(A)
    for j in 1:n, i in j:n
        inject!(L[i, j], A[i, j])  # like copyto! but L can be more general than A
    end
    for (j, λ) in enumerate(Λ)
        for i in j:n
            rmul!(L[i, j], λ)
        end
        for jj in 1:j
            lmul!(adjoint(λ), L[j, jj])
        end
        L[j, j] += I
    end
    for j in 1:n
        Ljj = L[j, j]
        cholUnblocked!(Ljj, Val{:L})
        Ljjlt = isa(Ljj, Diagonal) ? Ljj : LowerTriangular(Ljj)
        for i in (j + 1):n
            LinAlg.A_rdiv_Bc!(L[i, j], Ljjlt)
        end
        for i in (j + 1):n
            Lij = L[i, j]
            Lii = L[i, i]
            rankUpdate!(-one(T), Lij, isa(Lii, Diagonal) ? Lii : Hermitian(Lii, :L))
            for jj in (i + 1):n
                A_mul_Bc!(-one(T), L[jj, j], Lij, one(T), L[jj, i])
            end
        end
    end
    m
end
```
