Consistency

.Rhistory

@@ -174,3 +174,7 @@ install.packages("tidyverse")
 roxygen2::roxygenise()
 pkgdown::build_site()
 pkgdown::build_site()
+roxygen2::roxygenise()
+pkgdown::build_site()
+exit()
+q()

.Rproj.user/67B916F4/pcs/source-pane.pper

@@ -1,3 +1,3 @@
 {
-    "activeTab": 0
+    "activeTab": -1
 }

.Rproj.user/67B916F4/pcs/windowlayoutstate.pper

@@ -1,14 +1,14 @@
 {
     "left": {
-        "splitterpos": 600,
-        "topwindowstate": "NORMAL",
-        "panelheight": 1042,
-        "windowheight": 1080
+        "splitterpos": 477,
+        "topwindowstate": "HIDE",
+        "panelheight": 830,
+        "windowheight": 868
     },
     "right": {
-        "splitterpos": 653,
+        "splitterpos": 520,
         "topwindowstate": "NORMAL",
-        "panelheight": 1042,
-        "windowheight": 1080
+        "panelheight": 830,
+        "windowheight": 868
     }
 }

.Rproj.user/67B916F4/persistent-state

@@ -1,6 +1,6 @@
 build-last-errors="[]"
-build-last-errors-base-dir="~/Library/CloudStorage/Dropbox/artigo_cribari_marina/hcinfer/"
-build-last-outputs="[{\"type\":0,\"output\":\"==> devtools::document(roclets = c('rd', 'collate', 'namespace'))\\n\\n\"},{\"type\":1,\"output\":\"\\u001B[1m\\u001B[22m\\u001B[36mℹ\\u001B[39m Updating \\u001B[34mhcinfer\\u001B[39m documentation\\n\"},{\"type\":1,\"output\":\"\\u001B[1m\\u001B[22m\\u001B[36mℹ\\u001B[39m Loading \\u001B[34mhcinfer\\u001B[39m\\n\"},{\"type\":1,\"output\":\"Documentation completed\\n\\n\"},{\"type\":0,\"output\":\"==> R CMD INSTALL --preclean --no-multiarch --with-keep.source hcinfer\\n\\n\"},{\"type\":1,\"output\":\"* installing to library ‘/opt/homebrew/lib/R/4.5/site-library/_build’\\n\"},{\"type\":1,\"output\":\"* installing *source* package ‘hcinfer’ ...\\n\"},{\"type\":1,\"output\":\"** this is package ‘hcinfer’ version ‘0.1.1’\\n\"},{\"type\":1,\"output\":\"** using staged installation\\n\"},{\"type\":1,\"output\":\"** R\\n\"},{\"type\":1,\"output\":\"** data\\n\"},{\"type\":1,\"output\":\"*** moving datasets to lazyload DB\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** inst\\n\"},{\"type\":1,\"output\":\"** byte-compile and prepare package for lazy loading\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** help\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"*** installing help indices\\n\"},{\"type\":1,\"output\":\"*** copying figures\\n\"},{\"type\":1,\"output\":\"** building package indices\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** installing vignettes\\n\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from temporary location\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from final location\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package keeps a record of temporary installation path\\n\"},{\"type\":1,\"output\":\"* DONE (hcinfer)\\n\"},{\"type\":1,\"output\":\"\"}]"
+build-last-errors-base-dir="~/Library/CloudStorage/Dropbox/correcoes_artigo_hc_beta/hcinfer/"
+build-last-outputs="[{\"type\":0,\"output\":\"==> devtools::document(roclets = c('rd', 'collate', 'namespace'))\\n\\n\"},{\"type\":1,\"output\":\"\\u001B[1m\\u001B[22m\\u001B[36mℹ\\u001B[39m Updating \\u001B[34mhcinfer\\u001B[39m documentation\\n\"},{\"type\":1,\"output\":\"\\u001B[1m\\u001B[22m\\u001B[36mℹ\\u001B[39m Loading \\u001B[34mhcinfer\\u001B[39m\\n\"},{\"type\":1,\"output\":\"Documentation completed\\n\\n\"},{\"type\":0,\"output\":\"==> R CMD INSTALL --preclean --no-multiarch --with-keep.source hcinfer\\n\\n\"},{\"type\":1,\"output\":\"* installing to library ‘/opt/homebrew/lib/R/4.5/site-library/_build’\\n\"},{\"type\":1,\"output\":\"* installing *source* package ‘hcinfer’ ...\\n\"},{\"type\":1,\"output\":\"** this is package ‘hcinfer’ version ‘0.1.1.9000’\\n\"},{\"type\":1,\"output\":\"** using staged installation\\n\"},{\"type\":1,\"output\":\"** R\\n\"},{\"type\":1,\"output\":\"** data\\n\"},{\"type\":1,\"output\":\"*** moving datasets to lazyload DB\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** inst\\n\"},{\"type\":1,\"output\":\"** byte-compile and prepare package for lazy loading\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** help\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"*** installing help indices\\n\"},{\"type\":1,\"output\":\"*** copying figures\\n\"},{\"type\":1,\"output\":\"** building package indices\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** installing vignettes\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from temporary location\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package can be loaded from final location\\n\"},{\"type\":1,\"output\":\"\"},{\"type\":1,\"output\":\"** testing if installed package keeps a record of temporary installation path\\n\"},{\"type\":1,\"output\":\"* DONE (hcinfer)\\n\"},{\"type\":1,\"output\":\"\"}]"
 chat_resume_conversation="0"
 chat_suspended="0"
 compile_pdf_state="{\"tab_visible\":false,\"running\":false,\"target_file\":\"\",\"output\":\"\",\"errors\":[]}"

