Logistic regression under separation

The sex2 dataset from logistf illustrates logistic regression with separation (Heinze and Schemper, 2002). We compare finite bias-reduced estimates from fastglm with those from logistf. The maximum likelihood estimate, abbreviated MLE, refers to the ordinary fit without the bias adjustment.

data(sex2, package = "logistf")
X <- model.matrix(~ age + oc + vic + vicl + vis + dia, data = sex2)
y <- sex2$case

fit_f <- fastglm(X, y, family = binomial(), firth = TRUE)
fit_l <- logistf::logistf(case ~ age + oc + vic + vicl + vis + dia,
                          data = sex2, pl = FALSE, plconf = NULL,
                          control = logistf::logistf.control(
                              lconv = 1e-12, gconv = 1e-12, xconv = 1e-12,
                              maxit = 200L))

cbind(fastglm = unname(coef(fit_f)), logistf = unname(coef(fit_l)))
#>          fastglm     logistf
#> [1,]  0.12025405  0.12025405
#> [2,] -1.10598133 -1.10598133
#> [3,] -0.06881673 -0.06881673
#> [4,]  2.26887465  2.26887465
#> [5,] -2.11140819 -2.11140819
#> [6,] -0.78831695 -0.78831695
#> [7,]  3.09601183  3.09601183
max(abs(unname(coef(fit_f)) - unname(coef(fit_l))))
#> [1] 1.478671e-09

The maximum absolute coefficient difference above measures numerical agreement on this dataset. The following code compares elapsed time for the same bias-reduced regression, excluding profile-likelihood confidence intervals from both fits. The size of any runtime advantage depends on the problem and computing environment:

mb <- microbenchmark::microbenchmark(
    fastglm = fastglm(X, y, family = binomial(), firth = TRUE),
    logistf = logistf::logistf(case ~ age + oc + vic + vicl + vis + dia,
                               data = sex2, pl = FALSE, plconf = NULL),
    times = 3L
)
knitr::kable(summary(mb, unit = "ms")[, c("expr", "min", "median", "max")], digits = 2,
             col.names = c("Method", "Minimum (ms)", "Median (ms)", "Maximum (ms)"), row.names = FALSE)

General GLM families

The supported response distributions are binomial, Poisson, Gaussian, Gamma, inverse Gaussian, and negative binomial with theta supplied as a fixed value. The standard family definitions must be used because the bias adjustment depends on distribution derivatives beyond those needed for an ordinary GLM. Quasi families, Tweedie families, and user-modified family calculations are therefore excluded from this estimator.

The following examples compare coefficient estimates with the AS_mean option in brglm2::brglmFit(). Numerical agreement depends on the stopping tolerances and on how well the predictors distinguish the separate coefficient effects. The log, inverse, and other links specify how the response mean is related to the linear predictor.

library(brglm2)
set.seed(123)
n <- 500
x1 <- rnorm(n);  x2 <- rnorm(n)
X <- cbind(1, x1, x2)
eta_true <- 0.5 + 0.3 * x1 - 0.2 * x2

# Largest coefficient difference from brglm2, by link.
difference_table <- function(differences) {
    knitr::kable(data.frame(Link = names(differences),
                            Difference = formatC(differences, format = "e", digits = 1)),
                 col.names = c("Link", "Largest coefficient difference"), row.names = FALSE)
}
y_bin <- rbinom(n, 1, plogis(eta_true))
df <- data.frame(y = y_bin, x1 = x1, x2 = x2)

differences <- sapply(c("logit", "probit", "cloglog"), function(lnk) {
    fam <- binomial(lnk)
    f <- fastglm(X, y_bin, family = fam, firth = TRUE, tol = 1e-10)
    b <- glm(y ~ x1 + x2, data = df, family = fam, method = "brglmFit",
             type = "AS_mean", control = list(epsilon = 1e-10, maxit = 200))
    max(abs(unname(coef(f)) - unname(coef(b))))
})
difference_table(differences)

y_pois <- rpois(n, exp(eta_true))

differences <- sapply(c("log", "sqrt"), function(lnk) {
    fam <- poisson(lnk)
    f <- fastglm(X, y_pois, family = fam, firth = TRUE, tol = 1e-10)
    b <- glm(y ~ x1 + x2, data = data.frame(y = y_pois, x1 = x1, x2 = x2),
             family = fam, method = "brglmFit", type = "AS_mean",
             control = list(epsilon = 1e-10, maxit = 200))
    max(abs(unname(coef(f)) - unname(coef(b))))
})
difference_table(differences)

