Nothing
test_that("optim", {
set.seed(1289)
## parameters
ntips <- 20
rateTest <- 0.1
## Tree
tree <- rphylo(ntips, 0.1, 0)
tree$edge.length <- tree$edge.length / max(vcv(tree))
## root.value
mu <- 1
disp <- 0.1
## Sim
trait <- simulateTipsCauchy(tree, mu, disp)
## Opt
res <- fitCauchy(tree, trait, method = "fixed.root")
res
capture_output(print(res))
## Change init
expect_equal(res$logLik,
fitCauchy(tree, trait, method = "fixed.root", method.init.disp = "Sn")$logLik)
expect_equal(res$logLik,
fitCauchy(tree, trait, method = "fixed.root", method.init.disp = "MAD")$logLik)
expect_equal(res$logLik,
fitCauchy(tree, trait, method = "fixed.root", method.init.disp = "IQR")$logLik)
## Compare with true value
expect_true(logDensityTipsCauchy(tree, trait, mu, disp, method = "fixed.root") <= res$logLik)
## Global
resGlob <- fitCauchy(tree, trait, method = "fixed.root", optim = "global")
expect_equal(res$logLik, resGlob$logLik)
## vcov
expect_equal(colnames(vcov(res)), c("x0", "disp"))
expect_equal(unname(sqrt(diag(vcov(res)))), c(0.03243854, 0.01808950), tolerance = 1e-2)
## lambda
reslam <- fitCauchy(tree, trait, method = 'fixed.root', root.edge = root.edge, model = "lambda")
expect_true(reslam$logLik >= res$logLik)
reslamglob <- fitCauchy(tree, trait, method = 'fixed.root', root.edge = root.edge, model = "lambda", optim = "global")
expect_equal(reslam$logLik, reslamglob$logLik)
## vcov
expect_equal(colnames(vcov(reslam)), c("x0", "disp", "lambda"))
expect_equal(unname(diag(vcov(reslam))), c(1.457096e-03, 6.289735e-04, 1.100595e-05), tolerance = 1e-2)
})
test_that("optim, random root", {
set.seed(1289)
## parameters
ntips <- 20
rateTest <- 0.1
## Tree
tree <- rphylo(ntips, 0.1, 0)
tree$edge.length <- tree$edge.length / max(vcv(tree))
## params
disp <- 0.1
root.edge <- 100
## Sim
mu <- rcauchy(1, 0, disp * root.edge)
trait <- simulateTipsCauchy(tree, mu, disp)
## Opt
res <- fitCauchy(tree, trait, method = 'random.root', root.edge = root.edge)
res
capture_output(print(res))
## Change init
expect_equal(res$logLik,
fitCauchy(tree, trait, method = "random.root", method.init.disp = "Sn")$logLik)
expect_equal(res$logLik,
fitCauchy(tree, trait, method = "random.root", method.init.disp = "IQR")$logLik)
## Compare with true value
treerr <- res$phy
treerr$root.edge <- root.edge
expect_true(logDensityTipsCauchy(treerr, trait, 0, disp, method = 'random.root') <= res$logLik)
## Global
resGlob <- fitCauchy(tree, trait, method = 'random.root', root.edge = root.edge, optim = "global")
expect_equal(res$logLik, resGlob$logLik)
## vcov
expect_equal(colnames(vcov(res)), c("disp"))
expect_equal(unname(sqrt(diag(vcov(res)))), c(0.021), tolerance = 1e-1)
## lambda
reslam <- fitCauchy(tree, trait, method = 'random.root', root.edge = root.edge, model = "lambda")
capture_output(print(reslam))
expect_true(reslam$logLik > res$logLik)
## vcov
expect_equal(colnames(vcov(reslam)), c("disp", "lambda"))
expect_equal(unname(sqrt(diag(vcov(reslam)))), c(0.0285, 0.0415), tolerance = 1e-1)
## lambda fixed
reslam2 <- fitCauchy(tree, trait, method = 'random.root', root.edge = root.edge, model = "lambda",
starting.value = list(lambda = 1), lower.bound = list(lambda = 1), upper.bound = list(lambda = 1))
expect_equal(reslam2$lambda, 1.0, ignore_attr = TRUE)
expect_equal(reslam2$logLik, res$logLik)
})
test_that("optim, reml", {
set.seed(1289)
## parameters
ntips <- 20
rateTest <- 0.1
## Tree
tree <- rphylo(ntips, 0.1, 0)
tree$edge.length <- tree$edge.length / max(vcv(tree))
## params
disp <- 0.1
root.edge <- 1
## Sim
mu <- rcauchy(1, 0, disp * root.edge)
trait <- simulateTipsCauchy(tree, mu, disp)
## Opt
res <- fitCauchy(tree, trait, method = 'reml')
res
capture_output(print(res))
## Change init
expect_equal(res$logLik,
fitCauchy(tree, trait, method = "reml", method.init.disp = "Sn")$logLik)
expect_equal(res$logLik,
fitCauchy(tree, trait, method = "reml", method.init.disp = "MAD")$logLik)
expect_equal(res$logLik,
fitCauchy(tree, trait, method = "reml", method.init.disp = "IQR")$logLik)
## Compare with true value
expect_true(logDensityTipsCauchy(res$phy, trait, NULL, disp, method = "reml") <= res$logLik)
# logDensityTipsCauchy(res$phy, trait, NULL, res$disp, method = "reml")
# dd <- seq(0.05, 0.2, 0.005)
# plot(dd,
# sapply(dd, function(d) logDensityTipsCauchy(res$phy, trait, NULL, d, method = "reml")))
## Global
resGlob <- fitCauchy(tree, trait, method = 'reml', optim = "global")
expect_equal(res$logLik, resGlob$logLik)
## vcov
expect_equal(colnames(vcov(res)), c("disp"))
expect_equal(unname(sqrt(diag(vcov(res)))), c(0.02016588), tolerance = 1e-2)
## lambda
reslam <- fitCauchy(tree, trait, method = 'reml', model = "lambda")
expect_true(reslam$logLik > res$logLik)
## vcov
expect_equal(colnames(vcov(reslam)), c("disp", "lambda"))
expect_equal(unname(sqrt(diag(vcov(reslam)))), c(0.0283, 0.0574), tolerance = 1e-2)
# ## comparison with random root
# root.edge <- 1000000000000
# resran <- fitCauchy(tree, trait, method = 'random.root', root.edge = root.edge)
# expect_equal(resran$disp, res$disp)
})
test_that("cauphylm", {
set.seed(1289)
## parameters
ntips <- 20
rateTest <- 0.1
## Tree
tree <- rphylo(ntips, 0.1, 0)
tree$edge.length <- tree$edge.length / max(vcv(tree))
## root.value
mu <- 1
disp <- 0.1
## Sim
trait <- simulateTipsCauchy(tree, mu, disp)
## Reg
reg <- ape::rTraitCont(tree, model = "BM", sigma = 0.1, root.value = 0)
dat <- data.frame(tt = trait, rr = reg)
## Opt
reslm <- cauphylm(trait ~ 1, phy = tree)
reslmdat <- cauphylm(tt ~ 1, data = dat, phy = tree)
res <- fitCauchy(tree, trait, method = "fixed.root")
capture_output(print(res))
## Compare fits
expect_equal(res$logLik, reslm$logLik)
expect_equal(res$aic, reslm$aic)
expect_equal(res$x0, reslm$coefficients, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(res$disp, reslm$disp, ignore_attr = TRUE, tolerance = 1e-3)
expect_equal(reslmdat[!(names(reslmdat) %in% c("formula", "call"))],
reslm[!(names(reslm) %in% c("formula", "call"))])
## Compare vcov
expect_equal(unname(vcov(reslm)), unname(vcov(res)), tolerance = 1e-3)
## Regression
trait <- -3 + 2 * reg + simulateTipsCauchy(tree, 0, disp)
dat <- data.frame(trait = trait, reg = reg)
reslm <- cauphylm(trait ~ reg, data = dat, phy = tree)
reslmnull <- cauphylm(trait ~ 1, data = dat, phy = tree)
expect_true(reslmnull$logLik <= reslm$logLik)
})
test_that("cauphylm lambda", {
set.seed(1289)
## parameters
ntips <- 20
rateTest <- 0.1
## Tree
tree <- rphylo(ntips, 0.1, 0)
tree$edge.length <- tree$edge.length / max(vcv(tree))
tree_lambda <- phylolm::transf.branch.lengths(tree, model = "lambda", parameters = list(lambda = 0.6))$tree
## root.value
mu <- 1
disp <- 0.1
## Sim
trait <- simulateTipsCauchy(tree_lambda, mu, disp)
## Reg
reg <- ape::rTraitCont(tree, model = "BM", sigma = 0.1, root.value = 0)
dat <- data.frame(tt = trait, rr = reg)
## Opt
reslm <- cauphylm(trait ~ 1, phy = tree, model = "lambda")
reslmdat <- cauphylm(tt ~ 1, data = dat, phy = tree, model = "lambda")
res <- fitCauchy(tree, trait, method = "fixed.root", model = "lambda")
capture_output(print(reslm))
## Compare fits
expect_equal(reslmdat$logLik, reslm$logLik)
expect_equal(reslmdat$coefficient, reslm$coefficients, ignore_attr = TRUE)
expect_equal(reslmdat$disp, reslm$disp, ignore_attr = TRUE)
expect_equal(reslmdat$lambda, reslm$lambda, ignore_attr = TRUE)
expect_equal(res$logLik, reslm$logLik)
expect_equal(res$aic, reslm$aic)
expect_equal(res$x0, reslm$coefficients, ignore_attr = TRUE, tolerance = 1e-4)
expect_equal(res$disp, reslm$disp, ignore_attr = TRUE, tolerance = 1e-3)
expect_equal(res$lambda, reslm$lambda, ignore_attr = TRUE, tolerance = 1e-3)
expect_equal(reslmdat[!(names(reslmdat) %in% c("formula", "call"))],
reslm[!(names(reslm) %in% c("formula", "call"))])
## Compare without lambda
reslmcau <- cauphylm(trait ~ 1, phy = tree, model = "cauchy")
expect_true(reslmcau$logLik <= reslm$logLik)
## Compare vcov
expect_equal(unname(vcov(reslm)), unname(vcov(res)), tolerance = 1e-2)
## Regression
trait <- -3 + 2 * reg + rTraitCauchy(1, tree, "lambda", parameters = list(root.value = 0, disp = disp, lambda = 0.5))
dat <- data.frame(trait = trait, reg = reg)
reslmstart <- cauphylm(trait ~ reg, data = dat, phy = tree, model = "lambda", starting.value = list(lambda = 1))
reslm <- cauphylm(trait ~ reg, data = dat, phy = tree, model = "lambda")
reslmnull <- cauphylm(trait ~ 1, data = dat, phy = tree, model = "lambda")
reslmcau <- cauphylm(trait ~ reg, data = dat, phy = tree, model = "cauchy")
expect_true(reslmnull$logLik <= reslm$logLik)
expect_true(reslmcau$logLik <= reslmstart$logLik)
# expect_true(reslmstart$logLik <= reslm$logLik)
})
test_that("cauphylm helper functions", {
set.seed(1289)
## parameters
ntips <- 20
rateTest <- 0.1
## Tree
tree <- rphylo(ntips, 0.1, 0)
tree$edge.length <- tree$edge.length / max(vcv(tree))
tree_lambda <- phylolm::transf.branch.lengths(tree, model = "lambda", parameters = list(lambda = 0.6))$tree
## root.value
mu <- 1
disp <- 0.1
## Sim
trait <- simulateTipsCauchy(tree_lambda, mu, disp)
## Reg
reg <- ape::rTraitCont(tree, model = "BM", sigma = 0.1, root.value = 0)
dat <- data.frame(tt = trait, rr = reg)
## Opt
reslmdat <- cauphylm(tt ~ 1, data = dat, phy = tree, model = "lambda", hessian = TRUE)
## Functions
reslmdat
logLik(reslmdat)
expect_equal(logLik(reslmdat), 13.296526, ignore_attr = TRUE)
expect_equal(AIC(reslmdat), -20.5930517)
# expect_equal(extractAIC(reslmdat)[2], -20.5930517)
# expect_equal(nobs(reslmdat), ntips)
expect_equal(unique(predict(reslmdat)), unname(reslmdat$coefficients))
expect_error(predict(reslmdat, se.fit = TRUE))
expect_equal(unname(sqrt(diag(vcov(reslmdat)))), c(0.04988679, 0.06517500, 0.18057543), tolerance = 1e-2)
expect_message(cc <- confint(reslmdat, level = 0.9), "Approximated asymptotic confidence interval using the Hessian")
expect_true(cc[1, 1] <= mu && cc[1, 2] >= mu)
expect_true(cc[2, 1] <= disp && cc[2, 2] >= disp)
expect_true(cc[3, 1] <= 1.0 && cc[3, 2] >= 1.0)
expect_equal(coef(reslmdat), reslmdat$all_params)
## No Hessian
reslmdatbis <- cauphylm(tt ~ 1, data = dat, phy = tree, model = "lambda", hessian = FALSE)
expect_equal(vcov(compute_vcov(reslmdatbis)), vcov(reslmdat))
## Compare vcov with fitCauchy
resfitcau <- fitCauchy(tree, trait, model = "lambda", method = "fixed.root", hessian = TRUE)
expect_equal(vcov(resfitcau), vcov(reslmdat), ignore_attr = TRUE, tolerance = 1e-2)
expect_equal(vcov(compute_vcov(fitCauchy(tree, trait, model = "lambda", method = "fixed.root", hessian = FALSE))),
vcov(resfitcau))
expect_equal(coef(resfitcau), resfitcau$all_params)
## Functions fitCauchy
resfitcau
expect_equal(logLik(resfitcau), logLik(reslmdat), ignore_attr = TRUE)
expect_equal(AIC(resfitcau), AIC(reslmdat))
expect_equal(unname(sqrt(diag(vcov(resfitcau)))), unname(sqrt(diag(vcov(reslmdat)))), tolerance = 1e-2)
expect_message(cc <- confint(resfitcau, level = 0.9), "Approximated asymptotic confidence interval using the Hessian")
expect_true(cc[1, 1] <= mu && cc[1, 2] >= mu)
expect_true(cc[2, 1] <= disp && cc[2, 2] >= disp)
expect_true(cc[3, 1] <= 1.0 && cc[3, 2] >= 1.0)
})
test_that("profile likelihood", {
set.seed(1289)
## parameters
ntips <- 20
rateTest <- 0.1
## Tree
tree <- rphylo(ntips, 0.1, 0)
tree$edge.length <- tree$edge.length / max(vcv(tree))
tree_lambda <- phylolm::transf.branch.lengths(tree, model = "lambda", parameters = list(lambda = 0.6))$tree
## root.value
mu <- 1
disp <- 0.1
## Sim
trait <- simulateTipsCauchy(tree_lambda, mu, disp)
## Fit
fit <- fitCauchy(tree, trait, model = "lambda", method = "fixed.root")
## Profile
expect_message(pr <- profile(fit, which = c("x0", "disp", "blob")),
"Parameters blob are not in the fitted object, and will be ignored.")
expect_message(pr <- profile(fit, which = c(1, 2, 5)),
"Parameters with indexes 5 are not in the fitted object, and will be ignored.")
pr <- profile(fit)
expect_equal(max(pr$lambda$profLogLik), fit$logLik)
expect_true(mean(pr$x0$profLogLik) <= fit$logLik)
expect_true(mean(pr$disp$profLogLik) <= fit$logLik)
expect_true(mean(pr$lambda$profLogLik) <= fit$logLik)
plot(pr)
})
# test_that("testFitBig", {
#
# ## Parameters
# mu <- 0.8
# disp <- 0.1
#
# ## tree with three tips
# set.seed(1289)
# n <- 500
# tree <- rphylo(n, 0.1, 0)
# tree_height <- max(diag(vcv(tree)))
#
# ## data
# trait <- simulateTipsCauchy(tree, mu, disp)
#
# ## Opt
# reslm <- cauphylm(trait ~ 1, phy = tree)
# res <- fitCauchy(tree, trait, method = "fixed.root")
# expect_equal(res$logLik, reslm$logLik)
#
# })
test_that("Errors with species names", {
set.seed(12891026)
## Tree
ntips <- 20
tree <- ape::rphylo(ntips, 0.1, 0)
mat_tree <- ape::vcv(tree)
## data
# no names
y_data <- rnorm(ntips)
expect_error(checkTraitTree(y_data, tree, name = "trait"),
"`trait` and/or the tips of the phylogeny are not named.")
# wrong order
names(y_data) <- sample(tree$tip.label, ntips)
expect_warning(trait1 <- checkTraitTree(y_data, tree, name = "trait"),
"trait` was not sorted in the correct order")
y_data <- y_data[match(tree$tip.label, names(y_data))]
trait2 <- checkTraitTree(y_data, tree, name = "trait")
expect_equal(trait1, trait2)
# wrong names
names(y_data)[1] <- "moustache"
expect_error(checkTraitTree(y_data, tree, name = "trait"),
"Species 't1' are in the tree but not in trait.")
# Correct name and order
names(y_data) <- tree$tip.label
# wrong names tree
tree_wrong <- tree
tree_wrong$tip.label[1] <- "moustache"
expect_error(checkTraitTree(y_data, tree_wrong, name = "trait"),
"Species 'moustache' are in the tree but not in trait.")
## Design
design <- matrix(1, nrow = ntips, ncol = 2)
design[sample(1:ntips, floor(ntips / 2)), 2] <- 0
# wrong dimension
expect_error(checkTraitTree(t(design), tree, name = "design"),
"`design` should have as many rows as the number of taxa in the tree.")
# no names
expect_error(checkTraitTree(design, tree, name = "design"),
"`design` and/or the tips of the phylogeny are not named.")
# wrong order
rownames(design) <- sample(tree$tip.label, ntips)
expect_warning(res1 <- checkTraitTree(design, tree, name = "design"),
"`design` was not sorted in the correct order, when compared with the tips label.")
design <- design[match(tree$tip.label, rownames(design)), , drop = FALSE]
res2 <- checkTraitTree(design, tree, name = "design")
expect_equal(res1, res2)
# wrong names
rownames(design)[1] <- "moustache"
expect_error(checkTraitTree(design, tree, name = "design"),
"Species 't1' are in the tree but not in design.")
rownames(design) <- tree$tip.label
expect_error(checkTraitTree(design, tree_wrong, name = "design"),
"Species 'moustache' are in the tree but not in design")
## Duplicates
y_data[1] <- y_data[2]
expect_error(checkDuplicates(y_data, tree),
"The trait vector has duplicated entries on tips that are equidistant from the root")
tree$edge.length[tree$edge[, 2] == 1] <- 1.0
expect_no_error(checkDuplicates(y_data, tree))
})
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.