NEWS.md

@@ -2,6 +2,7 @@
 
 * Added `boot_pairs()` for pairs (case) bootstrap standard errors and confidence intervals of ordinary least squares coefficients. It resamples the observations with replacement, refits the model on each replicate, and summarizes the sampling distribution of the coefficients, providing an assumption-free empirical reference for the analytic heteroskedasticity-consistent standard errors from `hcinfer()` and `vcov_hc()`. Percentile, basic, and normal intervals are available, the resampling is reproducible through the `seed` argument, and the replicate fits can optionally run in parallel via `purrr::in_parallel()` and `mirai` without changing the numeric result.
 * Added `coef()`, `vcov()`, `confint()`, `print()`, and `plot()` methods for the `hcinfer_boot` objects returned by `boot_pairs()`. `vcov()` returns the bootstrap covariance matrix of the coefficients, `confint()` can recompute intervals at a different `level` or `type` directly from the stored replicates, and `plot()` draws the bootstrap confidence intervals, coloring each coefficient by whether its interval excludes or includes zero.
+* `hcinfer()` and `vcov_hc()` now accept independent HCbeta shape caps from 50 through 25000 inclusive, with defaults of 10000. HC0, HC1, and HCbeta also remain defined for an exact leverage value of one, while HC2, HC3, HC4, HC4m, HC5, and HC5m retain the positive leverage-complement requirement.
 
 # hcinfer 0.1.1
 

R/data.R

@@ -86,3 +86,78 @@
 #' table(PublicSchools2$south)
 #'
 "PublicSchools2"
