Merge pull request #39 from ikosmidis/develop

Develop
ikosmidis · Sep 12, 2024 · 444ccf9 · 444ccf9
2 parents f748bc9 + 2769231
commit 444ccf9
Show file tree

Hide file tree

Showing 18 changed files with 86 additions and 62 deletions.
diff --git a/.Rbuildignore b/.Rbuildignore
@@ -30,3 +30,5 @@ README_cache
 ^CRAN-RELEASE$
 ^code_of_conduct\.md$
 ^CRAN-SUBMISSION$
+^doc$
+^Meta$
diff --git a/.gitignore b/.gitignore
@@ -115,4 +115,6 @@ flycheck_*.el
 .dir-locals.el
 
 # network security
-/network-security.data
+/network-security.data
+/doc/
+/Meta/
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: brglm2
 Title: Bias Reduction in Generalized Linear Models
-Version: 0.9.1
+Version: 0.9.2
 Authors@R: c(person(given = "Ioannis", family = "Kosmidis", role = c("aut", "cre"), email = "ioannis.kosmidis@warwick.ac.uk", comment = c(ORCID = "0000-0003-1556-0302")),
 	     person(given = c("Euloge", "Clovis"), family = c("Kenne Pagui"), role = "aut", email = "kenne@stat.unipd.it"),
 	     person(given = "Kjell", family = "Konis", role = "ctb", email = "kjell.konis@me.com"),

diff --git a/NEWS.md b/NEWS.md
@@ -1,6 +1,17 @@
+# brglm2 0.9.2
+
+## Improvements, updates and additions
+
+* Convergence of the `brglm_fit` iterations is now determined if the L^Inf norm of the step size (rather than the L^1 as it was previously) of the quasi-Fisher scoring procedure is less than `epsilon` (see `?brglm_control` for the definition of `epsilon`). This is more natural as `epsilon` then determines directly the precision of the reported estimates and does not depend on their number.
+
+* `brglm_control()` now checks that the supplied value of `max_step_factor` is numeric and greater or equal to `1`. If not, then it is set to the default value of `12`.
+
+* Vignette updates
+
+
 # brglm2 0.9.1
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 
 * Added the `enzymes` and `hepatitis` data sets (from the [*pmlr*](https://cran.r-project.org/package=pmlr)) to support examples and tests.
 
@@ -14,7 +25,7 @@
 
 * Fixed a bug where the dispersion in the resulting object would not be transformed even if `transformation != "identity"` when `type` is `ML` or `AS_median` or `AS_mixed`.
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 
 * Moved unit tests to [**tinytest**](https://cran.r-project.org/package=tinytest).
 
@@ -26,7 +37,7 @@
 
 # brglm2 0.8.2
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 
 * Housekeeping.
 * Removed lpSolveAPI from imports.
@@ -37,7 +48,7 @@
 
 * Fixed a bug when predicting from `bracl` objects with non-identifiable parameters.
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 
 * Work on output consistently from `print()` methods for `summary.XYZ`
   objects; estimator type is now printed and other fixes.
@@ -72,7 +83,7 @@
 * `brglmFit()` iteration returns last estimates that worked if
   iteration fails.
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 
 * Documentation and example updates.
 
@@ -92,7 +103,7 @@
   (`check_aliasing = FALSE`) rank deficiency checks (through a QR
   decomposition of the model matrix), saving some computational effort.
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * updated DOI links in documentation and some http -> https fixes.
 
 # brglm2 0.7.0
@@ -104,15 +115,15 @@
 ## New functionality
 * `confint` method for `brmulitnom` objects
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * Updated reference to [Kenne Pagui et al (2017)](https://doi.org/10.1093/biomet/asx046).q
 * Updated reference to [Kosmidis and Firth (2020)](https://doi.org/10.1093/biomet/asaa052).
 * Fixed issues with references.
 * Updated documentation.
 
 # brglm2 0.6.2
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * `vcov.brglmFit()` now uses `vcov.summary.glm()` and supports the
   `complete` argument for controlling whether the variance covariance
   matrix should include rows and columns for aliased parameters.
@@ -130,7 +141,7 @@
   `bracl` objects.
 * `detect_separation()` now handles one-column model matrices correctly.
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * Documentation improvements and typo fixes.
 
 # brglm2 0.6
@@ -140,7 +151,7 @@
   supported generalized linear models. See the help files of
   `brglmControl()` and `brglmFit()` for details.
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * Documentation updates and improvements.
 * Updated vignettes to include maximum penalized likelihood with
   powers of the Jeffreys prior as penalty.
@@ -156,14 +167,14 @@
   for more fine-tuning of the starting values when `brglmFit()` is
   called with `start = NULL`.
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * Documentation updates and improvements.
 * Added Kosmidis et al (2019) in the description file.
 * Added tests for `brglmControl()`.
 
 # brglm2 0.5.1
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * Fixed typos in vignettes and documentation.
 * Added ORCHID for Ioannis Kosmidis in DESCRIPTION.
 
@@ -191,15 +202,15 @@ Added `residuals()` methods for `brmultinom` and `bracl` objects
   misclassification in binomial response models (Neuhaus, 1999,
   Biometrika).
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * Improved `summary()` method for `brmultinom` objects.
 * Better starting values for null fits.
 * Added references to [arxiv:1804.04085](https://arxiv.org/abs/1804.04085) in
   documentation.
 * Updated reference to [Kenne Pagui et al (2017)](https://doi.org/10.1093/biomet/asx046).
 
 # brglm2 0.1.8
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * Improved documentation examples.
 * Removed warning about observations with non-positive weights from brmultinom.
 * Updated email address for Ioannis Kosmidis in brglmFit.
@@ -212,15 +223,15 @@ Added `residuals()` methods for `brmultinom` and `bracl` objects
 
 # brglm2 0.1.7
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * Eliminated errors from markdown chunks in multinomial vignette.
 
 # brglm2 0.1.6
 
 ## Bug fixes
 * Compatibility with new version of enrichwith.
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * New email for Ioannis Kosmidis.
 
 # brglm2 0.1.5
@@ -237,7 +248,7 @@ Added `residuals()` methods for `brmultinom` and `bracl` objects
   `detect_separation()` methods in line with the update of
   `glm.fit()`.
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * less strict tolerance in `brglm_control()`.
 * Updates to help files.
 * Fixed typos in iteration vignette.
@@ -263,7 +274,7 @@ Added `residuals()` methods for `brmultinom` and `bracl` objects
   `detectSeparationControl()`, `detect_separation_control()`,
   `checkInfiniteEstimates()`, `check_infinite_estimates()`).
 
-## Other improvements, updates and additions
+## Improvements, updates and additions
 * Minor enhancements in the codebase.
 * The inverse expected information matrix is computed internally using
   `cho2inv()`.
@@ -272,12 +283,8 @@ Added `residuals()` methods for `brmultinom` and `bracl` objects
 
 # brglm2 0.1.3
 
-## Bug fixes
-
-## New functionality
-
-## Other improvements, updates and additions
-* Fixed typo in f_{Y_i}(y) in iteration vignette (thanks to Eugene
+## Improvements, updates and additions
+* Fixed typo in $f_{Y_i}(y)$ in iteration vignette (thanks to Eugene
   Clovis Kenne Pagui for spotting),
 
 # brglm2 0.1.2

diff --git a/R/brglm2-package.R b/R/brglm2-package.R
@@ -96,6 +96,7 @@
 #' 185-196. \doi{10.1002/wics.1296}.
 #'
 #' @docType package
+#' @aliases brglm2-package
 #' @name brglm2
 #' @import stats
 #' @import enrichwith

diff --git a/R/brglmControl.R b/R/brglmControl.R
@@ -201,6 +201,11 @@ brglmControl <- function(epsilon = 1e-06, maxit = 100,
     }
     if (!is.numeric(epsilon) || epsilon <= 0)
         stop("value of 'epsilon' must be > 0")
+
+    if (!is.numeric(max_step_factor) || max_step_factor < 1) {
+        warning("`max_step_factor = ", deparse(max_step_factor), "` is not a permissible value. Defaulting to 12")
+        max_step_factor <- 12
+    }
     list(epsilon = epsilon, maxit = maxit, trace = trace,
          check_aliasing = check_aliasing,
          response_adjustment = response_adjustment,

diff --git a/R/brglmFit.R b/R/brglmFit.R
@@ -876,29 +876,35 @@ brglmFit <- function(x, y, weights = rep(1, nobs), start = NULL, etastart = NULL
             for (iter in seq.int(control$maxit)) {
                 step_factor <- 0
                 testhalf <- TRUE
+
                 ## Inner iteration
                 while (testhalf & step_factor < control$max_step_factor) {
                     ## store previous values
-                    betas0 <- betas; dispersion0 <- dispersion
+                    ## betas0 <- betas
+                    ## dispersion0 <- dispersion
                     step_beta_previous <- step_beta
                     step_zeta_previous <- step_zeta
+
                     ## Update betas
                     betas <- betas + slowit * 2^(-step_factor) * step_beta
+
                     ## Update zetas
                     if (!no_dispersion & df_residual > 0) {
                         transformed_dispersion <- eval(control$Trans)
                         transformed_dispersion <- transformed_dispersion + 2^(-step_factor) * step_zeta
                         dispersion <- eval(control$inverseTrans)
                     }
+
                     ## Compute key quantities
                     theta <- c(betas, dispersion)
                     transformed_dispersion <- eval(control$Trans)
-                    ## Mean quantities
 
+                    ## Mean quantities
                     quantities <- try(key_quantities(theta, y = y, level = 2 * !no_dispersion, scale_totals = has_fixed_totals, qr = TRUE), silent = TRUE)
                     ## This is to capture qr failing and revering to previous estimates
                     if (failed_adjustment_beta <- inherits(quantities, "try-error")) {
-                        betas <- betas0; dispersion <- dispersion0
+                        ## betas <- betas0
+                        ## dispersion <- dispersion0
                         warning("failed to calculate score adjustment")
                         break
                     }
@@ -912,10 +918,9 @@ brglmFit <- function(x, y, weights = rep(1, nobs), start = NULL, etastart = NULL
                         warning("failed to calculate score adjustment")
                         break
                     }
-                    adjusted_grad_beta <- with(step_components_beta, {
-                        grad + adjustment
-                    })
+                    adjusted_grad_beta <- with(step_components_beta, grad + adjustment)
                     step_beta <- drop(step_components_beta$inverse_info %*% adjusted_grad_beta)
+
                     ## Dispersion quantities
                     if (no_dispersion) {
                         adjusted_grad_zeta <- step_zeta <- NA_real_
@@ -929,28 +934,28 @@ brglmFit <- function(x, y, weights = rep(1, nobs), start = NULL, etastart = NULL
                             warning("failed to calculate score adjustment")
                             break
                         }
-                        adjusted_grad_zeta <- with(step_components_zeta, {
-                            grad + adjustment
-                        })
+                        adjusted_grad_zeta <- with(step_components_zeta, grad + adjustment)
                         step_zeta <- as.vector(adjusted_grad_zeta * step_components_zeta$inverse_info)
                     }
+
+                    ## Convergence criteria
+                    linf_current <- max(abs(c(step_beta, step_zeta)), na.rm = TRUE)
+                    linf_previous <- max(abs(c(step_beta_previous, step_zeta_previous)), na.rm = TRUE)
+                    testhalf <- linf_current > linf_previous
+
                     ## Continue inner loop
-                    if (step_factor == 0 & iter == 1)  {
-                        testhalf <- TRUE
-                    } else {
-                        s2 <- c(abs(step_beta), abs(step_zeta))
-                        s1 <- c(abs(step_beta_previous), abs(step_zeta_previous))
-                        testhalf <- sum(s2, na.rm = TRUE) > sum(s1, na.rm = TRUE)
-                    }
+                    ## if (step_factor == 0 & iter == 1)  {
+                    ##     testhalf <- TRUE
+                    ## }
                     step_factor <- step_factor + 1
+
                     ##  Trace here
                     if (control$trace) {
-
                         trace_iteration()
                     }
                 }
                 failed <- failed_adjustment_beta | failed_inversion_beta | failed_adjustment_zeta | failed_inversion_zeta
-                if (failed | sum(abs(c(step_beta, step_zeta)), na.rm = TRUE) < control$epsilon) {
+                if (failed | linf_current < control$epsilon) {
                     break
                 }
             }
@@ -987,6 +992,7 @@ brglmFit <- function(x, y, weights = rep(1, nobs), start = NULL, etastart = NULL
         ## residuals and working_weights
 
         quantities <- key_quantities(c(betas, dispersion), y = y, level = 2 * !no_dispersion, scale_totals = has_fixed_totals, qr = TRUE)
+
         qr.Wx <- quantities$qr_decomposition
 
         mus <- quantities$mus

diff --git a/R/brmultinom.R b/R/brmultinom.R
@@ -142,7 +142,7 @@
 #' contrasts(hepat$type) <- contr.treatment(3, base = 1)
 #'
 #' # Maximum likelihood estimation fails to converge because some estimates are infinite
-#' hepML <- brmultinom(type ~ group * time, data = hepat, weights = counts, type = "ML")
+#' hepML <- brmultinom(type ~ group * time, data = hepat, weights = counts, type = "ML", slowit = 0.1)
 #'
 #' # Mean bias reduction returns finite estimates
 #' hep_meanBR <- brmultinom(type ~ group * time, data = hepat, weights = counts, type = "AS_mean")

diff --git a/inst/tinytest/test-binomial.R b/inst/tinytest/test-binomial.R
@@ -2,7 +2,7 @@ library("brglm")
 data("lizards", package = "brglm2")
 links <- lapply(c("logit", "probit", "cloglog", "cauchit"), make.link)
 
-tol <- 1e-10
+tol <- 1e-08
 for (l in seq_along(links)) {
     expect_warning(
         lizardsBRlegacy <- brglm(cbind(grahami, opalinus) ~ height + diameter +

diff --git a/inst/tinytest/test-bracl.R b/inst/tinytest/test-bracl.R
@@ -85,19 +85,19 @@ s2p <- summary(fit_vgam_p)
 
 tol <- 1e-06
 ## summary method for bracl returns the correct coef mat
-expect_equal(coef(s1), coef(s2), tolerance = tol, check.attributes = FALSE)
-expect_equal(coef(s1p), coef(s2p), tolerance = tol, check.attributes = FALSE)
+expect_equal(coef(s1), s2@coef3, tolerance = tol, check.attributes = FALSE)
+expect_equal(coef(s1p), s2p@coef3, tolerance = tol, check.attributes = FALSE)
 
 newdata <- expand.grid(gender = c("male", "female"),  religion = c("moderate", "fundamentalist"))
 ## predict.bracl works as expected
 pp <- predict(fit_bracl_p, newdata = stemcell, type = "probs")
 p <- predict(fit_bracl, newdata = stemcell, type = "probs")
 expect_equal(predict(fit_vgam_p, type = "response"),
              pp[19:24, ],
-             tolerance = 1e-08, check.attributes = FALSE)
+             tolerance = 1e-06, check.attributes = FALSE)
 expect_equal(predict(fit_vgam, type = "response"),
              p[19:24, ],
-             tolerance = 1e-08, check.attributes = FALSE)
+             tolerance = 1e-06, check.attributes = FALSE)
 
 ## no intercept returns warning
 expect_warning(

diff --git a/inst/tinytest/test-jeffreys.R b/inst/tinytest/test-jeffreys.R
@@ -41,7 +41,7 @@ data("lizards", package = "brglm2")
 links <- lapply(c("logit", "probit", "cloglog", "cauchit"), make.link)
 
 
-tol <- 1e-10
+tol <- 1e-08
 for (l in seq_along(links)) {
     expect_warning(
         lizardsBRlegacy <- brglm(cbind(grahami, opalinus) ~ height + diameter +

diff --git a/inst/tinytest/test-transformation.R b/inst/tinytest/test-transformation.R
@@ -24,7 +24,7 @@ c_sqrt[5] <- c_sqrt[5]^2
 c_inverse <- coef(mod_inverse, model = "full")
 c_inverse[5] <- 1/c_inverse[5]
 
-tol <- sqrt(.Machine$double.eps)
+tol <- 1e-08
 ## ML estimate of gamma dispersion from brglmFit is invariant to trasnformation
 expect_equal(c_identity, c_log, tolerance = tol, check.attributes = FALSE)
 expect_equal(c_identity, c_sqrt, tolerance = tol, check.attributes = FALSE)

diff --git a/man/brglm2.Rd b/man/brglm2.Rd