Adding a model class not covered by the default cv()
method
Independently sampled cases
Suppose that we want to cross-validate a multinomial logistic
regression model fit by the multinom() function in the
nnet package (Venables &
Ripley, 2002). We borrow an example from Fox (2016, sec. 14.2.1), with data from the
British Election Panel Study on vote choice in the 2001 British
election. Data for the example are in the BEPS data frame
in the carData package:
data("BEPS", package="carData")
head(BEPS)
#> vote age economic.cond.national economic.cond.household Blair
#> 1 Liberal Democrat 43 3 3 4
#> 2 Labour 36 4 4 4
#> 3 Labour 35 4 4 5
#> 4 Labour 24 4 2 2
#> 5 Labour 41 2 2 1
#> 6 Labour 47 3 4 4
#> Hague Kennedy Europe political.knowledge gender
#> 1 1 4 2 2 female
#> 2 4 4 5 2 male
#> 3 2 3 3 2 male
#> 4 1 3 4 0 female
#> 5 1 4 6 2 male
#> 6 4 2 4 2 maleThe polytomous (multi-category) response variable is
vote, a factor with levels "Conservative",
"Labour", and "Liberal Democrat". The
predictors of vote are:
age, in years;econ.cond.nationalandecon.cond.household, the respondent’s ratings of the state of the economy, on 1 to 5 scales.Blair,Hague, andKennedy, ratings of the leaders of the Labour, Conservative, and Liberal Democratic parties, on 1 to 5 scales.Europe, an 11-point scale on attitude towards European integration, with high scores representing “Euro-skepticism.”political.knowledge, knowledge of the parties’ positions on European integration, with scores from 0 to 3.gender,"female"or"male".
The model fit to the data includes an interaction between
Europe and political.knowledge; the other
predictors enter the model additively:
library("nnet")
m.beps <- multinom(
vote ~ age + gender + economic.cond.national +
economic.cond.household + Blair + Hague + Kennedy +
Europe * political.knowledge,
data = BEPS
)
#> # weights: 36 (22 variable)
#> initial value 1675.383740
#> iter 10 value 1240.047788
#> iter 20 value 1163.199642
#> iter 30 value 1116.519687
#> final value 1116.519666
#> converged
car::Anova(m.beps)
#> Analysis of Deviance Table (Type II tests)
#>
#> Response: vote
#> LR Chisq Df Pr(>Chisq)
#> age 13.9 2 0.00097 ***
#> gender 0.5 2 0.79726
#> economic.cond.national 30.6 2 2.3e-07 ***
#> economic.cond.household 5.7 2 0.05926 .
#> Blair 135.4 2 < 2e-16 ***
#> Hague 166.8 2 < 2e-16 ***
#> Kennedy 68.9 2 1.1e-15 ***
#> Europe 78.0 2 < 2e-16 ***
#> political.knowledge 55.6 2 8.6e-13 ***
#> Europe:political.knowledge 50.8 2 9.3e-12 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Most of the predictors, including the Europe \(\times\) political.knowledge
interaction, are associated with very small \(p\)-values; the Anova()
function is from the car package (Fox & Weisberg, 2019).
Here’s an “effect plot”, using the the effects
package (Fox & Weisberg, 2019) to
visualize the Europe \(\times\) political.knowledge
interaction in a “stacked-area” graph:
plot(
effects::Effect(
c("Europe", "political.knowledge"),
m.beps,
xlevels = list(Europe = 1:11, political.knowledge = 0:3),
fixed.predictors = list(given.values = c(gendermale = 0.5))
),
lines = list(col = c("blue", "red", "orange")),
axes = list(x = list(rug = FALSE), y = list(style = "stacked"))
)
To cross-validate this multinomial-logit model we need an appropriate
cost criterion. None of the criteria supplied by the cv
package—for example, neither mse(), which is appropriate
for a numeric response, nor BayesRule(), which is
appropriate for a binary response—will do. One possibility is to adapt
Bayes rule to a polytomous response:
head(BEPS$vote)
#> [1] Liberal Democrat Labour Labour Labour
#> [5] Labour Labour
#> Levels: Conservative Labour Liberal Democrat
yhat <- predict(m.beps, type = "class")
head(yhat)
#> [1] Labour Labour Labour Labour
#> [5] Liberal Democrat Labour
#> Levels: Conservative Labour Liberal Democrat
BayesRuleMulti <- function(y, yhat) {
result <- mean(y != yhat)
attr(result, "casewise loss") <- "y != yhat"
result
}
BayesRuleMulti(BEPS$vote, yhat)
#> [1] 0.31869
#> attr(,"casewise loss")
#> [1] "y != yhat"The predict() method for "multinom" models
called with argument type="class" reports the Bayes-rule
prediction for each case—that is, the response category with the highest
predicted probability. Our BayesRuleMulti() function
calculates the proportion of misclassified cases. Because this value is
the mean of casewise components, we attach a
"casewise loss" attribute to the result (as explained in
the preceding section).
The marginal proportions for the response categories are
xtabs(~ vote, data=BEPS)/nrow(BEPS)
#> vote
#> Conservative Labour Liberal Democrat
#> 0.30295 0.47213 0.22492and so the marginal Bayes-rule prediction, that everyone will vote Labour, produces an error rate of \(1 - 0.47213 = 0.52787\). The multinomial-logit model appears to do substantially better than that, but does its performance hold up to cross-validation?
We check first whether the default cv() method works
“out-of-the-box” for the "multinom" model:
cv(m.beps, seed=3465, criterion=BayesRuleMulti)
#> Error in GetResponse.default(model): non-vector responseThe default method of GetResponse() (a function supplied
by the cv package—see ?GetResponse) fails
for a "multinom" object. A straightforward solution is to
supply a GetResponse.multinom() method that returns the
factor response (using the
get_response() function from the insight
package, Lüdecke, Waggoner, & Makowski, 2019),
GetResponse.multinom <- function(model, ...) {
insight::get_response(model)
}
head(GetResponse(m.beps))
#> [1] Liberal Democrat Labour Labour Labour
#> [5] Labour Labour
#> Levels: Conservative Labour Liberal Democratand to try again:
cv(m.beps, seed=3465, criterion=BayesRuleMulti)
#> R RNG seed set to 3465
#> # weights: 36 (22 variable)
#> initial value 1507.296060
#> iter 10 value 1134.575036
#> iter 20 value 1037.413231
#> iter 30 value 1007.705242
#> iter 30 value 1007.705235
#> iter 30 value 1007.705235
#> final value 1007.705235
#> converged
#> Error in match.arg(type): 'arg' should be one of "class", "probs"A traceback() (not shown) reveals that the problem is
that the default method of cv() calls the
"multinom" method for predict() with the
argument type="response", when the correct argument should
be type="class". We therefore must write a
“multinom” method for cv(), but that proves to
be very simple:
cv.multinom <-
function (model, data, criterion = BayesRuleMulti, k, reps,
seed, ...) {
model <- update(model, trace = FALSE)
NextMethod(
type = "class",
criterion = criterion,
criterion.name = deparse(substitute(criterion))
)
}That is, we simply call the default cv() method with the
type argument properly set. In addition to supplying the
correct type argument, our method sets the default
criterion for the cv.multinom() method to
BayesRuleMulti. Adding the argument
criterion.name=deparse(substitute(criterion)) is
inessential, but it insures that printed output will include the name of
the criterion function that’s employed, whether it’s the default
BayesRuleMulti or something else. Prior to invoking
NextMethod(), we called update() with
trace=FALSE to suppress the iteration history reported by
default by multinom()—it would be tedious to see the
iteration history for each fold.
Then:
summary(cv(m.beps, seed=3465))
#> R RNG seed set to 3465
#> 10-Fold Cross Validation
#> criterion: BayesRuleMulti
#> cross-validation criterion = 0.32459
#> bias-adjusted cross-validation criterion = 0.32368
#> 95% CI for bias-adjusted CV criterion = (0.30017, 0.34718)
#> full-sample criterion = 0.31869The cross-validated polytomous Bayes-rule criterion confirms that the fitted model does substantially better than the marginal Bayes-rule prediction that everyone votes for Labour.
Calling cvCompute()
cv() methods for independently sampled cases, such as
cv.default(), cv.lm(), and
cv.glm(), work by setting up calls to the
cvCompute() function, which is exported from the
cv package to support development of cv()
methods for additional classes of regression models. In most cases,
however, such as the preceding cv.multinom() example, it
will suffice and be much simpler to set up a suitable call to
cv.default() via NextMethod().
To illustrate how to use cvCompute() directly, we write
an alternative, and necessarily more complicated, version of
cv.multinom().
cv.multinom <- function(model,
data = insight::get_data(model),
criterion = BayesRuleMulti,
k = 10,
reps = 1,
seed = NULL,
details = k <= 10,
confint = n >= 400,
level = 0.95,
ncores = 1,
start = FALSE,
...) {
f <- function(i) {
# helper function to compute to compute fitted values,
# etc., for each fold i
indices.i <- fold(folds, i)
model.i <- if (start) {
update(model,
data = data[-indices.i,],
start = b,
trace = FALSE)
} else {
update(model, data = data[-indices.i,], trace = FALSE)
}
fit.all.i <- predict(model.i, newdata = data, type = "class")
fit.i <- fit.all.i[indices.i]
# returns:
# fit.i: fitted values for the i-th fold
# crit.all.i: CV criterion for all cases based on model with
# i-th fold omitted
# coef.i: coefficients for the model with i-th fold omitted
list(
fit.i = fit.i,
crit.all.i = criterion(y, fit.all.i),
coef.i = coef(model.i)
)
}
fPara <- function(i, multinom, ...) {
# helper function for parallel computation
# argument multinom makes multinom() locally available
# ... is necessary but not used
indices.i <- fold(folds, i)
model.i <- if (start) {
update(model,
data = data[-indices.i,],
start = b,
trace = FALSE)
} else {
update(model, data = data[-indices.i,], trace = FALSE)
}
fit.all.i <- predict(model.i, newdata = data, type = "class")
fit.i <- fit.all.i[indices.i]
list(
fit.i = fit.i,
crit.all.i = criterion(y, fit.all.i),
coef.i = coef(model.i)
)
}
n <- nrow(data)
# see ?cvCompute for definitions of arguments
cvCompute(
model = model,
data = data,
criterion = criterion,
criterion.name = deparse(substitute(criterion)),
k = k,
reps = reps,
seed = seed,
details = details,
confint = confint,
level = level,
ncores = ncores,
type = "class",
start = start,
f = f,
fPara = fPara,
multinom = nnet::multinom
)
}Notice that separate “helper” functions are defined for non-parallel
and parallel computations.1 The new version of
cv.multinom() produces the same results as the version that
calls cv.default():2
summary(cv(m.beps, seed=3465))
#> R RNG seed set to 3465
#> 10-Fold Cross Validation
#> criterion: BayesRuleMulti
#> cross-validation criterion = 0.32459
#> bias-adjusted cross-validation criterion = 0.32368
#> 95% CI for bias-adjusted CV criterion = (0.30017, 0.34718)
#> full-sample criterion = 0.31869Mixed-effects models
Adding a cv() method for a mixed-model class is somewhat
more complicated. We provide the cvMixed() function to
facilitate this process, and to see how that works, consider the
"lme" method from the cv package:
cv:::cv.lme
#> function (model, data = insight::get_data(model), criterion = mse,
#> k = NULL, reps = 1L, seed, details = NULL, ncores = 1L, clusterVariables,
#> ...)
#> {
#> cvMixed(model, package = "nlme", data = data, criterion = criterion,
#> criterion.name = deparse(substitute(criterion)), k = k,
#> reps = reps, seed = seed, details = details, ncores = ncores,
#> clusterVariables = clusterVariables, predict.clusters.args = list(object = model,
#> newdata = data, level = 0), predict.cases.args = list(object = model,
#> newdata = data, level = 1), fixed.effects = nlme::fixef,
#> ...)
#> }
#> <bytecode: 0x113913378>
#> <environment: namespace:cv>Notice that cv.lme() sets up a call to
cvMixed(), which does the computational work.
Most of the arguments of cvMixed() are familiar:
modelis the mixed-model object, here of class"lme".packageis the name of the package in which the mixed-modeling function used to fit the model, herelme(), resides—i.e.,"nlme";cvMixed()uses this argument to retrieve the package namespace.datais the data set to which the model is fit, by default extracted by theget_data()function in the insight package.criterionis the CV criterion, defaulting to themse()function.kis the number of CV folds, defaulting to"loo"for CV by clusters and10for CV by cases.repsis the number of times the CV process is repeated, defaulting to1.seedis the seed for R’s random-number generator, defaulting to a randomly selected (and saved) value.ncoresis the number of cores to use for parallel computation; if1, the default, then the computation isn’t parallelized.clusterVariablesis a character vector of the names of variables defining clusters; if missing, then CV is based on cases rather than clusters.
The remaining arguments are unfamiliar:
predict.clusters.argsis a named list of arguments to be passed to thepredict()function to obtain predictions for the full data set from a model fit to a subset of the data for cluster-based CV. The first two arguments should beobjectandnewdata. It is typically necessary to tellcvMixed()how to base predictions only on fixed effects; in the case of"lme"models, this is done by settinglevel = 0.Similarly,
predict.cases.argsis a named list of arguments to be passed topredict()for case-based CV. Settinglevel = 1includes random effects in the predictions.fixed.effectsis used to compute detailed fold-based statistics.
Finally, any additional arguments, absorbed by ..., are
passed to update() when the model is refit with each fold
omitted. cvMixed() returns an object of class
"cv".
Now imagine that we want to support a new class of mixed-effects
models. To be concrete, we illustrate with the glmmPQL()
function in the MASS package (Venables & Ripley, 2002), which fits
generalized-linear mixed-effects models by penalized quasi-likelihood.3 Not
coincidentally, the arguments of glmmPQL() are similar to
those of lme() (with an additional family
argument), because the former iteratively invokes the latter; so
cv.glmmPQL() should resemble cv.lme().
As it turns out, neither the default method for
GetResponse() nor insight::get_data() work for
"glmmPQL" objects. These objects include a
"data" element, however, and so we can simply extract this
element as the default for the data argument of our
cv.glmmPQL() method.
To get the response variable is more complicated: We refit the fixed part of the model as a GLM with only the regression constant on the right-hand side, and extract the response from that; because all we need is the response variable, we limit the number of GLM iterations to 1 and suppress warning messages about non-convergence:
GetResponse.glmmPQL <- function(model, ...) {
f <- formula(model)
f[[3]] <- 1 # regression constant only on RHS
model <-
suppressWarnings(glm(
f,
data = model$data,
family = model$family,
control = list(maxit = 1)
))
cv::GetResponse(model)
}Writing the cv() method is then straightforward:
cv.glmmPQL <- function(model,
data = model$data,
criterion = mse,
k,
reps = 1,
seed,
ncores = 1,
clusterVariables,
...) {
cvMixed(
model,
package = "MASS",
data = data,
criterion = criterion,
k = k,
reps = reps,
seed = seed,
ncores = ncores,
clusterVariables = clusterVariables,
predict.clusters.args = list(
object = model,
newdata = data,
level = 0,
type = "response"
),
predict.cases.args = list(
object = model,
newdata = data,
level = 1,
type = "response"
),
fixed.effects = nlme::fixef,
verbose = FALSE,
...
)
}We set the argument verbose=FALSE to suppress
glmmPQL()’s iteration counter when cvMixed()
calls update().
Let’s apply our newly minted method to a logistic regression with a
random intercept in an example that appears in
?glmmPQL:
library("MASS")
m.pql <- glmmPQL(
y ~ trt + I(week > 2),
random = ~ 1 | ID,
family = binomial,
data = bacteria
)
#> iteration 1
#> iteration 2
#> iteration 3
#> iteration 4
#> iteration 5
#> iteration 6
summary(m.pql)
#> Linear mixed-effects model fit by maximum likelihood
#> Data: bacteria
#> AIC BIC logLik
#> NA NA NA
#>
#> Random effects:
#> Formula: ~1 | ID
#> (Intercept) Residual
#> StdDev: 1.4106 0.78005
#>
#> Variance function:
#> Structure: fixed weights
#> Formula: ~invwt
#> Fixed effects: y ~ trt + I(week > 2)
#> Value Std.Error DF t-value p-value
#> (Intercept) 3.4120 0.51850 169 6.5805 0.0000
#> trtdrug -1.2474 0.64406 47 -1.9367 0.0588
#> trtdrug+ -0.7543 0.64540 47 -1.1688 0.2484
#> I(week > 2)TRUE -1.6073 0.35834 169 -4.4853 0.0000
#> Correlation:
#> (Intr) trtdrg trtdr+
#> trtdrug -0.598
#> trtdrug+ -0.571 0.460
#> I(week > 2)TRUE -0.537 0.047 -0.001
#>
#> Standardized Within-Group Residuals:
#> Min Q1 Med Q3 Max
#> -5.19854 0.15723 0.35131 0.49495 1.74488
#>
#> Number of Observations: 220
#> Number of Groups: 50We compare this result to that obtained from glmer() in
the lme4 package:
library("lme4")
#> Loading required package: Matrix
m.glmer <- glmer(y ~ trt + I(week > 2) + (1 | ID),
family = binomial, data = bacteria)
summary(m.glmer)
#> Generalized linear mixed model fit by maximum likelihood (Laplace
#> Approximation) [glmerMod]
#> Family: binomial ( logit )
#> Formula: y ~ trt + I(week > 2) + (1 | ID)
#> Data: bacteria
#>
#> AIC BIC logLik deviance df.resid
#> 202.3 219.2 -96.1 192.3 215
#>
#> Scaled residuals:
#> Min 1Q Median 3Q Max
#> -4.561 0.136 0.302 0.422 1.128
#>
#> Random effects:
#> Groups Name Variance Std.Dev.
#> ID (Intercept) 1.54 1.24
#> Number of obs: 220, groups: ID, 50
#>
#> Fixed effects:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.548 0.696 5.10 3.4e-07 ***
#> trtdrug -1.367 0.677 -2.02 0.04352 *
#> trtdrug+ -0.783 0.683 -1.15 0.25193
#> I(week > 2)TRUE -1.598 0.476 -3.36 0.00078 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Correlation of Fixed Effects:
#> (Intr) trtdrg trtdr+
#> trtdrug -0.593
#> trtdrug+ -0.537 0.487
#> I(wk>2)TRUE -0.656 0.126 0.064
# comparison of fixed effects:
car::compareCoefs(m.pql, m.glmer)
#> Warning in car::compareCoefs(m.pql, m.glmer): models to be compared are of
#> different classes
#> Calls:
#> 1: glmmPQL(fixed = y ~ trt + I(week > 2), random = ~1 | ID, family =
#> binomial, data = bacteria)
#> 2: glmer(formula = y ~ trt + I(week > 2) + (1 | ID), data = bacteria,
#> family = binomial)
#>
#> Model 1 Model 2
#> (Intercept) 3.412 3.548
#> SE 0.514 0.696
#>
#> trtdrug -1.247 -1.367
#> SE 0.638 0.677
#>
#> trtdrug+ -0.754 -0.783
#> SE 0.640 0.683
#>
#> I(week > 2)TRUE -1.607 -1.598
#> SE 0.355 0.476
#> The two sets of estimates are similar, but not identical
Finally, we try out our cv.glmmPQL() method,
cross-validating both by clusters and by cases,
summary(cv(m.pql, clusterVariables="ID", criterion=BayesRule))
#> n-Fold Cross Validation based on 50 {ID} clusters
#> cross-validation criterion = 0.19545
#> bias-adjusted cross-validation criterion = 0.19545
#> full-sample criterion = 0.19545
summary(cv(m.pql, data=bacteria, criterion=BayesRule, seed=1490))
#> R RNG seed set to 1490
#> 10-Fold Cross Validation
#> cross-validation criterion = 0.20909
#> bias-adjusted cross-validation criterion = 0.20727
#> full-sample criterion = 0.14545and again compare to glmer():
summary(cv(m.glmer, clusterVariables="ID", criterion=BayesRule))
#> n-Fold Cross Validation based on 50 {ID} clusters
#> criterion: BayesRule
#> cross-validation criterion = 0.19545
#> bias-adjusted cross-validation criterion = 0.19545
#> full-sample criterion = 0.19545
summary(cv(m.glmer, data=bacteria, criterion=BayesRule, seed=1490))
#> R RNG seed set to 1490
#> 10-Fold Cross Validation
#> criterion: BayesRule
#> cross-validation criterion = 0.19545
#> bias-adjusted cross-validation criterion = 0.19364
#> full-sample criterion = 0.15