Think before you shrink: Alternatives to default shrinkage methods can improve
prediction accuracy, calibration and coverage

Mark A. van de Wiel1, Gwenaël G.R. Leday2, Jeroen Hoogland1, Martijn W.
Heymans1, Erik W. van Zwet3, Ailko H. Zwinderman1

1Dept of Epidemiology and Data Science, Amsterdam Public Health research
institute, Amsterdam University Medical Centers, Amsterdam, The Netherlands;
Biometris, Wageningen University and Research, Wageningen, The Netherlands;
3Dept of Biomedical Data Sciences, Leiden University Medical Center, Leiden,
The Netherlands

###### Abstract

While shrinkage is essential in high-dimensional settings, its use for low-
dimensional regression-based prediction has been debated. It reduces variance,
often leading to improved prediction accuracy. However, it also inevitably
introduces bias, which may harm two other measures of predictive performance:
calibration and coverage of confidence intervals. Much of the criticism stems
from the usage of standard shrinkage methods, such as lasso and ridge with a
single, cross-validated penalty. Our aim is to show that readily available
alternatives can strongly improve predictive performance, in terms of
accuracy, calibration or coverage. For linear regression, we use small sample
splits of a large, fairly typical epidemiological data set to illustrate this.
We show that usage of differential ridge penalties for covariate groups may
enhance prediction accuracy, while calibration and coverage benefit from
additional shrinkage of the penalties. In the logistic setting, we apply an
external simulation to demonstrate that local shrinkage improves calibration
with respect to global shrinkage, while providing better prediction accuracy
than other solutions, like Firth’s correction. The benefits of the alternative
shrinkage methods are easily accessible via example implementations using mgcv
and r-stan, including the estimation of multiple penalties. A synthetic copy
of the large data set is shared for reproducibility.

##  1 Introduction

Shrinkage has become a standard technique in statistics to counter over-
fitting in regression models. In particular in high-dimensional settings, with
the number of variables $p$ larger than sample size $n$, application of
shrinkage is necessary to obtain parameter estimates and predictions. In low-
dimensional settings, the benefit of shrinkage depends on the $n:p$ ratio, as
well as the main purpose of the analysis, e.g. prediction, parameter
estimation, variable selection or causal effect estimation. Here, we focus on
prediction. In general, the consensus seems to be that, on average, shrinkage
improves prediction accuracy on test data from the same population as the
training data Van Calster et al. (2020); Riley et al. (2021). This is no
surprise as the penalty parameter is usually tuned, e.g. by cross-validation,
to maximize out-of-bag prediction accuracy. Hence, it will adapt to the $n:p$
ratio: if large, little shrinkage is necessary; if small, more shrinkage is
required. The use of shrinkage, however, comes at a price. It inevitably
biases the parameter estimates, which in turn may lead to bad calibration of
the prediction Van Calster et al. (2020). Moreover, the penalty parameter can
be rather instable Riley et al. (2021), and such instability is often not
communicated when presenting resulting models. Most empirical results that
support these critiques on shrinkage are based on standard penalization
techniques like lasso and ridge. While these issues are to some extent
intrinsic to any shrinkage method, we show that they can be substantially
alleviated by using alternatives that are only slightly more advanced.
Therefore, our aim is to show a broad audience of applied statisticians and
epidemiologists that some ‘non-standard’ shrinkage methods may be very useful
for typical epidemiological studies that target multivariable prognostic
models.

In terms of shrinkage we focus mainly on several variants of ridge regression.
We do briefly discuss comparison of prediction accuracy with classical lasso
and stepwise selection (which could be regarded as an extremely bimodal type
of shrinkage) for the main data example, showing the latter two are inferior
to the ridge variations for this purpose. We emphasize that the latter two
cater mainly for variable selection, rendering the comparison with ridge-type
models somewhat unfair. Therefore, we focus on the latter.

Our setting is regression, where shrinkage is either effectuated by using a
penalty for the parameters or by using an informative prior on these
parameters in a Bayesian setting. In the latter case, the scale parameter of
the prior distribution acts effectively as a penalty parameter, as the mean of
the prior is usually set to zero. First, we explore differential shrinkage by
the use of different penalties for (groups of) variables. We argue that in
many settings, it is in fact unreasonable to assume the same penalty for all
variables. We show that both prediction accuracy and calibration may be
greatly enhanced by using differential penalties. In addition, the use of
differential penalties provides a more objective solution to down-weighting a
set of variables than simply discarding it prior to the analysis, as is
sometimes suggested as a (partial) solution for a low $n:p$ ratio Van Calster
et al. (2020). Importantly, these penalties can be estimated automatically
without (time-consuming) cross-validation in both the classical and full
Bayesian setting, using R-packages mgcv Wood (2011) and R-Stan Stan
Development Team (2022), respectively.

Global shrinkage methods, such as ridge, penalize all parameters equally. If
the primary aim is to optimize prediction accuracy, we show that these methods
perform well, in particular when one allows for differential shrinkage.
Calibration, however, is often bad, as the reduced variability is countered by
increased bias. While this is partly inevitable, we show that local shrinkage
Gelman (2006); Polson and Scott (2012), which penalizes (groups of) variables
differently, may improve calibration considerably. This comes at the cost of
some loss in prediction accuracy as compared to global shrinkage, but local
shrinkage still outperforms ordinary least squares (OLS) by a fair margin,
while being competitive in terms of calibration.

The reported instability of the penalty parameter(s) Riley et al. (2021) may
affect the uncertainty of the predictions. Unlike in settings with cross-
validated or fixed penalty parameter(s), uncertainty propagation of the
penalty parameter(s) is accomplished very naturally in a Bayesian hierarchical
shrinkage model. Hereby, it may aid in communicating appropriate prediction
uncertainties by improving (frequentist) coverage of their credibility
intervals. Moreover, the additional regularization of the penalty parameter(s)
may increase stability of their estimates, in particular in the presence of
multiple penalties. In the frequentist paradigm penalty parameter uncertainty
is not explicitly modeled. Marra and Wood (2012) show, however, that modifying
the classical intervals such that the covariance matrix of the regression
coefficients matches with the Bayesian one (which includes the additional
prior uncertainty of those coefficients) improves coverage. Here, we also
consider those ‘unconditional’ intervals, as implemented in mgcv .

Our main strategy to compare various methods is to divide a large
($N=21,570$), fairly typical epidemiological study into many subsets of sizes
$n=50$ and $n=200$. We study systolic blood pressure (SBP) as outcome, and use
twelve variables to predict SBP. We add five random noise covariates to make
sure that some covariates are not linked to the outcome at all, resulting in
$p=17$. The OLS results on the entire set are used as benchmark, as these
estimates are unbiased and very precise due to the large sample size. Then,
the subsets are analysed with various shrinkage methods and compared to the
benchmark to evaluate prediction accuracy (using out-of-bag test samples),
calibration and coverage of the prediction interval. We also analyse two
simulated data sets: an external one for the logistic setting Van Calster et
al. (2020), which focuses on calibration; and a second one serving as the
following introductory, motivating example.

###  1.1 Introductory example

Our motivating simulation example is loosely inspired by Model A in Riley et
al. (2021), which relates seven covariates to SBP. Treatment (yes/no) is one
of the covariates. The results of OLS in Riley et al. (2021) show a strong
treatment effect as compared to the other six parameters. Such a scenario may
be fairly common in clinical studies. Therefore, we simulate regression
coefficients that adhere to this scenario:
$\beta_{1}=\beta_{2}=\beta_{3}=-0.05,\beta_{4}=\beta_{5}=\beta_{6}=0.05$,
$\beta_{\text{Treat}}=-0.25$ and error variance $\sigma^{2}=1$. We simulate
1,000 data sets of size $n=50$ from the linear model. To each of these, we
apply standard ridge shrinkage, as in Riley et al. (2021), and two-penalty
ridge shrinkage (ridge_2) allowing differential shrinkage for the treatment
parameter. We argue the latter is reasonable, as treatment is an intervention
and therefore a different type of covariate than the others. In the Bayesian
setting (Bay_2) the square-roots of the two inverse penalties (i.e. standard
deviations) are endowed with a standard $\text{C}^{+}(0,1)$ (half-Cauchy)
prior Gelman (2006). Table 1 shows the evaluation results in terms of mean
squared error of the prediction (MSEp) and mean coverage of the 95% confidence
intervals of the predictions, as evaluated on a large test set
($n_{\text{test}}=1,000$).

Method | MSEp, mean | MSEp, 90% | Cover
---|---|---|---
ridge | 0.096 | 0.157 | 0.844
ridge_2 | 0.079 | 0.150 | 0.898
Bay_2 | 0.068 | 0.131 | 0.969
Table 1: Results from simulation example. MSEp: mean squared error of the
prediction; mean and 90% quantile across data sets; Cover: mean coverage of
the 95% confidence intervals of the predictions (target: 0.95)

The results clearly show the benefit of two-penalty shrinkage as compared to
standard ridge: ridge_2 and Bay_2 decrease the mean MSEp by 17.7% and 29.2%,
respectively. Most striking is the difference in interval coverage. While
Bay_2 is somewhat conservative it is much closer to the target than ridge and
ridge_2. We repeated the simulation for a setting that is more beneficial to
global ridge shrinkage: all regression parameters equal $\pm 0.1$. Then,
prediction accuracies of all three methods are very similar, implying the cost
of using one redundant penalty parameter, as for ridge_2 and Bay_2, is
limited. The interval coverage of Bay_2 is still superior to that of ridge
(ridge_2) as it averages to 95.7% versus 86.0% (87.8%).

Below we study various frequentist and Bayesian shrinkage methods in more
detail for real data in the linear regression setting. We evaluate methods by
considering prediction accuracy, coverage of the confidence intervals of the
predictions, and calibration of the predictions. As the latter is particularly
also a concern in binary prediction problems Van Calster et al. (2020) we
address this in an external simulation setting. To enhance reproducibility of
our results we provide i) a synthetic copy of our primary data set, for which
we show it renders qualitatively the same results as the real data; and ii)
all code to run and evaluate the models, including the simulations. The latter
should also assist researchers to use some of the ‘non-standard’ solutions we
propose for their own data. We end with a discussion, which includes
recommendations and several extensions.

###  1.2 Data

