--- a
+++ b/tests/testthat/test-projection.R
@@ -0,0 +1,208 @@
+test_that("projsel for a model with no penalized predictors",
+{
+    SW({
+        sel.base.00 <- projsel(hs.base, max.iters=0)
+        sel.base.02 <- projsel(hs.base, max.iters=0, start.from="X2")
+        sel.base.1 <- projsel(hs.base)
+        sel.base.2 <- projsel(hs.base, start.from="X2")
+    })
+
+    expect_equal(nrow(sel.base.00),
+                 1)
+    expect_true(is.na(sel.base.00$rel.kl))
+    expect_equal(sel.base.00$rel.kl.null,
+                 0)
+
+    expect_equal(nrow(sel.base.02),
+                 2)
+    expect_equal(sel.base.02$rel.kl[2],
+                 0)
+    expect_equal(sel.base.02$rel.kl.null[1],
+                 0)
+    expect_equal(sel.base.02$rel.kl.null[2],
+                 0.38193515, tolerance=tol)
+
+    expect_equal(nrow(sel.base.1),
+                 length(hs.base$betas$unpenalized))
+    expect_equal(attr(sel.base.1, "start.from"),
+                 character(0))
+
+    expect_equal(nrow(sel.base.2),
+                 length(hs.base$betas$unpenalized))
+    expect_equal(sel.base.2$var[2],
+                 "Initial submodel")
+    expect_equal(attr(sel.base.2, "start.from"),
+                 "X2")
+})
+
+test_that("projsel for gaussian family",
+{
+    SW({
+        sel.gauss <- projsel(hs.gauss, out.csv="out.csv")
+    })
+
+    expect_s3_class(sel.gauss,
+                    "projsel")
+    expect_equal(colnames(sel.gauss),
+                 c("var", "kl", "rel.kl.null", "rel.kl", "elpd", "delta.elpd"))
+    expect_equal(nrow(sel.gauss),
+                 length(pen) + 2)
+    expect_equal(attr(sel.gauss, "start.from"),
+                 hs.gauss$model.terms$unpenalized)
+    expect_equal(sel.gauss$var[1:2],
+                 c("Intercept only", "Initial submodel"))
+    expect_equal(sel.gauss$var[-c(1:2)],
+                 paste0("X", c(9, 4, 8, 7, 5, 6, 10)))
+    expect_equal(sel.gauss$kl[2],
+                 0.02948656, tolerance=tol)
+    expect_equal(tail(sel.gauss$kl, n=1),
+                 0)
+    expect_equal(sel.gauss$rel.kl.null[1],
+                 0)
+    expect_equal(tail(sel.gauss$rel.kl.null, n=1),
+                 1)
+    expect_true(is.na(sel.gauss$rel.kl[1]))
+    expect_equal(sel.gauss$rel.kl[2],
+                 0)
+    expect_equal(tail(sel.gauss$rel.kl, n=1),
+                 1)
+    expect_equal(sel.gauss$elpd[2],
+                 -111.370949, tolerance=tol)
+    expect_equal(tail(sel.gauss$elpd, n=1),
+                 -109.695740, tolerance=tol)
+    expect_equal(sel.gauss$delta.elpd[2],
+                 -1.67520810, tolerance=tol)
+    expect_true(all(diff(sel.gauss$kl) < 0))
+    expect_true(file.exists("out.csv"))
+    unlink("out.csv")
+
+    expect_s3_class(plot(sel.gauss, from.covariates=FALSE, max.points=3),
+                    "ggplot")
+    expect_silent(print(plot(sel.gauss, title="Test", max.labels=3)))
+    unlink("Rplots.pdf")
+})
+
+test_that("projsel for binomial family",
+{
+    num.sel <- 4
+    SW({
+        sel.binom <- projsel(hs.binom, max.iters=num.sel)
+    })
+
+    expect_equal(nrow(sel.binom),
+                 num.sel + 2)
+    expect_equal(sel.binom$var[-c(1:2)],
+                 paste0("X", c(6, 9, 5, 8)))
+    expect_equal(sel.binom$kl[2],
+                 0.03245156, tolerance=tol)
+    expect_equal(tail(sel.binom$kl, n=1),
+                 0.00478297, tolerance=tol)
+    expect_equal(sel.binom$rel.kl.null[1],
+                 0)
+    expect_true(is.na(sel.binom$rel.kl[1]))
+    expect_lt(tail(sel.binom$rel.kl, n=1),
+              1)
+    expect_lt(tail(sel.binom$rel.kl.null, n=1),
+              1)
+    expect_equal(sel.binom$elpd[2],
+                 -35.3793754, tolerance=tol)
+    expect_equal(tail(sel.binom$elpd, n=1),
+                 -34.0411588, tolerance=tol)
+    expect_equal(sel.binom$delta.elpd[2],
+                 -1.39197661, tolerance=tol)
+    expect_true(all(diff(sel.binom$kl) < 0))
+})
+
+test_that("projsel for a model with interaction terms",
+{
+    SW({
+        sel.inter <- projsel(hs.inter, start.from="X1:X3")
+    })
+    expect_equal(attr(sel.inter, "start.from"),
+                 c("X1", "X3", "X1:X3"))
+})
+
+test_that("projsel from the intercept-only model",
+{
+    SW({
+        sel.gauss.1 <- projsel(hs.gauss, start.from=character(0))
+        sel.gauss.2 <- projsel(hs.gauss, start.from="X2")
+    })
+
+    expect_equal(sel.gauss.1$var[1:2],
+                 c("Intercept only", "X2"))
+    expect_equal(sel.gauss.2$var[1:2],
+                 c("Intercept only", "Initial submodel"))
+    expect_equivalent(sel.gauss.1[-2, ],
+                      sel.gauss.2[-2, ])
+    expect_equal(nrow(sel.gauss.1),
+                 length(c(hs.gauss$betas$unpenalized, hs.gauss$betas$penalized)))
+    expect_equal(attr(sel.gauss.1, "start.from"),
+                 character(0))
+    expect_equal(attr(sel.gauss.2, "start.from"),
+                 "X2")
+})
+
+test_that("projsel from a non-default starting submodel",
+{
+    SW({
+        sel.gauss <- projsel(hs.gauss, start.from="X1")
+    })
+
+    expect_equal(nrow(sel.gauss),
+                 length(c(grep("X1", names(hs.gauss$beta$unpenalized), inv=TRUE),
+                          hs.gauss$betas$penalized)) + 1)
+})
+
+test_that("projsel from a submodel that includes all variables",
+{
+    SW({
+        sel.gauss <- projsel(hs.gauss, start.from=paste0("X", 1:P))
+    })
+
+    expect_equal(sel.gauss$var,
+                 c("Intercept only", "Initial submodel"))
+    expect_equal(sel.gauss$rel.kl.null,
+                 c(0, 1))
+    expect_equal(sel.gauss$rel.kl,
+                 c(NA, 1))
+    expect_equal(sel.gauss$delta.elpd[2],
+                 0)
+    expect_equal(attr(sel.gauss, "row.names"),
+                 1:2)
+    expect_equal(attr(sel.gauss, "start.from"),
+                 paste0("X", 1:P))
+})
+
+test_that("projsel for a model with few observations",
+{
+    SW({
+        hs.small <- hsstan(df[1:5, ], y.gauss ~ X2 + X3, pen,
+                           iter=iters,  chains=chains, refresh=0)
+        expect_message(sel <- projsel(hs.small),
+                      "Fully saturated model reached")
+    })
+    expect_equal(tail(sel$kl, 1),
+                 0)
+})
+
+test_that("projsel for a cross-validated object",
+{
+    SW({
+        sel.gauss <- projsel(cv.gauss$fits[[1]])
+    })
+
+    expect_equal(tail(sel.gauss$elpd, n=1),
+                 sum(colMeans(log_lik(cv.gauss$fits[[1]],
+                                      newdata=df[folds == 2, ]))))
+})
+
+test_that("projsel with invalid inputs",
+{
+    expect_error(projsel(hs.gauss, max.iters=-1),
+                 "must be a non-negative integer")
+    expect_error(projsel(hs.gauss, max.iters=2.3),
+                 "must be a non-negative integer")
+    expect_error(projsel(hs.gauss, max.iters=iris),
+                 "must be a non-negative integer")
+})