test_that("hsstan",

{

    expect_s3_class(hs.gauss,

                    "hsstan")

    expect_s4_class(hs.gauss$stanfit,

                    "stanfit")

    expect_named(hs.gauss,

                 c("stanfit", "betas", "call", "data", "model.terms", "family",

                   "hsstan.settings"))

    expect_false("lambda[1]" %in% names(hs.gauss$stanfit))

    expect_equal(hs.gauss$family,

                 gaussian())

    expect_named(hs.gauss$betas$unpenalized,

                 c("(Intercept)", "X1b", "X1c", "X2", "X3"))

    expect_named(hs.gauss$betas$penalized,

                 hs.gauss$model.terms$penalized)

    expect_length(hs.gauss$betas, 2)

    expect_named(hs.gauss$hsstan.settings,

                 c("adapt.delta", "qr", "seed", "scale.u", "regularized", "nu",

                   "par.ratio", "global.scale", "global.df", "slab.scale",

                   "slab.df"))

    expect_equal(hs.gauss$hsstan.settings$global.scale,

                 0.007071068)

    expect_equal(hs.gauss$hsstan.settings$adapt.delta,

                 0.99)

})

test_that("hsstan with no penalized predictors",

{

    expect_null(hs.base$betas$penalized)

    expect_length(hs.base$penalized,

                  0)

    expect_named(hs.base$hsstan.settings,

                 c("adapt.delta", "qr", "seed", "scale.u"))

    expect_equal(hs.base$hsstan.settings$adapt.delta,

                 0.95)

})

test_that("hsstan handles categorical variables in the penalized predictors",

{

    SW({

        hs.1 <- hsstan(df, y.gauss ~ X2 + X3, "X1", iter=250,

                       keep.hs.pars=TRUE, refresh=0)

    })

    expect_named(hs.1$betas$unpenalized,

                 c("(Intercept)", "X2", "X3"))

    expect_named(hs.1$betas$penalized,

                 c("X1b", "X1c"))

})

test_that("hsstan handles penalized predictors appearing in the formula",

{

    SW({

        hs.1 <- hsstan(df, y.gauss ~ X1 + X2 + X3, "X2", iter=250,

                       keep.hs.pars=TRUE, refresh=0)

        hs.2 <- hsstan(df, y.gauss ~ X1 + X3, "X2", iter=250,

                       keep.hs.pars=TRUE, refresh=0)

    })

    for (val in c("beta", "data", "model.terms"))

        expect_equal(hs.1[[val]],

                     hs.2[[val]])

})

test_that("hsstan handles interaction terms correctly",

{

    SW({

        hs.int.2 <- hs(y.gauss ~ X1 + X3 + X2 + X1b_X3 + X1c_X3 + X3_X2, gaussian)

    })

    expect_equal(names(hs.inter$betas$unpenalized),

                 gsub("_", ":", names(hs.int.2$betas$unpenalized)))

    expect_equivalent(hs.inter$betas$unpenalized,

                      hs.int.2$betas$unpenalized)

    expect_equal(hs.inter$betas$penalized,

                 hs.int.2$betas$penalized)

})

test_that("hsstan handles interaction term without main effects",

{

    SW({

        hs.int.0 <- hsstan(df, y.gauss ~ X1:X3, iter=200, refresh=0)

    })

    expect_named(hs.int.0$betas$unpenalized,

                 c("(Intercept)", "X1a:X3", "X1b:X3", "X1c:X3"))

})

test_that("hsstan doesn't use the QR decomposition if P > N",

{

    SW({

        hs.noqr <- hsstan(df[1:5, ], mod.gauss, pen, iter=100, qr=TRUE,

                          keep.hs.pars=TRUE, refresh=0)

    })

    expect_false(hs.noqr$hsstan.settings$qr)

    expect_match(names(hs.noqr$stanfit),

                 "lambda", all=FALSE)

    expect_match(names(hs.noqr$stanfit),

                 "tau", all=FALSE)

})

test_that("kfold",

{

    expect_s3_class(cv.gauss,

                    c("kfold", "loo"))

    expect_output(print(cv.gauss),

                  "Based on 2-fold cross-validation")

    expect_named(cv.gauss,

                 c("estimates", "pointwise", "fits", "data"))

    expect_equal(rownames(cv.gauss$estimates),

                 c("elpd_kfold", "p_kfold", "kfoldic"))

    expect_equal(colnames(cv.gauss$estimates),

                 c("Estimate", "SE"))

    expect_equal(nrow(cv.gauss$pointwise),

                 N)

    expect_equal(colnames(cv.gauss$pointwise),

                 c("elpd_kfold", "p_kfold", "kfoldic"))

    expect_true(all(is.na(cv.gauss$pointwise[, "p_kfold"])))

    expect_length(cv.gauss$fits[[1]]$stanfit@stan_args,

                  2)

    expect_named(cv.binom,

                 c("estimates", "pointwise", "fits", "data"))

    for (i in 1:max(folds))

        expect_s3_class(cv.binom$fits[[i]],

                        "hsstan")

    expect_silent(validate.samples(cv.binom$fits[[1]]))

    expect_equal(nrow(cv.binom$fits),

                 2)

    expect_length(cv.binom$fits[[1]]$stanfit@stan_args,

                  1)

    expect_length(cv.binom$fits,

                  max(folds) * 2)

    expect_equal(cv.binom$fits[[max(folds) + 1]],

                 which(folds == 1))

})

test_that("kfold with store.fits=FALSE",

{

    expect_named(cv.nofit,

                 c("estimates", "pointwise"))

})

test_that("hsstan with invalid inputs",

{

    expect_error(hsstan(df, mod.gauss, adapt.delta=1),

                 "'adapt.delta' must be less than 1")

    expect_error(hsstan(df, mod.gauss, iter=0),

                 "'iter' must be a positive integer")

    expect_error(hsstan(df, mod.gauss, iter=-1),

                 "'iter' must be a positive integer")

    expect_error(hsstan(df, mod.gauss, warmup=0),

                 "'warmup' must be a positive integer")

    expect_error(hsstan(df, mod.gauss, warmup=-1),

                 "'warmup' must be a positive integer")

    expect_error(hsstan(df, mod.gauss, iter=1000, warmup=1000),

                 "'warmup' must be smaller than 'iter'")

    expect_error(hsstan(df, mod.gauss, chains=0),

                 "rstan::sampling failed")

})

test_that("kfold with invalid inputs",

{

    expect_error(kfold(hs.gauss, folds, chains=0),

                 "'chains' must be a positive integer")

    expect_error(kfold(hs.gauss, folds, chains=4.4),

                 "'chains' must be a positive integer")

})