+
+#' State crime rates and socioeconomic indicators, 2009
+#'
+#' @description
+#' Violent-crime and murder rates together with socioeconomic indicators for the
+#' 50 U.S. states and the District of Columbia in 2009. The data are useful for
+#' illustrating heteroskedasticity-consistent inference in a cross-sectional
+#' design with influential observations.
+#'
+#' @format A tibble with 51 rows and 8 variables:
+#' \describe{
+#'   \item{state}{Name of one of the 50 U.S. states or the District of Columbia.}
+#'   \item{violent}{Violent-crime rate per 100,000 population.}
+#'   \item{murder}{Murder rate per 100,000 population.}
+#'   \item{hs_grad}{Percentage of the population that graduated from high school
+#'   or higher.}
+#'   \item{poverty}{Percentage of the population living below the poverty line.}
+#'   \item{single}{Percentage of households headed by a single parent.}
+#'   \item{white}{Percentage of the population that is white.}
+#'   \item{urban}{Percentage of the population living in urban areas.}
+#' }
+#'
+#' @source
+#' French, J. P. (2023). *api2lm: Functions and Data Sets for the Book 'A
+#' Progressive Introduction to Linear Models'*. R package version 0.2.
+#' \doi{10.32614/CRAN.package.api2lm}. The same data are distributed as the
+#' `statecrime` dataset in the Python `statsmodels` package (Seabold and
+#' Perktold, 2010, <https://www.statsmodels.org/>); the underlying figures come
+#' from the *Statistical Abstract of the United States* (2009) and are in the
+#' public domain.
+#'
+#' @examples
+#' data(Crime2009)
+#' Crime2009[Crime2009$state == "Alabama", ]
+#'
+#' fit <- lm(violent ~ poverty + single, data = Crime2009)
+#' hcinfer(fit, type = "hcbeta")
+#'
+"Crime2009"
+
+#' Boston-area home prices, 1990
+#'
+#' @description
+#' Sale prices, assessed values, and physical characteristics of 88 homes sold
+#' in the Boston, Massachusetts area in 1990. The data are widely used to
+#' illustrate regression and heteroskedasticity-consistent inference.
+#'
+#' @format A tibble with 88 rows and 10 variables:
+#' \describe{
+#'   \item{price}{House price, in thousands of U.S. dollars.}
+#'   \item{assess}{Assessed value, in thousands of U.S. dollars.}
+#'   \item{bdrms}{Number of bedrooms.}
+#'   \item{lotsize}{Size of the lot, in square feet.}
+#'   \item{sqrft}{Size of the house, in square feet.}
+#'   \item{colonial}{Indicator equal to 1 if the home is of colonial style.}
+#'   \item{lprice}{Natural logarithm of `price`.}
+#'   \item{lassess}{Natural logarithm of `assess`.}
+#'   \item{llotsize}{Natural logarithm of `lotsize`.}
+#'   \item{lsqrft}{Natural logarithm of `sqrft`.}
+#' }
+#'
+#' @source
+#' Wooldridge, J. M. (2020). *Introductory Econometrics: A Modern Approach*,
+#' 7th ed. Cengage Learning, Boston, MA. The `hprice1` data are distributed with
+#' the `wooldridge` R package and were originally collected from the real estate
+#' pages of the *Boston Globe*.
+#'
+#' @examples
+#' data(Hprice)
+#' head(Hprice)
+#'
+#' fit <- lm(price ~ lotsize + sqrft + bdrms, data = Hprice)
+#' hcinfer(fit, type = "hcbeta")
+#'
+"Hprice"

R/hc-weights.R

@@ -30,7 +30,7 @@ default_hc_args <- function(type, dots, call = rlang::caller_env()) {
       args
     },
     hcbeta = {
-      allowed <- c("c1", "c2", "lower", "upper")
+      allowed <- c("c1", "c2", "lower", "upper", "a_max", "b_max")
       check_allowed_dots(dots, allowed, call = call)
       args <- utils::modifyList(hc_default_arguments()$hcbeta, dots)
       args$c1 <- check_nonnegative_scalar(args$c1, "c1", call = call)
@@ -50,6 +50,12 @@ default_hc_args <- function(type, dots, call = rlang::caller_env()) {
           call = call
         )
       }
+      args$a_max <- check_scalar_number(args$a_max, "a_max",
+        lower = 50, upper = 25000, call = call
+      )
+      args$b_max <- check_scalar_number(args$b_max, "b_max",
+        lower = 50, upper = 25000, call = call
+      )
       args
     }
   )
@@ -62,22 +68,25 @@ hcbeta_components <- function(leverage, n, p, args, call = rlang::caller_env())
   mu_hat <- mean(w)
   s2_w <- sum((w - mu_hat)^2) / (n - 1)
 
