text:
# Bootstrapping Clustered Data in R using lmeresampler

Linear mixed-effects models are commonly used to analyze clustered data structures. There are numerous packages to fit these models in R and conduct likelihood-based inference. The implementation of resampling-based procedures for inference are more limited. In this paper, we introduce the lmeresampler package for bootstrapping nested linear mixed-effects models fit via lme4 or nlme. Bootstrap estimation allows for bias correction, adjusted standard errors and confidence intervals for small samples sizes and when distributional assumptions break down. We will also illustrate how bootstrap resampling can be used to diagnose this model class. In addition, lmeresampler makes it easy to construct interval estimates of functions of model parameters.

Adam Loy https://aloy.rbind.io/ (Carleton College) , Jenna Korobova (Carleton College)
2023-02-10

## Introduction

Clustered data structures occur in a wide range of studies. For example, students are organized within classrooms, schools, and districts, imposing a correlation structure that must be accounted for in the modeling process. Similarly, cholesterol measurements could be tracked across time for a number of subjects, resulting in measurements being grouped by subject. Other names for clustered data structures include grouped, nested, multilevel, hierarchical, longitudinal, repeated measurements, and blocked data. The covariance structure imposed by clustered data is commonly modeled using linear mixed-effects (LME) models, also referred to as hierarchical linear or multilevel linear models .

In this paper, we restrict attention to the Gaussian response LME model for clustered data structures. For cluster $$i=1, \ldots, g$$, this model is expressed as $$$\underset{(n_i \times 1)}{\boldsymbol{y}_i} = \underset{(n_i \times p)}{\boldsymbol{X}_i} \ \underset{(p \times 1)}{\boldsymbol{\beta}} + \underset{(n_i \times q)}{\boldsymbol{Z}_i} \ \underset{(q \times 1)}{\boldsymbol{b}_i} + \underset{(n_i \times 1)}{\boldsymbol{\varepsilon}_i}, \tag{1}$$$ where $$$\boldsymbol{\varepsilon}_i \sim \mathcal{N}(\boldsymbol{0}, \sigma^2 \boldsymbol{I}_{n_i}), \boldsymbol{b}_i \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{D}), \tag{2}$$$ where $$\boldsymbol{\varepsilon}_i$$ and $$\boldsymbol{b}_i$$ are independent for all $$i$$, $$\boldsymbol{\varepsilon}_i$$ is independent of $$\boldsymbol{\varepsilon}_j$$ for $$i\ne j$$, and $$\boldsymbol{b}_i$$ is independent of $$\boldsymbol{b}_j$$ for $$i\ne j$$. Here, $$\boldsymbol{0}$$ denotes a vector of zeros of length $$n_i$$ (the number of observations in group $$i$$), $$\boldsymbol{I}_{n_i}$$ denotes the $$n_i$$-dimensional identity matrix, and $$\boldsymbol{D}$$ is a $$q \times q$$ covariance matrix.

In R, the two most popular packages to fit LME models are nlme and lme4 . Both packages fit LME models using either maximum likelihood or restricted maximum likelihood methods. These methods rely on the distributional assumptions placed on the residual quantities, $$\boldsymbol{\varepsilon}_i$$ and $$\boldsymbol{b}_i$$, as well as large enough sample sizes. While some aspects of LME models are quite robust to model misspecification, others are more sensitive, leading to biased estimates and/or incorrect standard errors. Jacqmin-Gadda et al. (2007) found that inference for the fixed effects is robust if the distribution of the error terms, $$\boldsymbol{\varepsilon}_i$$, is non-normal or heteroscedastic, but that variance components and random effects are biased if the covariance structure for the error terms is misspecified. The fixed intercept is not robust to misspecification of the random effects . In addition, misspecification of the random effects distribution can lead to biased estimates of the variance components and undercoverage of the confidence intervals . In cases where distributional assumptions are suspect, bootstrapping provides an alternative inferential approach that leads to consistent, bias-corrected parameter estimates, standard errors, and confidence intervals. Standard errors and confidence intervals for functions of model parameters are also easily calculated using a bootstrap procedure, and are available even in situations where closed-form solutions are not.

A variety of bootstrap procedures for clustered data and the LME model have been proposed and investigated, including the cases (non-parametric) bootstrap, the residual bootstrap, the parametric bootstrap, the random effect block (REB) bootstrap, and the wild bootstrap . Van der Leeden et al. (2008) provide a thorough overview of bootstrapping LME models. Sánchez-Espigares and Ocaña (2009) developed a full-featured bootstrapping framework for LME models in R; however, this package is not available and there appears to be no active development. Consequently, R users must pick and choose packages based on what bootstrap procedure they wish to implement. The parametric bootstrap is available in lme4 via the bootMer() function, as is the semi-parametric (residual) bootstrap proposed by Morris (2002). The simulateY() function in nlmeU makes it easy to simulate values of the response variable for lme model objects; however, the user is required to implement the remainder of the parametric bootstrap. Chambers and Chandra (2013) made R code available to implement their REB bootstrap as well as the parametric bootstrap and the residual bootstrap proposed by Carpenter et al. (2003) for LME models fix via nlme::lme(). Unfortunately, these functions were not published on CRAN or extended to models fit via lme4::lmer(). Bootstrap procedures for specific inferential tasks have also been implemented. The parametric bootstrap has been implemented to carry out likelihood ratio tests for the presence of variance components in RLRsim . pbkrtest provides a Kenward-Roger approximation for performing F-tests using a parametric bootstrap approach for LME models and generalized LME models. rptR uses the parametric bootstrap to estimate repeatabilities for LME models . Finally, constrained inference for LMEs using the residual bootstrap is implemented in CLME .

In this paper, we introduce the lmeresampler package which implements a suite of bootstrap methods for LME models fit using either nlme or lme4 in a single package. lmeresampler provides a unified bootstrap() command to reduce cognitive load on the user and also provides access to procedures that were not previously available on CRAN. In the next section, we will clarify some notation and terminology for LME models. In Bootstrap procedures for clustered data, we provide an overview of the bootstrap methods implemented in lmeresampler. We then give an Overview of lmeresampler, discuss a variety of Example applications for LME models, and show how to Bootstrapping in parallel works in lmeresampler. We conclude with a Summary and highlight areas for future development.

## Bootstrap procedures for clustered data

In lmeresampler, we implement five bootstrap procedures for clustered data: a cases (i.e., nonparamteric) bootstrap, a residual bootstrap (i.e., semiparametric), a parametric bootstrap, a wild bootstrap, and a block bootstrap. In this section, we provide an overview of these bootstrap approaches. Our discussion focuses on two-level models, but procedures generalize to higher-level models unless otherwise noted.

### The cases bootstrap

The cases bootstrap is a fully non-parametric bootstrap that resamples the clusters in the data set to generate bootstrap resamples. Depending on the nature of the data, this resampling can be done only for the top-level cluster, only at the observation-level within a cluster, or at both levels. The choice of the exact resampling scheme should be dictated by the way the data were generated, since the cases bootstrap generates new data by mimicking this process. Van der Leeden et al. (2008) provide a cogent explanation of how to select a resampling scheme. To help ground the idea of resampling, consider a two-level hierarchical data set where students are organized into schools.

One version of the cases bootstrap is implemented by only resampling the clusters. This version of the bootstrap is what Field and Welsh (2007) term the cluster bootstrap and Goldstein (2011) term the non-parametric bootstrap. We would choose this resampling scheme, for example, if schools were chosen at random and then all students within each school were observed. In this case, the bootstrap proceeds as follows:

1. Draw a random sample, with replacement, of size $$g$$ from the clusters.
2. For each selected cluster, $$k$$, extract all of the cases to form the bootstrap sample $$\left(\boldsymbol{y}^*_k, \boldsymbol{X}^*_k, \boldsymbol{Z}^*_k \right)$$. Since entire clusters are sampled, the total sample size may differ from the that of the original data set.
3. Refit the model to the bootstrap sample and extract the parameter estimates of interest.
4. Repeat steps 1–3 $$B$$ times.

An alternative version of the cases bootstrap only resamples the observations within clusters, which makes sense in our example if the schools were fixed and students were randomly sampled within schools.