y_gam <- rgamma(n, shape = 5, rate = 5 / exp(eta_true))

differences <- sapply(c("log", "inverse"), function(lnk) {
    fam <- Gamma(lnk)
    f <- fastglm(X, y_gam, family = fam, firth = TRUE, tol = 1e-10)
    b <- glm(y ~ x1 + x2, data = data.frame(y = y_gam, x1 = x1, x2 = x2),
             family = fam, method = "brglmFit", type = "AS_mean",
             control = list(epsilon = 1e-10, maxit = 200))
    max(abs(unname(coef(f)) - unname(coef(b))))
})
difference_table(differences)

y_gauss <- eta_true + rnorm(n, 0, 0.5)

differences <- sapply(c("identity", "log"), function(lnk) {
    fam <- gaussian(lnk)
    y0 <- if (lnk == "log") pmax(y_gauss, 0.01) else y_gauss
    f <- fastglm(X, y0, family = fam, firth = TRUE, tol = 1e-10)
    b <- glm(y ~ x1 + x2, data = data.frame(y = y0, x1 = x1, x2 = x2),
             family = fam, method = "brglmFit", type = "AS_mean",
             control = list(epsilon = 1e-10, maxit = 200))
    max(abs(unname(coef(f)) - unname(coef(b))))
})
difference_table(differences)

mu_ig <- exp(eta_true)
y_ig <- statmod::rinvgauss(n, mean = mu_ig, shape = 5)
y_ig[y_ig <= 0] <- 0.01

f <- fastglm(X, y_ig, family = inverse.gaussian("log"), firth = TRUE, tol = 1e-10)
b <- glm(y ~ x1 + x2, data = data.frame(y = y_ig, x1 = x1, x2 = x2),
         family = inverse.gaussian("log"), method = "brglmFit", type = "AS_mean",
         control = list(epsilon = 1e-10, maxit = 200))
difference_table(c(log = max(abs(unname(coef(f)) - unname(coef(b))))))

Standard errors

Standard errors describe uncertainty around the fitted coefficients. The reported covariance uses the expected information, which measures how much information the fitted response distribution supplies about those coefficients, evaluated at the bias-reduced estimate. Its diagonal entries give squared standard errors. This is a model-based local approximation; it does not use the curvature of a penalized log likelihood.

Binomial, Poisson, and fixed-theta negative-binomial models have dispersion one. For Gaussian and inverse Gaussian models, the mean-adjusted dispersion estimate is \(D/(n-p)\), where \(D\) is the deviance, \(n\) counts observations with positive weight, and \(p\) is the design rank. The deviance measures discrepancy between fitted and observed outcomes. Gamma dispersion requires an additional adjustment to its estimating equation. The fitted dispersion is used consistently in the coefficient adjustment and the covariance calculation.

f <- fastglm(X, y_bin, family = binomial(), firth = TRUE, tol = 1e-10)
b <- glm(y ~ x1 + x2, data = df, family = binomial(), method = "brglmFit",
         type = "AS_mean", control = list(epsilon = 1e-10, maxit = 200))

cbind(fastglm_SE = f$se, brglm2_SE = summary(b)$coefficients[, "Std. Error"])
#>    fastglm_SE  brglm2_SE
#>    0.09423309 0.09423309
#> x1 0.09910290 0.09910290
#> x2 0.09527990 0.09527990
max(abs(f$se - summary(b)$coefficients[, "Std. Error"]))
#> [1] 4.145295e-14

Assessing convergence

A convergence flag records whether the numerical stopping rule was satisfied. For a bias-reduced fit, the adjusted scores provide a more direct assessment of the equations being solved. An ordinary likelihood score can remain appreciably different from zero at a correctly computed bias-reduced estimate, because the adjustment changes the estimating equation. fastglm_diagnostics() therefore reports the ordinary and adjusted scores separately.

fit_gamma <- fastglm(X, y_gam, family=Gamma("log"), firth=TRUE,
                    tol=1e-10)