The main data we use throughout the manuscript is obtained from the Helius
study Snijder et al. (2017). In a linear regression context, we study response
$Y$: systolic blood pressure (SBP), as a function of covariates $X$: age,
gender, BMI, ethnicity (5 levels), smoking (binary), packyears, coffee
(binary), glucose (log), cholesterol, rendering twelve covariates after dummy
coding the nominal covariate. We apply minimal preprocessing to the data: only
2.7% of the samples had at least one missing value; these samples were
removed, rendering a sample size $N=21,570$. Continuous covariates were
standardized, as this is common practice before applying shrinkage for the
penalty to have the same effect on all corresponding regressing parameters.
For binary covariates, the standardization itself may become unstable for
small data sets, in particular when the classes are unbalanced. This may
hamper generalization to test settings. Therefore, we opted for a default -1,
1 coding, as this standardizes a balanced binary covariate. For our data, we
compared results with those from complete standardization. Differences were
small, but marginally better for the proposed coding. Finally, we added five
independent standard normal noise covariates. Hereby, we are sure to include
some covariates completely unrelated to the outcome, SBP. Hence, the total
number of covariates equals 17.

We chose this data set for various reasons. First, it addresses a fairly
standard, and well-known prediction problem with an interesting mix of
covariates (binary, nominal and continuous). Second, its large sample size
allows to a) use the OLS estimates from the entire set as the benchmark,
because these have very small standard errors (typically $\approx 0.01$); and
b) to split the set in many independent subsets with sample sizes as often
encountered in clinical studies, such as $n=50,200$. This enables us to
evaluate various (shrinkage) methods on many real, relatively small sample
data sets. Third, from applying OLS to the entire data set, we have
$R^{2}=0.34$, which is neither trivially low nor high. The data cannot be
shared, but we provide a synthetic copy which qualitatively renders the same
results as those presented here.

##  2 Methods

All methods fit the linear model
$Y_{i}=\beta_{0}+X_{i}\bm{\beta}=\beta_{0}+\sum_{j=1}^{17}\beta_{j}X_{ij}+\epsilon_{i},\epsilon_{i}\sim
N(0,\sigma^{2})$, with $i=1,\ldots,n$,
$\bm{\beta}=(\beta_{1},\ldots,\beta_{17})$ and covariates
$X_{i}=(X_{i1},\ldots,X_{i17})$. Before discussing the methods, we state the
evaluation criteria, which all focus on prediction.

###  2.1 Evaluation criteria

First, note that the ‘true’ value of $\bm{\beta}$ is obtained from applying
OLS to the entire data set, as this rendered estimates with extremely small
standard errors. Training sets $T^{s}$ and their complementary test sets
$T^{\prime s}$ are indexed by $s$. The true value of the prediction for any
(test) sample $i$ with covariates $X_{i}$ then equals
$\eta_{i}=\beta_{0}+X_{i}\bm{\beta}$. We evaluate the methods on three
criteria:

  1. 1.

_Prediction accuracy_. We use the mean squared error of the predictions
(MSEp). For a given subset $s$ this is defined as

$\text{MSEp}^{s}=1/|T^{\prime s}|\sum_{i\in T^{\prime
s}}(\eta_{i}-\hat{\eta}^{s}_{i})^{2},$ (1)

with, for test sample $i$,
$\hat{\eta}^{s}_{i}=\hat{\beta}_{0}^{s}+X_{i}\hat{\bm{\beta}}^{s}$, where
$\hat{\beta}_{0}^{s}$ and $\bm{\beta}^{s}$ are estimated from samples in
training set $T_{s}$. As we have the ‘true’ predicted values $\eta_{i}$ in
this setting, we use those in (1) instead of $y_{i}$, which is often used when
$\eta_{i}$ is not avalaible, rendering the prediction squared error. Of note:
we checked that the latter leads to the same conclusions as the former.

  2. 2.

_Calibration of the predictions_. As we are in the setting of knowing the true
values, we use the term ‘calibration’ as it is used in measurement technology:
we quantify how close the predictions (which are derived from measurements)
are to an accurate benchmark, in this case the true values. For this, we
regress predictions $\hat{\eta}^{s}_{i}$ on their true values $\eta_{i}$ per
subset $s$, including an intercept as well. The slope of this regression,
which we term ‘cslope’, serves as a calibration metric as it should be close
to 1, and quantifies the possibly biasing effect of shrinkage across the range
of predictions. Note that this deviates from what has become known as the
‘calibration slope’ in the literature, which results from regressing the
observed values on the predictions. The latter provides an assessment of the
overall agreement of observed values and predictions, but has been claimed to
be a misnomer, partly because it is (also) a measure of spread of the
predictions Stevens and Poppe (2020).

  3. 3.

_Uncertainty of the predictions_. We evaluate this by the mean coverage of the
confidence intervals of the predictions
$\bm{\eta}=(\eta_{1},\ldots,\eta_{N})$. In the classical penalized ridge
regression setting we use the 95% Gaussian Wald intervals with standard errors
corrected for prior uncertainty of $\bm{\beta}$, as derived in Marra and Wood
(2012). For Bayesian methods, we simply use the (2.5%, 97.5%) quantiles of the
posteriors to create 95% credible intervals for the predictions $\eta_{i}$.
Note that confidence intervals of predictions $\bm{\eta}$ do not include the
measurement error $\epsilon_{ij}$ as prediction intervals do. For simplicity,
we decided to focus on the former as these more directly relate to the
shrinkage of regression coefficients $\bm{\beta}$, and less so to how the
error variance $\sigma^{2}$ is estimated.

###  2.2 Standard solutions

We study OLS, stepwise selection, lasso and ridge as standard solutions with
no or just one global penalty parameter. The intercept is not penalized. Note
that step is also a shrinkage method, but of an extreme nature: it shrinks
completely to 0 or not at all, hence somewhat similar to the use of a spike-
and-slab prior in a Bayesian model. Here, step is the standard step
implementation in R, using forward-backward selection and AIC to select a
model. We use lasso as implemented in glmnet, using lambda.min resulting from
10-fold cross-validation (CV). For ridge, we estimated penalties using
marginal likelihood optimization as implemented in mgcv Wood (2011), which, as
opposed to glmnet, provides confidence intervals Marra and Wood (2012).
Results from 10-fold CV were very similar, and hence not shown.

###  2.3 Multi-penalty solutions

Multi-penalty solutions allow the use of different penalties for groups of
covariates. We focus on ridge-type solutions for three reasons: a) our focus
lies on prediction and not on variable selection; b) standard software
implementations such as mgcv allow efficient estimation of the penalties; and
c) standard ridge performs relatively well (see Supp Fig 2). In some settings,
the use of covariate groups is very natural, e.g. when some covariates
represent similar entities (e.g. genes), or when main effects plus
interactions are included. In other cases, like ours, choice of the covariate
groups implies a level of subjectivity. Therefore, we also assess the
performance of multi-penalty approaches when the covariate groups are chosen
randomly; that is, when the prior information used to form the covariate
groups is useless. Besides differential penalization another option is not to
shrink one (group of) covariates, for example because there is substantial
evidence that these effect the outcome. We will explore this option as well.

We now define the covariate groups that were used:

  * •

Two groups, $G=2$. The first 3 covariates (age, gender, BMI) are one group,
supposedly because these are known to relate to SBP from previous studies. The
other 14 covariates are in the second group. Two options are explored: ridge_2
penalizes both groups, separately, whereas ridge_2un leaves the first group
unpenalized.

  * •

Three groups, $G=3$, ridge_3. As above, but the nominal covariate, ethnicity,
is a seperate group as this is a nominal covariate with five levels; here,
‘Dutch’ is chosen as the baseline level.

  * •

Two random groups. Three covariates are randomly picked to belong to group 1,
the other 14 belong to group 2. Randomization is repeated for each training
data set. Referred to as ridge_2r and ridge_2unr, the latter corresponding to
leaving the smallest group of covariates unpenalized.

Riley et al. (2021) and others have reported the instability of the penalty
parameters for standard solutions like ridge. This instability is likely to
increase when estimating multiple penalty parameters, as less data per penalty
parameter is available then. Bayesian solutions may be a worthwhile
alternative to shrink these penalties with an additional prior to counter such
instability. Moreover, the extreme case: one penalty per parameter, referred
to as local shrinkage, may benefit calibration, as illustrated further on.

###  2.4 Bayesian solutions

Shrinkage is a very natural concept in Bayesian statistics. Depending on the
type of shrinkage, frequentist and Bayesian solutions may be very similar. For
example, the frequentist ridge and lasso estimates of $\bm{\beta}$ are equal
to the posterior mode estimate of $\bm{\beta}$ when using a Gaussian or
Laplacian prior with a fixed precision parameter that relates proportionally
to the penalty parameter for fixed sample size.

Bayesian methodology, however, can model the variability of the penalty
parameter, which holds several promises. First, it counters potential
instability of estimates of penalty parameter(s) by using a (weakly)
informative prior for it. Second, it propagates uncertainty of penalty
parameter(s). This follows from the fact that for a single $\beta$ and random
penalty parameter $\lambda$ we have
$V(\beta)=E_{\lambda}[V(\beta|\lambda)]+V_{\lambda}[E(\beta|\lambda)].$
Finally, it allows local shrinkage, either by group, or even per covariate.
Note that, just as for OLS, an advantage of local shrinkage is that the
results are less sensitive to (differences in) scale and effect sizes, as the
local penalty can adapt.

Extensions to Bayesian ridge regression differ in how they allow the penalty
parameters to vary, a priori. Table 2 specifies the priors that we used.

Name | Prior $\beta$ | Prior penalty $\lambda$ | Type of shrinkage
---|---|---|---
Bay_EB | $\beta_{k}\sim N(0,\sigma^{2}\lambda^{-1})$ | $\lambda=\hat{\lambda}_{\text{EB}}$ | global, fixed
Bay_IG | $\beta_{k}\sim N(0,\sigma^{2}\lambda^{-1})$ | $\lambda^{-1}\sim\text{IG}(0.001,0.001)$ | global, vague
Bay_glo | $\beta_{k}\sim N(0,\sigma^{2}\lambda^{-1})$ | $\lambda^{-1/2}\sim\text{C}^{+}(0,1)$ | global, weakly informative
Bay_2 (3) | $\beta_{k}\sim N(0,\sigma^{2}\lambda^{-1}_{g(k))}$ | $\lambda^{-1/2}_{g}\sim\text{C}^{+}(0,1)$ | grouped, weakly informative
Bay_loc | $\beta_{k}\sim N(0,\sigma^{2}\lambda^{-1}_{k})$ | $\lambda^{-1/2}_{k}\sim\text{C}^{+}(0,1)$ | local, weakly informative
Table 2: Priors

