Example: Contrived hierarchical data
We introduce an artificial data set that exemplifies aspects of cross-validation particular to hierarchical models. Using this data set, 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.
Consider a researcher studying improvement in a skill, yodeling, for example, among students enrolled in a four-year yodeling program. The plan is to measure each student’s skill level at the beginning of the program and every year thereafter until the end of the program, resulting in five annual measurements for each student. It turns out that yodeling appeals to students of all ages, and students enrolling in the program range in age from 20 to 70. Moreover, participants’ untrained yodeling skill is similar at all ages, as is their rate of progress with training. All students complete the four-year program.
The researcher, who has more expertise in yodeling than in modeling, decides to model the response, \(y\), yodeling skill, as a function of age, \(x\), reasoning that students get older during their stay in the program, and (incorrectly) that age can serve as a proxy for elapsed time. The researcher knows that a mixed model should be used to account for clustering due to the expected similarity of measurements taken from each student.
We start by generating the data, using parameters consistent with the description above and meant to highlight the issues that arise in cross-validating mixed-effects models:4
# Parameters:
set.seed(9693)
Nb <- 100 # number of groups
Nw <- 5 # number of individuals within groups
Bb <- 0 # between-group regression coefficient on group mean
SDre <-
2.0 # between-group SD of random level relative to group mean of x
SDwithin <- 0.5 # within group SD
Bw <- 1 # within group effect of x
Ay <- 10 # intercept for response
Ax <- 20 # starting level of x
Nx <- Nw * 10 # number of distinct x values
Data <- data.frame(group = factor(rep(1:Nb, each = Nw)),
x = Ax + rep(1:Nx, length.out = Nw * Nb)) |>
within({
xm <- ave(x, group, FUN = mean) # within-group mean
y <- Ay +
Bb * xm + # contextual effect
Bw * (x - xm) + # within-group effect
rnorm(Nb, sd = SDre)[group] + # random level by group
rnorm(Nb * Nw, sd = SDwithin) # random error within groups
})Here is a scatterplot of the data for a representative group of 10 (without loss of generality, the first 10) of 100 students, showing the 95% concentration ellipse for each cluster:5
library("lattice")
library("latticeExtra")
plot <- xyplot(
y ~ x,
data = Data[1:Nx,],
group = group,
ylim = c(4, 16),
par.settings = list(superpose.symbol = list(pch = 1, cex =
0.7))
) +
layer(panel.ellipse(..., center.cex = 0))
plot # display graph
Hierarchical data set, showing the first 10 of 100 students.
The between-student effect of age is 0 but the within-student effect is 1. Due to the large variation in ages between students, the least-squares regression of yodeling skill on age (for the 500 observations among all 100 students) produces an estimated slope close to 0 (though with a small \(p\)-value), because the slope is heavily weighted toward the between-student effect:
summary(lm(y ~ x, data=Data))
#>
#> Call:
#> lm(formula = y ~ x, data = Data)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -5.771 -1.658 -0.089 1.552 7.624
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 9.05043 0.34719 26.07 <2e-16 ***
#> x 0.02091 0.00727 2.87 0.0042 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 2.35 on 498 degrees of freedom
#> Multiple R-squared: 0.0163, Adjusted R-squared: 0.0143
#> F-statistic: 8.26 on 1 and 498 DF, p-value: 0.00422The initial mixed-effects model that we fit to the data is a simple random-intercepts model:
# random intercept only:
mod.0 <- lmer(y ~ 1 + (1 | group), Data)
summary(mod.0)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: y ~ 1 + (1 | group)
#> Data: Data
#>
#> REML criterion at convergence: 2103.1
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.0351 -0.7264 -0.0117 0.7848 2.0438
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> group (Intercept) 2.90 1.70
#> Residual 2.71 1.65
#> Number of obs: 500, groups: group, 100
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 10.002 0.186 53.9We will shortly consider three other, more complex, mixed models; because of data-management considerations, it is convenient to fit them now, but we defer discussion of these models:
# effect of x and random intercept:
mod.1 <- lmer(y ~ x + (1 | group), Data)
# effect of x, contextual (student) mean of x, and random intercept:
mod.2 <- lmer(y ~ x + xm + (1 | group), Data)
# equivalent to y ~ I(x - xm) + xm + (1 | group)
# model generating the data (where Bb = 0)
mod.3 <- lmer(y ~ I(x - xm) + (1 | group), Data)We proceed to obtain predictions from the random-intercept model
(mod.0) and the other models (mod.1,
mod.2, and mod.3) based on fixed effects
alone, as would be used for cross-validation based on clusters (i.e.,
students), 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 students):
Data <- within(Data, {
fit_mod0.fe <- predict(mod.0, re.form = ~ 0) # fixed effects only
fit_mod0.re <- predict(mod.0) # fixed and random effects (BLUPs)
fit_mod1.fe <- predict(mod.1, re.form = ~ 0)
fit_mod1.re <- predict(mod.1)
fit_mod2.fe <- predict(mod.2, re.form = ~ 0)
fit_mod2.re <- predict(mod.2)
fit_mod3.fe <- predict(mod.3, re.form = ~ 0)
fit_mod3.re <- predict(mod.3)
})We then prepare the data for plotting:
Data_long <- reshape(Data[1:Nx, ], direction = "long", sep = ".",
timevar = "effect", varying = grep("\\.", names(Data[1:Nx, ])))
Data_long$id <- 1:nrow(Data_long)
Data_long <- reshape(Data_long, direction = "long", sep = "_",
timevar = "modelcode", varying = grep("_", names(Data_long)))
Data_long$model <- factor(
c("~ 1", "~ 1 + x", "~ 1 + x + xm", "~ 1 + I(x - xm)")
[match(Data_long$modelcode, c("mod0", "mod1", "mod2", "mod3"))]
)Predictions based on the random-intercept model mod.0
for the first 10 students are shown in the following graph:
(
plot +
xyplot(
fit ~ x,
subset(Data_long, modelcode == "mod0" & effect == "fe"),
groups = group,
type = "l",
lwd = 2
) +
xyplot(
fit ~ x,
subset(Data_long, modelcode == "mod0" & effect == "re"),
groups = group,
type = "l",
lwd = 2,
lty = 3
)
) |> update(
main="Model: y ~ 1 + (1 | group)",
key=list(
corner=c(0.05, 0.05),
text=list(c("fixed effects only","fixed and random")),
lines=list(lty=c(1, 3))))
Predictions from the random intercept model.
The fixed-effect predictions for the various individuals are identical—the estimated fixed-effects intercept or estimated general mean of \(y\)—while the BLUPs are the sums of the fixed-effects intercept and the random intercepts, and are only slightly shrunken towards the general mean. Because in our artificial data there is no population relationship between age and skill, the fixed-effect-only predictions and the BLUPs are not very different.
Our next model, mod.1, includes a fixed intercept and
fixed effect of x along with a random intercept:
summary(mod.1)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: y ~ x + (1 | group)
#> Data: Data
#>
#> REML criterion at convergence: 1564.5
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.9016 -0.6350 0.0188 0.5541 2.8293
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> group (Intercept) 192.941 13.890
#> Residual 0.257 0.507
#> Number of obs: 500, groups: group, 100
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) -33.9189 1.5645 -21.7
#> x 0.9653 0.0158 61.0
#>
#> Correlation of Fixed Effects:
#> (Intr)
#> x -0.460Predictions from this model appear in the following graph:
(
plot +
xyplot(
fit ~ x,
subset(Data_long, modelcode == "mod1" & effect == "fe"),
groups = group,
type = "l",
lwd = 2
) +
xyplot(
fit ~ x,
subset(Data_long, modelcode == "mod1" & effect == "re"),
groups = group,
type = "l",
lwd = 2,
lty = 3
)
) |> update(
main="Model: y ~ 1 + x + (1 | group)",
ylim=c(-15, 35),
key=list(
corner=c(0.95, 0.05),
text=list(c("fixed effects only","fixed and random")),
lines=list(lty=c(1, 3))))
Predictions from the model with random intercepts and \(x\) as a fixed-effect predictor.
The BLUPs fit the observed data very closely, but predictions based
on the fixed effects alone, 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,
mod.0. We therefore anticipate (and show later in this
section) that case-based cross-validation will prefer mod1
to mod0, but that cluster-based cross-validation will
prefer mod0 to mod1.
Our third model, mod.2, includes the contextual effect
of \(x\)—that is, the cluster mean
xm—along with \(x\) and
the intercept in the fixed-effect part of the model, and a random
intercept:
summary(mod.2)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: y ~ x + xm + (1 | group)
#> Data: Data
#>
#> REML criterion at convergence: 1169.2
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.9847 -0.6375 0.0019 0.5568 2.7325
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> group (Intercept) 3.399 1.844
#> Residual 0.255 0.505
#> Number of obs: 500, groups: group, 100
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 9.4787 0.6171 15.4
#> x 0.9915 0.0160 62.1
#> xm -0.9800 0.0206 -47.7
#>
#> Correlation of Fixed Effects:
#> (Intr) x
#> x 0.000
#> xm -0.600 -0.777This model is equivalent to fitting
y ~ I(x - xm) + xm + (1 | group), which is the model that
generated the data once the coefficient of the contextual predictor
xm is set to 0 (as it is in mod.3, discussed
below).
Predictions from model mod.2 appear in the following
graph:
(
plot +
xyplot(
fit ~ x,
subset(Data_long, modelcode == "mod2" & effect == "fe"),
groups = group,
type = "l",
lwd = 2
) +
xyplot(
fit ~ x,
subset(Data_long, modelcode == "mod2" & effect == "re"),
groups = group,
type = "l",
lwd = 2,
lty = 3
)
) |> update(
main="Model: y ~ 1 + x + xm + (1 | group)",
ylim=c(4, 16),
key=list(
corner=c(0.05, 0.05),
text=list(c("fixed effects only","fixed and random")),
lines=list(lty=c(1, 3))))
Predictors from the model with random intercepts, \(x\), and the group (student) mean of \(x\) as predictors.
Depending on the estimated variance parameters of the model, a mixed
model like mod.2 will apply varying degrees of shrinkage to
the random-intercept BLUPs that correspond to variation in the heights
of the parallel fitted lines for the individual students. In our
contrived data, the mod.2 applies little shrinkage,
allowing substantial variability in the heights of the fitted lines,
which closely approach the observed values for each student. The fit of
the mixed model mod.2 is consequently similar to that of a
fixed-effects model with age and a categorical predictor for individual
students (i.e., treating students as a factor, and not shown here).
The mixed model mod.2 therefore fits individual
observations well, and we anticipate a favorable assessment using
individual-based cross-validation. In contrast, the large variability in
the BLUPs results in larger residuals for predictions based on fixed
effects alone, and so we expect that cluster-based cross-validation
won’t show an advantage for model mod.2 compared to the
smaller model mod.0, which includes only fixed and random
intercepts.
Had the mixed model applied considerable shrinkage, then neither cluster-based nor case-based cross-validation would show much improvement over the random-intercept-only model. In our experience, the degree of shrinkage does not vary smoothly as parameters are changed but tends to be “all or nothing,” and near the tipping point, the behavior of estimates can be affected considerably by the choice of algorithm used to fit the model.
Finally, mod.3 directly estimates the model used to
generate the data. As mentioned, it is a constrained version of
mod.2, with the coefficient of xm set to 0,
and with x expressed as a deviation from the cluster mean
xm:
summary(mod.3)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: y ~ I(x - xm) + (1 | group)
#> Data: Data
#>
#> REML criterion at convergence: 1163.2
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -2.9770 -0.6320 0.0063 0.5603 2.7249
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> group (Intercept) 3.391 1.842
#> Residual 0.255 0.505
#> Number of obs: 500, groups: group, 100
#>
#> Fixed effects:
#> Estimate Std. Error t value
#> (Intercept) 10.002 0.185 53.9
#> I(x - xm) 0.992 0.016 62.1
#>
#> Correlation of Fixed Effects:
#> (Intr)
#> I(x - xm) 0.000The predictions from mod.3 are therefore similar to
those from mod.2:
(
plot +
xyplot(
fit ~ x,
subset(Data_long, modelcode == "mod3" & effect == "fe"),
groups = group,
type = "l",
lwd = 2
) +
xyplot(
fit ~ x,
subset(Data_long, modelcode == "mod3" & effect == "re"),
groups = group,
type = "l",
lwd = 2,
lty = 3
)
) |> update(
main="Model: y ~ 1 + I(x - xm) + (1 | group)",
ylim=c(4, 16),
key=list(
corner=c(0.05, 0.05),
text=list(c("fixed effects only","fixed and random")),
lines=list(lty=c(1, 3))))
Predictions from the estimated model generating the data.
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, 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:
modlist <- models(
"~ 1" = mod.0,
"~ 1 + x" = mod.1,
"~ 1 + x + xm" = mod.2,
"~ 1 + I(x - xm)" = mod.3
)
cvs_clusters <-
cv(
modlist,
data = Data,
cluster = "group",
k = 10,
seed = 6449
)
plot(cvs_clusters, main = "Model Comparison, Cluster-Based CV")
10-fold cluster-based cross-validation comparing random intercept models with varying fixed effects. The error bars show the 95% confidence interval around the CV estimate of the MSE for each model.
cvs_cases <- cv(modlist, data = Data, seed = 9693)
plot(cvs_cases, main = "Model Comparison, Case-Based CV")
10-fold case-based cross-validation comparing random intercept models with varying fixed effects.
In summary, model mod.1, 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.
Model mod.2, which includes both \(x\) and its contextual mean, produces
better results using both cluster-based and case-based CV. The
data-generating model, mod.3, which includes the fixed
effect of x - xm in place of separate terms in
x and xm, isn’t distinguishable from model
mod.2, which includes x and xm
separately, even though mod.2 has an unnecessary parameter
(recall that the population coefficient of xm is 0 when
x is expressed as deviations from the contextual mean).
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.