fastglm_diagnostics(fit_gamma)
#> Numerical diagnostics
#> Convergence: converged (7 iterations)
#> Solver: column-pivoted QR
#> Estimable coefficients (rank): 3
#> Boundary flag: no
#> 
#>                          Check     Value
#>    Standardized ordinary score   0.01401
#>    Standardized adjusted score 1.941e-11
#>  Standardized dispersion score 3.965e-14
#>    Covariance condition number     1.157
#> 
#> Note: Use the adjusted scores to assess bias-reduced convergence; the ordinary
#>   score need not be zero.

The standardized coefficient diagnostic divides each adjusted score by the square root of its expected information and reports the largest absolute value. This accounts for the units of the predictors when assessing how far the equations remain from zero. The dispersion diagnostic applies the same scaling to the dispersion equation. Both quantities should be close to zero at a numerical solution; they are numerical diagnostics rather than hypothesis-test statistics. Families with fixed dispersion, including binomial and Poisson models, have no dispersion equation to check.

Highly skewed Gamma responses can produce an initial coefficient step that overshoots the solution. For the Gamma log link, the fitting procedure holds the current dispersion fixed and reduces such a step until the coefficient objective decreases sufficiently. The final adjusted scores are then checked, so a small step caused by repeated halving does not by itself establish convergence. convergence.reason distinguishes a successful fit from an exhausted iteration limit or an unsuccessful step search. An unsuccessful fit requires examination before its coefficients or uncertainty summaries are used.

Binomial probit and complementary log-log regressions can also produce excessively large coefficient updates. These links transform a linear predictor into a probability, and a large update can move probabilities so close to zero or one that subsequent calculations become unstable. For these two links, an update is reduced when it would change any fitted linear predictor with positive weight by more than one unit. This changes the numerical path while preserving the adjusted equations. A reduced update alone cannot establish convergence.

The housing data provide an example with grouped binary outcomes. Each row represents a number of respondents in one satisfaction category. The first response column counts respondents reporting low satisfaction, and the second counts all other respondents in that row. Combining rows with the same predictors gives the same grouped-binomial model.

housing <- MASS::housing
housing_design <- model.matrix(~Infl+Type+Cont,housing)
housing_response <- with(housing,cbind(Freq*(Sat=="Low"),Freq*(Sat!="Low")))
fit_housing <- fastglm(housing_design,housing_response,
    family=binomial("cloglog"),firth=TRUE,tol=1e-10)
fastglm_diagnostics(fit_housing)
#> Numerical diagnostics
#> Convergence: converged (11 iterations)
#> Solver: column-pivoted QR
#> Reduced scoring steps: 3
#> Estimable coefficients (rank): 7
#> Boundary flag: no
#> 
#>                          Check       Value
#>    Standardized ordinary score     0.07924
#>    Standardized adjusted score   2.182e-11
#>  Standardized dispersion score Unavailable
#>    Covariance condition number       50.37
#> 
#> Note: Use the adjusted scores to assess bias-reduced convergence; the ordinary
#>   score need not be zero.

The display reports how many scoring steps were reduced and the adjusted score at the returned estimates. Sparse, memory-mapped, and streamed fits use the current information matrix for the bias adjustment with these links; using information from an earlier iteration can cause oscillation in difficult cases. Streamed fits consequently require additional reads of the source. These numerical safeguards improve stability but do not guarantee convergence for every dataset or starting value.

Accurate probability calculations are also necessary when outcomes are almost determined by the predictors. A probability can be so close to one that ordinary floating-point subtraction loses the complementary probability, even though that small quantity still affects the bias adjustment. Replacing small probabilities and derivatives separately by numerical lower bounds can therefore change the equations being solved. Probit and complementary log-log bias reduction evaluates these quantities using logarithms of probabilities and stable tail formulas. The likelihood, residuals, covariance, and score diagnostics follow the same calculations. Displayed probabilities can still round to zero or one.

These probability calculations also enter post-fit likelihood and residual summaries and prediction standard errors. The required summaries are computed together, and prediction derivatives are evaluated separately when only those derivatives are needed. This avoids repeated probability calculations and temporary vectors while preserving the adjusted equations and reported quantities. The computational benefit depends on the design and the requested calculations.

The following example separates the binary outcomes at zero exposure. The nuisance predictor, frequency weights, and offsets remain part of the model. The adjusted-score diagnostic assesses the numerical solution directly; agreement between two coefficient vectors would not establish accuracy if both calculations used the same probability approximation.

set.seed(112101)
separation_design <- cbind("(Intercept)"=1,
    exposure=seq(-2,2,length.out=240),nuisance=rnorm(240,sd=.25))