-  if (!is.finite(s2_w) || s2_w <= .Machine$double.eps) {
-    cli::cli_abort(
-      c(
-        "HCbeta cannot estimate beta parameters from these leverages.",
-        "x" = "The sample variance of truncated leverage complements is zero or too small."
-      ),
-      call = call
-    )
-  }
-
-  phi_hat <- mu_hat * (1 - mu_hat) / s2_w - 1
-  a_hat <- mu_hat * phi_hat
-  b_hat <- (1 - mu_hat) * phi_hat
   zeta <- n / (n + 50)
-  a_tilde <- (1 - zeta) + zeta * a_hat
-  b_tilde <- (1 - zeta) + zeta * b_hat
+
+  if (is.finite(s2_w) && s2_w > .Machine$double.eps) {
+    phi_hat <- mu_hat * (1 - mu_hat) / s2_w - 1
+    a_hat <- mu_hat * phi_hat
+    b_hat <- (1 - mu_hat) * phi_hat
+    a_tilde <- min((1 - zeta) + zeta * a_hat, args$a_max)
+    b_tilde <- min((1 - zeta) + zeta * b_hat, args$b_max)
+  } else {
+    # Degenerate leverage configuration: all truncated complements coincide, so
+    # the method-of-moments variance vanishes and the raw shape estimates
+    # diverge. Following the article, the caps define the limiting shapes, which
+    # keeps HCbeta well-defined and drives g_t to the HC1 scale n / (n - p).
+    phi_hat <- Inf
+    a_hat <- Inf
+    b_hat <- Inf
+    a_tilde <- args$a_max
+    b_tilde <- args$b_max
+  }
 
   if (a_tilde <= 0 || b_tilde <= 0 || !is.finite(a_tilde) || !is.finite(b_tilde)) {
     cli::cli_abort(
@@ -90,15 +99,10 @@ hcbeta_components <- function(leverage, n, p, args, call = rlang::caller_env())
     )
   }
 