1. For each cluster $$i=1,\ldots,g$$, draw a random sample of the rows for cluster $$i$$, with replacement, to form the bootstrap sample $$\left(\boldsymbol{y}^*_i, \boldsymbol{X}^*_i, \boldsymbol{Z}^*_i \right)$$.
2. Refit the model to the bootstrap sample and extract the parameter estimates of interest.
3. Repeat steps 1–2 $$B$$ times.

A third version of the cases bootstrap resamples both clusters and cases within clusters. This is what Field and Welsh (2007) term the two-state bootstrap. We would choose this resampling scheme if both schools and students were sampled during the data collection process.

All three versions of the case bootstrap are implemented in lmeresampler. We explain how the resampling is specified in our Overview of lmeresampler. Regardless of which version of the cases bootstrap you choose, it requires the weakest conditions: it only requires that the hierarchical structure in the data set is correctly specified.

### The parametric bootstrap

The parametric bootstrap simulates random effects and error terms from the fitted distributions to form bootstrap resamples. Consequently, it requires the strongest conditions—that is, the parametric bootstrap assumes that the model, as specified by Equations (1) and (2), is correctly specified.

Let $$\widehat{\boldsymbol \beta}$$, $$\widehat{\boldsymbol D}$$, and $$\widehat{\sigma}^2$$ be maximum likelihood or restricted maximum likelihood estimates of the fixed effects and variance components from the fitted model. The parametric bootstrap is then implemented through the following steps:

1. Simulate $$g$$ error term vectors, $$\boldsymbol{e}_i^*$$, of length $$n_i$$ from $$\mathcal{N}\left(\boldsymbol{0},\widehat{\sigma}^2 \boldsymbol{I}_{n_i} \right)$$.
2. Simulate $$g$$ random effects vectors, $$\boldsymbol{b}_i^*$$, from $$\mathcal{N}\left(\boldsymbol{0},\widehat{\boldsymbol{D}} \right)$$.
3. Generate bootstrap responses $$\boldsymbol{y}^*_i = \boldsymbol{X}_i \widehat{\boldsymbol{\beta}} + \boldsymbol{Z}_i \boldsymbol{b}_i^* + \boldsymbol{e}_i^*$$.
4. Refit the model to the bootstrap responses and extract the parameter estimates of interest.
5. Repeat steps 2–4 $$B$$ times.

### The residual bootstrap

The residual bootstrap resamples the residual quantities from the fitted LME model in order to generate bootstrap resamples. There are three general types of residuals for LME models . correspond to the errors made using the marginal predictions, $$\widehat{\boldsymbol r}_i = \boldsymbol{y_i} - \widehat{\boldsymbol{y}}_i = \boldsymbol{y}_i - \boldsymbol{X}_i \widehat{\boldsymbol{\beta}}$$. correspond to the errors made using predictions conditioned on the clusters, $$\widehat{\boldsymbol{e}}_i = \boldsymbol{y}_i - \widehat{\boldsymbol{y}}_i| \boldsymbol{b}_i = \boldsymbol{y}_i - \boldsymbol{X}_i \widehat{\boldsymbol{\beta}} - \boldsymbol{Z}_i \widehat{\boldsymbol{b}}_i$$. The are the last type of residual quantity and are defined as $$\widehat{\boldsymbol{b}}_i = \widehat{\boldsymbol{D}} \boldsymbol{Z}_i^\prime \widehat{\boldsymbol{V}}_i^{-1}\left(\boldsymbol{y}_i - \boldsymbol{X}_i \widehat{\boldsymbol{\beta}} \right)$$ where $$\widehat{\boldsymbol{V}}_i = \boldsymbol{Z}_i \widehat{\boldsymbol{D}}\boldsymbol{Z}_i^\prime + \widehat{\sigma}^2 \boldsymbol{I}_{n_i}$$.

A naive implementation of the residual bootstrap would draw random samples, with replacement, from the estimated conditional residuals and the best linear unbiased predictors (BLUPS); however, this will consistently underestimate the variability in the data because the residuals are shrunken toward zero . Carpenter et al. (2003) solve this problem by “reflating” the residual quantities so that the empirical covariance matrices match the estimated covariance matrices prior to resampling:

1. Fit the model and calculate the empirical BLUPs, $$\widehat{\boldsymbol{b}}_i$$, and the predicted conditional residuals, $$\widehat{\boldsymbol{e}}_i$$.

2. Mean center each residual quantity and reflate the centered residuals. Only the process to reflate the predicted random effects is discussed below, but the process is analogous for the conditional residuals.

1. Arrange the random effects into a $$g \times q$$ matrix, where each row contains the predicted random effects from a single group. Denote this matrix as $$\widehat{\boldsymbol{U}}$$. Define $$\boldsymbol{\widehat{\Gamma}} = \boldsymbol{I}_g \otimes \widehat{\boldsymbol{D}}$$, the block diagonal covariance matrix of $$\widehat{\boldsymbol{U}}$$.
2. Calculate the empirical covariance matrix as $$\boldsymbol{S} = \widehat{\boldsymbol{U}}^\prime \widehat{\boldsymbol{U}} / g$$. We follow the approach of Carpenter et al. (2003), dividing by the number of groups, $$g$$, rather than $$g-1$$.
3. Find a transformation of $$\widehat{\boldsymbol{U}}$$, $$\widehat{\boldsymbol{U}}^* = \widehat{\boldsymbol{U}} \boldsymbol{A}$$, such that $$\widehat{\boldsymbol{U}}^{*\prime} \widehat{\boldsymbol{U}}^* / g = \boldsymbol{\widehat{\Gamma}}$$. Specifically, we will find $$\boldsymbol{A}$$ such that $$\boldsymbol{A}^\prime \widehat{\boldsymbol{U}}^\prime \widehat{\boldsymbol{U}} \boldsymbol{A} / g = \boldsymbol{A}^\prime \boldsymbol{SA} = \boldsymbol{\widehat{\Gamma}}$$. The choice of $$\boldsymbol{A}$$ is not unique, so we use the recommendation given by Carpenter et al. (2003): $$\boldsymbol{A} = \left(\boldsymbol{L}_D \boldsymbol{L}_S^{-1} \right)^\prime$$ where $$\boldsymbol{L}_D$$ and $$\boldsymbol{L}_S$$ are the Cholesky factors of $$\boldsymbol{\widehat{\Gamma}}$$ and $$\boldsymbol{S}$$, respectively.
3. Draw a random sample, with replacement, from the set $$\lbrace \boldsymbol{u}_i^* \rbrace$$ of size $$g$$, where $$\boldsymbol{u}_i^*$$ is the $$i$$th row of the centered and reflated random effects matrix, $$\widehat{\boldsymbol{U}}^*$$.

4. Draw $$g$$ random samples, with replacement, of sizes $$n_i$$ from the set of the centered and reflated conditional residuals, $$\lbrace \boldsymbol{e}_i^* \rbrace$$.

5. Generate the bootstrap responses, $$\boldsymbol{y}^*_i$$, using the fitted model equation: $$\boldsymbol{y}^*_i = \boldsymbol{X}_i \widehat{\boldsymbol{\beta}} + \boldsymbol{Z}_i \widehat{\boldsymbol{u}}^*_i + \boldsymbol{e}_i^*$$.

6. Refit the model to the bootstrap responses and extract the parameter estimates of interest.

7. Repeat steps 3–6 B times.

Notice that the residual bootstrap is a semiparametric bootstrap, since it depends on the model structure (both the mean function and the covariance structure) but not the distributional conditions .

### The random-effects block bootstrap

Another semiparametric bootstrap is the random effect block (REB) bootstrap proposed by Chambers and Chandra (2013). The REB bootstrap can be viewed as a version of the residual bootstrap where conditional residuals are resampled from within clusters (i.e., blocks) to allow for weaker assumptions on the covariance structure of the residuals. The residual bootstrap requires that the conditional residuals are independent and identically distributed, whereas the REB bootstrap relaxes this to only require that the covariance structure of the error terms is similar across clusters. In addition, the REB bootstrap utilizes the marginal residuals to calculate non-parametric predicted random effects rather than relying on the model-based empirical best linear unbiased predictors (EBLUPS). Chambers and Chandra (2013) developed three versions of the REB bootstrap, all of which have been implemented in lmeresampler. We refer the reader to Chambers and Chandra (2013) for a discussion of when each should be used. It’s important to note that at the time of this writing, that the REB bootstrap has only been explored and implemented for use with two-level models.

