Nothing
context("Utilities: predict")
## run the example from predict.gMAP
suppressMessages(example("predict.gMAP", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE))
## check that we got for each input data item a prediction
test_that("correct # of predictions are generated", expect_equal(nrow(map$data), ncol(samp)))
## check that the predictive distribution has a variance which is
## larger in accordance to the betwee-trial heterogeniety (needs to be
## done on the link scale)
test_that("variances have correct ordering", {
pred_cov_link <- predict(map, type="link")
within_var <- (summary(pred_cov_link)[,"sd"])^2
pred_cov_link_pred <- predict(map, trans_cov, type="link")
pred_var_pred <- summary(pred_cov_link_pred)[,"sd"]
tau_est <- summary(map)$tau[,"mean"]
## the predictive must include between and within; as such it is
## larger than within
expect_true(all(pred_var_pred > tau_est))
## ensure that predictive has larger variance than the model estimate
expect_true(all(summary(pred_cov_link_pred)[,"sd"] > summary(pred_cov_link)[,"sd"]))
})
## new prediction was done for a single data item
test_that("correct # of new predictions are generated", expect_equal(ncol(pred_new), 1))
## must have larger sd than between-trial alone (on link scale)
test_that("predictive variances have correct ordering",{
pred_new_link <- predict(map, data.frame(country="CH", study=11), type="link")
tau_est <- summary(map)$tau[,"mean"]
expect_true(summary(pred_new_link)[,"sd"] > tau_est)
})
## whenever the same study/covariate combination is requested, then
## the MAP must be numerically exactly the same. This ensures that per
## study the random effect is sampled just once in each iteration.
test_that("predictive distributions for the same study & covariate must match exactly", {
trans_cov_new <- data.frame(study="new", n=50, r=0, country=levels(trans_cov$country)[c(1,1)])
post_trans <- as.matrix(predict(map, newdata=trans_cov_new))
expect_equal(post_trans[,1], post_trans[,2])
})
context("Utilities: (auto)mixfit")
test_that("automixfit attempts K=4 different models and returns best fitting", {
auto_map <- automixfit(map, Nc=1:4, k=6)
models <- attr(auto_map, "models")
expect_equal(length(models), 4)
perf <- sapply(models, AIC, k=6)
## ensure that performance is decreasing
expect_true(all(diff(perf) > 0))
expect_true("betaMix" %in% class(auto_map))
})
test_that("mixfit for prediction handles response and link scale", {
pred_map <- mixfit(pred_new, Nc=2)
expect_true(is.list(pred_map))
expect_true("betaMix" %in% class(pred_map[[1]]))
expect_equal(ncol(pred_map[[1]]), 2)
pred_new_link <- predict(map, data.frame(country="CH", study=11), type="link")
pred_map_link <- mixfit(pred_new_link, Nc=2)
expect_true(is.list(pred_map_link))
expect_true("normMix" %in% class(pred_map_link[[1]]))
expect_equal(ncol(pred_map_link[[1]]), 2)
})
context("Utilities: mixcombine")
example("mixcombine", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE)
test_that("combination of mixtures", {
m1 <- mixcombine(bm, unif, weight=c(9, 1))
m2 <- mixcombine(bm, unif, unif, weight=c(8, 1, 1))
expect_equivalent(m1[1,], c(bm[1,] - 0.1/2, 0.1))
expect_equivalent(m1[2:3,1:2], bm[2:3,1:2])
expect_equivalent(m2[2:3,1:2], bm[2:3,1:2])
})
test_that("throws an error if more weights than mixtures given", {
## giving 3 weights but only 2 mixtures must not work
expect_error(mixcombine(bm, unif, weight=c(8, 1, 1)), "length(weight) not equal to length(comp)", fixed=TRUE)
})
test_that("combination of normal mixtures without default sigma works", {
norm_ui <- mixnorm(c(1, 0, 2))
norm_ui_mix <- mixcombine(norm_ui, norm_ui, weight=c(0.5,0.5))
expect_true(ncol(norm_ui_mix) == 2)
})
context("Utilities: robustify")
example("robustify", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE)
test_that("beta mixture is robustified with Beta(1,1)", {
expect_equal(ncol(bmix)+1, ncol(rbmix))
expect_equivalent(rbmix[,ncol(rbmix)], c(0.1, 1, 1))
})
test_that("beta mixture is robustified with Beta(0.5,0.5)", {
rbmix2 <- robustify(bmix, w=0.1, n=0)
expect_equal(ncol(bmix)+1, ncol(rbmix2))
expect_equivalent(rbmix2[,ncol(rbmix2)], c(0.1, 0.5, 0.5))
})
test_that("gamma mixture is robustified with n=1 equivalent prior", {
m <- summary(gmnMix)["mean"]
nr <- ncol(rgmnMix)
expect_equivalent(rgmnMix[[nr, rescale=TRUE]], mixgamma(c(1, m, 1), param="mn"))
expect_equal(rgmnMix[1,nr], 0.1)
})
test_that("gamma mixture is robustified with n=5 equivalent prior", {
m <- summary(gmnMix)["mean"]
rgmnMix2 <- robustify(gmnMix, w=0.1, n=5)
nr <- ncol(rgmnMix2)
expect_equivalent(rgmnMix2[[nr, rescale=TRUE]], mixgamma(c(1, m, 5), param="mn"))
expect_equal(rgmnMix2[1,nr], 0.1)
})
test_that("normal mixture is robustified with n=1 equivalent prior", {
nr <- ncol(rnMix)
expect_equivalent(rnMix[[nr, rescale=TRUE]], mixnorm(c(1, 0, 1), param="mn", sigma=sigma(nm)))
expect_equal(rnMix[1,nr], 0.1)
})
test_that("normal mixture is robustified with n=5 equivalent prior", {
rnMix2 <- robustify(nm, w=0.1, mean=0, n=5, sigma=sigma(nm))
nr <- ncol(rnMix2)
expect_equivalent(rnMix2[[nr, rescale=TRUE]], mixnorm(c(1, 0, 5), param="mn", sigma=sigma(nm)))
expect_equal(rnMix2[1,nr], 0.1)
})
context("Utilities: Plotting of Mixtures")
test_that("plotting of normal mixtures without default sigma works", {
norm_ui <- mixnorm(c(1, 0, 2))
norm_mix_ui <- mixcombine(norm_ui, norm_ui, weight=c(0.5,0.5))
pl <- plot(norm_mix_ui)
expect_true(inherits(pl, "ggplot"))
})
context("Utilities: Mixture Effective Sample Size")
example("ess", package="RBesT", echo=FALSE, ask=FALSE, verbose=FALSE)
test_that("conjugate beta case matches canonical formula", {
expect_equal(a+b, ess(prior, "moment"))
expect_equal(a+b, round(ess(prior, "morita")))
expect_equal(a+b, ess(prior, "elir"))
})
test_that("ess elir for beta mixtures gives a warning for a<1 & b<1 densities", {
unconstrain1 <- mixbeta(c(0.95, 10, 5), c(0.05, 0.9, 2))
unconstrain2 <- mixbeta(c(0.95, 10, 5), c(0.05, 2, 0.9))
expect_error(ess(unconstrain1, "elir"), "At least one parameter of the beta mixtures is less than 1")
expect_error(ess(unconstrain2, "elir"), "At least one parameter of the beta mixtures is less than 1")
## this one can trigger errors if the integration is not setup properly
constrained <- mixbeta(c(0.48, 1, 11), c(0.34, 6.9, 173), c(0.18, 1.0, 1.13))
expect_numeric(ess(constrained, "elir"), lower=0, finite=TRUE, any.missing=FALSE, len=1)
})
test_that("ess elir for normal mixtures returns correct values", {
mix <- mixnorm( inf1=c(0.5026,-191.1869,127.4207),inf2=c(0.2647,-187.5895,31.6130),inf3=c(0.2326,-184.7445,345.3849), sigma=270.4877)
expect_gt(ess(mix), 0)
})
test_that("moment matching for beta mixtures is correct", {
expect_equal(ess(bmix, method="moment"), sum(ab_matched))
})
test_that("normal mixtures have reference scale used correctly", {
nmix_sigma_small <- nmix
sigma_large <- RBesT::sigma(nmix)
sigma(nmix_sigma_small) <- sigma_large/sqrt(2)
suppressMessages(e1m <- ess(nmix, "moment"))
suppressMessages(e2m <- ess(nmix_sigma_small, "moment"))
expect_gt(e1m, e2m)
expect_equal(floor(abs(e2m - e1m/2)), 0)
suppressMessages(e1b <- ess(nmix, "morita"))
suppressMessages(e2b <- ess(nmix_sigma_small, "morita"))
expect_gt(e1b, e2b)
expect_equal(floor(abs(e2b - e1b/2)), 0)
suppressMessages(e1r <- ess(nmix, "elir"))
suppressMessages(e2r <- ess(nmix_sigma_small, "elir"))
expect_gt(e1r, e2r)
expect_equal(floor(abs(e2r - e1r/2)), 0)
})
test_that("gamma mixtures have likelihood property respected", {
gmix1 <- gmix
likelihood(gmix1) <- "poisson"
gmix2 <- gmix
likelihood(gmix2) <- "exp"
e1m <- ess(gmix1, "moment")
e2m <- ess(gmix2, "moment")
expect_true(e1m != e2m)
e1b <- ess(gmix1, "morita")
e2b <- ess(gmix2, "morita")
expect_true(e1b != e2b)
e1r <- ess(gmix1, "morita")
e2r <- ess(gmix2, "morita")
expect_true(e1r != e2r)
})
test_that("gamma 1-component density gives canonical results", {
guni1 <- gmix[[1, rescale=TRUE]]
likelihood(guni1) <- "poisson"
guni2 <- gmix[[1, rescale=TRUE]]
likelihood(guni2) <- "exp"
e1m <- ess(guni1, "moment")
e2m <- ess(guni2, "moment")
expect_true(e1m != e2m)
expect_equal(guni1[3,1], e1m)
expect_equal(guni2[2,1], e2m)
e1b <- round(ess(guni1, "morita"))
e2b <- round(ess(guni2, "morita"))
expect_true(e1b != e2b)
expect_equal(guni1[3,1], e1b)
expect_equal(guni2[2,1], e2b)
e1r <- ess(guni1, "elir")
e2r <- ess(guni2, "elir")
expect_true(e1r != e2r)
expect_true(abs(guni1[3,1] - e1r) < 1E-4)
## ELIR gives a-1 as ESS
expect_true(abs(guni2[2,1] - (e2r+1)) < 1E-4)
})
## check predictive consistency of ELIR
elir_predictive_consistent <- function(dens, m, Nsim, seed, stat, ...) {
## simulated from predictve which is m events equivalent to
suppressMessages(pdens <- preddist(dens, n=m))
set.seed(seed)
psamp <- rmix(pdens, Nsim)
if(inherits(dens, "gammaMix"))
psamp <- psamp / m
posterior_ess <- function(mix, method, stat, ...) {
args <- c(list(priormix=mix, stat=0), list(...))
names(args)[2] <- stat
fn <- function(x) {
args[[stat]] <- x
suppressMessages(res <- ess(do.call(postmix, args), method=method))
res
}
Vectorize(fn)
}
## obtain ess of each posterior
pred_ess <- posterior_ess(dens, "elir", stat, ...)
ess_psamp <- pred_ess(psamp)
suppressMessages(elir_prior <- ess(dens, "elir"))
## the average over the predicitve of the posterior ESS must match
## the the elir value taken directly (when m is subtracted, of
## course)
elir_pred <- mean(ess_psamp) - m
expect_true(abs(elir_prior - elir_pred) < 0.75)
}
test_that("ESS elir is predictively consistent for normal mixtures", {
skip_on_cran()
nmix <- mixnorm(rob=c(0.5, 0, 2), inf=c(0.5, 3, 4), sigma=10)
elir_predictive_consistent(nmix, m=3E2, Nsim=1E3, seed=3435, stat="m", se=10/sqrt(3E2))
})
test_that("ESS elir is predictively consistent for beta mixtures", {
skip_on_cran()
bmix <- mixbeta(rob=c(0.2, 1, 1), inf=c(0.8, 10, 2))
elir_predictive_consistent(bmix, m=1E2, Nsim=1E3, seed=355435, stat="r", n=1E2)
})
test_that("ESS elir is predictively consistent for gamma mixtures (Poisson likelihood)", {
skip_on_cran()
gmixP <- mixgamma(rob=c(0.3, 20, 4), inf=c(0.7, 50, 10), likelihood="poisson")
elir_predictive_consistent(gmixP, m=1E2, Nsim=1E3, seed=355435, stat="m", n=1E2)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.