Stating the optimization problem

Let \(q = (q_1, \ldots, q_n)\) be the resampling distribution. Here, \(n\) is the number of rows in the original dataset, and \(q_i\) is the probability that row \(i\) will be chosen for a given row of the bootstrap dataset. Let \(p = (p_1, \ldots, p_n)\) be uniform (i.e., \(p_i = 1/n\), the resampling distribution of the ordinary bootstrap).

We want to find a \(q\) that is as close to \(p\) as possible. We consider three possible measures of distance. First, the Euclidean distance is defined as: \[ D_{Euclidian}= \sum_{i = 1}^n (p_i - q_i)^2 \] Second, the Backwards Kullback-Leibler (KL) divergence (which is technically not considered a distance metric) from \(p\) to \(q\) is defined as: \[ \begin{aligned} D_{KL}(q||p) &= \sum_{i = 1}^n q_i \log \frac{q_i}{p_i} \\ &= \sum_{i = 1}^n q_i \log nq_i \end{aligned} \]

Third, the Forwards Kullback-Leibler (KL) divergence from \(q\) to \(p\) is defined as: \[ D_{KL}(p||q) = \sum_{i = 1}^n p_i \log \frac{p_i}{q_i} \\ \]

Note that minimizing \(D_{KL}(p||q)\) is the same as minimizing \(-\sum_{i = 1}^n \log q_i\).

The goal of determining the weights for the ‘tilted bootstrap’ may be formulated as a constrained optimization problem: \[ q=\underset{q \in S}{\operatorname{argmin}} {D(p,q)} \] where,

\[ S=\{q|\space q_i>0 \space for \space i \in {1...n}, \space and \space \sum_{i = 1}^n q_i=1, \space and \space X^Tq = y \space \}. \] Here \(y\) is a vector of target means for the columns.

Solving the optimization problem

In the case of the Euclidean distance, the problem is a standard “Quadratic Programming” problem which we solve using the ‘solve.QP’ function from the ‘quadprog’ package.

To solve the problem for the Backwards KL distance, we use the Lagrangian approach. Define \(A=[X|1]\) (i.e. \(A\) is the concatenation of \(X\) with a conformal column of all ones). Also, let \(b=[y^T|1]^T\) (i.e. \(b\) is the \(y\) vector with a “one” element concatenated to the end). The Lagrangian for this optimization problem is