separated_fit <- fastglm(separation_design,
    as.numeric(separation_design[,"exposure"]>0),
    family=binomial("cloglog"),firth=TRUE,
    weights=rep(c(0,1,2,1),60),offset=cos(seq_len(240))/8,
    tol=1e-10,maxit=500)
summary(separated_fit)
#> Generalized linear model
#> 
#> Call:
#> fastglm.default(x = separation_design, y = as.numeric(separation_design[, 
#>     "exposure"] > 0), family = binomial("cloglog"), weights = rep(c(0, 
#>     1, 2, 1), 60), offset = cos(seq_len(240))/8, tol = 1e-10, 
#>     maxit = 500, firth = TRUE)
#> Family: binomial (cloglog link)
#> Estimation: mean bias reduction
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)   
#> (Intercept)  -0.7601     0.5135  -1.480  0.13877   
#> exposure     18.3399     5.8189   3.152  0.00162 **
#> nuisance     -2.7106     1.7560  -1.544  0.12268   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>  Dispersion Null deviance Null df Deviance Residual df  AIC
#>           1         333.7     179    6.298         177 12.3
#> Convergence: converged (65 iterations)
fastglm_diagnostics(separated_fit)
#> Numerical diagnostics
#> Convergence: converged (65 iterations)
#> Solver: column-pivoted QR
#> Reduced scoring steps: 36
#> Estimable coefficients (rank): 3
#> Boundary flag: yes
#> 
#>                          Check       Value
#>    Standardized ordinary score       1.149
#>    Standardized adjusted score   5.174e-12
#>  Standardized dispersion score Unavailable
#>    Covariance condition number       200.4
#> 
#> Note: Use the adjusted scores to assess bias-reduced convergence; the ordinary
#>   score need not be zero.

Solving the adjusted equations is a numerical requirement. It does not make Wald intervals or hypothesis tests exact in small samples, and this example does not assess their coverage or rejection rates.

The adjusted diagnostics are retained with supported fits, including prepared Gamma-log fits that omit their data or covariance. Ordinary scores require the original observations and predictors, and older fitted objects may not contain the stored adjusted diagnostics. An unavailable diagnostic is reported as missing rather than as a zero score.

For streaming fits using method="qr", convergence also depends on the change in fitted linear predictors implied by the remaining adjusted equations. This check is useful when nearly redundant predictors have large, opposite coefficients: their individual updates may remain appreciable while their combined effect changes very little. A joint score diagnostic accounts for combinations of coefficients that separate marginal score checks can miss. The following reader supplies the same Gamma observations in blocks.

gamma_blocks <- split(seq_len(nrow(X)),ceiling(seq_len(nrow(X))/250L))
gamma_reader <- local({
    design <- X
    response <- y_gam
    blocks <- gamma_blocks
    function(k) {
        i <- blocks[[k]]
        list(X=design[i,,drop=FALSE],y=response[i])
    }
})
qr_gamma <- fastglm_streaming(gamma_reader,length(gamma_blocks),Gamma("log"),
                              method="qr",firth=TRUE,tol=1e-8)
fastglm_diagnostics(qr_gamma)
#> Numerical diagnostics
#> Convergence: converged (8 iterations)
#> Solver: column-pivoted QR
#> Estimable coefficients (rank): 3
#> Boundary flag: no
#> 
#>                                 Check      Value
#>           Standardized ordinary score    0.01401
#>           Standardized adjusted score  2.632e-09
#>         Standardized dispersion score -7.194e-09
#>  Remaining change in linear predictor   2.22e-10
#>          Joint adjusted score per row  1.589e-10
#>  Remaining relative dispersion change    4.4e-10
#>        Joint dispersion score per row  3.217e-10
#>           Covariance condition number      1.157
#> 
#> Note: Use the adjusted scores to assess bias-reduced convergence; the ordinary
#>   score need not be zero.

remaining.predictor.step is the largest proposed change in a fitted linear predictor among observations with positive weight. The joint score is scaled by expected information and the number of positive-weight rows; the dispersion diagnostics measure the remaining relative dispersion update and its scaled score. All checks use the reported coefficients and current information. Verification requires additional block reads, and difficult designs can require a further factorization to assess the equations accurately. These costs can matter when reading the source is slow. A failed step search or insufficient numerical precision still produces an unsuccessful fit.

