Skip to contents

The HCbeta algorithm and its consistency

Source: vignettes/hcinfer-hcbeta-algorithm.Rmd
hcinfer-hcbeta-algorithm.Rmd

This vignette presents the HCbeta adjustment factor as an explicit algorithm and records the argument that it yields a consistent covariance estimator. It complements vignette("hcinfer-methodology", package = "hcinfer"), which places HCbeta within the broader family of HC estimators, and vignette("hcinfer-hcbeta", package = "hcinfer"), which is usage oriented.

Motivation

The classical leverage corrections write the adjustment factor as a power of the leverage complement u_t = 1 - h_t. Because u_t \in [0, 1], it can be read as the value at 1 - h_t of the cumulative distribution function of a Uniform(0, 1) random variable. HCbeta keeps this reading but replaces the uniform cdf with a Beta cdf evaluated at a truncated complement w_t. The two shape parameters give the correction curve an extra degree of flexibility, so the adjustment can adapt to the observed leverage configuration instead of following a fixed exponent.

This flexibility aims to temper the overshooting that HC3, HC4, and HC4m can display under strong leverage, where small values of 1 - h_t raised to a fixed negative exponent produce very large factors. HCbeta calibrates the Beta shape parameters to the first two sample moments of w_t, so the resulting factor reflects the actual spread of the leverages rather than a worst case exponent.

HCbeta does not assume that the leverage complements follow a Beta distribution. The Beta cdf is a calibration device whose shape parameters are matched to the sample mean and variance of the truncated complements, and the construction would be unchanged if the leverages arose from another mechanism. The Beta family is used because its support is [0, 1] and its two parameters are enough to bend the uniform baseline in either direction.

The algorithm

