context("compare coefficients to glmnet when no external data (gaussian)")
##### Code used to generate data files for all tests #####
# set.seed(123)
# n <- 100
# p <- 50
# q <- 5
# meanx_true <- rep(0, p)
# covx_true <- matrix(NA, nrow = p, ncol = p)
# for (i in 1:p) {
# for (j in 1:p) {
# covx_true[i, j] <- 0.5^abs(i-j)
# }
# }
# meanz_true <- rep(0, q)
# covz_true <- matrix(NA, nrow = q, ncol = q)
# for (i in 1:q) {
# for (j in 1:q) {
# covz_true[i, j] <- 0.5^abs(i-j)
# }
# }
# a0 <- 0.01
# a <- c(0.1, -0.1, rep(0, q - 2))
# z <- mvrnorm(n = p, mu = meanz_true, Sigma = covz_true)
# b <- drop(a0 + z %*% a + 0.2*rnorm(p))
# x <- mvrnorm(n = n, mu = meanx_true, Sigma = covx_true)
# y <- drop(x %*% b + rnorm(n))
# mean_x <- colMeans(x)
# var_x <- apply(x, 2, var) * (n-1) / n
# sd_x <- sqrt(var_x)
# xscaled <- matrix(NA, nrow = n, ncol = p)
# for (i in 1:p) {
# xscaled[, i] <- (x[, i] - mean_x[i]) / sd_x[i]
# }
# mean_z <- colMeans(z)
# sd_z <- sqrt(apply(z, 2, var) * (p-1) / p)
# zscaled <- matrix(NA, nrow = p, ncol = q)
# for (i in 1:q) {
# zscaled[, i] <- (z[, i] - mean_z[i]) / sd_z[i]
# }
# Ridge Regression #
test_that("x standardized, intercept", {
# coef_glmnet <- glmnet(x = x, y = y, family = "gaussian", alpha = 0, thresh = 1e-15)
penalty <- define_penalty(penalty_type = 0, num_penalty = 100)
expect_equal(
betas_glmnet[, 1],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
expect_equal(
b0_glmnet[1],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
penalty_main = penalty,
control = test_control
)$beta0[10, 1] * sd_y,
tolerance = 1e-5
)
})
test_that("x NOT standardized, intercept", {
# coef_glmnet <- glmnet(x = x, y = y, family = "gaussian", standardize = FALSE, alpha = 0, thresh = 1e-15)
penalty <- define_penalty(penalty_type = 0, num_penalty = 100)
expect_equal(
betas_glmnet[, 2],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
expect_equal(
b0_glmnet[2],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$beta0[10, 1] * sd_y,
tolerance = 1e-5
)
})
test_that("x standardized, NO intercept", {
# coef_glmnet <- glmnet(x = x, y = y, family = "gaussian", intercept = FALSE, alpha = 0, thresh = 1e-15)
penalty <- define_penalty(penalty_type = 0, num_penalty = 100)
expect_equal(
betas_glmnet[, 3],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(F, F),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
})
test_that("x NOT standardized, NO intercept", {
# coef_glmnet <- glmnet(x = x, y = y, family = "gaussian", standardize = FALSE, intercept = FALSE, alpha = 0, thresh = 1e-15)
penalty <- define_penalty(penalty_type = 0, num_penalty = 100)
expect_equal(
betas_glmnet[, 4],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(F, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
expect_equal(
b0_glmnet[4],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(F, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$beta0[10, 1] * sd_y,
tolerance = 1e-5
)
})
# Lasso Regression #
test_that("x standardized, intercept", {
# coef_glmnet <- glmnet(x = x, y = y, family = "gaussian", alpha = 1, thresh = 1e-15)
penalty <- define_penalty(penalty_type = 1, num_penalty = 100)
expect_equal(
betas_glmnet[, 5],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
expect_equal(
b0_glmnet[5],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
penalty_main = penalty,
control = test_control
)$beta0[10, 1] * sd_y,
tolerance = 1e-5
)
})
test_that("x NOT standardized, intercept", {
# coef_glmnet <- glmnet(x = x, y = y, family = "gaussian", standardize = FALSE, alpha = 1, thresh = 1e-15)
penalty <- define_penalty(penalty_type = 1, num_penalty = 100)
expect_equal(
betas_glmnet[, 6],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
expect_equal(
b0_glmnet[6],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$beta0[10, 1] * sd_y,
tolerance = 1e-5
)
})
test_that("x standardized, NO intercept", {
# coef_glmnet <- glmnet(x = x, y = y, family = "gaussian", intercept = FALSE, alpha = 1, thresh = 1e-15)
penalty <- define_penalty(penalty_type = 1, num_penalty = 100)
expect_equal(
betas_glmnet[, 7],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(F, F),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
})
test_that("x NOT standardized, NO intercept", {
# coef_glmnet <- glmnet::glmnet(x = x, y = y, family = "gaussian", standardize = FALSE, intercept = FALSE, alpha = 1, thresh = 1e-15)
penalty <- define_penalty(penalty_type = 1, num_penalty = 100)
expect_equal(
betas_glmnet[, 8],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(F, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
})
# Elastic Net Regression #
test_that("x standardized, intercept", {
# coef_glmnet <- glmnet::glmnet(x = x, y = y, family = "gaussian", alpha = 0.5, thresh = 1e-15)
# betas_glmnet9 <- drop(unname(as.matrix(coef_glmnet$beta[,10])))
# b0_glmnet9 <- unname(coef_glmnet$a0[10])
penalty <- define_penalty(penalty_type = 0.5, num_penalty = 100)
expect_equal(
betas_glmnet[, 9],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
expect_equal(
b0_glmnet[9],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
penalty_main = penalty,
control = test_control
)$beta0[10, 1] * sd_y,
tolerance = 1e-5
)
})
test_that("x NOT standardized, intercept", {
# coef_glmnet <- glmnet::glmnet(x = x, y = y, family = "gaussian", standardize = FALSE, alpha = 0.5, thresh = 1e-15)
# betas_glmnet10 <- drop(unname(as.matrix(coef_glmnet$beta[,10])))
# b0_glmnet10 <- unname(coef_glmnet$a0[10])
penalty <- define_penalty(penalty_type = 0.5, num_penalty = 100)
expect_equal(
betas_glmnet[, 10],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
expect_equal(
b0_glmnet[10],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(T, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$beta0[10, 1] * sd_y,
tolerance = 1e-5
)
})
test_that("x standardized, NO intercept", {
# coef_glmnet <- glmnet::glmnet(x = x, y = y, family = "gaussian", intercept = FALSE, alpha = 0.5, thresh = 1e-15)
# betas_glmnet11 <- drop(unname(as.matrix(coef_glmnet$beta[,10])))
# b0_glmnet11 <- unname(coef_glmnet$a0[10])
penalty <- define_penalty(penalty_type = 0.5, num_penalty = 100)
expect_equal(
betas_glmnet[, 11],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(F, F),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
})
test_that("x NOT standardized, NO intercept", {
# coef_glmnet <- glmnet::glmnet(x = x, y = y, family = "gaussian", standardize = FALSE, intercept = FALSE, alpha = 0.5, thresh = 1e-15)
# betas_glmnet12 <- drop(unname(as.matrix(coef_glmnet$beta[,10])))
# b0_glmnet12 <- unname(coef_glmnet$a0[10])
penalty <- define_penalty(penalty_type = 0.5, num_penalty = 100)
expect_equal(
betas_glmnet[, 12],
xrnet(
x = xtest,
y = ytest_scaled,
family = "gaussian",
intercept = c(F, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1] * sd_y,
tolerance = 1e-5
)
})
# Ridge - No Penalty on 1st two variables #
test_that("x NOT standardized, intercept", {
pf <- c(rep(0, 2), rep(1, NCOL(xtest) - 2))
# coef_glmnet <- glmnet(x = xtest, y = ytest, family = "gaussian", standardize = FALSE, alpha = 0, thresh = 1e-15, penalty.factor = pf)
penalty <- define_penalty(
penalty_type = 0,
user_penalty = path_ridge,
custom_multiplier = NCOL(xtest) * pf[-c(1, 2)] / sum(pf)
)
expect_equal(
betas_glmnet[-c(1, 2), 13],
xrnet(
x = xtest[, -c(1, 2)],
y = ytest,
unpen = xtest[, c(1, 2)],
family = "gaussian",
intercept = c(T, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1],
tolerance = 1e-5
)
})
# Lasso - No Penalty on 1st two variables #
test_that("x NOT standardized, intercept", {
pf <- c(rep(0, 2), rep(1, NCOL(xtest) - 2))
# coef_glmnet <- glmnet(x = xtest, y = ytest, family = "gaussian", standardize = FALSE, alpha = 1, thresh = 1e-15,penalty.factor = pf)
penalty <- define_penalty(
penalty_type = 1,
user_penalty = path_lasso,
custom_multiplier = NCOL(xtest) * pf[-c(1, 2)] / sum(pf)
)
expect_equal(
betas_glmnet[-c(1, 2), 14],
xrnet(
x = xtest[, -c(1, 2)],
y = ytest,
unpen = xtest[, c(1, 2)],
family = "gaussian",
intercept = c(T, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1],
tolerance = 1e-5
)
})
# Elastic Net - No Penalty on 1st two variables #
test_that("x NOT standardized, intercept", {
pf <- c(rep(0, 2), rep(1, NCOL(xtest) - 2))
# coef_glmnet <- glmnet(x = xtest, y = ytest, family = "gaussian", standardize = FALSE, alpha = 0.5, thresh = 1e-15, penalty.factor = pf)
penalty <- define_penalty(
penalty_type = 0.5,
user_penalty = path_en,
custom_multiplier = NCOL(xtest) * pf[-c(1, 2)] / sum(pf)
)
expect_equal(
betas_glmnet[-c(1, 2), 15],
xrnet(
x = xtest[, -c(1, 2)],
y = ytest,
unpen = xtest[, c(1, 2)],
family = "gaussian",
intercept = c(T, F),
standardize = c(F, T),
penalty_main = penalty,
control = test_control
)$betas[, 10, 1],
tolerance = 1e-5
)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.