solve_se: Solve the MDYPL state evolution equations with or without...

In brglm2: Bias Reduction in Generalized Linear Models

Solve the MDYPL state evolution equations with or without intercept, with signal strength or contaminated signal strength

Description

Solve the MDYPL state evolution equations with or without intercept, with signal strength or contaminated signal strength

Usage

solve_se(
  kappa,
  ss,
  alpha,
  intercept = NULL,
  start,
  corrupted = FALSE,
  gh = NULL,
  prox_tol = 1e-10,
  transform = TRUE,
  init_method = "Nelder-Mead",
  init_iter = 50,
  ...
)

Arguments

kappa	asymptotic ratio of columns/rows of the design matrix. kappa should be in ⁠(0, 1)⁠.
ss	signal strength or corrupted signal strength, depending on whether corrupted = TRUE or not. See Details.
alpha	the shrinkage parameter of the MDYPL estimator. alpha should be in ⁠(0, 1)⁠.
intercept	if NULL (default) then the MDYPL state evolution equations for the model with no intercept parameter are solved. If a real then the equations for the models with intercept parameter equal to intercept are solved. See Details.
start	a vector with starting values for mu, b,sigma (and iota if intercept is numeric).
corrupted	if FALSE (default) then ss is signal strength and intercept, if numeric, is the oracle intercept value. If TRUE, then ss is the corrupted signal strength, and intercept, if numeric, is the limit of the estimator computed by mdyplFit() with shrinkage parameter alpha. See Details.
gh	A list with the Gauss-Hermite quadrature nodes and weights, as returned from statmod::gauss.quad() with kind = "hermite". Default is NULL, in which case gh is set to statmod::gauss.quad(200, kind = "hermite").
prox_tol	tolerance for the computation of the proximal operator; default is 1e-10.
transform	if TRUE (default), the optimization is with respect to log(mu), log(b),log(sigma), (and iota if intercept is numeric). If FALSE, then it is over mu, b, sigma (and iota if intercept is numeric). The solution is returned in terms of the latter set, regardless of how optimization took place.
init_method	The method to be passed to optim(). Default is "Nelder-Mead".
init_iter	how many iterations of optim() should we make to get starting values for nleqslv::nleqslv()? Default is 50, but can also be 0 in which case start is directly passed to nleqslv:nleqslv(). init_iter = "only" results in only optim() being used. See Details.
...	further arguments to be passed to nleqslv::nleqslv(), unless init_iter = "only", in which case ... is further arguments to be passed to optim().

Details

init_iter iterations of optim() with method = init_method are used towards minimizing sum(se)^2, where se is a vector of the state evolution functions. The solution is then passed to nleqslv::nleqslv() for a more aggressive iteration. The state evolution equations are given in expressions (8) (model without intercept) and expression (15) (model with intercept) in Sterzinger & Kosmidis (2024).

If corrupted = FALSE (default), then ss is the square root of the signal strength, which is the limit \gamma^2 of var(X \beta). If corrupted = TRUE, then ss is the square root of the corrupted signal strength which is the limit \nu^2 of var(X hat(beta)(\alpha)), where hat(\beta)(\alpha) is the maximum Diaconis-Ylvisaker prior penalized likelihood (MDYPL) estimator as computed by mdyplFit() with shrinkage parameter alpha.

If intercept = NULL, then the state evolution equations are solved for the model without intercept. If intercept is a real number, then the state evolution equations for the model with intercept are solved (i.e. with predictor \eta_i = \theta + x_i^T \beta). In that case, what intercept represents depends on the value of corrupted. If corrupted = FALSE, intercept represents the oracle value of \theta, otherwise it represents the limit iota of the MDYPL estimator of \theta as computed by mdyplFit() with shrinkage parameter alpha.

Note that start is always for mu, b,sigma, as is the result, regardless whether transform = TRUE or not. Transformations during optimization are done internally.

Value

If intercept = NULL, a vector with the values of mu, b,sigma. Otherwise, a vector with the values of mu, b,sigma, and iota, if corrupted = FALSE, or the value of the intercept otherwise. The vector has attributes the state evolution functions at the solution ("funcs"), the number of iterations used by the last optimization method ("iter"), any messages from the last optimization method ("message"), and information on the optimization methods used ("optimization-chain").

Author(s)

Ioannis Kosmidis ⁠[aut, cre]⁠ [e-mail], Federico Boiocchi ⁠[ctb]⁠ [e-mail], Philipp Sterzinger ⁠[ctb, earlier Julia code by]⁠ [e-mail]

References

Zhao Q, Sur P, Cand\'es E J (2022). The asymptotic distribution of the MLE in high-dimensional logistic models: Arbitrary covariance. Bernoulli, 28, 1835–1861. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.3150/21-BEJ1401")}.

Sterzinger P, Kosmidis I (2024). Diaconis-Ylvisaker prior penalized likelihood for p/n \to \kappa \in (0,1) logistic regression. arXiv:2311.07419v2, https://arxiv.org/abs/2311.07419.

Examples

## Reproducing Table 13 of Zhao et al. (2022, DOI: 10.3150/21-BEJ1401)
## Not run: 

thetas <- c(0, 0.5, 1, 2, 2.5)
gamma0 <- 5
pars3 <- matrix(NA, length(thetas), 3)
pars4 <- matrix(NA, length(thetas), 4)
colnames(pars4) <- c("I_mu", "I_b", "I_sigma", "I_iota")
colnames(pars3) <- c("II_mu", "II_b", "II_sigma")
for (i in seq_along(thetas)) {
    start3 <- c(0.5, 1, 1)
    pars3[i, ] <- solve_se(kappa = 0.2, ss = sqrt(5 + thetas[i]^2),
                           alpha = 1, start = start3, init_iter = 0)
    start4 <- c(pars3[i, ], thetas[i])
    pars4[i, ] <- solve_se(kappa = 0.2, ss = sqrt(5), intercept = thetas[i],
                           alpha = 1, start = start4, init_iter = 0)
}

cbind(pars3, pars4)


## End(Not run)

brglm2 documentation built on Aug. 29, 2025, 5:25 p.m.