Fits a linear mixed model by REML or maximum likelihood with the
profiled-deviance method of lme4: the fixed effects and the residual
variance are profiled out, and the deviance is minimised over the
relative covariance parameters theta only. Formulas use
lme4's syntax – (1 | g), (x | g), uncorrelated
(x || g), nested (1 | a/b), crossed
(1 | a) + (1 | b) – or nlme's, a fixed formula with
random = ~ x | g.
Usage
morie_lmm(
formula,
data,
REML = TRUE,
random = NULL,
weights = NULL,
offset = NULL,
start = NULL,
control = list(),
...
)Arguments
- formula
A two-sided formula: fixed effects plus random terms
(terms | group)or(terms || group). With nlme's interface this is the fixed-effects formula.- data
A data frame.
- REML
Logical;
TRUE(default) for REML,FALSEfor maximum likelihood.- random
Optional nlme-style random effects: a one-sided formula
~ x | g(also~ 1 | a/b) or a named list of one-sided formulaslist(a = ~ 1, b = ~ 1)(b nested in a).- weights
Optional prior weights (residual variance \(\sigma^2 / w_i\)).
- offset
Optional offset added to the linear predictor.
- start
Optional starting value for
theta.- control
Optional list passed to
stats::nlminb'scontrol.- ...
Ignored (for compatibility with calls written for
lmer).
Value
An object of class morie_lmm: coefficients,
vcov, fixef_table (estimate, SE, t, and nlme's
containment df and p-value when the grouping factors are nested),
sigma, theta, varcorr (per term: SDs,
correlations, covariance), varcorr_table (as
as.data.frame(VarCorr())), ranef (conditional modes
per grouping factor), fitted, residuals,
logLik, df, AIC, BIC,
deviance (ML deviance at the estimates, or NA for
REML), REMLcrit, ngroups, nobs,
convergence (nlminb code and message, iterations,
evaluations, largest free-parameter gradient, singular flag) and
devfun, the profiled deviance as a function of theta (this
object's theta order), for checking against other fitters.
Details
The computation never forms an \(n \times n\) (or \(n \times q\)) matrix. The grouping factor with the most random effects is handled level by level in vectorised form from per-level crossproducts accumulated in one pass over the data, so a single grouping factor costs \(O(n q^2 + m q^3)\) per deviance evaluation. Other grouping factors enter through a Schur complement and a dense Cholesky factor of their total number of random effects \(q_2\) (the update is accumulated sparsely for crossed designs); that is fast up to a few thousand (about 0.5 s per evaluation at \(q_2 = 1000\)), warns above 5000, and needs the dense coupling block (\(m_1 q_1 \times q_2\)) to fit in memory.
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.
Examples
set.seed(1)
d <- data.frame(g = factor(rep(1:20, each = 8)), x = rep(0:7, 20))
d$y <- 1 + 0.5 * d$x + rnorm(20)[d$g] + rnorm(20, 0, 0.2)[d$g] * d$x +
rnorm(160)
fit <- morie_lmm(y ~ x + (x | g), d)
fit
#> Linear mixed model fit by REML
#> Formula: y ~ x + (x | g)
#> REML criterion
#> 496.0721
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> g (Intercept) 0.959890 0.97974
#> g x 0.014059 0.11857
#> Residual 0.890450 0.94363
#> Correlations:
#> Groups var1 var2 Corr
#> g (Intercept) x -0.028
#> Number of obs: 160, groups: g, 20
#>
#> Fixed effects:
#> Estimate Std.Error t value df Pr(>|t|)
#> (Intercept) 1.17940 0.257960 4.5719 139 1.0582e-05
#> x 0.50801 0.041988 12.0990 139 1.8162e-23
fit2 <- morie_lmm(y ~ x, d, random = ~ x | g)