Here, $g(k)$ denotes the map of covariates to groups:
$\\{1,\ldots,17\\}\rightarrow\\{1,\ldots,G\\}$, with $G=2,3$. EB refers to
‘empirical Bayes’, meaning that the penalty parameter is considered fixed and
estimated by maximizing marginal likelihood. IG is the conjugate inverse-gamma
prior, and $\text{C}^{+}(0,1)$ is the half-Cauchy prior (see Supp Fig 1) for
the standard deviation(s), which has fairly heavy tails and was advocated as a
good default prior by Gelman (2006); Polson and Scott (2012). Moreover,
$\sigma^{2}$, the error variance in the linear model (set to 1 in above’s
prior for the logistic), is endowed with either a Jeffrey’s prior or a vague
inverse gamma (results very similar). It is standard practice to include
$\sigma^{2}$ in the prior of $\beta_{k}$ in Bayesian linear regression. Since
we observed that the global versions of Bay_EB and Bay_IG performed very
similarly to their frequentist counterpart (ridge; see Supp Fig 3), we
primarily focus on the C+ prior for the grouped and local settings.

##  3 Results

###  3.1 Prediction accuracy

We first consider prediction accuracy, as measured by the mean squared error
of the predictions, MSEp (1). In all figures, results are based on 400
training subsets of sizes $n=50,200$, which due to the size of the entire data
set show no or little overlap in samples. Predictions are evaluated on
complementary test sets. Below, we discuss the results for several
comparisons.

Standard solutions: We first compare the standard methods. Suppl. Fig. 2
depicts the performances. For $n=50$, we observe that ridge performs better
than lasso, which is substantially better than stepwise selection and OLS. For
$n=200$ the gap between ridge and lasso becomes much smaller, and so do the
other gaps.

Multi-penalty solutions: Fig 1 shows that use of multiple penalties can
improve prediction accuracy (lower MSEp). For n=50, we observe a substantial
improvement for $G=2$ penalties, likely caused by the large difference in
penalties for the two groups (Supp Fig 4; first and third boxplot). For n=200,
the improvement for $G=2$ is marginal, but the larger sample size benefits the
estimation of one more penalty rendering improved accuracy for $G=3$. In this
setting ridge_2un, which does not penalize the first three covariates,
performs on par with ridge_2. This is reasonable, because these covariates are
relatively strong, so need little penalization. The figure also shows,
however, that when the two groups are assigned at random ridge_2unr performs
inferior to ridge_2r (with suffix ‘r’ denoting the random group setting) as it
cannot adapt when the unpenalized group is relatively weak, and would hence
have benefitted from penalization. Importantly, also observe that even when
the groups are random, ridge_2r performs - on average - (nearly) on par with
ridge, again due to the adaptive penalties. The improvement with the multi-
penalty approach becomes very tangible in Fig. 3. It shows that OLS and ridge
need many more samples to achieve the same median MSEp as ridge_2, f.e.
approximately 135 and 85, respectively, versus 50.

Adding Bayesian solutions: We first compare the standard (empirical) Bayesian
ridge regression solutions, Bay_EB and Bay_IG, with their frequentist
counterpart, ridge. Suppl. Fig. 3 shows very little difference in MSEp between
those methods. This is expected as for both methods the posterior mode
estimate coincides with the frequentist estimate. Small differences are due to
the use of posterior mean instead of posterior mode as a summary for the
predictions.

Fig. 2 compares each of the three Bayesian methods with C+ priors to their
frequentist counterparts. Here, Bay_loc is contrasted with OLS as both methods
do not imply any grouping on the covariates. Likewise, Bay_2 is contrasted
with ridge_2 (2 groups) and Bay_glo is contrasted with ridge (1 group). For
the latter two comparisons we observe that differences in MSEp are minor. A
more striking difference is observed for OLS vs Bay_loc. While the latter
cannot compete with the global shrinkage methods (for $n=50$), it does perform
much better than OLS. Apparently, the default shrinkage of the normal-C+ prior
substantially stabilizes the $\beta$ estimates in such small sample size
settings. The effect of the local normal-C+ prior for $\beta$ is depicted in
Supp. Fig. 5 for two $n=50$ subsets with extreme OLS estimates of one of the
most important coefficients, $\beta_{\text{BMI}}$. We observe that the
normal-C+ prior of Bay_loc shrinks the extreme OLS estimates in the right
direction.

###  3.2 Calibration

Fig. 4 shows the ‘cslopes’ resulting from regressing test sample predictions
$\hat{\eta}^{s}_{i}$ against true $\eta_{i}$ for all subsets $s$ for OLS,
Bay_loc, ridge, Bay_glo, ridge_2 and Bay_2. We observe that OLS and Bay_loc
outperform the other methods on calibration. OLS and Bay_loc are competitive:
the first being unbiased, but with larger variability of the slopes; in fact,
the root MSE of the slope estimates (as compared to 1) are comparable: 0.249
(OLS) vs 0.232 (Bay_loc). The bias of ridge and Bay_glo is evident; shrinkage
compresses the slopes. Differential penalization using two covariate groups
improves results as compared to global regularization, with a small edge for
Bay_glo2 w.r.t. ridge_2 in terms of decreased variability likely due to the
extra regularization of the penalties. Differences in calibration performance
are somewhat smaller for $n=200$ (Supp. Fig. 6).

To elaborate on the inferior calibration of ridge and Bay_glo compared to OLS
and Bay_loc consider the estimated coefficients of a strong covariate (in the
entire study), BMI, from the $n=50$ subsets. Supp. Fig. 7 clearly shows that
ridge and Bay_glo over-shrink $\beta_{\text{BMI}}$ due to the common shrinkage
factor shared with the weaker covariates, whereas Bay_loc strikes a good
balance between a small amount of bias and reducing variability compared to
the OLS estimates.

###  3.3 Quantifying uncertainty: confidence interval of predictions

Uncertainty of the predictions $\bm{\eta}$ results from the variability of
$\bm{\beta}$, which includes two components: one conditional on fixed penalty
parameter(s) $\bm{\lambda}$ and one reflecting the variability of
$\bm{\lambda}$. Riley et al. (2021) show that the latter, the ridge penalties,
may be very variable across small subsets. Indeed, the left-hand side of
Figure 5 confirms this for our data: the penalty estimates may differ 2-3
natural logs in magnitude from one subset to another ($n=50$), with standard
ridge and Bay_glo alike. Hence, the extra regularization of Bay_glo has little
effect here. When multiple penalties are estimated, however, extra
regularization can improve stability of those, as illustrated by Supp. Fig. 4:
Bay_2 compares favourably to ridge_2 for stability of $\lambda_{2}$.
Nevertheless, between subset variability of $\bm{\lambda}$ remains
considerably large. Therefore, we argue that it is important to propagate the
uncertainty of $\bm{\lambda}$ into estimation of $\bm{\beta}$ when analysing
_one_ subset, provided that this reflects the _between_ subset variability, as
the latter is usually not available. The right-hand side of Figure 5 shows
this: when using Bay_glo the posterior variability of $\bm{\lambda}$ for ten
random subsets approximates the between subset variability of the ridge
penalty (left-hand side) fairly well, in particular in order of magnitude. The
penalties are usually not the primary estimands of interest; what matters most
is to what extent their variability propagates towards $\bm{\beta}$ and,
eventually, the predictions $\bm{\eta}$. Therefore, we finalize the comparison
by evaluating coverages of the confidence intervals of $\bm{\eta}$.

As before, we focus on OLS, ridge and ridge_2 and contrast these with Bay_loc,
Bay_glo and Bay_2 in Table 3. Both ridge and ridge_2 are used with uncertainty
propagation of the penalties, as implemented in mgcv’s predict.gam function
with argument unconditional = TRUE. First, while we observe that OLS and
Bay_loc are competitive in terms of coverage, the latter results in narrower,
and hence slightly more useful, intervals. Second, unlike for the introductory
example, ridge and Bay_glo are competitive in terms of coverage, and also in
terms of mean width of the intervals. Finally, for the grouped penalties the
extra shrinkage imposed by Bay_2 pays off as it maintains good coverage,
unlike ridge_2, which seems to render too narrow intervals. Note that in the
latter case, the 10% quantile of the coverage of the predictions (as computed
from the test individuals) drops to only 0.772 for ridge_2, whereas it equals
0.925 for Bay_2. Note that the mean coverages of ridge and ridge_2 drop
slightly (from 0.944 to 0.920 and from 0.896 to 0.864, respectively) when not
adjusting the standard errors (unconditional = FALSE in predict.gam). For
$n=200$, results are qualitatively similar, although the intervals of OLS are
now somewhat narrower than those of Bay_loc (Supp Table 4).

| Coverage | Width
---|---|---
Methods | Classical | Bayes | Classical | Bayes
OLS, Bay_loc | 0.939 | 0.971 | 2.334 | 2.149
ridge, Bay_glo | 0.944 | 0.961 | 1.521 | 1.607
ridge_2, Bay_2 | 0.896 | 0.966 | 1.216 | 1.470
Table 3: Mean coverage (target: 0.95) and width of confidence intervals of
predictions for $n=50$, contrasting classical method with Bayesian counterpart

###  3.4 Conclusions, linear case

Below we list the conclusions from the linear regression case.

  1. 1.

Shrinkage benefits prediction accuracy compared to OLS and stepwise selection.

  2. 2.

Differential shrinkage for groups of covariates may substantially benefit
prediction accuracy compared to global shrinkage. Therefore, the same
prediction accuracy can be achieved with a smaller sample size.

  3. 3.

Penalty estimates vary substantially between small subsets of samples. Hence,
it is relevant to propagate this variability when quantifying uncertainty of
the predictions.

  4. 4.

Shrinkage may indeed lead to some level of miscalibration, although minimally
for Bayesian shrinkage with a local C+ prior. The latter is competitive to OLS
on this matter and on coverage, while rendering better prediction accuracy.
Prediction accuracy is worse, though, than that of the global regularization
methods.

  5. 5.

Bayesian shrinkage with a grouped C+ prior outperforms its frequentist
counterparts in terms of interval estimation. In case of global shrinkage, the
two are very competitive, with the latter slightly less conservative.