#### REB/0

The base algorithm for the REB bootstrap (also known as REB/0) is as follows:

1. Calculate non-parametric residual quantities for the model.
1. Calculate the marginal residuals for each group, $$\widehat{\boldsymbol{r}}_i = \boldsymbol{y}_i - \boldsymbol{X}_i \widehat{\boldsymbol{\beta}}$$.
2. Calculate the non-parametric predicted random effects, $$\tilde{\boldsymbol{b}}_i = \left( \boldsymbol{Z}_i^\prime \boldsymbol{Z}_i \right)^{-1} \boldsymbol{Z}^\prime_i \boldsymbol{r}_i$$.
3. Calculate the non-parametric conditional residuals using the residuals quantities obtained in the previous two steps, $$\boldsymbol{\tilde{e}}_i = \widehat{\boldsymbol{r}}_i - \boldsymbol{Z}_i \boldsymbol{\tilde{b}}_i$$.
2. Take a random sample, with replacement, of size $$g$$ from the set $$\lbrace \tilde{\boldsymbol{b}}_i \rbrace$$. Denote these resampled random effects as $$\boldsymbol{\tilde{b}}^*_i$$.
3. Take a random sample, with replacement, of size $$g$$ from the cluster ids. For each sampled cluster, draw a random sample, with replacement, of size $$n_i$$ from that cluster’s vector of error terms, $$\boldsymbol{\tilde{e}}_i$$.
4. Generate bootstrap responses, $$\boldsymbol{y}^*_i$$, using the fitted model equation: $$\boldsymbol{y}^*_i = \boldsymbol{X}_i \widehat{\beta} + \boldsymbol{Z}_i \boldsymbol{\tilde{b}}^*_i + \boldsymbol{\tilde{e}}^*_i$$.
5. Refit the model to the bootstrap sample and extract the parameter estimates of interest.
6. Repeat steps 2–5 $$B$$ times.

#### REB/1

The first variation of the REB bootstrap (REB/1) zero centers and reflates the residual quantities prior to resampling in order to satisfy the conditions for consistency . This is the same process outlined in Step 2 of the residual bootstrap outlined above.

#### REB/2

The second variation of the REB bootstrap (REB/2 or postscaled REB) addresses two issues: potential non-zero covariances in the joint bootstrap distribution of the variance components and potential bias in the parameter estimates. After the REB/0 algorithm is run, the following post processing is performed:

##### Uncorrelate the variance components.

To uncorrelate the bootstrap estimates of the variance components produced by REB/0, Chambers and Chandra (2013) propose the following procedure:

1. Apply natural logarithms to the bootstrap distribution of each variance component and form the following matrices. Note that we use $$\nu$$ to denote the total number of variance components.

• $$\boldsymbol{S}^*$$: a $$B \times \nu$$ matrix formed by binding the columns of these distributions together
• $$\boldsymbol{M}^*$$: a $$B \times \nu$$ matrix where each column contains the column mean from the corresponding column in $$\boldsymbol{S}^*$$
• $$\boldsymbol{D}^*$$: a $$B \times \nu$$ matrix where each column contains the column standard deviation from the corresponding column in $$\boldsymbol{S}^*$$
2. Calculate $$\boldsymbol{C}^*=\mathrm{cov} \left( \boldsymbol{S}^* \right)$$.

3. Calculate $$\boldsymbol{L}^* = \boldsymbol{M}^* + \lbrace \left(\boldsymbol{S}^* - \boldsymbol{M}^*\right)\boldsymbol{C}^{*-1/2} \rbrace \circ \boldsymbol{D}^*$$, where $$\circ$$ denotes the Hadamard (elementwise) product.

4. Exponentiate the elements of $$\boldsymbol{L}^*$$. The columns of $$\boldsymbol{L}^*$$ are then uncorrelated versions of the bootstrap variance components.

##### Center the bootstrap estimates at the original estimate values.

To correct bias in the estimation of the fixed effects, apply a mean correction. For each parameter estimate, $$\widehat{\beta}_k$$, adjust the bootstrapped estimates, $$\widehat{\boldsymbol{\beta}}_k^*$$ as follows: $$\widetilde{\boldsymbol{\beta}}_k^{*} = \widehat{\beta}_k \boldsymbol{1}_{B} + \widehat{\boldsymbol{\beta}}_k^* -{\rm avg} \left( \widehat{ \boldsymbol{\beta} }_k^* \right)$$. To correct bias in the estimation of the variance components, apply a ratio correction. For each estimated variance component , $$\widehat{\sigma}^2_v$$, adjust the uncorrelated bootstrapped estimates, $$\widehat{\boldsymbol{\sigma}}_v^{2*}$$ as follows: $$\widetilde{\boldsymbol{\sigma}}_{v}^{2*} = \widehat{\boldsymbol{\sigma}}_v^{2*} \circ \lbrace \widehat{\sigma}^2_v / {\rm avg} \left( \widehat{\boldsymbol{\sigma}}_v^{2*} \right) \rbrace$$

### The wild bootstrap

The wild bootstrap also relaxes the assumptions made on the error terms of the model, allowing heteroscedasticity both within and across groups. The wild bootstrap is well developed for the ordinary regression model and Modugno and Giannerini (2015) adapt it for the nested LME model.

To begin, we can re-express model (1) as

$$$\boldsymbol{y}_i = \boldsymbol{X}_i \boldsymbol{\beta} + \boldsymbol{v}_i, \text{ where } \boldsymbol{v}_i = \boldsymbol{Z}_i \boldsymbol{b}_i + \boldsymbol{\varepsilon}_i. \tag{3}$$$

The wild bootstrap proceeds as follows:

1. Draw a random sample, $$w_1, w_2, \ldots, w_g$$, from an auxiliary distribution with mean zero and unit variance.

2. Generate bootstrap responses using the re-expressed model equation (3): $$\boldsymbol{y}^*_i = \boldsymbol{X}_i \widehat{\boldsymbol{\beta}} + \tilde{\boldsymbol{v}}_i w_j$$, where $$\tilde{\boldsymbol{v}}_i$$ is a heteroscedasticity consistent covariance matrix estimator. Modugno and Giannerini (2015) suggest using what Flachaire (2005) calls $${{\rm HC}_2}$$ or $${{\rm HC}_3}$$ in the regression context:

\begin{align} {{\rm HC}_2}&: \tilde{\boldsymbol{v}}_i = {\rm diag} \left( \boldsymbol{I}_{n_i} - \boldsymbol{H}_i \right)^{-1/2} \circ \boldsymbol{r}_i\\ {{\rm HC}_3}&: \tilde{\boldsymbol{v}}_i = {\rm diag} \left( \boldsymbol{I}_{n_i} - \boldsymbol{H}_i \right) \circ \boldsymbol{r}_i, \end{align}

where $$\boldsymbol{H}_i = \boldsymbol{X}_i \left(\boldsymbol{X}_i^\prime \boldsymbol{X}_i \right)^\prime \boldsymbol{X}_i^\prime$$, the $$i$$th diagonal block of the orthogonal projection matrix, and $$\boldsymbol{r}_i$$ is the vector of marginal residuals for group $$i$$.

3. Refit the model to the bootstrap sample and extract the parameter estimates of interest.

4. Repeat steps 1–3 $$B$$ times.

Five options for the auxiliary distribution are implemented in lmeresampler:

• The two-point distribution proposed by Mammen (1993) takes value $$w_i=-(\sqrt{5} - 1)/2$$ with probability $$p=(\sqrt{5}+1) / (2 \sqrt{5})$$ and $$w_i = (\sqrt{5} + 1)/2$$ with probability $$1-p$$.

• The two-point Rademacher distribution places equal probability on $$w_i = \pm 1$$.

• The six-point distribution proposed by Webb (2013) places equal probability on $$w_i = \pm \sqrt{1/2}, \pm1, \pm \sqrt{3/2}$$.

• The standard normal distribution.

• The gamma distribution with shape parameter 4 and scale parameter $$1/2$$ (Liu 1988) that is centered to have mean zero. $$w_i^*$$ are randomly sampled from the gamma distribution, and then centered via $$w_i = w_i^* - 2$$.

Modugno and Giannerini (2015) found the Mammen distribution to be preferred based on a Monte Carlo study, but only compared it to the Rademacher distribution.