Adjusted scores and penalized likelihood

Mean bias reduction modifies the score equations. For logistic regression, this corresponds to adding half the log determinant of the expected information to the log likelihood. This equivalence does not extend to every family and link. penalized.deviance is therefore reported only for the supported canonical cases; ordinary deviance remains available for all supported families.

f <- fastglm(X, y_bin, family=binomial(), firth=TRUE)
c(deviance=f$deviance, penalized.deviance=f$penalized.deviance)
#>           deviance penalized.deviance 
#>           644.0851           630.0293

The implemented adjustment targets mean bias, which concerns the average estimation error over repeated samples. An adjustment targeting median bias and user-selected strengths of the Jeffreys penalty are outside the current interface.

Profile intervals and tests for logistic regression

An interval based on a coefficient estimate plus or minus a multiple of its standard error is called a Wald interval. This calculation uses a local quadratic approximation to the likelihood, yielding an interval that is symmetric around the estimate. A profile interval instead varies the coefficient of interest and re-estimates the remaining coefficients at each value. The resulting interval can be asymmetric because it follows the shape of the fitted likelihood over a wider range.

For Firth logistic regression with a logit link, confint(..., method="profile") uses the penalized log likelihood described above. Each constrained fit retains all columns of the original information matrix in the penalty. Dropping the tested predictor would change the penalty and therefore the comparison. A candidate coefficient value is included in the interval when twice the decrease in the maximized penalized log likelihood is below the appropriate chi-square cutoff with one degree of freedom.

fit_profile <- fastglm(case ~ age + oc + vic + vicl + vis + dia,
                       data=sex2, family=binomial(), firth=TRUE,
                       tol=1e-10, maxit=200)
profile_ci <- confint(fit_profile, parm=c("age","dia"), method="profile")
cbind(estimate=coef(fit_profile)[rownames(profile_ci)], profile_ci)
#>      estimate      2.5 %     97.5 %
#> age -1.105981 -1.9737884 -0.3074251
#> dia  3.096012  0.7745682  8.0302935
attr(profile_ci,"diagnostics")[,c("parameter","side","converged","lr.error")]
#>   parameter  side converged      lr.error
#> 1       age lower      TRUE -2.331574e-07
#> 2       age upper      TRUE -2.764065e-09
#> 3       dia lower      TRUE  3.265113e-10
#> 4       dia upper      TRUE -6.932233e-12

The parm argument selects the coefficients to profile, limiting the additional fitting work. The diagnostic column lr.error records the difference between the likelihood-ratio statistic at an endpoint and the target cutoff. A failed endpoint is reported as missing with a warning and a diagnostic message. These diagnostics help distinguish an interval obtained successfully from a numerical calculation that requires further examination.

fastglm_profile_test() applies the same comparison to a specified hypothesis. With several coefficients in parm, it tests their hypothesized values jointly while re-estimating the remaining coefficients. The following example tests whether the coefficients for oc and vic are both zero, conditional on the other predictors:

fastglm_profile_test(fit_profile, parm=c("oc","vic"), null=0)
#> Firth penalized likelihood ratio test (asymptotic chi-square reference)
#> Null fit: Firth penalized likelihood
#> 
#>  Statistic df p-value
#>      28.46  2 < 1e-04
#> 
#> Tested coefficients:
#>     Estimate     Null
#> oc  -0.06882  0.00000
#> vic  2.26887  0.00000

The reference chi-square distribution has one degree of freedom for each constrained coefficient. A small p-value indicates that the constrained model fits less well relative to the unconstrained model under this reference distribution. Both the tests and profile intervals use large-sample approximations under independent binomial sampling; profiling does not make their finite-sample calibration exact. Grouped binomial counts and frequency weights are supported when the effective successes and failures are integer counts. Fractional survey or power-likelihood weights require a different inferential justification and are excluded from these procedures.

Profile inference currently requires a retained dense design and response from a converged, full-rank binomial-logit Firth fit. Sparse, mapped, and streaming designs remain useful for fitting, but the profile functions do not convert them to dense form. The default confint() calculation remains the Wald interval, and failed profile calculations return their diagnostics without replacing the result by another interval.

Speed comparison across families