##  4 Calibration for logistic regression: external simulation

Here, we study the small sample logistic regression setting. This differs from
the linear one in three essential ways: a) binary outcomes are much less
information-rich than continuous ones; b) the logistic function flattens the
impact of large $\beta$’s; and c) the maximum likelihood estimator (MLE) is
biased Firth (1993). We focus on prediction accuracy and calibration, as
standard shrinkage methods like lasso and ridge with cross-validated penalties
generally improve the first, but may deteriorate the latter. This is intrinsic
due to the introduced bias, but may be worsened by the conventional tuning of
the penalty parameter: minimization of cross-validated prediction error, which
may be far from optimal for calibration. Hence, we discuss some alternatives.
Firth’s solution to adjust bias Firth (1993) uses a fixed penalty, defined on
the level of the information matrix, and was shown to perform favourably to
the tuning methods Van Calster et al. (2020). Sullivan and Greenland (2013)
propose a similar, alternative solution, which also accounts for the
aforementioned issues a) and b): simply fix the prior variance of the
$\beta$’s to 0.5, equivalent to setting the ridge penalty $\lambda=1/0.5=2$.
Their argument is that for the small sample logistic setting one has to ‘dare’
to be fairly informative, as the outcomes are information-poor. Based on this,
we propose a somewhat more objective alternative that does not require any
tuning either: use Bay_glo or Bay_loc as defined in Table 2, but with a
$C^{+}(0,\sqrt{0.5})$ instead of a $C^{+}(0,1)$ prior for the standard
deviation, as the former matches to a median prior variance of 0.5. These
priors are depicted in Supp. Fig. 1. So, we allow for some variation in the
prior variances, thereby enabling adaptation to cases in which (some of the)
$\beta$’s are more (or less) extreme.

For evaluation we follow one of the simulation set-ups by Van Calster et al.
(2020): $n=50,\bm{\beta}=(\beta_{1},\ldots,\beta_{5})=(0.2,0.2,0.2,0.5,0.8)$
and intercept $\beta_{0}$ tuned such that the average event probability equals
0.5, rendering a challenging setting of 25/5 = 5 events per covariate.
Covariates $X$ are simulated from a multivariate normal with means 0,
variances 1, and correlations 0.5. Response $Y$ is then simulated from a
Bernoulli with success probabilities $1/(1+\exp(-\beta_{0}-X\bm{\beta})).$ We
refer to this setting as ‘Moderate signal’. Additionally, to study how results
adapt to the signal strength, we simulate the ‘Weak and Strong signal’
settings by using $\bm{\beta}_{\text{weak}}=\bm{\beta}/3$, and
$\bm{\beta}_{\text{strong}}=3\bm{\beta}$. Calibration is assessed by the
‘cslope’, which results from regressing the estimated linear predictors
against the true ones (using the real $\bm{\beta}$) for a large test set. A
slope close to 1 indicates good calibration. Note that Van Calster et al.
(2020) report the classical calibration slope, which results from regressing
observations on predictions. The latter measures in fact a mix of calibration
and spread of the predictions Stevens and Poppe (2020), which is why we opt
for the cslope. In our case, too strong shrinkage will render cslope $\ll 1$.
The MSE of the predictions (MSEp) is assessed by the average squared error of
those predicted probabilities w.r.t. the true ones. Figure 6 displays the
results based on 50 training sets for ML (maximum likelihood); Firth; ridgeCV:
ridge with cross-validated penalty; ridge05: ridge with fixed penalty 1/0.5;
Bay_glo05 and Bay_loc05: global and local regularization with
$C^{+}(0,\sqrt{0.5})$ priors for the standard deviation(s). Below we list the
conclusions.

  * •

We confirm that ridgeCV performs well in terms of prediction accuracy, but
poor on calibration, with clopes $\ll 1$. Firth calibrates better than both ML
and ridgeCV.

  * •

The solution by Sullivan and Greenland (2013), ridge05, is competitive to
Firth in terms of calibration, when the signal is weak or moderate, but
overshrinking when the signal is strong. It is generally superior to ML and
Firth for prediction.

  * •

Bay_glo05 is generally fairly competitive to ridgeCV in terms of prediction
accuracy. On average it calibrates better, but shows more variability. It
calibrates less well than Firth, and Bay_loc05.

  * •

Bay_loc05 is competitive to Firth in terms of calibration, but is superior for
prediction in particular in the weak signal setting.

  * •

Bay_loc05 is competitive to ridge05 in terms of calibration. It adapts better
to the change in signal at the price of somewhat more variation. It is
competitive in terms of prediction accuracy, taking into account that
variation towards low MSEp is actually desirable.

Supp Fig 9 and 10 present the results of two alternative scenarios (larger
sample size: $n=100$, and adding five $\beta$’s equalling 0). Results are
qualitatively similar to those presented here. For completeness, we have also
co-plotted the (inverted) cslope and the classical calibration slope for the
original simulation scenario in Supp Fig 11. These two agree fairly well for
the strong signal setting, but less so for the weak and moderate signal
settings, as the classical calibration slopes are more affected by the
variability of the predictions.

All-in-all, we conclude that for calibration Firth, ridge05, and Bay_loc05 are
competitive to one another and superior to ML and the global penalization
methods with adaptive penalty, ridgeCV and Bay_glo05. The latter two, though,
are superior in terms of prediction accuracy when the signal is weak. In that
case Bay_loc05 outperforms Firth in terms of prediction accuracy. Unlike
ridge05, Bay_loc05 adapts to the signal strength, but shows somewhat more
variability due to the non-fixed penalty. Note that an important advantage of
Bay_loc05 is that unlike Firth and ridge05, it straightforwardly provides
uncertainty estimation of the predictions, including the uncertainties of the
penalties.

##  5 Software and reproducibility

###  5.1 Software

Below, we list the software packages used for the various models.

  1. 1.

OLS and stepwise selection: R’s lm and step functions.

  2. 2.

Linear ridge, global and grouped penalty: gam function in R-package mgcv Wood
(2011). Unlike glmnet this allows for uncertainty computations.

  3. 3.

Logistic ridge: glmnet, as this was also used by Van Calster et al. (2020).
Firth’s correction: logistf Heinze and Schemper (2002).

  4. 4.

Lasso: glmnet Friedman et al. (2010).

  5. 5.

Bayes, linear: shrinkage software Leday (2022), which is optimized for
computational efficiency, hence convenient for the repeated calculations on
data splits. Results highly agree with those of r-stan.

  6. 6.

Bayes, logistic: r-stan Stan Development Team (2022).

###  5.2 Data sharing

As is the case for many cohort studies, our primary data source, the Helius
data, cannot be shared publicly. This is impeding methodological studies like
this, because such large $N$ studies are very useful to check results of
various methods on small subsets as we did. Therefore, we share a synthetic
copy of our data set which may be used to a) qualitatively reproduce our
results; or b) evaluate other methods than those proposed here. We verified
whether results from the synthetic copy agreed with those from the real data,
which is the case (see Supp Fig 8 for MSEp). The imputation-based method to
generate the synthetic data is described in the Supp Mat.

###  5.3 Do it yourself

One may certainly argue that results and conclusions might differ somewhat for
other data than those used here. The following scheme may be useful to assess
what shrinkage method works best for the data at hand. Prediction accuracy can
be assessed by cross-validation principles. To evaluate calibration and
coverage, one needs to know ‘real’ predictions in a setting that mimics one’s
own. For that, run an OLS or, preferrably, the more stable Bay_loc on the data
and obtain $\hat{\bm{\beta}}$ as well as an estimate of the noise. Then
simulate response $y$ many times from the regression model using the real
design matrix $X$ and those estimated parameters. Then, our scripts (or other
software for the proposed models) can be used to assess the performances of
the methods in one’s own data-centric simulation setting.

###  5.4 Availability

Annotated software scripts and R-markdown files to reproduce our results are
available from github: https://github.com/markvdwiel/ThinkBeforeShrink. This
repository also contains a synthetic copy of the Helius data and an example
script to guide users.

##  6 Discussion

We observed that grouping of covariates to allow for differential penalties
may improve predictive performance. Here, we only studied the low-dimensional
linear regression setting, but we have observed similar results for the high-
dimensional logistic ridge setting Van de Wiel et al. (2016); van Nee et al.
(2021). One could argue that the grouping of covariates is subjective.
However, we demonstrated that the results are robust against misspecification
of the groups. Moreover, it is much less subjective than a commonly used
strategy, as stated in Van Calster et al. (2020): “Alternatively, a less
complicated model can be considered, for example by discarding many predictors
a priori.”. The latter strategy basically assigns an infinite penalty to the
excluded covariates a priori. We believe it is better to leave those in, but
allow differential shrinkage for those as compared to the others.

If calibration is key, global shrinkage methods with flexible penalties, both
frequentist and Bayesian ones, are usually not very suitable. This is not
surprising as these penalty parameter(s) adapt to predict optimally in terms
of accuracy, which is a different goal than calibration. We showed that
Bayesian local shrinkage is a good alternative which competes with
conventional methods that do not tune penalties, such as OLS, ML, Firth and
ridge with a fixed penalty. Note that the importance of having calibrated
predictions depends on a) the possibility of recalibrating predictions in a
later stage of the study; and b) the need or desire to interpret the
predictions on an absolute scale. For the latter, the design of the study is
also very relevant; often recalibration will be needed anyhow when applying
the predictor to a population with somewhat different characteristics than the
one represented by the study, a very common case for medical studies.

While the grouping of covariates can aid in improving prediction accuracy, it
may deteriorate coverage of the confidence intervals of predictions in the
classical setting. Here, the additional regularization by the $C^{+}(0,1)$
prior as invoked by the Bayesian procedure helps to improve coverage, while
maintaining competitiveness on prediction accuracy and calibration. The
importance of quantifying uncertainty of the predictions differs from study to
study. However, in any study it will be useful to know it, also to assess
whether the sample size should be increased to lower this uncertainty to an
acceptable level.