## Overview of lmeresampler

The lmeresampler package implements the five bootstrap procedures outlined in the previous section for Gaussian response LME models for clustered data structures fit using either nlme or lme4. The workhorse function in lmeresampler is

bootstrap(model, .f, type, B, resample = NULL, reb_type = NULL, hccme, aux.dist,
orig_data, .refit)

The four required parameters to are:

• model, an lme or merMod fitted model object.

• .f, a function defining the parameter(s) of interest that should be extracted/calculated for each bootstrap iteration.

• type, a character string specifying the type of bootstrap to run. Possible values include: "parametric", "residual", "reb", "wild", and "case".

• B, the number of bootstrap resamples to generate.

• .refit, whether the model should be refit to the bootstrap sample. This defaults to TRUE.

There are also five optional parameters: resample, reb_type, hccme, aux.dist, and orig_data.

• If the user sets type = "case", then they must also set resample to specify what cluster resampling scheme to use. resample requires a logical vector of length equal to the number of levels in the model. A value of TRUE in the $$i$$th position indicates that cases/clusters at that level should be resampled. For example, to only resample the clusters (i.e., level 2 units) in a two-level model, the user would specify resample = c(FALSE, TRUE).

• If the user sets type = "reb", then they must also set reb_type to indicate which version of the REB bootstrap to run. reb_type accepts the integers 0, 1, and 2 to indicate REB/0, REB/1, and REB/2, respectively.

• If the user sets type = "wild", then they must specify both hccme and aux.dist. Currently, hccme can be set to "hc2" or "hc3" and aux.dist can be set to "mammen", "rademacher", "norm", "webb", or "gamma".

• If variables are transformed within the model formula and lmer() is used to fit the LME model, then orig_data should be set to the original data frame, since this cannot be recovered from the fitted model object.

The bootstrap() function is a generic function that calls functions for each type of bootstrap. The user can call the specific bootstrap function directly if preferred. An overview of the specific bootstrap functions is given in Table 1.

Table 1: Summary of the specific bootstrap functions called by bootstrap() and their required arguments.
Bootstrap Function name Required arguments
Cases case_bootstrap model, .f, type, B, resample, orig_data, .refit
Residual resid_bootstrap model, .f, type, B, orig_data, .refit
REB reb_bootstrap model, .f, type, B, reb_type, orig_data, .refit
Wild wild_bootstrap model, .f, type, B, hccme, aux.dist, orig_data, .refit
Parametric parametric_bootstrap model, .f, type, B, orig_data, .refit

Each of the specific bootstrap functions performs four general steps:

1. Setup. Key information (parameter estimates, design matrices, etc.) is extracted from the fitted model to eliminate repetitive actions during the resampling process.

2. Resampling. The setup information is passed to an internal resampler function to generate the B bootstrap samples.

3. Refitting. The model is refit for each of the bootstrap samples and the specified parameters are extracted/calculated.

4. Clean up. An internal completion function formats the original and bootstrapped quantities to return a list to the user.

Each function returns an object of class lmeresamp, which is a list with elements outlined in Table 2. print(), summary(), plot(), and confint() methods are available for lmeresamp objects.

Table 2: Summary of the specific bootstrap functions called by bootstrap() and their required arguments.
Element Description
observed values for the original model parameter estimates.
model the original fitted model object.
.f the function call defining the parameters of interest.
replicates a $$B \times p$$ tibble containing the bootstrapped quantities. Each column contains a single bootstrap distribution.
stats a tibble containing the observed, rep.mean (bootstrap mean), se (bootstrap standard error), and bias values for each parameter.
B the number of bootstrap resamples performed.
data the original data set.
seed a vector of randomly generated seeds that are used by the bootstrap.
type a character string specifying the type of bootstrap performed.
call the user’s call to the bootstrap function.
message a list of length B giving any messages generated during refitting. An entry will be NULL if no message was generated.
warning a list of length B giving any warnings generated during refitting. An entry will be NULL if no warning was generated.
error a list of length B giving any errors encountered during refitting. An entry will be NULL if no error was encountered.

lmeresampler also provides the extract_parameters() helper function to extract the fixed effects and variance components from merMod and lme objects as a named vector.

## Example applications

### A two-level example: JSP data

As a first application of the bootstrap for nested LME models, consider the junior school project (JSP) data that is stored as jsp728 in lmeresampler. The data set is comprised of measurements taken on 728 elementary school students across 48 schools in London.

library(lmeresampler)
tibble::as_tibble(jsp728)
# A tibble: 728 × 9
mathA…¹ mathA…² gender class school normA…³ normA…⁴ schoo…⁵ mathA…⁶
<dbl>   <dbl> <fct>  <fct> <fct>    <dbl>   <dbl>   <dbl>   <dbl>
1      39      36 M      nonm… 1       1.80    1.55      22.4  13.6
2      11      19 F      manu… 1      -2.29   -0.980     22.4  -3.42
3      32      31 F      manu… 1      -0.0413  0.638     22.4   8.58
4      27      23 F      nonm… 1      -0.750  -0.460     22.4   0.579
5      36      39 F      nonm… 1       0.743   2.15      22.4  16.6
6      33      25 M      manu… 1       0.163  -0.182     22.4   2.58
7      30      27 M      manu… 1      -0.372   0.0724    22.4   4.58
8      17      14 M      manu… 1      -1.63   -1.52      22.4  -8.42
9      33      30 M      manu… 1       0.163   0.454     22.4   7.58
10      20      19 M      manu… 1      -1.40   -0.980     22.4  -3.42
# … with 718 more rows, and abbreviated variable names ¹mathAge11,
#   ²mathAge8, ³normAge11, ⁴normAge8, ⁵schoolMathAge8, ⁶mathAge8c

Suppose we wish to fit a model using the math score at age 8, gender, and the father’s social class to describe math scores at age 11, including a random intercept for school . This LME model can be fit using the lmer() function in lme4.

library(lme4)
jsp_mod <- lmer(mathAge11 ~ mathAge8 + gender + class + (1 | school), data = jsp728)

To implement the residual bootstrap to estimate the fixed effects, we can use the bootstrap() function and set type = "residual".

(jsp_boot <- bootstrap(jsp_mod, .f = fixef, type = "residual", B = 2000))
Bootstrap type: residual

Number of resamples: 2000

term   observed   rep.mean        se          bias
1    (Intercept) 14.1577509 14.1697871 0.6702835  0.0120361440
2       mathAge8  0.6388895  0.6391441 0.0248334  0.0002546297
3        genderM -0.3571922 -0.3654058 0.3391162 -0.0082135863
4 classnonmanual  0.7200815  0.7154916 0.3783149 -0.0045898872

There were 0 messages, 0 warnings, and 0 errors.

We can then calculate normal, percentile, and basic bootstrap confidence intervals via confint().

confint(jsp_boot)
# A tibble: 12 × 6
term           estimate     lower  upper type  level
<chr>             <dbl>     <dbl>  <dbl> <chr> <dbl>
1 (Intercept)      14.2   12.8      15.5   norm   0.95
2 mathAge8          0.639  0.590     0.687 norm   0.95
3 genderM          -0.357 -1.01      0.316 norm   0.95
4 classnonmanual    0.720 -0.0168    1.47  norm   0.95
5 (Intercept)      14.2   12.8      15.5   basic  0.95
6 mathAge8          0.639  0.588     0.687 basic  0.95
7 genderM          -0.357 -1.00      0.323 basic  0.95
8 classnonmanual    0.720 -0.000162  1.46  basic  0.95
9 (Intercept)      14.2   12.9      15.5   perc   0.95
10 mathAge8          0.639  0.590     0.690 perc   0.95
11 genderM          -0.357 -1.04      0.287 perc   0.95
12 classnonmanual    0.720 -0.0187    1.44  perc   0.95

The default setting is to calculate all three intervals, but this can be restricted by setting the type parameter to "norm", "basic", or "perc".

### User-specified statistics: Estimating repeatability/intraclass correlation

The beauty of the bootstrap is its flexibility. Interval estimates can be constructed for functions of model parameters that would otherwise require more complex derivations. For example, the bootstrap can be used to estimate the intraclass correlation. The intraclass correlation measures the proportion of the total variance in the response accounted for by groups, and is an important measure of repeatability in ecology and evolutionary biology . As a simple example, we’ll consider the BeetlesBody data set in rptR . This simulated data set contains information on body length (BodyL) and the Population from which the beetles were sampled. A simple Guassian-response LME model of the form

