Linear Mixed Models (REML/ML) or MixedModelFit

This package fits linear mixed models in Igor Pro, by restricted maximum

likelihood (REML) or by maximum likelihood. It is pure Igor procedure code: no

XOP and no external dependencies, nothing beyond `MatrixOP` and the `Optimize`

operation.

What is a linear mixed model?

Ordinary regression assumes that your observations are independent of one

another. A great deal of real experimental data is not. If you measure the same

subject ten times, those ten numbers are more alike than ten numbers taken from

ten different subjects. If you run samples in batches, on plates, in sessions,

at different sites, or over several days, observations sharing a batch tend to

share whatever that batch did differently. The grouping is part of how the data

were produced, and ignoring it does not make it go away.

A linear mixed model is ordinary regression plus an explicit description of

that grouping structure. It splits the effects in your model into two kinds.

Fixed effects are the quantities you set out to measure, and which you want

a single number for: the slope of response against dose, the difference between

treated and control, the overall mean. They are assumed to be the same

underlying constants for everyone, and a mixed model reports an estimate and a

standard error for each, just as ordinary regression does.

Random effects are the departures of individual groups from those common

values. Each subject, batch or site gets its own offset, and possibly its own

slope. Crucially, you do not estimate those offsets as free parameters, one per

group, the way you would with a dummy variable for every subject. Instead you

assume they are drawn from a normal distribution with mean zero, and you

estimate the *variance* of that distribution. One number, the variance

component, stands in for the whole population of groups, including the ones you

did not sample.

That shift in viewpoint is what makes the model useful.

You get the standard errors right.  Treating correlated observations as

independent inflates your effective sample size. Ordinary least squares then

reports standard errors that are too small and p values that are too

optimistic, sometimes dramatically so. Measuring one subject a hundred times

does not give you a hundred subjects' worth of information about the

population, and a mixed model knows that.

Variation between groups becomes a result, not a nuisance.  The fit reports

the between group standard deviation alongside the within group standard

deviation. In a method validation study that is reproducibility versus

repeatability. In an animal study it is how much of the spread is biological

and how much is measurement. The question "where is my variability coming

from?" is answered directly, with an estimate and not a hand wave.

Unbalanced data are handled naturally.  Groups do not need equal numbers of

observations, and nothing needs to be discarded to make the design balanced.

Groups with more data simply contribute more.

Estimates for individual groups are shrunk sensibly.  The model also

predicts each group's own random effect, the BLUP (Best Linear Unbiased Predictor). These predictions are pulled

toward the overall mean by an amount that depends on how much data that group

has and how large the variance components are. A subject measured twice gets

pulled a long way toward the average; a subject measured fifty times barely

moves. This is borrowing strength across groups, and it gives better individual

predictions than fitting each group alone.

The simplest and most common case is a random intercept: every group sits

at its own level, scattered about the overall mean. The next step is a random

slope: the effect of a covariate, time or dose or concentration, differs from

group to group, and the model estimates how much it varies and whether groups

that start high also change fastest.

Formally the model is

    y = X*beta + Z*b + eps,   b ~ N(0, sigma^2 Lam Lam'),   eps ~ N(0, sigma^2 I)

where `X` holds the fixed effect predictors, `Z` is built from your grouping

variable, and the covariance of the random effects is what gets estimated. The

algorithm follows Bates et al. (2015), the paper behind the R package `lme4`:

`beta` and `sigma^2` are profiled out analytically, so the numerical

optimization runs only over the relative covariance factor. That is what makes

the fit fast and well behaved.

REML, the default, is the right choice for estimating variance components,

because it corrects for the degrees of freedom consumed by the fixed effects.

Maximum likelihood is also available, and is what you need if you want to

compare models that differ in their fixed effects.

What the package does

- Random intercepts, random slopes, or both, with a full covariance between

  them

- REML (default) or ML

- Fixed effect estimates with standard errors, z values and asymptotic p values

- Variance components, reported as variances, standard deviations and

  correlations

- BLUPs for every group

- Fitted values and residuals

- Results left in a data folder as ordinary waves and global variables, ready

  for further analysis

Getting started

Put MixedModelFit.ipf in your User Procedures folder and add

#include "MixedModelFit" to your own procedure file, or simply open the file

and drag it into Igor.

MixedModelGuide.md is the user guide: how to lay out your response, design

matrix and grouping wave, worked examples for random intercepts and random

slopes, what each output wave contains, a function reference and a

troubleshooting section. MixedModelDemo.pxp is a demo experiment holding two

classic data sets, the Dyestuff yields and the sleepstudy reaction times,

set up so that a single command line fits each of them.

Validation

I would rather be explicit about this than offer the usual boilerplate, because

the numbers a mixed model produces are easy to trust and hard to check by eye.

Verified against exact closed form results, with no external reference:

- the REML and ML criterion, by two independent implementations, agreeing to

  1e-13

- random intercept variance components on a balanced design against the ANOVA

  closed form, to 1e-6

- BLUPs against the closed form shrinkage estimator, to 1e-12

- the fixed effect standard error against its closed form

- optimizer convergence, checked against golden section search, at the

  numerical precision floor

- for random slopes, the covariance factor and the Z column layout against a

  third implementation that shares no code with them, to 1e-13

Verified against an external reference: the `Dyestuff` REML criterion

reproduces the published `lme4` value, 319.6543.

Not yet verified: a complete random slopes fit against `lme4` or SAS

`PROC MIXED` on real data. The machinery underneath random slopes is validated

and simulated data is recovered correctly, but the fitted values themselves

have not been compared against an independent implementation. Treat random

slope results as provisional until that check is done. If anyone has `lme4` or

SAS handy and wants to run the comparison, I would be glad of the help.

The validation suite ships with the package as `MixedModelTests.ipf`.

Limitations

Please read these before relying on the output.

- No denominator degrees of freedom. There is no Satterthwaite or

  Kenward-Roger approximation. The p values are asymptotic, from the normal

  distribution. With fewer than about 30 groups they are anticonservative, that

  is they make effects look more significant than they are. Treat them as a

  rough guide.

- No AIC or BIC. Deliberately omitted rather than shipped with an unchecked

  convention. The log likelihood is available if you want to compute them.

- No likelihood ratio tests or profile confidence intervals.

- One grouping factor only. No crossed random effects. Nested factors can

  be handled by combining the codes, for example `plot*1000 + subplot`.

- Gaussian responses only. No generalized linear mixed models.

- Independent residuals only.  No AR(1) or spatial structures.

- You build the dummy columns for categorical predictors yourself.

- Dense linear algebra. A few hundred groups is comfortable; several

  thousand will be slow and memory hungry.

Requirements

Igor Pro 9.00 or later.

What is in the zip:

- MixedModelFit.ipf -- the package

- MixedModelGuide.md -- the user guide

- MixedModelDemo.pxp -- the demo experiment

- MixedModelTests.ipf -- the validation suite, optional

- README.md -- short install and demo instructions

References

Bates, D., Maechler, M., Bolker, B., and Walker, S. (2015). "Fitting Linear

Mixed-Effects Models Using lme4." *Journal of Statistical Software* 67(1),

1-48.

Pinheiro, J.C., and Bates, D.M. (2000). Mixed-Effects Models in S and S-PLUS.

Springer.

---

Comments, bug reports and especially cross-checks against other mixed model

software are welcome.