We discussed the linear and logistic case, which could be regarded as two
extremes: the former being information-rich with unbiased ML (= OLS)
estimates, the latter being information-poor with biased ML estimates. This
means that in the latter case mild shrinkage is also beneficial for
calibration. Moreover, it may be relatively more beneficial to a use a
subjective prior (or penalty) in the logistic setting. We discussed two
examples that performed well, prior variance equal to 1/2 Sullivan and
Greenland (2013), or less subjective, a hyperprior for the variance with
median equal to 1/2. A good alternative may be to use a historical prior,
which is tuned to available data from similar studies Neuenschwander et al.
(2010); MacLehose and Hamra (2014). Then, the solution with a hyperprior may
be preferable over one that fixes the variance, as the former allows the
current study to deviate somewhat more when it would not behave similarly.

For our applications, computational time was not an issue at all: all methods
ran within seconds for given data sets. Nevertheless, the classical methods
fitted substantially faster than the Bayesian ones, rendering the former an
edge for analysing larger scale studies. Whether this balances against the
demonstrated benefit of higher-level shrinkage that Bayesian methods can
provide, will depend on the main aim(s) of the study (prediction accuracy /
calibration / uncertainty quantification). In general, Bayesian methods gain
popularity in epidemiological research MacLehose and Hamra (2014), which
assists their acceptance for general use. Note that, on purpose, we did not
tweak the Bayesian methods in terms of convergence checks or use of other
hyper-priors, as these would imply further tuning, rendering the comparison
with less flexible methods unfair. In terms of software, we recommend to use
either mgcv, or R-stan with the $C^{+}$ priors, as both incorporate estimation
of multiple penalty parameters and provide intervals with, in most settings,
fairly good coverage. The latter has an edge when calibration is key as it
allows for local penalties, or when stabilization of grouped penalties is
relevant, e.g. for uncertainty quantification. Supp Fig 12 shows a flowchart
with our recommendations on when to use which method.

Our work is limited in scope, so we discuss several important extensions. The
first one is variable selection. As this is a different goal than prediction,
other priors and penalties should be discussed. In the classical setting many
variants of the lasso (adaptive, group, hierarchical) become relevant,
compared to more traditional stepwise selection techniques. In the Bayesian
setting, one may wish to include spike-and-slab priors and/or Zellner’s
g-prior, as these are targeting for variable selection George and Foster
(2000). Alternatively, posterior selection techniques that apply to the fairly
dense ridge-type methods used here may be very competitive to those more
sparse formulations Bondell and Reich (2012).

A second extension is to study the effect of shrinkage in the context of
causal inference. Many causal inference frameworks make use of prediction
methods to account for confounding or non-random treatment allocation. Unlike
in predictive modelling, the emphasis in causal inference is primarily placed
on bias, and much less so on variance reduction. This makes the use of
shrinkage less natural, in particular for the treatment effect. Nevertheless,
some (local) shrinkage of the other covariates may be beneficial when $p$ is
relatively large compared to $n$, possibly in combination with double-robust
estimation Avagyan and Vansteelandt (2021) to counter misspecification of the
model.

A third extension is the multi-regression setting, as often encountered in
high-dimensional multiple testing setting, e.g. when relating gene expression
to phenotype, correcting for confounders like age and gender. In such a
setting, shrinking effects across similar features (e.g. genes) may improve
effect size estimates and multiple testing properties Van de Wiel et al.
(2012). Moreover, a multi-regression setting allows tuning of hyper-parameters
of the penalty’s prior across features using empirical Bayes Leday et al.
(2017).

Shrinkage is only a partial solution for underpowered studies. In the end,
increase of sample size will often be needed to draw firm conclusions Riley et
al. (2020). However, we do believe shrinkage may play an important role in the
research cycle, which for many practical reasons often starts with a fairly
small study. Then, as demonstrated, well-thought global or grouped shrinkage
methods can help to better assess the predictive potential of the study and
quantify uncertainty of the predictions, whereas local shrinkage methods are
able to reduce bias and improve calibration. Therefore, these shrinkage
methods aid in deciding whether and how to extend the study.

## References

  * Avagyan and Vansteelandt (2021) Avagyan, V. and Vansteelandt, S. (2021).  High-dimensional inference for the average treatment effect under model misspecification using penalized bias-reduced double-robust estimation.  _Biostatistics & Epidemiology_, pages 1–18.
  * Bondell and Reich (2012) Bondell, H. and Reich, B. (2012).  Consistent high-dimensional bayesian variable selection via penalized credible regions.  _J Amer Statist Assoc_ , 107, 1610–1624.
  * Breiman (2001) Breiman, L. (2001).  Random Forests.  _Machine Learning_ , 45, 5–32.  doi:10.1023/A:1010933404324.
URL http://link.springer.com/10.1023/A:1010933404324

  * Buuren (2018) Buuren, S.v. (2018).  _Flexible imputation of missing data_.  Chapman and Hall/CRC interdisciplinary statistics series. CRC Press, Taylor & Francis Group, Boca Raton, second edition.
  * Buuren and Groothuis-Oudshoorn (2011) Buuren, S.v. and Groothuis-Oudshoorn, K. (2011).  mice : Multivariate Imputation by Chained Equations in R.  _Journal of Statistical Software_ , 45, 1–67.  doi:10.18637/jss.v045.i03.
  * Firth (1993) Firth, D. (1993).  Bias reduction of maximum likelihood estimates.  _Biometrika_ , 80, 27–38.
  * Friedman et al. (2010) Friedman, J. et al. (2010).  Regularization paths for generalized linear models via coordinate descent.  _J Statist Soft_ , 33, 1.
  * Gelman (2006) Gelman, A. (2006).  Prior distributions for variance parameters in hierarchical models.  _Bayesian analysis_ , 1, 515–534.
  * George and Foster (2000) George, E.I. and Foster, D.P. (2000).  Calibration and empirical Bayes variable selection.  _Biometrika_ , 87, 731–747.
  * Heinze and Schemper (2002) Heinze, G. and Schemper, M. (2002).  A solution to the problem of separation in logistic regression.  _Statist Med_ , 21, 2409–2419.
  * Leday et al. (2017) Leday, G.G.R. et al. (2017).  Gene network reconstruction using global-local shrinkage priors.  _Ann Appl Statist_ , 11, 41–68.
  * Leday (2022) Leday, G. (2022).  shrinkage: An R package to fit Bayesian regression models with shrinkage priors.  https://github.com/gleday/shrinkage.  R package version 1.2.
URL https://github.com/gleday/shrinkage

  * Liaw and Wiener (2002) Liaw, A. and Wiener, M. (2002).  Classification and Regression by randomForest.  _R News_ , 2, 18–22.
URL https://CRAN.R-project.org/doc/Rnews/

  * MacLehose and Hamra (2014) MacLehose, R.F. and Hamra, G.B. (2014).  Applications of Bayesian methods to epidemiologic research.  _Curr Epidemiol Rep_ , 1, 103–109.
  * Marra and Wood (2012) Marra, G. and Wood, S.N. (2012).  Coverage properties of confidence intervals for generalized additive model components.  _Scand J Statist_ , 39, 53–74.
  * Neuenschwander et al. (2010) Neuenschwander, B. et al. (2010).  Summarizing historical information on controls in clinical trials.  _Clin Trials_ , 7, 5–18.
  * Polson and Scott (2012) Polson, N.G. and Scott, J.G. (2012).  On the half-cauchy prior for a global scale parameter.  _Bayesian Analysis_ , 7, 887–902.
  * Riley et al. (2020) Riley, R.D. et al. (2020).  Calculating the sample size required for developing a clinical prediction model.  _BMJ_ , 368.
  * Riley et al. (2021) Riley, R.D. et al. (2021).  Penalization and shrinkage methods produced unreliable clinical prediction models especially when sample size was small.  _J Clin Epidemiol_ , 132, 88–96.
  * Snijder et al. (2017) Snijder, M.B. et al. (2017).  Cohort profile: the healthy life in an urban setting (HELIUS) study in Amsterdam, the Netherlands.  _BMJ open_ , 7, e017873.
  * Stan Development Team (2022) Stan Development Team (2022).  RStan: the R interface to Stan.  https://mc-stan.org/.  R package version 2.21.7.
URL https://mc-stan.org/

  * Stevens and Poppe (2020) Stevens, R.J. and Poppe, K.K. (2020).  Validation of clinical prediction models: what does the “calibration slope” really measure?  _J Clin Epidemiol_ , 118, 93–99.
  * Stronks et al. (2013) Stronks, K. et al. (2013).  Unravelling the impact of ethnicity on health in Europe: the HELIUS study.  _BMC Public Health_ , 13, 402.  doi:10.1186/1471-2458-13-402.
URL http://bmcpublichealth.biomedcentral.com/articles/10.1186/1471-2458-13-402

  * Sullivan and Greenland (2013) Sullivan, S.G. and Greenland, S. (2013).  Bayesian regression in SAS software.  _Int J Epidemiol_ , 42, 308–317.
  * Van Calster et al. (2020) Van Calster, B. et al. (2020).  Regression shrinkage methods for clinical prediction models do not guarantee improved performance: simulation study.  _Statist Meth Med Res_ , 29, 3166–3178.
  * Van de Wiel et al. (2012) Van de Wiel, M.A. et al. (2012).  Bayesian analysis of RNA sequencing data by estimating multiple shrinkage priors.  _Biostatistics_ , 14, 113–128.
  * Van de Wiel et al. (2016) Van de Wiel, M.A. et al. (2016).  Better prediction by use of co-data: adaptive group-regularized ridge regression.  _Statist Med_ , 35, 368–381.
  * van Nee et al. (2021) van Nee, M.M. et al. (2021).  Flexible co-data learning for high-dimensional prediction.  _Statist Med_ , 40, 5910–5925.
  * Volker and Vink (2021) Volker, T.B. and Vink, G. (2021).  Anonymiced Shareable Data: Using mice to Create and Analyze Multiply Imputed Synthetic Datasets.  _Psych_ , 3, 703–716.  doi:10.3390/psych3040045.
URL https://www.mdpi.com/2624-8611/3/4/45

  * White et al. (2011) White, I.R. et al. (2011).  Multiple imputation using chained equations: Issues and guidance for practice.  _Statistics in Medicine_ , 30, 377–399.  doi:10.1002/sim.4067.
URL http://doi.wiley.com/10.1002/sim.4067

  * Wood (2011) Wood, S.N. (2011).  Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models.  _J Roy Statist Soc B_ , 73, 3–36.

##  7 Figures

