brglm2 (version 1.0.1)

solve_se: 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,
  ...
)

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").

Arguments

kappa

asymptotic ratio of columns/rows of the design matrix. kappa should be in (0, 1).

ss

square root of signal strength or of 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 the square root of the signal strength and intercept, if numeric, is the oracle intercept value. If TRUE, then ss is the square root of 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().

Author

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

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.

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. 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)
# \donttest{

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)

# }