Example: Contrasting cluster-based and case-based CV
We introduce four artificial data sets that exemplify aspects of cross-validation particular to hierarchical models. Using these data sets, we show that model comparisons employing cluster-based and those employing case-based cross-validation may not agree on a “best” model. Furthermore, commonly used measures of fit, such as mean-squared error, do not necessarily become smaller as models become larger, even when the models are nested, and even when the measure of fit is computed for the whole data set.
The four datasets differ in the relative magnitude of between-cluster variance compared with within-cluster variance. They serve to illustrate how fitting mixed models, and, consequently, the cross-validation of mixed models, is sensitive to relative variance, which determines the degree of shrinkage of within-cluster estimates of effects towards between-cluster estimates.
For these analyses, we will use the glmmTMB() function
in the glmmTMB package (Brooks
et al., 2017) because, in our experience, it is more likely to
converge than functions in the nlme and the
lme4 packages for models with low between-cluster
variance.
Consider a researcher studying the effect of the dosage of a drug on the severity of symptoms for a hypothetical disease. The researcher has longitudinal data on 20 patients, each of whom was observed on five occasions in which patients received different dosages of the drug. The data are observational, with dosages prescribed by the patients’ physicians, so that patients who were more severely affected by the disease received higher dosages of the drug.
Our four contrived data sets (see below) illustrate possible results
for data obtained in such a scenario. The relative configuration of
dosages x and symptoms y are identical within
patients in each of the four datasets. Within patients, higher dosages
are generally associated with a reduction in symptoms.
Between patients, however, higher dosages are associated with higher levels of symptoms. A plausible mechanism is a reversal of causality: Within patients, higher dosages alleviate symptoms, but between patients higher morbidity causes the prescription of higher dosages.
The four data sets differ in the between-patient variance of patient centroids from a common between-patients regression line. The data sets exhibit a progression from low to high variance around the common regression line.
We start by generating a data set to serve as a common template for the four sample data sets. We then apply different multipliers to the between-patient variation using parameters consistent with the description above to highlight the issues that arise in cross-validating mixed-effects models:4
# Parameters:
Nb <- 20 # number of patients
Nw <- 5 # number of occasions for each patient
Bb <- 1.0 # between-patient regression coefficient on patient means
Bw <- -0.5 # within-patient effect of x
SD_between <- c(0, 5, 6, 8) # SD between patients
SD_within <- rep(2.5, length(SD_between)) # SD within patients
Nv <- length(SD_within) # number of variance profiles
SD_ratio <- paste0('SD ratio = ', SD_between,' / ',SD_within)
SD_ratio <- factor(SD_ratio, levels = SD_ratio)
set.seed(833885)
Data_template <- expand.grid(patient = 1:Nb, obs = 1:Nw) |>
within({
xw <- seq(-2, 2, length.out = Nw)[obs]
x <- patient + xw
xm <- ave(x, patient) # within-patient mean
# Scaled random error within each SD_ratio_i group
re_std <- scale(resid(lm(rnorm(Nb*Nw) ~ x)))
re_between <- ave(re_std, patient)
re_within <- re_std - re_between
re_between <- scale(re_between)/sqrt(Nw)
re_within <- scale(re_within)
})
Data <- do.call(
rbind,
lapply(
1:Nv,
function(i) {
cbind(Data_template, SD_ratio_i = i)
}
)
)
Data <- within(
Data,
{
SD_within_ <- SD_within[SD_ratio_i]
SD_between_ <- SD_between[SD_ratio_i]
SD_ratio <- SD_ratio[SD_ratio_i]
y <- 10 +
Bb * xm + # contextual effect
Bw * (x - xm) + # within-patient effect
SD_within_ * re_within + # within patient random effect
SD_between_ * re_between # adjustment to between patient random effect
}
)Here is a scatterplot of the data sets showing estimated 50% concentration ellipses for each cluster:5
library("lattice")
library("latticeExtra")
plot <- xyplot(y ~ x | SD_ratio, data = Data, group = patient,
layout = c(Nv, 1),
par.settings = list(superpose.symbol = list(pch = 1, cex = 0.7))) +
layer(panel.ellipse(..., center.pch = 16, center.cex = 1.5,
level = 0.5),
panel.abline(a = 10, b = 1))
plot # display graph
Data sets showing identical within-patient configurations with increasing between-patient variance. The labels above each panel show the between-patient SD/within-patient SD. The ellipses are estimated 50% concentration ellipses for each patient. The population between-patient regression line, \(E(y) = 10 + x\), is shown in each panel.
The population from which these datasets is generated has a
between-patient effect of dosage of 1 and a within-patient effect of
\(-0.5\). The researcher may attempt to
obtain an estimate of the within-patient effect of dosage through the
use of mixed models with a random intercept. We will illustrate how this
approach results in estimates that are highly sensitive to the relative
variance at the two levels of the mixed model by considering four
models, all of which include a random intercept. The first three models
have sequentially nested fixed effects: (1)~ 1, intercept
only; (2)~ 1 + x, intercept and effect of x;
and (3) ~ 1 + x + xm, intercept, effect of x,
and a contextual variable, xm, consisting of the
within-patient mean of x. The fourth model,
~ 1 + I(x - xm), uses an intercept and the
centered-within-group variable x - xm. We thus fit four
models to four datasets for a total of 16 models.
model.formulas <- c(
' ~ 1' = y ~ 1 + (1 | patient),
'~ 1 + x' = y ~ 1 + x + (1 | patient),
'~ 1 + x + xm' = y ~ 1 + x + xm + (1 | patient),
'~ 1 + I(x - xm)' = y ~ 1 + I(x - xm) + (1 | patient)
)
fits <- lapply(split(Data, ~ SD_ratio),
function(d) {
lapply(model.formulas, function(form) {
glmmTMB(form, data = d)
})
})We proceed to obtain predictions from each model based on fixed effects alone, as would be used for cross-validation based on clusters (i.e., patients), and for fixed and random effects—so-called best linear unbiased predictions or “BLUPs”—as would be used for cross-validation based on cases (i.e., occasions within patients):
# predicted fixed and random effects:
pred.BLUPs <- lapply(fits, lapply, predict)
# predicted fixed effects:
pred.fixed <- lapply(fits, lapply, predict, re.form = ~0) We then prepare the data for plotting:
Dataf <- lapply(split(Data, ~ SD_ratio),
function(d) {
lapply(names(model.formulas),
function(form) cbind(d, formula = form))
}) |>
lapply(function(dlist) do.call(rbind, dlist)) |>
do.call(rbind, args = _)
Dataf <- within(
Dataf,
{
pred.fixed <- unlist(pred.fixed)
pred.BLUPs <- unlist(pred.BLUPs)
panel <- factor(formula, levels = c(names(model.formulas), 'data'))
}
)
Data$panel <- factor('data', levels = c(names(model.formulas), 'data'))The fixed-effects predictions from these models are shown in the following graph:
{
xyplot(y ~ x |SD_ratio * panel, Data,
groups = patient, type = 'n',
par.strip.text = list(cex = 0.7),
drop.unused.levels = FALSE) +
glayer(panel.ellipse(..., center.pch = 16, center.cex = 0.5,
level = 0.5),
panel.abline(a = 10, b = 1)) +
xyplot(pred.fixed ~ x |SD_ratio * panel, Dataf, type = 'l',
groups = patient,
drop.unused.levels = F,
ylab = 'fixed-effect predictions')
}|>
useOuterStrips() |> print()
Fixed-effect predictions using each model applied to data sets with varying between- and within-patient variance ratios. The top row shows summaries of the within-patient data using estimated 50% concentration ellipses for each patient.
The BLUPs from these models are shown in the following graph:
{
xyplot(y ~ x | SD_ratio * panel, Data,
groups = patient, type = 'n',
drop.unused.levels = F,
par.strip.text = list(cex = 0.7),
ylab = 'fixed- and random-effect predictions (BLUPS)') +
glayer(panel.ellipse(..., center.pch = 16, center.cex = 0.5,
level = 0.5),
panel.abline(a = 10, b = 1)) +
xyplot(pred.BLUPs ~ x | SD_ratio * panel, Dataf, type = 'l',
groups = patient,
drop.unused.levels = F)
}|>
useOuterStrips() |> print()
Fixed- and random-effect predictions (BLUPs) using each model applied to data sets with varying between- and within-patient variance ratios. The top row shows summaries of the within-patient data using estimated 50% concentration ellipses for each patient.
Data sets with relatively low between-patient variance result in
strong shrinkage of fixed-effects predictions and also of BLUPs towards
the between-patient relationship between y and
x in the model with an intercept and x as
fixed-effect predictors, ~ 1 + x. The inclusion of a
contextual variable in the model corrects this problem.
Although the BLUPs fit the observed data more closely than
predictions based on fixed effects alone, the slopes of within-patient
BLUPs do not conform with the within-patient slopes for the
~ 1 + x model in the two datasets with the smallest
between-patient variances.
For data with a small between-patient variance, fixed-effects
predictions for the ~ 1 + x model have a slope that is
close to the between-patient slope but provide better overall
predictions than the fixed-effect predictions for datasets with larger
between-subject variance. With data whose between-patient variance is
relatively large, predictions based on the model with a common intercept
and slope for all clusters, are very poor—indeed, much worse than the
fixed-effects-only predictions based on the simpler random-intercept
model.
We therefore anticipate (and show later in this section) that
case-based cross-validation may prefer the intercept-only model,
~ 1 to the larger ~ 1 + x model when the
between-cluster variance is relatively small, but that cluster-based
cross-validation will prefer the latter to the former.
We will discover that case-based cross-validation prefers the
~ 1 + x model to the ~ 1 model for the ‘5 /
2.5’ dataset, but cluster-based cross-validation prefers the latter
model to the former. The situation is entirely reversed with the ‘8 /
2.5’ dataset.
The third model, ~ 1 + x + xm, includes a contextual
effect of x—that is, the cluster mean xm—along
with x and the intercept in the fixed-effect part of the
model, along with a random intercept. This model is equivalent to
fitting y ~ I(x - xm) + xm + (1 | patient), which is the
model that generated the data. The fit of the mixed model
~ 1 + x + xm is consequently similar to that of a
fixed-effects only model with x and a categorical predictor
for individual patients (i.e., y ~ factor(patient) + x,
treating patients as a factor, and not shown here).
We next carry out case-based cross-validation, which, as we have
explained, is based on both fixed and predicted random effects (i.e.,
BLUPs), and cluster-based cross-validation, which is based on fixed
effects only. In order to reduce between-model random variability in
comparisons of models on the same dataset, we apply cv() to
the list of models created by the models() function
(introduced previously), performing cross-validation with the same folds
for each model.
model_lists <- lapply(fits, function(fitlist) do.call(models, fitlist))
cvs_cases <-
lapply(1:Nv,
function(i){
cv(model_lists[[i]], k = 10,
data = split(Data, ~ SD_ratio)[[i]])
})
#> R RNG seed set to 982470
#> R RNG seed set to 903914
#> R RNG seed set to 246168
#> R RNG seed set to 948407
cvs_clusters <-
lapply(1:Nv,
function(i){
cv(model_lists[[i]], k = 10,
data = split(Data, ~SD_ratio)[[i]],
clusterVariables = 'patient')
})
#> R RNG seed set to 382198
#> R RNG seed set to 533366
#> R RNG seed set to 637459
#> R RNG seed set to 542289For a given dataset, we can plot the results for a list of models
using the plot method for "cvList" objects. For example
10-fold cluster-based cross-validation comparing random intercept mixed models with varying fixed effects.
We assemble the results for all datasets to show them in a common figure.
names(cvs_clusters) <- names(cvs_cases) <- SD_ratio
dsummary <- expand.grid(SD_ratio_i = names(cvs_cases), model = names(cvs_cases[[1]]))
dsummary$cases <-
sapply(1:nrow(dsummary), function(i){
with(dsummary[i,], cvs_cases[[SD_ratio_i]][[model]][['CV crit']])
})
dsummary$clusters <-
sapply(1:nrow(dsummary), function(i){
with(dsummary[i,], cvs_clusters[[SD_ratio_i]][[model]][['CV crit']])
})xyplot(cases + clusters ~ model|SD_ratio_i, dsummary,
auto.key = list(space = 'top', reverse.rows = T, columns = 2), type = 'b',
xlab = "Fixed Effects",
ylab = 'CV criterion (MSE)',
layout= c(Nv,1),
par.settings =
list(superpose.line=list(lty = c(2, 3), lwd = 3),
superpose.symbol=list(pch = 15:16, cex = 1.5)),
scales = list(y = list(log = TRUE), x = list(alternating = F, rot = 60))) |> print()
10-fold cluster- and case-based cross-validation comparing mixed models with random intercepts and various fixed effects.
In summary, when between-cluster variance is relatively large, the
model ~ 1 + x, with x alone and without the
contextual mean of x, is assessed as fitting very poorly by
cluster-based CV, but relatively much better by case-based CV. In all
our examples, the model ~ 1 + x + xm, which includes both
x and its contextual mean, produces better results using
both cluster-based and case-based CV. These conclusions are consistent
with our observations based on graphing predictions from the various
models, and they illustrate the desirability of assessing mixed-effect
models at different hierarchical levels.