Figure 1: MSEp (y-axis) for grouped penalization across 400 subsets. Digital
suffix denotes number of covariate groups, suffix ‘un’ denotes one unpenalized
group, suffix ‘r’ denotes random groups Figure 2: MSEp (y-axis) for classical
methods and their Bayesian counterpart across 400 subsets. Digital suffix
denotes number of covariate groups Figure 3: Median MSEp (y-axis) across 400
subsets (y-axis) for various sample sizes (x-axis) for OLS, ridge and 2-group
ridge Figure 4: Calibration: cslopes (y-axis) across all 400 subsets of size
$n=50$ Figure 5: Variability of ridge penalties. Y-axis: ridge penalties on
natural log-scale. Left panel (black): ridge penalties estimated by mgcv
(ridge) and shrinkage (Bay_glo) from _all_ 400 subsets: 2.5%, 50% (dot), 97.5%
quantiles. Right panel (grey): posterior estimates of ridge penalty for 10
random subsets using Bay_glo Figure 6: Calibration (top; cslope: winsorized at
$(1/3,3)$, log-scale) and prediction accuracy (bottom; MSEp) for logistic
regression simulations

##  8 Acknowledgements

The Amsterdam University Medical Centres and the Public Health Service of
Amsterdam (GGD Amsterdam) provided core financial support for HELIUS. The
HELIUS study is also funded by research grants of the Dutch Heart Foundation
(Hartstichting; grant no. 2010T084), the Netherlands Organization for Health
Research and Development (ZonMw; grant no. 200500003), the European
Integration Fund (EIF; grant no. 2013EIF013) and the European Union (Seventh
Framework Programme, FP-7; grant no. 278901). Moreover, we thank Hans Berkhof
for discussions on several aspects of this manuscript.

##  9 Supplement, incl. Figures

###  9.1 Synthetic data

Data from the Helius study Stronks et al. (2013) were used as a running
example throughout the manuscript, but they cannot be shared due to privacy
regulations. To enhance reproducibility of our results, we provide a synthetic
copy of the Helius study data to render qualitatively similar results as
compared to the real data. This supplement describes construction of the
synthetic data set.

In essence, the problem was approached from a missing data perspective. After
applying standardization and the removal of rows with missing data (as
described in the main manuscript), the original data set of size $n\times m$
was augmented with an empty data set of size $n\times m$. Subsequently, the
empty rows were imputed using imputation by chained equations White et al.
(2011); Buuren (2018) as implemented in the mice package Buuren and Groothuis-
Oudshoorn (2011) in R. This procedure makes use of univariate imputation
models for every incomplete variable, and uses an iterative updating scheme
that iterates between imputation and updating of the imputations models. The
specific application of (multiple) imputation by chained equations for the
creation of synthetic data was described by Volker and Vink (2021).

For this particular application, the variables systolic blood pressure, age,
gender, BMI, ethnicity (5 levels), smoking (binary), packyears, coffee
(binary), glucose (log), and cholesterol, had to be imputed. Predictive mean
matching Buuren (2018) was used for all continuous variables, logistic
regression imputation was used for binary variables, and a multinomial
imputation model was used to impute ethnicity. All imputation models were
linear additive models conditional on all of the remaining variables. For
example, the imputation models for age was a linear additive model with main
effects of systolic blood pressure, age, gender, BMI, ethnicity (dummy coded,
4 df), smoking (binary), packyears, coffee (binary), log glucose, and
cholesterol level. A single imputed data set was created to represent a
synthetic version of the helius data. Traceplots were stable after the chosen
number of 25 iterations for the mice algorithm. The synthetic data set was
found to closely resemble the original data with respect to mean structure,
variance, and covariance structure, with a mean absolute deviation from the
original data $<0.01$ on the standardized data for all three. The synthetic
data set did not contain any duplicated rows from the original data. As a
further check, a random forest was fitted with the aim to distinguish the
original and synthetic data Breiman (2001); Liaw and Wiener (2002). Half of
the data was used for training, and half for testing, and 500 trees were
fitted. The test Brier score was 0.23 (with random providing 0.25), indicating
that the forest could not accurately separate real and synthetic data.

It is worth noting that the linear additive nature of the imputation models
suited our purposes well, but that it does not provide a general solution for
synthetic data generations, because it does not reflect possible non-linear
relations and interactions. Other possibilities are available, such as
multiple imputation using classification regression trees Volker and Vink
(2021). Also, it should be noted that use of a single imputed data set does
not reflect the uncertainty of the imputed (synthetic) values. For synthetic
data to be used to replicate inference results, where this uncertainty is key,
multiple imputed/synthetic sets would be required Volker and Vink (2021).

###  9.2 Figures

Figure 1: Half-Cauchy priors: $C^{+}(0,1)$ (solid), $C^{+}(0,\sqrt{0.5})$ (dashed) Figure 2: MSEp (y-axis) for several standard methods across 400 subsets Figure 3: MSEp (y-axis) for standard ridge variants across 400 subsets Figure 4: Variability of multi-group (2) ridge penalties. Y-axis: ridge penalties on natural log-scale. Ridge penalties (lam1, lam2) estimated by mgcv (ridge_2) and shrinkage (Bay_2) from all 400 subsets: 2.5%, 50% (dot), 97.5% quantiles Figure 5: Normal-$C^{+}(0,1)$ prior and estimates of $\beta_{\text{BMI}}$ for two subsets of size $n=50$. Vertical lines: true value of $\beta_{\text{BMI}}$ (solid), OLS estimates (dotted) and Bay_loc estimates (dashed). Arrows indicate direction of shrinkage by Bay_loc. Figure 6: Calibration: cslopes across all 400 subsets of size $n=200$. Figure 7: Estimates of $\beta_{\text{BMI}}$ for subsets of size $n=50$. Horizontal line: true value of $\beta_{\text{BMI}}$.  | Coverage | Width
---|---|---
Methods | Classical | Bayes | Classical | Bayes
OLS, Bay_loc | 0.949 | 0.984 | 0.980 | 1.186
ridge, Bay_glo | 0.953 | 0.973 | 0.864 | 1.000
ridge_2, Bay_2 | 0.905 | 0.952 | 0.746 | 0.870
Table 4: Mean coverage (target: 0.95) and width of confidence intervals of
predictions for $n=200$ Figure 8: MSE of the predictions for real (left) and
synthetic (right) data. Upper row: $n=50$, bottom row: $n=200$ Figure 9:
Calibration (top; cslope: winsorized at $(1/3,3)$, log-scale) and prediction
accuracy (bottom; MSEp) for logistic regression simulations. Simulation
setting as described, but with sample size $n$ doubled to 100 Figure 10:
Calibration (top; cslope: winsorized at $(1/3,3)$, log-scale) and prediction
accuracy (bottom; MSEp) for logistic regression simulations. Simulation
setting as described, but with five $\beta_{j}$’s added for which
$\beta_{j}=0$ Figure 11: Inverted cslope (top) and calibration slope (bottom)
for logistic regression simulations; both winsorized at $(1/10,10)$ and on
log-scale. Simulation setting as in main document. As the predictions are used
as dependent variables for the cslope and as independent ones for the
calibration slope, the inverted cslope is displayed Figure 12: Flow chart for
recommended use of various methods, depending on primary goal (Calibration /
Accuracy) and secondary goal (+ Accuracy / + Coverage).
∗Coverage and effect of covariate groups have not been studied for the
logistic setting, so these are not included in this flow chart

*


# Sloshing in coffee as a pumped pendulum

Nicolás Guarín Zapata
<EMAIL_ADDRESS>
###### Abstract

We present a study of the dynamics of a cup of coffee while walking with it in
terms of a lumped model. The model considered is a planar (nonlinear) pendulum
with parametric and direct excitation. The position of the pivot is changed
with time, which leads to parametric excitation in the vertical direction; and
parametric/direct excitation in the horizontal one. For the former case, the
Method of Averaging is used to determine the regions of stability for two
steady-state solutions. In the vertical/horizontal excitation, we determined
the resonances of the system using the Method of Multiple Scales and show the
results computed with _exact_ numerical integration.

## Introduction

Spilling coffee while walking is a common problem. The sloshing dynamics of
the coffee, while walking, can be understood in terms of an oscillator with
parametric excitation caused by the motion of the cup along the path. The main
objective of this project was to study the slosh dynamics of a cup of coffee
(while walking) in terms of a lumped mechanical model [8].

According to [3], the motion of a free-liquid-surface has three regimes
(linear, weakly nonlinear, and strongly nonlinear). Depending on the regime
that we want to study the system, we can use a different —lumped— mechanical
model to understand its dynamics. A common way to model the sloshing of
liquids in a cylindrical container is to consider it as a pendulum, and obtain
the parameters (length, mass, dissipation coefficient) from the original
system [3, 5, 8]. The model considered is a pendulum with a (parametric)
excitation that comes from the _human walking_ , this excitation is composed
of two degrees of freedom: one in the plane (back-and-forth and lateral) and
one vertical. Due to the main frequencies in the _human walking_ , the model
will consist of a single pendulum111In [3], the author use several coupled
pendulums to take into account different modes of vibration in the model..

We present a study of the dynamics of a cup of coffee while walking with it in
terms of a lumped model. The model considered is a planar (nonlinear) pendulum
with parametric and direct excitation. The position of the pivot is changed
with time, which leads to parametric excitation in the vertical direction; and
parametric/direct excitation in the horizontal one. For the former case, the
Method of Averaging is used to determine the regions of stability for two
steady-state solutions. In the vertical/horizontal excitation, we determined
the resonances of the system using the Method of Multiple Scales and show the
results computed with _exact_ numerical integration. An interesting outcome of
an analytical solution is the capability of relating the parameters of the
system to obtain some insight into the design process.222Something that can be
useful in the development of devices like the one described in this patent
[11]..

##  1 Modeling

A cup of coffee is, essentially, a cylindrical-like container filled with
liquid —that can be assumed to behave like water. And this liquid will present
some motion when a person walks with a cup of coffee at hand. The behavior of
this system (for linear and weakly nonlinear regimes) can be interpreted in
terms of the natural frequencies of the fluid. The natural frequencies of
oscillation of a frictionless, vorticity-free, and incompressible liquid in a
cylindrical container (a cup of coffee in this case) with a free liquid
surface are given by [3]

$\omega_{nm}^{2}=\frac{g\epsilon_{mn}}{R}\tanh\left(\epsilon_{mn}\frac{H}{R}\right)\left[1+\frac{\sigma}{\rho
g}\left(\frac{\epsilon_{mn}}{R}\right)^{2}\right]\enspace,$