Let X be the n \times p design matrix of an OLS fit with residuals \hat e_t and leverages h_t = x_t'(X'X)^{-1}x_t. HCbeta computes the adjustment factor g_t in four steps.

  1. Truncated leverage complements. Form u_t = 1 - h_t and truncate it to w_t = \max\{\text{lower}, \min(u_t, \text{upper})\}, with defaults \text{lower} = 0.01 and \text{upper} = 0.99. The lower bound prevents the adjustment from exploding when h_t approaches one, and the upper bound keeps the values away from the boundary when leverages are small.

  2. Method of moments. Estimate the mean and variance of the truncated complements, \hat\mu = n^{-1}\sum_{t=1}^{n} w_t and s_w^2 = (n - 1)^{-1}\sum_{t=1}^{n}(w_t - \hat\mu)^2, and map them to Beta shape parameters through the mean precision parameterization, \hat\phi = \hat\mu(1 - \hat\mu)/s_w^2 - 1, \hat a = \hat\mu\hat\phi, and \hat b = (1 - \hat\mu)\hat\phi.

  3. Shrinkage with a fixed floor and caps. Shrink the estimates toward the uniform case a = b = 1 using \zeta = n/(n + 50), then set \tilde a = \min\{\max\{(1 - \zeta) + \zeta\hat a,\, \varepsilon\},\, A_{\max}\} and \tilde b = \min\{\max\{(1 - \zeta) + \zeta\hat b,\, \varepsilon\},\, B_{\max}\}. The floor is fixed at \varepsilon = 0.01, while the defaults A_{\max} = B_{\max} = 10000 are set through a_max and b_max, each of which accepts finite values in the inclusive range [50, 25000]. The floor is applied after shrinkage and before the caps. It is distinct from lower, which truncates w_t, and it is not a method argument.

  4. Adjustment factor and sandwich. Combine the Beta cdf with a decaying exponent, g_t = \frac{n}{n - p}\left\{1/F_B(w_t; \tilde a, \tilde b)\right\}^{c_1/n^{c_2}}, with defaults c_1 = 7 and c_2 = 0.75. The covariance estimator is the sandwich \widehat\Psi_{HC\beta} = (X'X)^{-1} X'\widehat\Omega_\beta X (X'X)^{-1} with \widehat\Omega_\beta = \operatorname{diag}(\hat e_t^2 g_t).

The degrees of freedom factor n/(n - p) enters g_t in the same way as in HC1. When \tilde a = \tilde b = 1 the Beta cdf returns w_t, and with c_1 = 0 the factor reduces to n/(n - p), so HCbeta coincides with HC1. This identity is a convenient implementation check and is exercised in the package tests.

Consistency

The consistency of HCbeta follows from the standard argument for HC0 given by White (1980) for fixed regressors, extended to stochastic regressors in the usual way. The role of the HCbeta construction is only to show that its adjustment factor does not disturb that argument asymptotically. The reasoning below is presented at an accessible level and does not introduce a new theorem.

Two features keep the adjustment under control. First, the truncation to [0.01, 0.99] and the bounds on the shape parameters place \tilde a and \tilde b in a compact interval [\varepsilon, \max(A_{\max}, B_{\max})] with the fixed value \varepsilon = 0.01. On that compact set, the Beta cdf evaluated at w_t is bounded away from zero, so that \inf_t F_B(w_t; \tilde a, \tilde b) \geq \delta for some \delta > 0, and therefore \sup_t |\log F_B(w_t; \tilde a, \tilde b)| \leq -\log\delta < \infty. The caps make the parameter space compact, keeping the log adjustments uniformly bounded.

The caps are essential when the leverage variance vanishes. Under the standard fixed-design condition \max_t h_t \to 0, the complements eventually reach the common upper truncation, so s_w^2 = 0 and the uncapped method-of-moments shape estimates diverge. HCbeta then uses A_{\max} and B_{\max} as the limiting shape parameters instead of aborting. This keeps the Beta cdf and g_t well-defined throughout the sequence; with the default caps, the common cdf is numerically one and g_t equals the HC1 scale n/(n - p).

Second, a first-order Taylor expansion of g_t around one gives \sup_t |g_t - 1| = O(n^{-c_2}) = o(1), because the exponent c_1/n^{c_2} tends to zero while it multiplies the bounded log cdf. It follows that \widehat\Omega_\beta and the HC0 residual matrix \operatorname{diag}(\hat e_t^2) are asymptotically equivalent, so \widehat\Psi_{HC\beta} and \widehat\Psi_{HC0} share the same limit. HCbeta is asymptotically equivalent to HC0, so it inherits the same consistency, while the finite sample factor g_t aims to temper the overshooting that fixed exponent methods can produce under strong leverage.

Numerical implementation

Four safeguards keep the computation of g_t numerically stable. First, the post-shrinkage shape parameters are clamped to the fixed floor 0.01 and their respective caps. Second, the moment branch is used only when is.finite(s2_w) && s2_w > .Machine$double.eps; otherwise, the shapes are set directly to their caps. Third, the Beta cdf is evaluated with pbeta(w, a_tilde, b_tilde, log.p = TRUE), which returns \log F_B directly and avoids the \log(0) = -\infty produced when a cdf underflows. Finally, the exponent -(c_1/n^{c_2})\log F_B(w_t; \tilde a, \tilde b) is capped at 700 before exponentiation so that exp() cannot overflow to infinity.

In the public-schools model shown here, the default caps A_{\max} = B_{\max} = 10000 are far from binding.

library(hcinfer)

schools <- transform(
  PublicSchools,
  income_scaled = income / 10000,
  income_scaled_sq = (income / 10000)^2
)
fit <- lm(expenditure ~ income_scaled + income_scaled_sq, data = schools)
cov <- vcov_hc(fit, type = "hcbeta")

cov$method_params[c("a_tilde", "b_tilde", "a_max", "b_max")]
#> $a_tilde
#> [1] 3.147193
#> 
#> $b_tilde
#> [1] 0.6689698
#> 
#> $a_max
#> [1] 10000
#> 
#> $b_max
#> [1] 10000
max(cov$weights)
#> [1] 4.580723

The adjusted shape parameters are well below the caps, and the largest adjustment factor matches the value reported for this model. By construction each factor is at least the HC1 scale n/(n - p), because the capped exponent is nonnegative and the leading factor is exactly n/(n - p). The code below verifies this lower-bound property.

n <- cov$n
p <- cov$p

all(cov$weights >= n / (n - p))
#> [1] TRUE