probe_one {probe}

Fitting PaRtitiOned empirical Bayes Ecm (PROBE) algorithm to sparse high-dimensional linear models.

Description

A wrapper function for the one-at-a-time variant of the PROBE algorithm.

Usage

probe_one(Y, X, ep = 0.001, maxit = 10000, Y_test = NULL, X_test = NULL, 
verbose = FALSE, signal = NULL, eta_i = NULL, alpha = 0.05, plot_ind = FALSE, 
order.method = "lasso", adj = 10, delta = 0.4, update_order= NULL, beta_start= NULL,
seed = NULL)

Arguments

Y

The outcome variable.

X

An n x M matrix of sparse predictors variables.

ep

Value against which to compare convergence criterion (default = 0.001).

maxit

Maximum number of iterations the algorithm will run for (default = 10000).

Y_test

(optional) Test Y data used plotting purposes only (doesn't impact results).

X_test

(optional) Test X data used plotting purposes only (doesn't impact results).

verbose

A logical (true/false) value whether to print algorithm iteration progress and summary quantities (default = FALSE).

signal

(optional) A vector of indicies of the true non-null coefficients. This is used to calculate the true and false discovery rates by iteration for simulated data. Used plotting purposes only (doesn't impact results).

eta_i

(optional) A vector of the true signal. This is used to calculate the MSE by iteration for simulated data. Used plotting purposes only (doesn't impact results).

alpha

(optional) significance level

plot_ind

A logical values (True/False) for whether to output plots on algorithm results and progress (default = FALSE)

order.method

Updating order and initial values of the algorithm. For lasso (default) or ridge, a lasso or a ridge regression model (fit with 10-fold CV) will be fitted and used. The update_order is defined by the absolute values of the coefficient and beta_start is the coefficient values. When using none, update_order and beta_start must be given. random will randomly select the updating order and use very small values for beta_start.

adj

Bandwidth parameter for empirical Bayes E-step. The bandwidth will be equal to adj times Silverman's 'rule of thumb' (default = 10).

delta

Learning rate for iteration t is (1 + t)^(-1 + delta) (default delta = 0.4).

update_order

Manual value for the updating order for when order.method = "none" is used.

beta_start

Manual value for the starting beta coefficients for when order.method = "none" is used.

seed

Seed value to ensure reproducibility when order.method = "lasso", order.method = "ridge", or order.method = "random".

Value

A list including

beta_ast_hat MAP estimates of the regression coefficients (\beta^\ast),

beta_hat, beta_hat_var MAP estimates of the posterior expectation (beta_hat) and variance (beta_hat_var) of the prior mean (\beta) of the regression coefficients assuming \gamma=1,

gamma_hat the posterior expectation of the latent \gamma variables,

sigma2_est MAP estimate of the residual variance,

E_step full results of the final E_step,

count the total number of iterations before convergence.

References

McLain, A. C., Zgodic, A., & Bondell, H. (2022). Sparse high-dimensional linear regression with a partitioned empirical Bayes ECM algorithm. arXiv preprint arXiv:2209.08139..

See Also

predict_probe_func to obtain predictions.

Examples

### Example
data(Sim_data)
data(Sim_data_test)
attach(Sim_data)
attach(Sim_data_test)
plot_ind <- TRUE
adj <- 10

# Run the analysis. Y_test and X_test are included for plotting purposes only
full_res <- probe_one( Y = Y, X = X, Y_test = Y_test, order.method = "random",
X_test = X_test, plot_ind = plot_ind, adj = adj)

# Predicting for test data
pred_res <- predict_probe_func(full_res, X = X_test)
sqrt(mean((Y_test - pred_res$Pred)^2))

# Estimate of the residual variance and true value
full_res$sigma2_est
sigma2_tr

# RMSE of estimated beta coefficients
beta_ast_est <- c(full_res$beta_ast_hat)
sqrt(mean((beta_ast_est - beta_tr)^2))

# Posterior expectation of gamma by true
gamma_est <- full_res$E_step$gamma
table(gamma_est > 0.5, beta_tr > 0)
sum(gamma_est)
sum(gamma_est[beta_tr>0])