$y_{ij} = \beta_0 + b_i + \varepsilon_{ij}, \qquad b_i \sim \mathcal{N}(0, \sigma^2_b), \qquad \varepsilon_{ij} \sim \mathcal{N}(0, \sigma^2),$

can be used to describe the body length of beetle $$j$$ from population $$i$$. The repeatability is then calculated as $$R = \sigma^2_b / (\sigma^2_b + \sigma^2)$$. Below we fit this model using lmer():

data("BeetlesBody", package = "rptR")
(beetle_mod <- lmer(BodyL ~ (1 | Population), data = BeetlesBody))
Linear mixed model fit by REML ['lmerMod']
Formula: BodyL ~ (1 | Population)
Data: BeetlesBody
REML criterion at convergence: 3893.268
Random effects:
Groups     Name        Std.Dev.
Population (Intercept) 1.173
Residual               1.798
Number of obs: 960, groups:  Population, 12
Fixed Effects:
(Intercept)
14.08

To construct a bootstrap confidence interval for the repeatability, we first must write a function to calculate it from the fitted model. Below we write a one-off function for this model to demonstrate a “typical” workflow rather than trying to be overly general.

repeatability <- function(object) {
vc <- as.data.frame(VarCorr(object))
vc$vcov[1] / (sum(vc$vcov))
}

The original estimate of repeatability can then be quickly calculated:

repeatability(beetle_mod)
[1] 0.2985548

To construct a bootstrap confidence interval for the repeatability, we run the desired bootstrap procedure, specifying .f = repeatability and then pass the results to confint().

(beetle_boot <- bootstrap(beetle_mod, .f = repeatability, type = "parametric", B = 2000))
Bootstrap type: parametric

Number of resamples: 2000

observed  rep.mean         se        bias
1 0.2985548 0.2874126 0.08787679 -0.01114219

There were 0 messages, 0 warnings, and 0 errors.
(beetle_ci <- confint(beetle_boot, type = "basic"))
# A tibble: 1 × 6
term  estimate lower upper type  level
<chr>    <dbl> <dbl> <dbl> <chr> <dbl>
1 ""       0.299 0.139 0.471 basic  0.95

Notice that the column of beetle_ci is an empty character string since we did not have return a named vector.

Alternatively, we can plot() the results, as shown in Figure 1. The plot method for lmeresamp objects uses stat_halfeye() from ggdist (Kay 2021) to render a density plot with associated 66% and 95% percentile intervals.

plot(beetle_boot, .width = c(.5, .9)) +
ggplot2::labs(
title = "Bootstrap repeatabilities",
y = "density",
x = "repeatability"
)

### Bootstrap tests for a single parameter

While lmeresampler was designed with a focus on estimation, the bootstrap functions can be used to conduct bootstrap tests on individual parameters. For example, returning to the JSP example, we might be interested in generating approximate $$p$$-values for the fixed effects:

summary(jsp_mod)$coefficients  Estimate Std. Error t value (Intercept) 14.1577509 0.73344165 19.303173 mathAge8 0.6388895 0.02544478 25.108865 genderM -0.3571922 0.34009284 -1.050279 classnonmanual 0.7200815 0.38696812 1.860829 To generate a bootstrap $$p$$-value for a fixed effect, we must generate $$B$$ bootstrap resamples under the reduced model (i.e., the null hypothesis), refit the full model, and then calculate the $$t$$-statistic for each resample, denoted $$t^*_i$$. The bootstrap $$p$$-value is then calculated as $$(n_{\rm extreme} + 1) / (B + 1)$$ . Using the bootstrap() function, you can implement this procedure for a single parameter. For example, if we wish to calculate the bootstrap $$p$$-value for the class fixed effect we first generate $$B=1000$$ bootstrap samples from the reduced model without class. In this example, we use the Wild bootstrap to illustrate that we are not restricted to a parametric bootstrap. reduced_model <- update(jsp_mod, . ~ . - class) reduced_boot <- bootstrap(reduced_model, type = "wild", B = 1000, hccme = "hc2", aux.dist = "mammen", .refit = FALSE) Next, we refit the full model, jsp_mod, to each simulation and extract the $$t$$-statistic for the class variable. Note that the function extract_t() extracts the specified $$t$$-statistic from the coefficient table from model summary. extract_t <- function(model, term) { coef(summary(model))[term, "t value"] } tstats <- purrr::map_dbl( reduced_boot, ~refit(jsp_mod, .x) %>% extract_t(., term = "classnonmanual") ) With the bootstrap $$t$$-statistics in hand, we can approximate the $$p$$-value using basic logical and arithmetic operators (sum(abs(tstats) >= extract_t(jsp_mod)) + 1) / (1000 + 1) [1] 0.2687313 While the above process is not particularly difficult to implement using the tools provided by lmeresampler, things get tedious if multiple parameters are of interest in this summary table. To help reduce this burden on the user, we have provided the bootstrap_pvals() function that will add bootstrap $$p$$-values for each term in the coefficient summary table. For jsp_mod, this is achieved in the below code chunk: bootstrap_pvals(jsp_mod, type = "wild", B = 1000, hccme = "hc2", aux.dist = "mammen") Bootstrap type: wild Number of resamples: 1000 # A tibble: 4 × 5 term Estimate Std. Error t value p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) 14.2 0.733 19.3 0.000999 2 mathAge8 0.639 0.0254 25.1 0.000999 3 genderM -0.357 0.340 -1.05 0.273 4 classnonmanual 0.720 0.387 1.86 0.0659  It’s important to note that running bootstraps for each term in the model is computationally demanding. To speed up the computation, you can run the command in parallel, as we discuss below in Bootstrapping in parallel. ### Model comparison The bootstrap can be useful during model selection. For example, if you are comparing a full and reduced model where the reduced model has fewer random effects, a 50:50 mixture of $$\chi^2$$ distributions is often used ; however, Pinheiro and Bates (2000) point out that this approximation is not always optimal. In this example, we explore the Machine data set discussed by Pinheiro and Bates (2000), which consists of productivity scores for six workers on three brands of machine. This data set can be loaded from nlme: data("Machines", package = "nlme") Pinheiro and Bates (2000) consider two LME models for these data. The first model has a fixed effect for the machine and a random intercept for the worker. reduced_mod <- lmer(score ~ Machine + (1 | Worker), data = Machines, REML = FALSE) The second model has the same fixed effects structure, but adds an additional random effect for the machine within the worker. full_mod <- lmer(score ~ Machine + (1 | Worker/Machine), data = Machines, REML = FALSE) Pinheiro and Bates (2000) note that the approximate null distribution given by $$0.5\chi^2_0 + 0.5 \chi^2_1$$ is not successful when the models are fit via maximum likelihood, and that the mixture is closer to $$0.65\chi^2_0 + 0.35 \chi^2_1$$. Instead of relying on the conventional approximation, a bootstrap test can be conducted using bootstrap() to simulate the responses from the reduced model. To conduct this bootstrap test, we first extract the observed statistic obtained via anova() and then generate $$B=1000$$ bootstrap responses from the reduced model, fm1_machine. Recall that specifying .refit = FALSE returns a data frame of the simulated responses. Here, we use a residual bootstrap for illustration. observed <- anova(full_mod, reduced_mod)$Chisq[2]

reduced_boot <- bootstrap(reduced_mod, type = "residual", B = 1000, .refit = FALSE)

Next, we must fit both the full and reduced models to the bootstrap responses and calculate the test statistic. The user-written compare_models() function performs this refit and calculation for given models and bootstrap responses. The control argument for the full model was set to reduce the number of convergence warnings, since the null model had a variance component of 0 for machines within workers, so we expect warnings as we fit an expanded model.

