Run the SWAG algorithm for generalized linear models

Load package

library(swaglm)

Define covariates and parameters

n <- 2000
p <- 100

# create design matrix and vector of coefficients
Sigma <- diag(rep(1 / p, p))
X <- MASS::mvrnorm(n = n, mu = rep(0, p), Sigma = Sigma)
beta <- c(-15, -10, 5, 10, 15, rep(0, p - 5))

Logistic regression

# --------------------- generate from logistic regression with an intercept of one
z <- 1 + X %*% beta
pr <- 1 / (1 + exp(-z))
set.seed(12345)
y <- as.factor(rbinom(n, 1, pr))
y <- as.numeric(y) - 1

Define SWAG hyperparameters

# define swag parameters
quantile_alpha <- .15
p_max <- 20

Run SWAG algorithm

swag_obj <- swaglm::swaglm(
  X = X, y = y, p_max = p_max, family = stats::binomial(),
  alpha = quantile_alpha, verbose = TRUE, seed = 123
)
#> Completed models of dimension 1
#> Completed models of dimension 2
#> Completed models of dimension 3
#> Completed models of dimension 4
#> Completed models of dimension 5
#> Completed models of dimension 6
#> Completed models of dimension 7
#> Completed models of dimension 8
#> Completed models of dimension 9
#> Completed models of dimension 10
#> Completed models of dimension 11
#> Completed models of dimension 12
#> Completed models of dimension 13
#> Completed models of dimension 14
#> Completed models of dimension 15
print(swag_obj)
#> SWAGLM results :
#> -----------------------------------------
#> Input matrix dimension:  2000 100 
#> Number of explored models:  881 
#> Number of dimensions explored:  15

Compute and plot network

swag_network <- compute_network(swag_obj)
plot(swag_network, scale_vertex = 0.05)

Perform test

B <- 10
res_test <- swaglm_test(swag_obj, B = B)
print(res_test)
#> SWAGLM Test Results:
#> ----------------------
#> p-value (Eigen): 1 
#> p-value (Freq): 0

Linear regression

sigma2 <- 4
set.seed(12345)
y <- 1 + X %*% beta + rnorm(n = n, mean = 0, sd = sqrt(sigma2))

Run SWAG algorithm

swag_obj <- swaglm::swaglm(
  X = X, y = y, p_max = p_max, family = stats::gaussian(),
  alpha = quantile_alpha, verbose = TRUE, seed = 123
)
#> Completed models of dimension 1
#> Completed models of dimension 2
#> Completed models of dimension 3
#> Completed models of dimension 4
#> Completed models of dimension 5
#> Completed models of dimension 6
#> Completed models of dimension 7
#> Completed models of dimension 8
#> Completed models of dimension 9
#> Completed models of dimension 10
#> Completed models of dimension 11
#> Completed models of dimension 12
#> Completed models of dimension 13
#> Completed models of dimension 14
#> Completed models of dimension 15
print(swag_obj)
#> SWAGLM results :
#> -----------------------------------------
#> Input matrix dimension:  2000 100 
#> Number of explored models:  841 
#> Number of dimensions explored:  15

Compute and plot network

swag_network <- compute_network(swag_obj)
plot(swag_network, scale_vertex = 0.05)

Perform test

res_test <- swaglm_test(swag_obj, B = B)
print(res_test)
#> SWAGLM Test Results:
#> ----------------------
#> p-value (Eigen): 0.7279 
#> p-value (Freq): 0

Poisson regression

eta <- 1 + X %*% beta
lambda <- exp(eta)
set.seed(12345)
y <- rpois(n = n, lambda = lambda)

Run SWAG algorithm

# Run swag procedure
swag_obj <- swaglm::swaglm(
  X = X, y = y, p_max = p_max, family = stats::poisson(),
  alpha = quantile_alpha, verbose = TRUE, seed = 123
)
#> Completed models of dimension 1
#> Completed models of dimension 2
#> Completed models of dimension 3
#> Completed models of dimension 4
#> Completed models of dimension 5
#> Completed models of dimension 6
#> Completed models of dimension 7
#> Completed models of dimension 8
#> Completed models of dimension 9
#> Completed models of dimension 10
#> Completed models of dimension 11
#> Completed models of dimension 12
#> Completed models of dimension 13
#> Completed models of dimension 14
#> Completed models of dimension 15
print(swag_obj)
#> SWAGLM results :
#> -----------------------------------------
#> Input matrix dimension:  2000 100 
#> Number of explored models:  830 
#> Number of dimensions explored:  15

Compute and plot network

swag_network <- compute_network(swag_obj)
plot(swag_network, scale_vertex = 0.05)

Perform test

res_test <- swaglm_test(swag_obj, B = B)
print(res_test)
#> SWAGLM Test Results:
#> ----------------------
#> p-value (Eigen): 0.9995 
#> p-value (Freq): 0.0015