\[ \begin{aligned} L(q, \lambda) &= \sum_{i = 1}^n q_i \log nq_i + \lambda (A'q - b) \\ &= \sum_{i = 1}^n \left [ q_i \log nq_i + \sum_{k = 1}^{K+1} \lambda_k (A_{ik} q_i - b_k) \right ] \end{aligned} \]

where

• The vector \(\lambda\) of length \(K+1\) represents the Lagrange multipliers.
• \(K\) is the number of columns of \(X\).
• \(A_{ik}\) is the element of \(A\) in row \(i\) and column \(k\).

To optimize, we start by setting the gradient of \(L\) equal to zero and solve for a locally optimal \(q\) in terms of \(\lambda\). \[ \begin{aligned} 0 &= \left . \frac{\partial L}{\partial q_i} \right |_{(q, \ \lambda)} \\ &= \log nq_i + \frac{1}{n} + A_i' \lambda \\ \log nq_i &= -\frac{1}{n} - A_i' \lambda \\ q_i &= e^{-A_i' \lambda^*} \end{aligned} \] Here, \(\lambda^*\) is scaled and recentered to absorb the constants in the equation.

Here, \(X_i\) is the \(i\)’th row of \(X\) (expressed as a column vector). Hence, the locally optimal \(q\) in terms of the vector \(\lambda^*\) is

\[q^*(\lambda^*) = \left (e^{-A_1' \lambda^*}, \ldots, e^{-A_n' \lambda^*} \right ) \] To find \(\lambda^*\), we solve the equation \(A'q^*(\lambda^*) = b\) by using the Newton-Raphson algorithm. Thus, it turns out that the backwards version of the KL distance results in what is called “exponential tilting.” Some of the elegant properties of exponential tilting are described in a different vignette.

Similarly, it may be shown for the forward version of the KL, that the locally optimal \(q\) in terms of \(\lambda\) is

\[q^*(\lambda) = \left (\frac{1}{A_1^T \lambda}, \ldots, \frac{1}{ A_n^T \lambda} \right ) \]

However, tboot’s current implementation of forward KL distance is not very stable. A warning will be issued if the user selects the forward KL distance in the tboot package.

Not full rank

If the data matrix used to determine the weights is not full rank, tweights will likely fail. Consider reducing the number of columns being constrained to deal with this issue.

Constraint not achievable

Sometimes the desired constraint may not be achievable. This often occurs when the desired target differs from the data by a large amount. If the target is not achievable, tweights will give a warning, then it will find the “closest” achievable constraint. “Closest” is measured by a weighted Euclidean distance. Thus the target (\(y\)) is replaced with: \[ y_{new}=X'q^* \] where \[ q^*=\underset{q^* \in S^*}{\operatorname{argmin}} {(X'q - y)'W(X'q - y)} \] and \(S^*={q^*|\space \sum_{i = 1}^n q_i^* = 1}\). The weights (\(W\)) is a diagonal matrix with diagonal entries \(W_{kk}=1/s^2_k\) where \(s^2_k\) is the sample variance of column \(k\) from the data matrix (\(X\)). We solve this problem using the ‘ipop’ function in the ‘kernlap’ package. Note that ‘kernlab’ is used instead of ‘quadprog’ because the problem is not easily formulated using a full rank matrix for the quadratic form in the problem.

Unable to find an optimum

The algorithm may fail to converge properly, and an error or warning will be issued.

Some samples have a very high weight

As previously mentioned, if any sample is too highly weighted, the tilted bootstrap approach is not recommended. Somtimes transformations or the removal of outliers can help remediate the issues. See Example 1 above for more discussion about this issue.

Augmented Sampling to Improved Stability

In many cases, the problems described above can be avoided by the use of the ‘Nindependent’ option in the ‘tweights’ function. We suggest that typically, if this option is used, the user set Nindependent=1. This means that one one additional “special” sample is assumed to exist in the original data. If ‘tboot’ draws this special sample, it will proceed to bootstrap each variable independently so that the variables are independent for that patient. The weights for each independent variable bootstrap are set so that (if possible) the “special” sample would average to the target. In effect, the Nindependent option makes the correlation structure slightly less dependent. In exchange for the small bias created, the user will be able to better achieve a target even when the dataset available is small or ill-conditioned. This method is used for the Bayesian marginal reconstruction method as it helps eliminate errors in corners of the parameters space that do not fit well with the data but still occasionally occur. In these cases, the independent sample may occur frequently, so the simulation defaults towards greater independence when the data is not consistent with the parameter. Here is an example of using the Nindependent option where correlation is slightly smaller with augmentation:

x1=rnorm(1000)
x2=rnorm(1000)*sqrt(.1) + x1*sqrt(.9)
dataset2=data.frame(x1=x1, x2=x2)

weights_no_augmentation <- tweights(dataset = dataset2, 
                    target = c(x1=0.2, x2=-0.2))

## ----------------------------------------------------------------
## Optimization was successful. The weights have a sampling
## distribution with means close to the attempted target:
##                      x1         x2
## Achieved Mean 0.1999999 -0.2000001
## Target Mean   0.2000000 -0.2000000
## Maximum weight was:  0.04102889 
## ----------------------------------------------------------------

weights_augmentation <- tweights(dataset = dataset2, 
                    target = c(x1=0.2, x2=-0.2), Nindependent=1)

## ----------------------------------------------------------------
## Optimization was successful. The weights have a sampling
## distribution with means close to the attempted target:
##                      x1         x2
## Achieved Mean 0.1999999 -0.2000001
## Target Mean   0.2000000 -0.2000000
## Maximum weight was:  0.0409359 
## Data augmented with 1 sample(s) with independent variables. 
## The final weight of the indpendent sample(s) was:  0.002266406 
## ----------------------------------------------------------------

boot_no_augmentation <- tboot(weights_no_augmentation, nrow = 1e5)
boot_augmentation <- tboot(weights_augmentation, nrow = 1e5)
cor(boot_no_augmentation)

##           x1        x2
## x1 1.0000000 0.9400113
## x2 0.9400113 1.0000000

cor(boot_augmentation)

##           x1        x2
## x1 1.0000000 0.9376946
## x2 0.9376946 1.0000000