# Singular covariance estimates in random regression models

This notebook explores the occurrence of singularity in the estimated covariance matrix of random regression models.
These are mixed-effects models with vector-valued random effects.

First, fit a model to the `sleepstudy` data from [`lme4`](https://github.com/lme4/lme4).

## Fitting a linear mixed-model to the sleepstudy data

Load the required packages

```{julia;term=true}
using DataFrames, FreqTables, LinearAlgebra, MixedModels, Random, RData
using Gadfly
using Gadfly.Geom: density, histogram, point
using Gadfly.Guide: xlabel, ylabel
const dat = Dict(Symbol(k)=>v for (k,v) in 
    load(joinpath(dirname(pathof(MixedModels)), "..", "test", "dat.rda")));
```

However, the `LinearMixedModel` constructor only creates a model structure but does not fit it.
An explicit call to `fit!` is required to fit the model.
As is customary (though not required) in Julia, a function whose name ends in `!` is a _mutating function_ that modifies one or more of its arguments.

An optional second argument of `true` in the call to `fit!` produces verbose output from the optimization.

```{julia;term=true}
sleepm = fit!(LinearMixedModel(@formula(Y ~ 1 + U + (1+U|G)), dat[:sleepstudy]), true)
```

The variables in the optimization are the elements of a lower triangular matrix, $\Lambda$, which is the relative covariance factor of the random effects.
The corresponding parameter vector is called $\theta$.

```{julia;term=true}
Λ = sleepm.λ[1]
```

The matrix $\Lambda$ is the left (or lower) Cholesky factor of the covariance matrix of the unconditional distribution of the vector-valued random effects, relative to the variance, $\sigma^2$, of the per-observation noise.
That is
\begin{equation}
    \Sigma = \sigma^2\Lambda\Lambda'
\end{equation}

In terms of the estimates,

```{julia;term=true}
s² = varest(sleepm)    # estimate of the residual variance
```

```{julia;term=true}
s² * Λ * Λ'   # unconditional covariance matrix of the random effects
```

The estimated correlation of the random effects can, of course, be evaluated from the covariance matrix.
Writing out the expressions for the elements of the covariance matrix in terms of the elements of Λ shows that many terms cancel in the evaluation of the correlation, resulting in the simpler formula.

```{julia;term=true}
Λ[2, 1] / sqrt(Λ[2, 1]^2 + Λ[2, 2]^2)
```

For a $2\times 2$ covariance matrix it is not terribly important to perform this calculation in an efficient and numerically stable way.
However, it is a good idea to pay attention to stability and efficiency in a calculation that can be repeated tens of thousands of times in a simulation or a parametric bootstrap.
The `norm` function evaluates the (geometric) length of a vector in a way that controls round-off better than the naive calculation.
The `view` function provides access to a subarray, such as the second row of $\Lambda$, without generating a copy.
Thus the estimated correlation can be written

```{julia;term=true}
Λ[2, 1] / norm(view(Λ, 2, :))
```

### Optimization with respect to θ

As described in section 3 of the 2015 _Journal of Statistical Software_ [paper](https://www.jstatsoft.org/index.php/jss/article/view/v067i01/v67i01.pdf) by Bates, Maechler, Bolker and Walker, using the relative covariance factor, $\Lambda$, in the formulation of mixed-effects models in the [`lme4`](https://github.com/lme4/lme4) and [`MixedModels`](https://github.com/dmbates/MixedModels) packages and using the vector $\theta$ as the optimization variable was a conscious choice.
Indeed, a great deal of effort went into creating this form so that the profiled log-likelihood can be easily evaluated and so that the constraints on the parameters, $\theta$, are simple "box" constraints.
In fact, the constraints are simple lower bounds.

```{julia;term=true}
show(sleepm.lowerbd)
```

In contrast, trying to optimize the log-likelihood with respect to standard deviations and correlations of the random effects
would be quite difficult because the constraints on the correlations when the covariance matrix is larger than $2\times 2$ are quite complicated.
Also, the correlation itself can be unstable.
Consider what happens to the expression for the correlation if both $\Lambda_{2,1}$ and $\Lambda_{2,2}$ are small in magnitude.
Small perturbations in $\Lambda_{2,1}$ that result in sign changes can move the correlation from near $-1$ to near $+1$ or vice-versa.

Some details on the optimization process are available in an `OptSummary` object stored as the `optsum` field of the model.

```{julia;term=true}
sleepm.optsum
```

### Convergence on the boundary

Determining if an estimated covariance matrix is singular is easy when using the $\theta$  parameters because singularity corresponds to points on the boundary of the allowable parameter space.
In other words, if the optimization converges to a vector in which either or both of $\theta_1$ or $\theta_3$ are zero, the covariance matrix is singular.
Otherwise it is non-singular.

The $\theta_1$ parameter is the estimated relative standard deviation of the random intercepts.
If this is zero then the correlation is undefined and reported as `NaN`.
If $\theta_3$ is zero and $\theta_2$ is non-zero then the estimated correlation is $\pm 1$ with the sign determined by the sign of $\theta_2$.
If both $\theta_2$ and $\theta_3$ are zero the correlation is `NaN` because the standard deviation of the random slopes will be zero.

Singular covariance matrices larger than $2\times 2$ do not necessarily result in particular values, like ±1, for the correlations.

Users of `lmer` or `lmm` are sometimes taken aback by convergence on the boundary if this produces correlations of `NaN` or $\pm 1$.
Some feel that this is a sign of model failure.
Others consider such estimates as a sign that Bayesian methods with priors that pull singular covariance matrices away from the boundary should be used.

This type of value judgement seems peculiar.
An important property of maximum likelihood estimates is that these estimates are well-defined once the probability model for the data has been specified.
It may be difficult to determine numerical values of the estimates but the definition itself is straightforward.
If there is a direct method of evaluating the log-likelihood at a particular value of the parameters, then, by definition, the mle's are the parameter values that maximize this log-likelihood.
Bates et al. (2015) provide such a method of evaluating the log-likelihood for a linear mixed-effects model.
Indeed they go further and describe how the fixed-effects parameters and one of the variance components can be profiled out of the log-likelihood evaluation, thereby reducing the dimension of the nonlinear, constrained optimization problem to be solved.

If the mle's correspond to a singular covariance matrix, this is a property of the model and the data.
It is not a mistake in some way.
It is just the way things are.
It reflects the fact that often the distribution of the estimator of a covariance matrix is diffuse.
It is difficult to estimate variances and covariances precisely.
A search for papers or books on "covariance estimation" will produce many results, often describing ways of regularizing the estimating process because the data themselves do not provide precise estimates.

For the example at hand a parametric bootstrap is one way of evaluating the precision of the estimates.

## The bootstrap function

The `MixedModels` package provides a `bootstrap` method to create a parametric bootstrap sample from a fitted model.

For reproducibility, set the random number seed to some arbitrary value.

```{julia;term=true}
Random.seed!(1234321);
```

Arguments to the `bootstrap` function are the number of samples to generate and the model from which to generate them.
By default the converged parameter estimates are those used to generate the samples.
Addition, named arguments can be used to override these parameter values, allowing `bootstrap` to be used for simulation.

`bootstrap` returns a `DataFrame` with columns:
- `obj`: the objective (-2 loglikelihood)
- `σ`: the standard deviation of the per-observation noise
- `β₁` to `βₚ`: the fixed-effects coefficients
- `θ₁` to `θₖ`: the covariance parameter elements
- `σ₁` to `σₛ`: the estimates standard deviations of random effects.
- `ρ₁` to `ρₜ`: the estimated correlations of random effects

The `ρᵢ` and `σᵢ` values are derived from the `θᵢ` and `σ` values.

```{julia;term=true}
sleepmbstrp = bootstrap(10000, sleepm);
show(names(sleepmbstrp))
```

Recall that the constrained parameters are $\theta_1$ and $\theta_3$ which both must be non-negative.
If either or both of these are zero (in practice the property to check is if they are "very small", which here is arbitrarily defined as less than 0.00001) then the covariance matrix is singular.

```{julia;term=true}
issmall(x) = x < 0.00001   # defines a one-liner function in Julia
```

```{julia;term=true}
freqtable(issmall.(sleepmbstrp[:θ₁]), issmall.(sleepmbstrp[:θ₃]))
```

Here the covariance matrix estimate is non-singular in 9,686 of the 10,000 samples, has an zero estimated intercept variance in 6 samples and is otherwise singular (i.e. correlation estimate of $\pm 1$) in 308 samples.

Empirical densities of the θ components are:

```{julia;echo=false;fig_width=8}
plot(sleepmbstrp, x = :θ₁, density, 
    xlabel("Relative standard deviation of intercept, θ₁"))
```

```{julia;echo=false;fig_width=8}
plot(sleepmbstrp, x = :θ₃, density, xlabel("θ₃"))
```

A density plot is typically a good way to visualize such a large sample.
However, when there is a spike such as the spike at zero here, a histogram provides a more informative plot.

```{julia;echo=false;fig_width=8}
plot(sleepmbstrp, x = :θ₃, histogram, xlabel("θ₃"))
```

```{julia;echo=false;fig_width=8}
plot(x = filter(isfinite, sleepmbstrp[:ρ₁]), histogram, xlabel("ρ₁₂"))
```

### Reciprocal condition number

The definitve way to assess singularity of the estimated covariance matrix is by its _condition number_ or, alternatively, its _reciprocal condition number_.
In general the condition number, $\kappa$, of a matrix is the ratio of the largest singular value to the smallest.
For singular matrices it is $\infty$, which is why it is often more convenient to evaluate and plot $\kappa^{-1}$.
Because $\kappa$ is a ratio of singular values it is unaffected by nonzero scale factors.
Thus
\begin{equation}
\kappa^{-1}(s^2\Lambda\Lambda') = \kappa^{-1}(\Lambda\Lambda') =
[\kappa^{-1}(\Lambda)]^2
\end{equation}

```{julia}
function recipcond(bstrp::DataFrame)
    T = eltype(bstrp[:θ₁])
    val = sizehint!(T[], size(bstrp, 1))
    d = Matrix{T}(undef, 2, 2)
    for (t1, t2, t3) in zip(bstrp[:θ₁], bstrp[:θ₂], bstrp[:θ₃])
        d[1, 1] = t1
        d[1, 2] = t2
        d[2, 2] = t3
        v = svdvals!(d)
        push!(val, v[2] / v[1])
    end
    val
end
rc = recipcond(sleepmbstrp)
extrema(rc)
```

```{julia;echo=false;fig_width=8}
plot(x = rc, histogram, xlabel("κ⁻¹"))
```

$\kappa^{-1}$ is small if either or both of $\theta_1$ or $\theta_3$ is small.

```{julia;term=true}
sum(issmall, rc)
```

The density of the estimated correlation

```{julia;echo=false;fig_width=8}
plot(x = filter(isfinite, sleepmbstrp[:ρ₁]), histogram, xlabel("ρ₁₂"))
```

```{julia;term=true}
sum(isfinite, sleepmbstrp[:ρ₁])  # recall that ρ = NaN in 7 cases
```

```{julia;term=true}
sum(x -> x == -1, sleepmbstrp[:ρ₁])  # number of cases of rho == -1
```

```{julia;term=true}
sum(x -> x == +1, sleepmbstrp[:ρ₁])  # number of cases of rho == +1
```

In this case the bootstrap simulations that resulted in $\rho = -1$ were not close to being indeterminant with respect to sign.
That is, the values of $\theta_2$ were definitely negative.

```{julia;term=true}
sleepmbstrp[:θ₂][findall(x -> x == -1, sleepmbstrp[:ρ₁])]
```

## The Oxboys data

In the `nlme` package for R there are several data sets to which random regression models are fit.
The `RCall` package for Julia provides the ability to run an embedded R process and communicate with it.
The simplest form of writing R code within Julia is to use character strings prepended with `R`.
In Julia strings are delimited by `"` or by `"""`.
With `"""` multi-line strings are allowed.

```{julia;term=true}
using RCall
R"""
library(nlme)
plot(Oxboys)
"""
```

```{julia;term=true}
oxboys = rcopy(R"Oxboys");
show(names(oxboys))
```

```{julia;term=true}
oxboysm = fit(LinearMixedModel, @formula(height ~ 1 + age + (1+age | Subject)), oxboys)
```

```{julia;term=true}
show(getθ(oxboysm))
```

As seen in the plot and by the estimate $\widehat{\theta_1} = 12.0$, the estimated standard deviation of the random intercepts is much greater than the residual standard deviation.
It is unlikely that bootstrap samples will include singular covariance estimates.

```{julia;term=true}
Random.seed!(4321234);
oxboysmbtstrp = bootstrap(10000, oxboysm);
```

In this bootstrap sample, there are no singular estimated covariance matrices for the random effects.

```{julia;term=true}
freqtable(issmall.(oxboysmbtstrp[:θ₁]), issmall.(oxboysmbtstrp[:θ₃]))
```

The empirical density of the correlation estimates shows that even in this case the correlation is not precisely estimated.

```{julia;echo=false;fig_width=8}
plot(oxboysmbtstrp, x = :ρ₁, histogram, xlabel("ρ₁₂"))
```

```{julia;term=true}
extrema(oxboysmbtstrp[:ρ₁])
```

The reciprocal condition number

```{julia;term=true}
rc = recipcond(oxboysmbtstrp);
extrema(rc)
```

does not get very close to zero.

```{julia;echo=false;fig_width=8}
plot(x = rc, histogram, xlabel("κ⁻¹"))
```

## The Orthodont data

```{julia;term=true}
R"plot(Orthodont)"
```

The subject labels distinguish between the male and the female subjects.  Consider first the female subjects only.

```{julia;term=true}
orthfemale = rcopy(R"subset(Orthodont, Sex == 'Female', -Sex)");
orthfm = fit!(LinearMixedModel(@formula(distance ~ 1 + age + (1 + age | Subject)), orthfemale))
```

```{julia;term=true}
srand(1234123)
orthfmbtstrp = bootstrap(10000, orthfm);
```

```{julia;term=true}
freqtable(issmall.(orthfmbtstrp[:θ₁]), issmall.(orthfmbtstrp[:θ₃]))
```

For this model almost 1/3 of the bootstrap samples converge to singular covariance estimates for the vector-valued random effects.
A histogram of the estimated correlations of the random effects is dominated by the boundary values.

```{julia;echo=false;fig_width=8}
plot(orthfmbtstrp, x = :ρ₁, histogram, xlabel("ρ₁₂"))
```

Even though the estimated correlation in the model is -0.30, more of the boundary values are at +1 than at -1.
This may be an artifact of the optimization routine used.
In some cases there may be multiple optima on the boundary.
It is difficult to determine the global optimum in these cases.

A histogram of the reciprocal condition number is also dominated by the boundary values.

```{julia;echo=false}
#plot(x = recipcond(orthfmbtstrp), density, xlabel("κ⁻¹"))
```

## Early childhood cognitive study

This example from Singer and Willett (2003), *Applied Longitudinal Data Analysis* was the motivation for reformulating the estimation methods to allow for singular covariance matrices. Cognitive scores (`cog`) were recorded at `age` 1, 1.5 and 2 years on 103 infants, of whom 58 were in the treatment group and 45 in the control group.  The treatment began at age 6 months (0.5 years).  The data are available as the `Early` data set in the `mlmRev` package for R.  In the model, time on study (`tos`) is used instead of age because the zero point on the time scale should be when the treatment begins.

```{julia;term=true}
R"""
suppressMessages(library(mlmRev))
library(lattice)
Early$tos <- Early$age - 0.5
Early$trttos <- Early$tos * (Early$trt == "Y")
xyplot(cog ~ tos | reorder(id, cog, min), Early, 
    type = c("p","l","g"), aspect="xy")
"""
```

Notice that most of these plots within subjects have a negative slope and that the scores at 1 year of age (`tos = 0.5`) are frequently greater than would be expected on an age-adjusted scale.

```{julia;term=true}
R"print(xtabs(cog ~ age + trt, Early) / xtabs(~ age + trt, Early))";
```

When the time origin is the beginning of the treatment there is not generally a "main effect" for the treatment but only an interaction of `trt` and `tos`.

```{julia;term=true}
early = rcopy(R"subset(Early, select = c(cog, tos, id, trt, trttos))");
earlym = fit(LinearMixedModel, @formula(cog ~ 1 + tos + trttos + (1 + tos | id)), early)
```

The model converges to a singular covariance matrix for the random effects.

```{julia;term=true}
getθ(earlym)
```

The conditional (on the observed responses) means of the random effects fall along a line.

```{julia;echo=false;fig_width=5}
rr = ranef(earlym)[1];
plot(x = view(rr, 1, :), y = view(rr, 2, :), point,
    xlabel("Conditional mean of intercept random effect"),
    ylabel("Conditional mean of slope random effect"))
```
