## Testing internal functions and testing functions in other packages
## are both considered to be poor form. So this is something that we
## might re-address before release, but it will be good for now.
context("diversitree (internals)")
data(geospiza)
dat <- suppressWarnings(treedata(geospiza$phy, geospiza$dat))
phy <- dat$phy
states <- dat$dat[,"wingL"]
lik.bm <- make.bm(phy, states)
lik.ou <- make.ou(phy, states)
lik.eb <- make.eb(phy, states)
lik.ou.opt <- make.ou(phy, states, with.optimum=TRUE,
control=list(method="pruning"))
fit.bm <- find.mle(lik.bm, .1)
fit.ou <- find.mle(lik.ou, c(.1, .1))
fit.eb <- find.mle(lik.eb, c(.1, 0))
## TODO: Potential problem when dealing with renamed names; won't
## translate properly. Potentially, find.mle should have dealt with
## this for us.
fit.ou.opt <- find.mle(lik.ou.opt, c(coef(fit.ou), theta=mean(states)))
set.seed(1)
samples.bm <- mcmc(lik.bm, coef(fit.bm), 100, w=0.1, print.every=0)
set.seed(1)
samples.ou <- mcmc(lik.ou, coef(fit.ou), 30,
w=c(0.1, 10),
lower=c(0, 0),
upper=c(Inf, 100),
print.every=0)
set.seed(1)
samples.eb <- mcmc(lik.eb, coef(fit.eb), 30,
w=c(.2, 10), lower=c(0, -Inf), upper=c(Inf, 0),
print.every=0)
set.seed(1)
samples.ou.opt <- mcmc(lik.ou.opt, coef(fit.ou.opt), 30,
w=c(0.1, 10, 10),
lower=c(0, 0, -10),
upper=c(Inf, 100, 10),
print.every=0)
test_that("Model types are correct", {
expect_that(model_type(fit.bm), is_identical_to("BM"))
expect_that(model_type(fit.ou), is_identical_to("OU"))
expect_that(model_type(fit.eb), is_identical_to("EB"))
expect_that(model_type(fit.ou.opt), is_identical_to("OU"))
expect_that(model_type(samples.bm), is_identical_to("BM"))
expect_that(model_type(samples.ou), is_identical_to("OU"))
expect_that(model_type(samples.eb), is_identical_to("EB"))
expect_that(model_type(samples.ou.opt), is_identical_to("OU"))
})
test_that("Models return their source data", {
cmp <- list(phy=phy, data=states[phy$tip.label])
expect_that(model_data(fit.bm), equals(cmp))
expect_that(model_data(fit.ou), equals(cmp))
expect_that(model_data(fit.eb), equals(cmp))
expect_that(model_data(fit.ou.opt), equals(cmp))
expect_that(model_data(samples.bm), equals(cmp))
expect_that(model_data(samples.ou), equals(cmp))
expect_that(model_data(samples.eb), equals(cmp))
expect_that(model_data(samples.ou.opt), equals(cmp))
})
test_that("Processed coefficient names are as expected", {
expect_that(names(model_pars(fit.bm)),
is_identical_to(parnames.bm()))
expect_that(names(model_pars(fit.ou)),
is_identical_to(parnames.ou()))
expect_that(names(model_pars(fit.eb)),
is_identical_to(parnames.eb()))
expect_that(names(model_pars(fit.ou.opt)),
is_identical_to(parnames.ou()))
expect_that(model_pars(fit.bm)$SE, equals(0))
expect_that(model_pars(fit.ou)$SE, equals(0))
expect_that(model_pars(fit.eb)$SE, equals(0))
## MCMC:
expect_that(names(model_pars(samples.bm)),
is_identical_to(parnames.bm()))
expect_that(names(model_pars(samples.ou)),
is_identical_to(parnames.ou()))
expect_that(names(model_pars(samples.eb)),
is_identical_to(parnames.eb()))
expect_that(names(model_pars(samples.ou.opt)),
is_identical_to(parnames.ou()))
expect_that(model_pars(samples.bm)$SE,
equals(rep(0, nrow(samples.bm))))
expect_that(model_pars(samples.ou)$SE,
equals(rep(0, nrow(samples.ou))))
expect_that(model_pars(samples.eb)$SE,
equals(rep(0, nrow(samples.eb))))
})
test_that("MCMC extraction options are passed from model_pars", {
n <- nrow(samples.bm)
expect_that(nrow(model_pars(samples.bm)), equals(n))
nb <- 10
expect_that(nrow(model_pars(samples.bm, burnin=nb)),
equals(nrow(coef(samples.bm, burnin=nb))))
pb <- 0.14234
expect_that(nrow(model_pars(samples.bm, burnin=pb)),
equals(nrow(coef(samples.bm, burnin=pb))))
thin <- 7
expect_that(nrow(model_pars(samples.bm, thin=thin)),
equals(nrow(coef(samples.bm, thin=thin))))
sample <- 10
set.seed(1)
p1 <- model_pars(samples.bm, burnin=nb, thin=thin, sample=sample)
set.seed(1)
p2 <- coef(samples.bm, burnin=nb, thin=thin, sample=sample)
expect_that(nrow(p1), equals(nrow(p2)))
expect_that(p1$sigsq, is_identical_to(as.numeric(p2[,"s2"])))
})
test_that("Overall processed object looks legit", {
obj.bm <- model_info(fit.bm)
obj.ou <- model_info(fit.ou)
obj.eb <- model_info(fit.eb)
obj.names <- c("data", "pars", "type")
expect_that(names(obj.bm), is_identical_to(obj.names))
expect_that(names(obj.ou), is_identical_to(obj.names))
expect_that(names(obj.eb), is_identical_to(obj.names))
obj.bm <- model_info(samples.bm)
obj.ou <- model_info(samples.ou)
obj.eb <- model_info(samples.eb)
expect_that(names(obj.bm), is_identical_to(obj.names))
expect_that(names(obj.ou), is_identical_to(obj.names))
expect_that(names(obj.eb), is_identical_to(obj.names))
})
test_that("MCMC extraction options are passed from model_info", {
n <- nrow(samples.bm)
expect_that(nrow(model_info(samples.bm)$pars), equals(n))
nb <- 10
pb <- 0.14234
thin <- 7
sample <- 10
expect_that(nrow(model_info(samples.bm, burnin=nb)$pars),
equals(nrow(coef(samples.bm, burnin=nb))))
expect_that(nrow(model_info(samples.bm, burnin=pb)$pars),
equals(nrow(coef(samples.bm, burnin=pb))))
expect_that(nrow(model_info(samples.bm, thin=thin)$pars),
equals(nrow(coef(samples.bm, thin=thin))))
set.seed(1)
p1 <- model_info(samples.bm, burnin=nb, thin=thin, sample=sample)$pars
set.seed(1)
p2 <- coef(samples.bm, burnin=nb, thin=thin, sample=sample)
expect_that(nrow(p1), equals(nrow(p2)))
expect_that(p1$sigsq, is_identical_to(as.numeric(p2[,"s2"])))
})
## Constrained functions
test_that("Constrained functions work", {
lik.ou.c <- constrain(lik.ou.opt, theta ~ 0)
set.seed(1)
fit.ou.c <- find.mle(lik.ou.c, c(.1, .1))
samples.ou.c <- mcmc(lik.ou.c, coef(fit.ou.c), 30,
w=c(0.1, 10),
lower=c(0, 0),
upper=c(Inf, 100),
print.every=0)
## Can get the type out:
expect_that(model_type(fit.ou.c), is_identical_to("OU"))
expect_that(model_type(samples.ou.c),
is_identical_to("OU"))
## Underlying data:
cmp <- list(phy=phy, data=states[phy$tip.label])
expect_that(model_data(fit.ou.c), equals(cmp))
expect_that(model_data(samples.ou.c), equals(cmp))
## Parameters -- these need to include all parameters:
p <- model_pars(fit.ou.c)
expect_that(names(p), is_identical_to(parnames.ou()))
expect_that(p$SE, equals(0))
expect_that(p$z0, equals(0)) # z0 is theta, translated
p <- model_pars(samples.ou.c)
expect_that(names(p), is_identical_to(parnames.ou()))
expect_that(p$SE, equals(rep(0, nrow(samples.ou))))
expect_that(p$z0, equals(rep(0, nrow(samples.ou))))
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.