The following comparison uses 10,000 observations and five predictors plus an intercept. Both functions receive the same preconstructed model matrix, so formula processing is excluded from both timings. Each fit is checked for convergence and coefficient agreement before timing. The methods retain their own starting-value procedures, and iteration counts are reported because these procedures and stopping rules affect the work required. These repeated timings describe this dataset and machine; they do not establish a uniform speed advantage across families or designs.

set.seed(123)
n <- 10000;  p <- 5
x1 <- rnorm(n);  x2 <- rnorm(n); x3 <- rnorm(n)
x4 <- rnorm(n);  x5 <- rnorm(n)
X <- cbind("(Intercept)"=1, x1, x2, x3, x4, x5)
eta <- 0.5 + 0.3 * x1 - 0.2 * x2 + 0.1 * x3

families <- list(
    `binomial logit`   = list(fam = binomial("logit"),
                              y = rbinom(n, 1, plogis(eta))),
    `binomial probit`  = list(fam = binomial("probit"),
                              y = rbinom(n, 1, plogis(eta))),
    `binomial cloglog` = list(fam = binomial("cloglog"),
                              y = rbinom(n, 1, plogis(eta))),
    `poisson log`      = list(fam = poisson("log"),
                              y = rpois(n, exp(eta))),
    `Gamma log`        = list(fam = Gamma("log"),
                              y = rgamma(n, shape = 5, rate = 5 / exp(eta))),
    `gaussian identity`= list(fam = gaussian("identity"),
                              y = eta + rnorm(n, 0, 0.5))
)

results <- list()
for (nm in names(families)) {
    fi <- families[[nm]]
    fit_fast <- function() fastglm(X,fi$y,family=fi$fam,firth=TRUE,
                                   tol=1e-8,maxit=100)
    fit_reference <- function() brglm2::brglmFit(X,fi$y,family=fi$fam,
        control=brglm2::brglmControl(type="AS_mean",epsilon=1e-8,maxit=100))
    f <- fit_fast(); b <- fit_reference()
    error <- max(abs(unname(coef(f))-unname(coef(b))))
    stopifnot(f$converged,b$converged,error<1e-5)
    mb <- microbenchmark::microbenchmark(
        fastglm=fit_fast(),brglm2=fit_reference(),times=5L)
    timing <- summary(mb,unit="ms")
    ms <- setNames(timing$median,as.character(timing$expr))
    results[[nm]] <- data.frame(family=nm,fastglm_ms=ms[["fastglm"]],
        brglm2_ms=ms[["brglm2"]],ratio=ms[["brglm2"]]/ms[["fastglm"]],
        fastglm_iterations=f$iter,brglm2_iterations=b$iter,
        max_coefficient_difference=error)
}
timing_table <- do.call(rbind,results)
timing_table$max_coefficient_difference <- formatC(
    timing_table$max_coefficient_difference,format="e",digits=2)
knitr::kable(timing_table,row.names=FALSE,digits=3,
    col.names=c("Family and link","fastglm (ms)","brglm2 (ms)",
                "Time ratio","fastglm iterations","brglm2 iterations",
                "Largest coefficient difference"))

The time ratio divides the brglm2 time by the fastglm time, so values above one favor fastglm. Gamma bias reduction requires an additional calculation for dispersion at each fitting iteration. When observations have equal prior weights, their contributions to part of this calculation are identical. Combining those contributions avoids repeating the same special-function evaluations for every observation. A numerical root search then solves the same dispersion equation with fewer evaluations. Gamma models with distinct prior weights also benefit from the improved search, but cannot reuse identical contributions to the same extent.

References

• Firth, D. (1993). Bias reduction of maximum likelihood estimates. Biometrika, 80(1), 27–38. doi:10.1093/biomet/80.1.27

• Heinze, G. and Schemper, M. (2002). A solution to the problem of separation in logistic regression. Statistics in Medicine, 21(16), 2409–2419. doi:10.1002/sim.1047

• Kosmidis, I. and Firth, D. (2009). Bias reduction in exponential family nonlinear models. Biometrika, 96(4), 793–804. doi:10.1093/biomet/asp055

• Kosmidis, I. and Firth, D. (2021). Jeffreys-prior penalty, finiteness and shrinkage in binomial-response generalized linear models. Biometrika, 108(1), 71–82. doi:10.1093/biomet/asaa052

• Kosmidis, I. (2014). Bias in parametric estimation: reduction and useful side-effects. WIREs Computational Statistics, 6(3), 185–196. doi:10.1002/wics.1296