context("test-plnnetworkfamily")
data(trichoptera)
## use a subset t save some time
trichoptera <- prepare_data(trichoptera$Abundance[1:20, 1:5], trichoptera$Covariate[1:20, ])
models <- ZIPLNnetwork(Abundance ~ 1, data = trichoptera)
test_that("ZIPLNnetwork: main function, fields access and methods", {
expect_equal(getBestModel(models), getBestModel(models, "BIC"))
Y <- as.matrix(trichoptera$Abundance)
data = list(
Y = Y,
X = model.matrix(Abundance ~ 1, data = trichoptera),
X0 = matrix(NA,0,0),
O = matrix(0, nrow(Y),ncol(Y)),
w = rep(1, nrow(Y)),
formula = Abundance ~ 1
)
## extract the data matrices and weights
ctrl <- PLNmodels:::ZIPLNnetwork_param(trace = 0)
ctrl$ziparam <- "single"
## instantiate
myPLN <- PLNmodels:::ZIPLNnetworkfamily$new(NULL, data, ctrl)
## optimize
myPLN$optimize(data, ctrl$config_optim)
expect_equal(myPLN$criteria$BIC, models$criteria$BIC)
## S3 methods
expect_true(PLNmodels:::isNetworkfamily(myPLN))
expect_is(plot(myPLN), "ggplot")
expect_is(plot(myPLN, reverse = TRUE), "ggplot")
expect_is(plot(myPLN, type = "diagnostic"), "ggplot")
expect_is(getBestModel(myPLN), "ZIPLNfit_sparse")
expect_is(getModel(myPLN, myPLN$penalties[1]), "ZIPLNfit_sparse")
## Field access
expect_true(all(myPLN$penalties > 0))
expect_null(myPLN$stability_path)
expect_true(anyNA(myPLN$stability))
## Other R6 methods
expect_true(is.data.frame(myPLN$coefficient_path()))
subs <- replicate(2,
sample.int(nrow(trichoptera), size = nrow(trichoptera)/2),
simplify = FALSE)
myPLN$stability_selection(subsamples = subs)
expect_is(plot(myPLN, type = "stability"), "ggplot")
expect_true(!is.null(myPLN$stability_path))
expect_true(inherits(myPLN$plot(), "ggplot"))
expect_true(inherits(myPLN$plot_objective(), "ggplot"))
expect_true(inherits(myPLN$plot_stars(), "ggplot"))
})
test_that("ZIPLNnetwork computes the stability path only once.", {
## extract_probs fails if stability selection has not been performed.
expect_error(extract_probs(models),
"Please perform stability selection using stability_selection(Robject) first", fixed = TRUE)
set.seed(1077)
subs <- replicate(2,
sample.int(nrow(trichoptera), size = nrow(trichoptera)/2),
simplify = FALSE)
stability_selection(models, subsamples = subs, force = TRUE)
models$stability_selection(subsamples = subs)
## Stability_path has correct dimensions
p <- getModel(models, index = 1)$p
expect_equal(dim(models$stability_path),
c(length(models$penalties) * p*(p-1)/2L, 5))
## try to compute it again
expect_message(stability_selection(models),
"Previous stability selection detected. Use \"force = TRUE\" to recompute it.")
## extracts the inclusion frequencies
expect_equal(dim(extract_probs(models, index = 1, format = "matrix")),
c(p, p))
expect_length(extract_probs(models, index = 1, format = "vector"),
p*(p-1)/2)
})
test_that("ZIPLNnetwork: matrix of penalties work", {
p <- ncol(trichoptera$Abundance)
W <- diag(1, p, p)
W[upper.tri(W)] <- runif(p*(p-1)/2, min = 1, max = 5)
W[lower.tri(W)] <- t(W)[lower.tri(W)]
myPLN <- ZIPLNnetwork(Abundance ~ 1, data = trichoptera, control = ZIPLNnetwork_param(penalty_weights = W))
## S3 methods
expect_true(PLNmodels:::isNetworkfamily(myPLN))
expect_is(plot(myPLN), "ggplot")
expect_is(plot(myPLN, reverse = TRUE), "ggplot")
expect_is(plot(myPLN, type = "diagnostic"), "ggplot")
expect_is(getBestModel(myPLN), "ZIPLNfit_sparse")
expect_is(getModel(myPLN, myPLN$penalties[1]), "ZIPLNfit_sparse")
## Field access
expect_true(all(myPLN$penalties > 0))
expect_null(myPLN$stability_path)
expect_true(anyNA(myPLN$stability))
## Other R6 methods
expect_true(is.data.frame(myPLN$coefficient_path()))
subs <- replicate(2,
sample.int(nrow(trichoptera), size = nrow(trichoptera)/2),
simplify = FALSE)
myPLN$stability_selection(subsamples = subs, control = ZIPLNnetwork_param(penalty_weights = W))
expect_is(plot(myPLN, type = "stability"), "ggplot")
expect_true(!is.null(myPLN$stability_path))
expect_true(inherits(myPLN$plot(), "ggplot"))
expect_true(inherits(myPLN$plot_objective(), "ggplot"))
expect_true(inherits(myPLN$plot_stars(), "ggplot"))
## misspecification of penalty weights should induce errors
## not symmetric
W <- diag(1, p, p)
W[upper.tri(W)] <- runif(p*(p-1)/2, min = 1, max = 5)
expect_error(PLNnetwork(Abundance ~ 1, data = trichoptera, control = PLNnetwork_param(penalty_weights = W)))
## not square
W <- matrix(1, p + 1, p)
expect_error(PLNnetwork(Abundance ~ 1, data = trichoptera, control = PLNnetwork_param(penalty_weights = W)))
## nonpositive entries
W <- matrix(0, p, p)
expect_error(PLNnetwork(Abundance ~ 1, data = trichoptera, control = PLNnetwork_param(penalty_weights = W)))
})
test_that("PLNnetwork: list of matrices of penalties work", {
p <- ncol(trichoptera$Abundance)
W <- diag(1, p, p)
W[upper.tri(W)] <- runif(p*(p-1)/2, min = 1, max = 5)
W[lower.tri(W)] <- t(W)[lower.tri(W)]
list_W <- lapply(seq(1, 1e-2, len = 30), function(rho) rho * W)
myPLN <- ZIPLNnetwork(Abundance ~ 1, data = trichoptera, control = ZIPLNnetwork_param(penalty_weights = list_W))
## S3 methods
expect_true(PLNmodels:::isNetworkfamily(myPLN))
expect_is(plot(myPLN), "ggplot")
expect_is(plot(myPLN, reverse = TRUE), "ggplot")
expect_is(plot(myPLN, type = "diagnostic"), "ggplot")
expect_is(getBestModel(myPLN), "ZIPLNfit_sparse")
expect_is(getModel(myPLN, myPLN$penalties[1]), "ZIPLNfit_sparse")
## Field access
expect_true(all(myPLN$penalties > 0))
expect_null(myPLN$stability_path)
expect_true(anyNA(myPLN$stability))
## Other R6 methods
expect_true(is.data.frame(myPLN$coefficient_path()))
subs <- replicate(2,
sample.int(nrow(trichoptera), size = nrow(trichoptera)/2),
simplify = FALSE)
myPLN$stability_selection(subsamples = subs, control = ZIPLNnetwork_param(penalty_weights = W))
expect_is(plot(myPLN, type = "stability"), "ggplot")
expect_true(!is.null(myPLN$stability_path))
expect_true(inherits(myPLN$plot(), "ggplot"))
expect_true(inherits(myPLN$plot_objective(), "ggplot"))
expect_true(inherits(myPLN$plot_stars(), "ggplot"))
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.