compare_models <- function(full, reduced, newdata) {
full_mod <- refit(full, newdata,
control = lmerControl(check.conv.singular = "ignore",
reduced_mod <- refit(reduced, newdata)
anova(full_mod, reduced_mod)$Chisq[2] } chisq_stats <- purrr::map_dbl(reduced_boot, ~compare_models(full_mod, reduced_mod, newdata = .x)) With the test statistics in hand, we can quickly calculate the $$p$$-value (sum(chisq_stats >= observed) + 1) / (1000 + 1) [1] 0.000999001 ### Simulation-based model diagnostics Our final example illustrates how the bootstrap() function can be used for model diagnosis. Loy et al. (2017) propose using the lineup protocol to diagnose LME models, since artificial structures often appear in conventional residual plots for this model class that are not indicative of a model deficiency. In this example, we consider the Dialyzer data set provided by nlme. The data arise from a study characterizing the water transportation characteristics of 20 high flux membrane dialyzers, which were introduced to reduce the time a patient spends on hemodialysis . The dialyzers were studied in vitro using bovine blood at flow rates of either 200 or 300 ml/min. The study measured the the ultrafiltration rate (ml/hr) at even transmembrane pressures (in mmHg). Pinheiro and Bates (2000) discuss modeling these data. Here, we explore how to create a lineup of residual plots to investigate the adequacy of the initial homoscedastic LME model fit by Pinheiro and Bates (2000). library(nlme) dialyzer_mod <- lme( rate ~ (pressure + I(pressure^2) + I(pressure^3) + I(pressure^4)) * QB, data = Dialyzer, random = ~ pressure + I(pressure^2) ) Pinheiro and Bates (2000) construct a residual plot of the conditional residuals plotted against the transmembrane pressure to explore the adequacy of the fitted model (Figure 2). There appears to be increasing spread of the conditional residuals, which would indicate that the homoscedastic model is not sufficient. To check if this pattern is actually indicative of a problem, we construct a lineup of residual plots. To do this, we must generate data from a number of properly specified models, say 19, to serve as decoy residual plots. Then, we create a faceted set of residual plots where the observed residual plot (Figure 2) is randomly assigned to a facet. To generate the residuals from properly specified models, we use the parametric bootstrap via bootstrap() and extract a data frame containing the residuals from each bootstrap sample using hlm_resid() from HLMdiag . set.seed(1234) library(HLMdiag) sim_resids <- bootstrap(dialyzer_mod, .f = hlm_resid, type = "parametric", B = 19) The simulated residuals are stored in the replicates element of the sim_resids list. sim_resids$replicates is a tibble containing the 19 bootstrap samples, with the the replicate number stored in the .n column.

dplyr::glimpse(sim_resids$replicates) Rows: 2,660 Columns: 15$ id              <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1…
$rate <dbl> -0.8579226, 17.0489029, 34.3658738, 44.84955…$ pressure        <dbl> 0.240, 0.505, 0.995, 1.485, 2.020, 2.495, 2.…
$I(pressure^2) <I<dbl>> 0.0576, 0.255025, 0.990025, 2.205225, …$ I(pressure^3) <I<dbl>>     0.013824,  0.128787625,  0.985074875,…
$I(pressure^4) <I<dbl>> 0.00331776, 0.065037...., 0.980149....,…$ QB              <fct> 200, 200, 200, 200, 200, 200, 200, 200, 200,…
$Subject <ord> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3,…$ .resid          <dbl> -1.19758632, 0.65600425, -0.90104247, 1.3913…
$.fitted <dbl> 0.3396637, 16.3928987, 35.2669162, 43.458244…$ .ls.resid       <dbl> -0.51757746, 1.23968783, -1.32231163, 0.5907…
$.ls.fitted <dbl> -0.3403452, 15.8092151, 35.6881854, 44.25884…$ .mar.resid      <dbl> -2.1236410, -0.3100360, -2.1298323, -0.34531…
$.mar.fitted <dbl> 1.265718, 17.358939, 36.495706, 45.194867, 4…$ .n              <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…

Now, we use the lineup() function from nullabor to generate the lineup data. lineup() will randomly insert the observed (true) data into the samples data, “encrypt” the position of the observed data, and print a message that you can later decrypt in the console.

library(nullabor)
lineup_data <- lineup(true = hlm_resid(dialyzer_mod), n = 19, samples = sim_resids$replicates) dplyr::glimpse(lineup_data) Rows: 2,800 Columns: 15$ id              <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1…
$rate <dbl> -0.8579226, 17.0489029, 34.3658738, 44.84955…$ pressure        <dbl> 0.240, 0.505, 0.995, 1.485, 2.020, 2.495, 2.…
$I(pressure^2) <I<dbl>> 0.0576, 0.255025, 0.990025, 2.205225, …$ I(pressure^3) <I<dbl>>     0.013824,  0.128787625,  0.985074875,…
$I(pressure^4) <I<dbl>> 0.00331776, 0.065037...., 0.980149....,…$ QB              <fct> 200, 200, 200, 200, 200, 200, 200, 200, 200,…
$Subject <ord> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3,…$ .resid          <dbl> -1.19758632, 0.65600425, -0.90104247, 1.3913…
$.fitted <dbl> 0.3396637, 16.3928987, 35.2669162, 43.458244…$ .ls.resid       <dbl> -0.51757746, 1.23968783, -1.32231163, 0.5907…
$.ls.fitted <dbl> -0.3403452, 15.8092151, 35.6881854, 44.25884…$ .mar.resid      <dbl> -2.1236410, -0.3100360, -2.1298323, -0.34531…
$.mar.fitted <dbl> 1.265718, 17.358939, 36.495706, 45.194867, 4…$ .sample         <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…

With the lineup data in hand, we can create a lineup of residual plots using facet_wrap():

ggplot(lineup_data, aes(x = pressure, y = .resid)) +
geom_hline(yintercept = 0, color = "gray60") +
geom_point(shape = 1) +
facet_wrap(~.sample) +
theme_bw() +
labs(x = "Transmembrane pressure (dmHg)", y = "Residuals (ml/hr)")


text:
# When is it worth using Hexadecimal over Decimal?

How many bits does a number need to have, to make it worth using hexadecimal notation of the form 0xHHH...H over its decimal representation of the form DDD...D?

Here, "being worth" is defined as "giving a shorter average length for uniformly distributed values in the range $0$ to $2^N-1$", where $N$ is the number of bits. Negative numbers are not considered.

For example, clearly for all $N \ge 64$ it's worth using hexadecimal over decimal, because 18446744073709551615 has length 20, while 0xFFFFFFFFFFFFFFFF has length 18, and there are many more decimal numbers within the range that are longer than their hexadecimal counterparts than shorter, making it pretty obvious that the decimal lengths in the range where the decimal numbers are shorter don't compensate for the hex gain with greater lengths.

Extra credit for giving solutions to the same problem but with these representations for hexadecimal:

• HHH...Hh (e.g. E5F3h, 3FEh)
• DHHH...Hh (i.e. add leading zero if it starts with a letter, e.g. 0E5F3h, 3FEh)

For $n$-bit numbers, $2^n$ times the average number of base-$b$ digits required -- equivalently, the total number of digits required to write down all the $n$-bit numbers -- equals (probability of needing at least 1 digit) + (probability of needing at least 2 digits) + etc., which equals $\sum_{k>0} (2^n-b^{k-1})$ with the provisos that (1) when $k=0$ we take $b^{k-1}=0$ and (2) when the summand is negative we replace it with zero. (Each summand is then the number of $n$-bit numbers that are not representable using fewer than $k$ base-$b$ digits.)

This is of course evaluable in closed form without much pain but I am lazy and just wrote a program to do it. The result is that (with the extra two characters for 0x) hex wins once

$n>=40$ (where $2^n$ times the average number of characters is 13182540049978 for decimal and 13120838758128 for hex).

With the more economical (but to my mind uglier) suffix-h notation, hex wins once

$n>=20$ (where $2^n$ times the average number of characters is 6228922 for decimal and 6221552 for hex).

The third notation (with a leading 0 when, and only when, the number's hex representation starts with a letter) is a little hairier. We need at least one character always; at least two when the number is at least 10; at least three when it is at least 160; etc. Again, there is certainly a not-too-painful closed form; again, I am lazy; if my code isn't made of bugs then hex wins once

$n>=24$, at which point $2^n$ times the average number of characters is 123106618 for decimal and 123032918 for hex.

• ̶h̶a̶i̶r̶i̶e̶r̶.̶.̶.̶l̶a̶t̶e̶r̶ ̶ too long to fit in the margin. (fixed that for you.) – Rubio Dec 30 '16 at 21:28
• I found a bit more margin space. – Gareth McCaughan Dec 30 '16 at 21:35
• A perfectly fe̶ormatted comment, wasted. Alas. :) – Rubio Dec 30 '16 at 21:45
• These were my numbers as well. Congrats! – Pedro Gimeno Dec 31 '16 at 1:35

The solution for the first problem is

41 bits.

Reasoning:

(counting from 1, as that makes the math easier)

DECIMAL
For the first 9 * 1 = 9 numbers, you need 9 * 1 = 9 digits.
For the next 9 * 10 = 90 numbers, you need 90 * 2 = 180 digits.
For the next 9 * 100 = 900 numbers, you need 900 * 3 = 2700 digits.
In general, for all numbers with $n$ digits (those between $10^{n-1}$ and $10^n$) you need $9n \cdot 10^{n-1}$ digits.
Summing this (it's a harmonic series) gives a formula for the number of digits needed to write all numbers with up to $n$ digits: $(n-\frac{1}{9}) \cdot 10^n + \frac{1}{9}$.
The average for all $10^n$ numbers is thus approximately $n-\frac{1}{9}$.

For the first 15 numbers, you need 15 * 3 = 45 digits.
For the next 15 * 16 numbers, you need 15 * 16 * 4 = 960 digits.
For the next 15 * 256 numbers, you need 15 * 256 * 5 = 19200 digits.
In general, for all numbers with $x$ digits (those between $16^{x-1}$ and $16^x$) you need $15(x+2) \cdot 16^{x-1}$ digits.
Summing this gives a formula for the number of digits needed to write all numbers with up to $x$ digits: $(x+2-\frac{1}{15}) \cdot 16^x - \frac{29}{15}$.
The average for all $16^x$ numbers is thus approximately $x+2-\frac{1}{15}$.

While these formulas are only valid for integer values of $n$ and $x$ (so for numbers up to $10^n$ for the total number of decimal digits, and up to $16^x$ for the hexadecimal case), it's interesting to see where they intersect, because from that point on, the hexadecimal numbers are bound to be shorter on average. We know $10^n = 16^x$, so $x = n \cdot \log_{16}10$. Plugging this into the equation $n - \frac{1}{9} = x + 2 - \frac{1}{15}$ gives:

$n - \frac{1}{9} = n \cdot \log_{16}10 + 2 - \frac{1}{15}$

$n (1 - \log_{16}10) = \frac{29}{15} + \frac{1}{9} = \frac{92}{45}$

$n = \frac{92}{45 (1 - \log_{16}10)}$

The number of bits is a factor $\log_2 10$ higher than $n$, so $\frac{92 \log_2 10}{45 (1 - \log_{16}10)} = \frac{92}{45(\log_{10}2 - \frac14)} \approx 40.06$. So while at 40 bits, decimal is shorter, at 41 bits, we expect the hexadecimal notation to win.

For the second problem, the main equation becomes

$n - \frac{1}{9} = x + 1 - \frac{1}{15}$

which leads to a number of bits equal to

$\frac{47}{45(\log_{10}2 - \frac14)} \approx 20.47$

so the cutoff is here already 21 bits.

• While your analysis is certainly on the right track, the solution is not correct. I can only say that you were close without giving too many spoilers. (Just to be sure: no, it's not me being picky about "41 or more", it's that the cutoff is not at 41 for the first problem, nor at 21 for the second problem). – Pedro Gimeno Dec 30 '16 at 19:33

I would have to say

16 bits or 15 bits technically

This is because

at 14 bits 0x3FFE = 16,382. At this point the DEC is larger than the HEX representation. So at this point if the binary was increased by another bit, more than 50% of the conversion from DEC to HEX would be shorter than the DEC representation. Thus, making it "worth it".

format: plain


text:
RUS  ENG JOURNALS   PEOPLE   ORGANISATIONS   CONFERENCES   SEMINARS   VIDEO LIBRARY   PACKAGE AMSBIB
 General information Latest issue Archive Impact factor Search papers Search references RSS Latest issue Current issues Archive issues What is RSS

 Vladikavkaz. Mat. Zh.: Year: Volume: Issue: Page: Find

 Vladikavkaz. Mat. Zh., 2018, Volume 20, Number 4, Pages 35–42 (Mi vmj674)

$L_p-L_q$-estimates for potential-type operators with oscillating kernels

M. N. Gurova, V. A. Noginb

a ISC «Lyceum of Classical Elite Education», 166A Pushkinskaya str., Rostov-on-don 344006, Russia
b Southern Mathematical Institute VSC RAS, 22 Marcus st., Vladikavkaz 362027, Russia

Abstract: We consider a class of multidimensional potential-type operators whose kernels are oscillating at infinity. The characteristics of these operators are from a wide class of functions including the product of a homogeneous function infinitely differentiable in $\Bbb R^n\setminus\{0\}$ and any function from $C^{m,\gamma}(\dot{R}^1_{+})$. We describe convex sets in the $(1/p;1/q)$-plane for which these operators are bounded from $L_p$ into $L_q$ and indicate the domains where they are not bounded. In some cases, the accuracy of the estimates obtained is proved. In particular, necessary and sufficient conditions for the boundedness of the operators under considered in $L_p$ are obtained. Currently, there is a number of papers on $L_p-L_q$-estimates for convolution operators with oscillating kernels, in particular, for the Bochner–Riesz operators and acoustic potentials arising in various problems of analysis and mathematical physics. These papers cover kernels containing only the radial characteristic $b(r)$, which stabilized at infinity as a Hëlder function. Due to this property, the derivation of estimates for the indicated operators was reduced to the case of an operator with the characteristic $b(r)\equiv1$. Such a reduction is impossible when the Riesz potential kernel contains a homogeneous characteristic $a(t')$. To receive the results we use new method which based on special representation of the symbols multidimensional potential-type operators. To these representations of the symbols we apply the technique of Fourier-multipliers, which degenerate or have singularities on the unit sphere in $\mathbb{R}^n$.

Key words: potential-type operators, oscillating kernel, $L_p-L_q$-estimates, $\mathcal{L}$-characteristics.

DOI: https://doi.org/10.23671/VNC.2018.4.23385

Full text: PDF file (281 kB)
References: PDF file   HTML file

UDC: 517.983.2
MSC: 46E35, 26A33

Citation: M. N. Gurov, V. A. Nogin, “$L_p-L_q$-estimates for potential-type operators with oscillating kernels”, Vladikavkaz. Mat. Zh., 20:4 (2018), 35–42

Citation in format AMSBIB
\Bibitem{GurNog18} \by M.~N.~Gurov, V.~A.~Nogin \paper $L_p-L_q$-estimates for potential-type operators with oscillating kernels \jour Vladikavkaz. Mat. Zh. \yr 2018 \vol 20 \issue 4 \pages 35--42 \mathnet{http://mi.mathnet.ru/vmj674} \crossref{https://doi.org/10.23671/VNC.2018.4.23385}

• http://mi.mathnet.ru/eng/vmj674
• http://mi.mathnet.ru/eng/vmj/v20/i4/p35

 SHARE:

Citing articles on Google Scholar: Russian citations, English citations
Related articles on Google Scholar: Russian articles, English articles
•  Number of views: This page: 48 Full text: 25 References: 6


text:
Lemma 15.9.7. Let $A$ be a ring, let $I \subset A$ be an ideal. Let $f \in A[x]$ be a monic polynomial. Let $\overline{f} = \overline{g} \overline{h}$ be a factorization of $f$ in $A/I[x]$ and assume that $\overline{g}$, $\overline{h}$ generate the unit ideal in $A/I[x]$. Then there exists an étale ring map $A \to A'$ which induces an isomorphism $A/I \to A'/IA'$ and a factorization $f = g' h'$ in $A'[x]$ lifting the given factorization over $A/I$.

Proof. Say $f = x^ d + a_1 x^{d - 1} + \ldots + a_ d$ has degree $d$. Write $\overline{g} = \sum \overline{b}_ j x^ j$ and $\overline{h} = \sum \overline{c}_ j x^ j$. Then we see that $1 = \sum \overline{b}_ j \overline{c}_{d - j}$. It follows that $\mathop{\mathrm{Spec}}(A/I)$ is covered by the standard opens $D(\overline{b}_ j \overline{c}_{d - j})$. However, each point $\mathfrak p$ of $\mathop{\mathrm{Spec}}(A/I)$ is contained in at most one of these as by looking at the induced factorization of $f$ over the field $\kappa (\mathfrak p)$ we see that $\deg (\overline{g} \bmod \mathfrak p) + \deg (\overline{h} \bmod \mathfrak p) = d$. Hence our open covering is a disjoint open covering. Applying Lemma 15.9.3 (and replacing $A$ by $A'$) we see that we may assume there is a corresponding disjoint open covering of $\mathop{\mathrm{Spec}}(A)$. This disjoint open covering corresponds to a product decomposition of $A$, see Algebra, Lemma 10.24.3. It follows that

$A = A_0 \times \ldots \times A_ d, \quad I = I_0 \times \ldots \times I_ d,$

where the image of $\overline{g}$, resp. $\overline{h}$ in $A_ j/I_ j$ has degree $j$, resp. $d - j$ with invertible leading coefficient. Clearly, it suffices to prove the result for each factor $A_ j$ separatedly. Hence the lemma follows from Lemma 15.9.6. $\square$

In your comment you can use Markdown and LaTeX style mathematics (enclose it like $\pi$). A preview option is available if you wish to see how it works out (just click on the eye in the toolbar).


text:
NIPS 2018
Sun Dec 2nd through Sat the 8th, 2018 at Palais des Congrès de Montréal
Paper ID: 5289 Generalized Cross Entropy Loss for Training Deep Neural Networks with Noisy Labels

Reviewer 1

I acknowledge that I read the author's response, and I feel that the revised version of the manuscript will be even stronger. As a result, I am raising my score to an 8. This paper is definitely in the Top 50%. ----- Summary: This paper shows that using the negative Box-Cox transformation as a loss function ($L_q$ or Lq loss) generalizes both categorical cross-entropy (CCE) and mean absolute error (MAE) losses and offers an elegant way to trade off their respective virtues when training on noisy labels. The key insight comes from analyzing the loss function gradients: they are equivalent, except that CCE includes a term that implicitly assigns higher weights to incorrect predictions. This makes training with CCE faster than with MAE but also makes it more susceptible to overfitting label noise. Like CCE, the gradient of Lq loss yields a weighting term but with an exponent parameter that we can choose. When q=0, we get CCE, and when q=1, the weighting term disappears, which is equivalent to MAE. The paper shows that a variant of a known risk bound for MAE under uniform label noise applies to Lq loss as q approaches 1. Experimental results are noticeably strong: Lq consistently outperforms CCE and MAE both and is competitive with several alternative strong baselines. My current feeling is that this paper is a pretty clear accept, but I look forward to reading the other reviewers' comments and authors' responses. Quality: This submission appears sound to me -- the bounds quoted from the Ghosh work are correct and the gradients for the CCE, MAE, and Lq losses are correct. The analyses of Lq loss (lines 147-153) and truncated Lq loss (lines 165-178) appear correct to me, but I can't guarantee their correctness. For what it's worth, I'd recommend including a more detailed proof in either an appendix or an extended version of the paper on arXiv. The experiments provide convincing evidence that Lq loss is superior to CCE and MAE in the presence of uniform and class conditional label noise. Lq loss is competitive with Sukhbaatar and Fergus algorithm, even when it has access to the noise rate. I do have two questions for the authors: (1) How do these results apply to binary classification (presumably they do since binary classification is a special case)? However, Ghosh showed in an earlier paper [4] that there for binary classification, there are other noise robust losses, e.g., ramp loss, that work quite well and might provide a stronger baseline than MAE. (2) Ghosh's results apply to more general label noise scenarios, including non-uniform noise. Does the analysis here also apply to this more challenging (and realistic) scenario? The experiments do have one weakness: there has been a lot of research on noisy labels in the last few years, so I think that a comparison with additional baselines is warranted. Examples include [1][2][3]. Frankly, I'm not sure ANY baseline comparisions is strictly necessary since the paper's central claim is that Lq loss is superior to CCE and MAE -- not that it's superior to all noisy label learning techniques. However, the inclusion of the "Forward" method opens the door to the "why didn't you include this other baselinee" critique. Clarity: This paper is clearly written for the most part, though as noted, the theoretical analysis of Lq loss (lines 147-153, 165-178) would benefit from greater detail, even in a separate document. The analysis of the gradients could use further clarification -- this sentence on lines 122-23 is particularly challenging to understand: "...less congruent with provided labels are implicitly weighed more in the gradient update." Staring at Equation 5, it's not immediately clear to me what this means. The discussion of truncated Lq loss could benefit from further discourse. For example, what is the motivation behind the introduction of the threshold k in Equation 9? How can we intuitively understand what it's doing, and how does it connect to Equation 11? Also, the details of the pruning step aren't entirely clear: since neural net training typically involves unconstrained SGD, how do we optimize v_i bounded between 0 and 1? Originality: This paper has a marginal amount of originality -- it leans heavily on the Ghosh and Ferrari/Yang papers, synthesizing their respective ideas (risk bounds for noise robust symmetric loss functions and Lq loss), rather than proposing entirely new ones. Nonetheless, I think it demonstrates innovtative thinking, especially when looking at the gradients of the loss functions and considering how they will impact the dynamics of neural net learning (wereas Ghosh focuses primarily on the statistical properties). Significance: I'm really uncertain about the significance of the paper. But my takeaway from the paper is "use Lq loss instead of CCE and MAE, in most cases." Given that Lq loss appears pretty easy to implement and differentiate, there seems to be no reason to replace CCE with Lq in any experiment where there may be nontrivial amounts of label noise. References [1] Reed, et al. Training Deep Neural Networks on Noisy Labels with Bootstrapping. ICLR 2015 Workshop Track. [2] Natarajan, et al. Learning with Noisy Labels Nagarajan. ICML 2013. [3] Ma, et al. Dimensionality-Driven Learning with Noisy Labels. ICML 2018. [4] Ghosh, et al. Making Risk Minimization Tolerant to Label Noise. Neurocomputing 2015.

Reviewer 2

Summary A method for learning a classifier robust to label noise is proposed. A new class of robust loss functions is considered to tackle the problem. The contribution is based a truncated and Box-Cox transformed version of the mean absolute loss, which was shown to be theoretically robust to uniform label noise. The proposed method (unlike most previous work) works for both closed and open-set noisy labels. Several experiments support the approach. Detailed comments I liked the insight given in formula (5), which can be considered the starting point of the paper argument. The theoretically-backed derivation of the truncated loss is nice -- in contrast with several heuristics proposed in this area. The straightforward application to open set noise is also interesting. I think this is a solid contribution and deserves publication. A weakness is the need of two new hyperparameters, i.e. k and q. Although, it seems that a default choice worked well throughout all experiments. It would be helpful to know more about the additional computational cost due to the pruning step. How slower is training? Moreover, I think that some details about the proposed alternating optimisation are omitted. How do you guarantee that the v_i are bounded in [0,1] by using SGD? Do you use a proximal method? Minors * line 113: why the fact that CCE is unbounded matter here? And isn’t MAE unbounded as well? * line 143: could you at least show on a footnote the derivation of CCE from Lq when q -> 0 ?

Reviewer 3

Update: after reading the authors' response, I have increased my overall score to "Top 50% of accepted NIPS papers". This paper proposes a new loss function to deal with label noise, which overcomes the drawbacks of standard loss functions (non-robustness to label errors) as well as of existing robust loss functions (they lead to very slow training and a drop in accuracy). Moreover, the proposed family of loss functions contains the cross-entropy loss and the mean absolute error loss as particular (extreme) cases. The paper is well written (although there are many missing determiners, such as "the" or "a") and the ideas are well presented, motivated, and justified. The experiments clearly illustrate the behavior of the proposed loss, confirming that it performs as expected. Important questions remain about the usefulness and applicability of this type of robust loss function. Namely, in the presence of some dataset, how to infer its level of noise (if any) and how to adjust the parameters (q and k) of the loss? If the dataset is noisy (contains wrong labels), arguably cross-validation will not be very informative.


text:
HOSTING A TOTAL OF 318 FORMULAS WITH CALCULATORS

## Volume of cube

The number of cubic units that will exactly fill a cube.The volume of a cube is found by multiplying the length of any edge by itself twice

## ${s}^{3}$

here s=side.

ENTER THE VARIABLES TO BE USED IN THE FORMULA

Similar formulas which you may find interesting.