where $m=0,1,2,\cdots$ and $n=1,2,\cdots$.333In our case the surface tension
is negligible and the only relevant parameters are the geometrical ones. In
this equation $H$ is the height of the cup, $R$ its radius, $g$ the gravity,
$\rho$ fluid density, and $\sigma$ is the surface tension. $\epsilon_{mn}$ are
the roots of the first derivative of the $m$-th order Bessel function —
$J^{\prime}_{m}(\epsilon)=0$.

One can match each of the modes of vibration of the liquid with a single
degree of freedom oscillator, e.g, a mass-spring system or a pendulum.
According to Ibrahim in [3]:

> A realistic representation of the liquid dynamics inside closed containers
> can be approximated by an equivalent mechanical system… For linear planar
> liquid otion, one can develop equivalent mechanical models in the form of a
> series of mass-spring dashpot systems or a set of simple pendulums. For
> nonlinear sloshing phenomena, other equivalent models such as spherical or
> compound pendulum may be developed to emulate rotational and chaotic
> sloshing.

Since the sloshing dynamics of the coffee (while walking) can be understood in
terms of an oscillator with parametric excitation caused by the motion of the
cup along the path [8]. A first approach is to consider the system as a planar
pendulum with parametric pumping. Taking into account the parameters that we
use parametric excitation for this model can come from: the change of length,
the motion of the point of gyration, the motion of the center of mass444This
is the case of a playground swing, where the center of mass oscillates around
a given position. When the oscillation is symmetric it can be proved to be
equivalent to the varying length pendulum [1]., and the change of
mass555Although this is less common..

The system can be considered as a pendulum with varying length, which is a
system that has been studied before [1, 12], and with moving —both, vertically
and horizontally— pivot. The parameters of the system are: $r_{0}$ original
length of the rod, $r$ length of the rod, $m$ mass of the bob, $g$
gravitational acceleration, $x_{0}$ is a function of time that describes the
horizontal position of the pivot, and $z_{0}$ is a function of time that
describes the vertical position of the pivot. The present work restricts the
excitation of the system to the motion of the pivot, which leads to _direct_
and _parametric_ excitations. 666Appendix A shows the parameters of the lumped
model in terms of the original _cup of coffee_.

Figure 1: Considered system: A pendulum with the point of gyration moving in
the plane. This motion cause direct and parametric excitation in the system.

The Lagrangian of the system is given by [4]

$L=T-V=\frac{m}{2}[(\dot{x}-\dot{x}_{0})^{2}+(\dot{z}+\dot{z}_{0})^{2}]+mg[z+z_{0}]$
(1)

taking

$x=r(t)\sin\theta,\quad z=-r(t)\cos\theta$

and using the Euler-Lagrange equation we get the differential equation, we get

$r^{2}\ddot{\theta}+r[g+\ddot{z}_{0}]\sin\theta+r\ddot{x}_{0}\cos\theta+2r\dot{r}\dot{\theta}=0\enspace.$
(2)

In the general case $x_{0}$, $z_{0}$ and $r$ are functions of time, and we can
talk about having _modulated stiffness_ , _damping_ , and _inertia_. Where the
$r$ accompanies the second derivative of the angle $\theta$ and the last term
that appears due to the time dependence of the length 777This term can be
compared with a damping term due to its dependence on the first derivative of
$\theta$. [7]. If we expand Eq. (2) in a two terms Taylor series in $\theta$
around 0, we get

$r^{2}\ddot{\theta}+r[g+\ddot{z}_{0}]\left\\{\theta-\frac{\theta^{3}}{6}\right\\}+r\ddot{x}_{0}\frac{\theta^{2}}{2}+2r\dot{r}\dot{\theta}=-r\ddot{x}_{0}\
\enspace,$ (3)

that is a Mathieu equation with quadratic an cubic nonlinearities and a direct
excitation term.

##  2 Results

In [8] the authors already presented an analysis based on numerical and
experimental results for the sloshing dynamics of the coffee. Here, we
emphasize our effort in approximated solutions via perturbation methods,
namely: the Method of Averaging and the Method of Multiple Scales [2].

###  2.1 Vertical motion of the pivot: Method of Averaging

If we consider just the vertical motion of the pivot, the (nondimensional)
differential equations turn to be

$u^{\prime\prime}+[1+\varepsilon\lambda\Omega^{2}\cos\Omega\tau]\left(u-\frac{\varepsilon^{2}u^{3}}{6}\right)=0\enspace,$
(4)

where

$\displaystyle\tau=$ $\displaystyle\omega t$ $\displaystyle\omega_{0}^{2}=$
$\displaystyle\frac{g}{r_{0}}$ $\displaystyle\varepsilon\lambda=$
$\displaystyle-\frac{\Delta z}{r_{0}}$ $\displaystyle\Omega=$
$\displaystyle\frac{\omega}{\omega_{0}}\enspace,$

being $\lambda$ the ratio between the _effective length_ and amplitude of
oscillation, and $\varepsilon$ a bookkeeping parameter.

We can study this nonlinear differential equation using the Method of
Averaging. In Cartesian coordinates, we propose a solution of the form

$\displaystyle u(\tau)=$ $\displaystyle
X(\tau)\cos\left(\frac{\Omega}{2}\tau\right)+Y(\tau)\sin\left(\frac{\Omega}{2}\tau\right)$
(5) $\displaystyle u^{\prime}(\tau)=$
$\displaystyle-\frac{\Omega}{2}X(\tau)\sin\left(\frac{\Omega}{2}\tau\right)+\frac{\Omega}{2}Y(\tau)\cos\left(\frac{\Omega}{2}\tau\right)$
(6)

that is a constrained coordinate transformation. From the restriction in the
first derivative, the substitution in the differential equation and averaging
over one period we obtain two (autonomous) differential equations for
$X^{\prime}$ and $Y^{\prime}$ and introducing a detuning parameter
($\varepsilon\sigma=\Omega-2$), namely

$\displaystyle X^{\prime}=$ $\displaystyle\frac{\pi\varepsilon
Y}{12(\varepsilon\sigma+2)^{2}}[2\lambda\varepsilon^{2}Y^{2}-3\varepsilon
Y^{2}-3\varepsilon X^{2}-6\varepsilon\sigma^{2}-24\sigma-12\lambda]$ (7)
$\displaystyle Y^{\prime}=$ $\displaystyle\frac{\pi\varepsilon
X}{12(\varepsilon\sigma+2)^{2}}[3\varepsilon
Y^{2}+2\lambda\varepsilon^{2}X^{2}+3\varepsilon
X^{2}+6\varepsilon\sigma^{2}+24\sigma-12\lambda]\enspace,$ (8)

we can solve this differential equation for the transient response, or equate
$(X^{\prime},Y^{\prime})=(0,0)$ and solve for $X$ and $Y$ to get the steady-
state solution. The solutions for the steady-state case are

$X=Y=0$

or

$\displaystyle X^{2}=$
$\displaystyle-\frac{3}{\lambda\varepsilon^{3}}(\varepsilon^{2}\sigma^{2}+4\varepsilon\sigma-2\lambda\varepsilon+6)$
$\displaystyle Y^{2}=$
$\displaystyle\frac{3}{\lambda\varepsilon^{3}}(\varepsilon^{2}\sigma^{2}+4\varepsilon\sigma+2\lambda\varepsilon+6)\enspace.$

We can compute the Jacobian to obtain

$J(X,Y)=\begin{bmatrix}-\frac{\pi\,{\varepsilon}^{2}\,X\,Y}{2\,{\left(\varepsilon\,\sigma+2\right)}^{2}}&\frac{\pi\,\varepsilon\,\left(2\,\lambda\,{\varepsilon}^{2}\,Y^{2}-3\,\varepsilon\,Y^{2}-\varepsilon\,X^{2}-2\,\varepsilon\,\sigma^{2}-8\,\sigma-4\,\lambda\right)}{4\,{\left(\varepsilon\,\sigma+2\right)}^{2}}\cr\frac{\pi\,\varepsilon\,\left(\varepsilon\,Y^{2}+2\,\lambda\,{\varepsilon}^{2}\,X^{2}+3\,\varepsilon\,X^{2}+2\,\varepsilon\,\sigma^{2}+8\,\sigma-4\,\lambda\right)}{4\,{\left(\varepsilon\,\sigma+2\right)}^{2}}&\frac{\pi\,{\varepsilon}^{2}\,X\,Y}{2\,{\left(\varepsilon\,\sigma+2\right)}^{2}}\end{bmatrix}$

around $(0,0)$ gives

$J(0,0)=\frac{\pi\varepsilon}{2(\varepsilon\sigma+2)^{2}}\begin{bmatrix}0&-(\varepsilon\sigma^{2}+4\sigma+2\lambda)\\\
\varepsilon\sigma^{2}+4\sigma+2\lambda&0\end{bmatrix}\enspace,$

the stability of this fixed point is given by the determinant (since the trace
is equal to zero)

$\Delta=\frac{\pi^{2}\varepsilon^{2}}{4(\varepsilon+2)^{4}}(\varepsilon\sigma^{2}+4\sigma-2\lambda)(\varepsilon\sigma^{2}+4\sigma+2\lambda)\enspace,$

this gives the critical values

$\lambda=-\frac{1}{2}(\varepsilon\sigma_{1}^{2}+4\sigma_{1}),\quad\lambda=\frac{1}{2}(\varepsilon\sigma_{2}^{2}+4\sigma_{2})\enspace.$

The stability of the system is depicted in figure 2.

Figure 2: Stability of the system around the point (0,0). Regions in blue
refer to stable behavior, while regions in red to unstable behavior.

For the other fixed point the determinant is

$\Delta=\frac{{\pi}^{2}\,\left({\varepsilon}^{2}\,{\sigma}^{2}+4\,\varepsilon\,\sigma-2\,\lambda\,\varepsilon+6\right)\,\left({\varepsilon}^{2}\,{\sigma}^{2}+4\,\varepsilon\,\sigma+2\,\lambda\,\varepsilon+6\right)}{{\left(\varepsilon\,\sigma+2\right)}^{4}}$

this gives the critical values

$\lambda=-\frac{1}{2\varepsilon}(\varepsilon^{2}\sigma_{1}^{2}+4\varepsilon\sigma_{1}+6),\quad\lambda=\frac{1}{2\varepsilon}(\varepsilon^{2}\sigma_{2}^{2}+4\varepsilon\sigma_{2}+6)\enspace.$