-  beta_cdf <- stats::pbeta(w, a_tilde, b_tilde)
-  if (any(!is.finite(beta_cdf)) || any(beta_cdf <= 0)) {
-    cli::cli_abort(
-      "HCbeta produced invalid beta distribution probabilities.",
-      call = call
-    )
-  }
-
-  weights <- (n / (n - p)) * (1 / beta_cdf)^(args$c1 / n^args$c2)
+  n_over_np <- n / (n - p)
+  log_beta_cdf <- stats::pbeta(w, a_tilde, b_tilde, log.p = TRUE)
+  exponent <- pmin(-(args$c1 / n^args$c2) * log_beta_cdf, 700)
+  weights <- n_over_np * exp(exponent)
 
   list(
     weights = weights,
@@ -127,6 +131,19 @@ compute_hc_weights <- function(type, leverage, n, p, dots,
   ratio <- leverage / h_bar
   u <- 1 - leverage
 
+  if (
+    type %in% c("hc2", "hc3", "hc4", "hc4m", "hc5", "hc5m") &&
+      any(u <= .Machine$double.eps)
+  ) {
+    cli::cli_abort(
+      c(
+        "At least one leverage value is too close to 1.",
+        "i" = "HC leverage corrections require positive {.code 1 - h_t}."
+      ),
+      call = call
+    )
+  }
+
   out <- switch(
     type,
     hc0 = list(weights = rep(1, n), params = list()),

R/hcinfer.R

@@ -30,8 +30,10 @@
 #' @param null Null values for the coefficient tests. Use a scalar to test all
 #'   coefficients against the same value, or a numeric vector with one value per
 #'   coefficient.
-#' @param ... Method-specific constants passed to [vcov_hc()]. Defaults are
-#'   documented in [vcov_hc()] and can be inspected with [hc_methods()].
+#' @param ... Method-specific constants passed to [vcov_hc()]. For HCbeta,
+#'   `a_max` and `b_max` default to 10000, may be set independently, and must
+#'   each be finite and lie in `[50, 25000]`. See [vcov_hc()] for all other
+#'   method-specific defaults and parameter domains.
 #'
 #' @return
 #' An object of class `hcinfer` containing the fitted HC covariance estimator,
@@ -87,7 +89,8 @@
 #' summary(result)
 #' confint(result)
 #'
-#' hcinfer(fit, type = "hcbeta", c1 = 7, c2 = 0.75, lower = 0.01, upper = 0.99)
+#' # Sensitivity analysis with nondefault HCbeta caps
+#' hcinfer(fit, type = "hcbeta", a_max = 20000, b_max = 20000)
 #' hcinfer(fit, type = "hc5", k = 0.7)
 #' hcinfer(fit, type = "hc5m", k = 0.7, k1 = 1, k2 = 0, k3 = 1)
 #'

R/methods.R

@@ -36,7 +36,7 @@ hc_default_arguments <- function() {
     hc4m = list(),
     hc5 = list(k = 0.7),
     hc5m = list(k = 0.7, k1 = 1, k2 = 0, k3 = 1, gamma1 = 1, gamma2 = 1.5),
-    hcbeta = list(c1 = 7, c2 = 0.75, lower = 0.01, upper = 0.99)
+    hcbeta = list(c1 = 7, c2 = 0.75, lower = 0.01, upper = 0.99, a_max = 10000, b_max = 10000)
   )
 }
 

R/model-info.R

@@ -64,15 +64,6 @@ model_info_lm <- function(object, call = rlang::caller_env()) {
 
   leverage <- pmin(pmax(leverage, 0), 1)
 
-  if (any(1 - leverage <= .Machine$double.eps)) {
-    cli::cli_abort(
-      c(
-        "At least one leverage value is too close to 1.",
-        "i" = "HC leverage corrections require positive {.code 1 - h_t}."
-      ),
-      call = call
-    )
-  }
 
   terms <- names(coefficients)
   if (is.null(terms)) {

R/vcov-hc.R

@@ -24,13 +24,20 @@
 #' * `"hc5"`: `k = 0.7`.
 #' * `"hc5m"`: `k = 0.7`, `k1 = 1`, `k2 = 0`, `k3 = 1`, `gamma1 = 1`, and
 #'   `gamma2 = 1.5`.
-#' * `"hcbeta"`: `c1 = 7`, `c2 = 0.75`, `lower = 0.01`, and `upper = 0.99`.
+#' * `"hcbeta"`: `c1 = 7`, `c2 = 0.75`, `lower = 0.01`, `upper = 0.99`,
+#'   `a_max = 10000`, and `b_max = 10000`.
 #'
 #' For `"hc5"` and `"hc5m"`, `k`, `k1`, `k2`, and `k3` must be nonnegative,
 #' while `gamma1` and `gamma2` must be positive. For `"hcbeta"`, `c1` must be
 #' nonnegative, `c2` must be positive, and `lower` and `upper` must lie in
 #' `(0, 1)` with `lower < upper`. The HCbeta truncation is
 #' \eqn{w_t = max(lower, min(1 - h_t, upper))}.
+#' Both `a_max` and `b_max` must be finite and lie in `[50, 25000]`; they cap
+#' the shrunk Beta shape parameters \eqn{\tilde a} and \eqn{\tilde b} and
+#' default to 10000. The two caps can be set independently through `...` for
+#' sensitivity analysis. HCbeta remains defined when \eqn{h_t = 1} because it
+#' truncates \eqn{1 - h_t} before evaluating the Beta CDF. In contrast, HC2
+#' through HC5m require a strictly positive leverage complement.
 #'
 #' @param object An ordinary least squares model fitted by [stats::lm()].
 #' @param type A character string specifying the HC estimator. The default is
@@ -91,7 +98,8 @@
 #' vcov(cov)
 #' plot(cov)
 #'
-#' vcov_hc(fit, type = "hcbeta", c1 = 7, c2 = 0.75, lower = 0.01, upper = 0.99)
+#' # Sensitivity analysis with nondefault HCbeta caps
+#' vcov_hc(fit, type = "hcbeta", a_max = 20000, b_max = 20000)
 #' vcov_hc(fit, type = "hc5", k = 0.7)
 #' vcov_hc(fit, type = "hc5m", k = 0.7, k1 = 1, k2 = 0, k3 = 1)
 #'

README.Rmd

@@ -69,6 +69,15 @@ result <- hcinfer(fit)
 The default estimator is HCbeta. Use `tests()` and `confint()` to extract
 the main inferential quantities as tibbles.
 
+HCbeta caps the two shrunk Beta shape parameters with `a_max` and `b_max`.
+Both arguments default to 10000 and accept finite values in the inclusive
+interval `[50, 25000]`. They can be set independently through `...` and should
+normally be changed only for a sensitivity analysis.
+
+```{r}
+hcinfer(fit, type = "hcbeta", a_max = 20000, b_max = 20000)
+```
+
 ```{r}
 tests(result)
 confint(result)

README.md

@@ -35,7 +35,7 @@ matrix estimators currently implemented in `hcinfer`.
 | hc4m | HC4m | Modified HC4 correction by Cribari-Neto and da Silva. | none |
 | hc5 | HC5 | High-leverage correction by Cribari-Neto, Souza, and Vasconcellos. | k = 0.7 |
 | hc5m | HC5m | Modified HC5 correction by Li, Zhang, Zhang, and Wang. | k = 0.7, k1 = 1, k2 = 0, k3 = 1, gamma1 = 1, gamma2 = 1.5 |
-| hcbeta | HCbeta | Beta-distribution leverage correction. | c1 = 7, c2 = 0.75, lower = 0.01, upper = 0.99 |
+| hcbeta | HCbeta | Beta-distribution leverage correction. | c1 = 7, c2 = 0.75, lower = 0.01, upper = 0.99, a_max = 10000, b_max = 10000 |
 
 ## Installation
 
@@ -67,6 +67,23 @@ result <- hcinfer(fit)
 The default estimator is HCbeta. Use `tests()` and `confint()` to
 extract the main inferential quantities as tibbles.
 
+HCbeta caps the two shrunk Beta shape parameters with `a_max` and
+`b_max`. Both arguments default to 10000 and accept finite values in the
+inclusive interval `[50, 25000]`. They can be set independently through
+`...` and should normally be changed only for a sensitivity analysis.
+
+``` r
+hcinfer(fit, type = "hcbeta", a_max = 20000, b_max = 20000)
+#> 
+#> ── HCbeta robust inference ─────────────────────────────────────────────────────
+#> Model: `expenditure ~ income_scaled + income_scaled_sq`
+#> Observations: 50 | Parameters: 3
+#> Robust covariance: HCbeta
+#> Confidence level: 95.0% | Normal critical value: 1.9600
+#> Use `summary()` for p-values, test results, confidence intervals, and
+#> diagnostics.
+```
+
 ``` r
 tests(result)
 #> # A tibble: 3 × 8
@@ -93,7 +110,7 @@ the null value used in the tests.
 plot(result)
 ```
 
-<img src="man/figures/README-unnamed-chunk-5-1.png" alt="Robust confidence intervals for the public-schools regression coefficients." width="672" />
+<img src="man/figures/README-unnamed-chunk-6-1.png" alt="Robust confidence intervals for the public-schools regression coefficients." width="672" />
 
 ## Diagnostics
 
@@ -106,7 +123,7 @@ cov_hcbeta <- vcov_hc(fit)
 plot(cov_hcbeta)
 ```
 
-<img src="man/figures/README-unnamed-chunk-6-1.png" alt="HCbeta adjustment factors plotted against leverage values for the public-schools regression." width="672" />
+<img src="man/figures/README-unnamed-chunk-7-1.png" alt="HCbeta adjustment factors plotted against leverage values for the public-schools regression." width="672" />
 
 ## Main Functions
 

_pkgdown.yml

@@ -20,6 +20,7 @@ articles:
     navbar: ~
     contents:
       - hcinfer-methodology
+      - hcinfer-hcbeta-algorithm
       - hcinfer-comparison
       - hcinfer-hcbeta
       - hcinfer-bootstrap
@@ -52,6 +53,8 @@ reference:
     contents:
       - PublicSchools
       - PublicSchools2
+      - Crime2009
+      - Hprice
   - title: S3 methods
     contents:
       - coef.hcinfer