This curves present an extremum at $\sigma=-2\epsilon$. The stability of the
system is depicted in figure 3.

Figure 3: Stability of the system around the second fixed point. Regions in
blue refer to stable behavior, while regions in red to unstable behavior.

For parametric pumping, the optimum effect is achieved by pumping at twice the
natural frequency of the pendulum. This is constant for the linear pendulum,
due to isochrony condition. When the pendulum increases the amplitude the
period increase as

$T=4\sqrt{\frac{r_{0}}{2g}}\int\limits_{0}^{\theta_{0}}\frac{1}{\sqrt{\cos\theta-\cos\theta_{0}}}d\theta\enspace,$

being $\theta_{0}$ the initial amplitude, and the angular frequency
diminishes.888Hence, as the amplitude increases the optimal frequency of
parametric excitation decreases. This gives us an explanation about the wider
stability for the second (non-trivial) solution. Around (0,0) the linear
analysis of stability resembles the linear case, while in the second fixed
point we are not around this point and the amplitude of the oscillation would
affect the (parametric) resonances of the system.

###  2.2 Vertical and horizontal motion of the pivot: Method of Multiple
Scales

If we keep both, the vertical and horizontal motion of the pivot we get
parametric and direct excitation in the differential equation. After
normalization reads

$u^{\prime\prime}+[1+\varepsilon\lambda_{1}\Omega^{2}\cos\Omega\tau]\left\\{u-\frac{\varepsilon^{2}u^{3}}{6}\right\\}-\varepsilon^{2}\lambda_{2}\Omega^{2}\cos\Omega\tau\frac{u^{2}}{2}=-\lambda_{2}\Omega^{2}\cos\Omega\tau$

with

$\displaystyle\tau=\omega t{\omega_{0}}$
$\displaystyle\varepsilon\lambda_{1}=-\frac{\Delta z}{r_{0}}$
$\displaystyle\omega_{0}^{2}=\frac{g}{r_{0}}$
$\displaystyle\Omega=\frac{\omega}{\omega_{0}}$
$\displaystyle\varepsilon\lambda_{1}=-\frac{\Delta z}{r_{0}}$
$\displaystyle\enspace,$

when, one more time $\lambda_{1}$ and $\lambda_{2}$ refers to the ratios
between the amplitude of the parametric excitations and the length of the
pendulum.

In this case, it is interesting to find the resonances that the system can
achieve, that will differ from the single resonance from the linear case.
Using the method of multiple scales we propose a solution

$u=u_{0}+\varepsilon u_{1}+\varepsilon^{2}u_{2}+\cdots\enspace,$

and

$t=T_{0}+\varepsilon T_{1}+\varepsilon^{2}T_{2}+\cdots\enspace.$

Replacing, and grouping by powers of $\varepsilon$ we get

$\displaystyle\varepsilon^{0}:D_{0}^{2}u_{0}+u_{0}=$
$\displaystyle=-\lambda_{2}\Omega^{2}\cos\Omega T_{0}$
$\displaystyle\varepsilon^{1}:D_{0}^{2}u_{1}+u_{1}=$
$\displaystyle-2D_{0}D_{1}u_{0}-u_{0}\lambda_{1}\cos\Omega
T_{0}+\frac{1}{2}\lambda_{2}\Omega^{2}u_{0}^{2}\cos\Omega T_{0}$
$\displaystyle\varepsilon^{2}:D_{0}^{2}u_{2}+u_{2}=$
$\displaystyle-\lambda_{1}u_{1}\cos\Omega T_{0}-2D_{0}D_{1}u_{1}$
$\displaystyle-
D_{1}^{2}u_{0}+\frac{u_{0}^{3}}{6}+\lambda_{2}u_{0}u_{1}\Omega^{2}\cos\Omega
T_{0}$

Solving sequentially the set of equations, and enforcing secular terms to be
zero we get as resonances

$\Omega\in\left\\{\frac{1}{2},\ 1,\ 2,\ 3\right\\}$

The appearence of three of the resonances was obtained via numerical
integration and are presented in Figure 4. Figure 5 show solutions for a range
of frequencies around the direct and primary parametric resonances, due to
nonlinearities the higher values are not present exactly at $\Omega=1,\,2$.

(a) $\varepsilon=0.5$, $\lambda_{1}=0.2$ and $\lambda_{2}=0.2$.

(b) $\varepsilon=0.2$, $\lambda_{1}=0.1$ and $\lambda_{2}=0.3$.

Figure 4: Time response for two sets of different parameters of the equation;
in both cases $\theta(0)=\pi/4$ and $\dot{\theta}(0)=0$.

(a) Solution around $\Omega=1$.

(b) Solution around $\Omega=2$.

Figure 5: Time response around the primary and direct resonances. The
parameters used were: $\lambda_{1}=0.2$, $\lambda_{2}=0.1$, $\epsilon_{0}.2$,
$\theta(0)=\pi/30$, and $\dot{\theta}(0)=0$.

###  2.3 Numerical Results

We computed the numerical solution using odeint from Scipy999odeint solve a
system of ordinary differential equations using lsoda from the Fortran library
odepack.
http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html.

(a) $\Omega=0.00$.

(b) $\Omega=1.01$.

(c) $\Omega=1.96$.

(d) $\Omega=2.60$.

(e) $\Omega=2.82$.

(f) $\Omega=2.99$.

(g) $\Omega=3.54$.

(h) $\Omega=4.00$.

Figure 5: Phase portrait for different excitation frequencies. In all the
cases $\lambda_{1}=0.2$, $\lambda_{2}=0.1$, $\varepsilon=0.2$, $\theta(0)=0.1$
and $\dot{\theta}(0)=0$. The red square is the initial point in the phase
space, while the blue one is the finishing one after $\tau=500$

##  3 Conclusions

A model consisting of a planar pendulum was used to model the sloshing
dynamics of coffee. Rotating sloshing was neglected since this will need two
degrees of freedom (a spherical pendulum model).

The vertically excited pendulum presents a parametric pumping with an
effective stiffness that varies with time. This is similar to the Mathieu
equation with Duffing-like nonlinearity. When excited both, vertically and
horizontally the system presents parametric and direct excitation. Once again,
the equation reads similar to Mathieu equation with quadratic and cubic
nonlinearities, which leads to multiple resonances.

If we consider the aspect ratio of cups _constant_ we can have as design
parameter the radius of the cup 101010See Appendix A.. Currently, the
frequency ranges for the cups and regular walking are overlapped, but we can
change this fact by _tuning_ the size of the cup. An increase in the size will
lead to a decrease in the natural frequency and will cease the overlapping.
From a practical point of view, this is not feasible, since one doesn’t want
to have a huge cup of coffee.

Small values of $\lambda$ will translate in longer walking with coffee without
spilling, ehich can be achieved with higher sizes of cup or smaller amplitudes
of oscillation—smoother walking. The latter can be obtained _focusing_ in not
spilling while walking, but this can be considered as a controlled system.
This is one of the facts discussed in [8].

A single harmonic excitation was considered due to its main importance [8],
but the real walking is a complex signal with a broadband frequency content
and several resonances will appear in that case.

## Appendix A: Physical Parameters

Figure 6 shows a schematic of the cup of coffee. According to [8], the
parameters for common cups are

$\displaystyle R\in[2.5,\ 6.7]\,\mathrm{cm}$ $\displaystyle H\in[5.7,\
8.9]\,\mathrm{cm}$ $\displaystyle h\in[5.10,\ 10]\,\mathrm{cm}$
$\displaystyle\alpha\in[8^{\circ},\ 16^{\circ}]$

Figure 6: Geometrical description of the cup of coffee.

This combination of parameters gives a range of frequencies for the cup

$f_{\text{cup}}\in[2.6,\ 4.3]\,\mathrm{Hz}\enspace,$

and the frequency of walking is [8]

$f_{\text{step}}\in[1,\ 2.5]\,\mathrm{Hz}\enspace.$

According to [3] the lumped parameters for a single degree of freedom system
are

$\displaystyle m=$ $\displaystyle\frac{c_{1}m_{0}}{1-c_{1}}$ $\displaystyle
m_{0}=$ $\displaystyle\rho\pi(H-h)R^{2}$ $\displaystyle c_{1}=$
$\displaystyle\frac{R}{2.2}\tanh\left(1.84\frac{H}{R}\right)$ $\displaystyle
r_{0}=$
$\displaystyle\frac{R}{1.84}\tanh\left(1.84\frac{H}{R}\right)\enspace,$

for simple calculation we can take $H/R\approx 1.5$, what yields

$r_{0}\approx 0.54R\enspace.$

## References

  * [1] Baker, G. L., and Blackburn, J. A. (2005). The pendulum: a case study in physics (Vol. 8). Oxford: Oxford University Press.
  * [2] Holmes, Mark H. Introduction to perturbation methods (2012). Vol. 20. Springer Science & Business Media.
  * [3] Ibrahim, R. A. (2005). Liquid sloshing dynamics: theory and applications. Cambridge University Press.
  * [4] Poole, C. P., Safko, J. L., and Goldstein, H. (2001). Classical Mechanics.
  * [5] Kana, D. D. (1989). Validated spherical pendulum model for rotary liquid slosh. Journal of Spacecraft and Rockets, 26(3), 188-195
  * [6] Kana, D. D., and Fox, D. J. (1995). Distinguishing the transition to chaos in a spherical pendulum. Chaos: An Interdisciplinary Journal of Nonlinear Science, 5(1), 298-310.
  * [7] Lanczos, Cornelius. The variational principles of mechanics. Vol. 4. Courier Dover Publications, 1970.
  * [8] Mayer, H. C., and Krechetnikov, R. (2012). Walking with coffee: Why does it spill?. Physical Review E, 85(4), 046117.
  * [9] NASA (1969), Slosh suppression, May 1969, 36p.
  * [10] Miles, J. (1984). Resonant motion of a spherical pendulum. Physica D: Nonlinear Phenomena, 11(3), 309-323.
  * [11] Rucci, C. D. (1979). U.S. Patent No. 4,141,462. Washington, DC: U.S. Patent and Trademark Office.
  * [12] Sanmartin, J. R. (1984). O Botafumeiro: Parametric pumping in the middle ages. American Journal of Physics, 52(10), 937-945.