context("Function simulate with shifts")
###############################################################################
test_that("Mean of the BM", {
set.seed(586)
ntaxa <- 20
tree <- rtree(ntaxa)
## Simulate Process
p <- 3
variance <- matrix(0.2, p, p) + diag(0.3, p, p)
root.state <- list(random = FALSE,
value.root = c(1, -1, 2),
exp.root = NA,
var.root = NA)
shifts = list(edges = c(18, 32),
values=cbind(c(4, -10, 3),
c(-5, 5, 0)),
relativeTimes = 0)
X1 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "BM",
variance = variance,
shifts = shifts)
X1.tips.exp <- extract_simulate_internal(X1, where = "tips", what = "exp")
## Compute expectations with tree matrix
T_tree <- incidence.matrix(tree)
Delta <- shifts.list_to_matrix(tree, shifts)
X1.tips.exp.mat <- tcrossprod(Delta, T_tree) + root.state$value.root
expect_that(X1.tips.exp, equals(X1.tips.exp.mat))
## Compute expectations of internal nodes
U_tree <- incidence.matrix.full(tree)
Delta <- shifts.list_to_matrix(tree, shifts)
X1.all.exp.mat <- tcrossprod(Delta, U_tree) + root.state$value.root
X1.nodes.exp <- extract_simulate_internal(X1, where = "nodes", what = "exp")
X1.all.exp <- cbind(X1.tips.exp, X1.nodes.exp)
expect_that(X1.all.exp.mat, equals(X1.all.exp))
## Without simulating
X2 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "BM",
variance = variance,
shifts = shifts,
simulate_random = FALSE)
X2.tips.exp <- extract_simulate_internal(X2, where = "tips", what = "exp")
expect_that(X1.tips.exp, equals(X2.tips.exp))
})
###############################################################################
test_that("Mean of the BM - random root", {
set.seed(586)
ntaxa <- 20
tree <- rtree(ntaxa)
## Simulate Process
p <- 3
variance <- matrix(0.2, p, p) + diag(0.3, p, p)
root.state <- list(random = TRUE,
value.root = NA,
exp.root = c(1, -1, 2),
var.root = matrix(0.1, p, p) + diag(0.1, p, p))
shifts = list(edges = c(18, 32),
values=cbind(c(4, -10, 3),
c(-5, 5, 0)),
relativeTimes = 0)
X1 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "BM",
variance = variance,
shifts = shifts)
X1.tips.exp <- extract_simulate_internal(X1, where = "tips", what = "exp")
## Does adding a root edge change anything ?
tree$root.edge <- 1
X2 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "BM",
variance = variance,
shifts = shifts)
X2.tips.exp <- extract_simulate_internal(X2, where = "tips", what = "exp")
expect_that(X1.tips.exp, equals(X2.tips.exp))
## Compute expectations with tree matrix
T_tree <- incidence.matrix(tree)
Delta <- shifts.list_to_matrix(tree, shifts)
X1.tips.exp.mat <- tcrossprod(Delta, T_tree) + root.state$exp.root
expect_that(X1.tips.exp, equals(X1.tips.exp.mat))
## Compute expectations of internal nodes
U_tree <- incidence.matrix.full(tree)
Delta <- shifts.list_to_matrix(tree, shifts)
X1.all.exp.mat <- tcrossprod(Delta, U_tree) + root.state$exp.root
X1.nodes.exp <- extract_simulate_internal(X1, where = "nodes", what = "exp")
X1.all.exp <- cbind(X1.tips.exp, X1.nodes.exp)
expect_that(X1.all.exp.mat, equals(X1.all.exp))
## Without simulating
X2 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "BM",
variance = variance,
shifts = shifts,
simulate_random = FALSE,
U_tree = U_tree)
X2.tips.exp <- extract_simulate_internal(X2, where = "tips", what = "exp")
expect_that(X1.tips.exp, equals(X2.tips.exp))
})
###############################################################################
test_that("Mean of the OU", {
set.seed(1899)
ntaxa <- 32
tree <- rcoal(ntaxa)
## Simulate Process
p <- 3
variance <- matrix(0.2, p, p) + diag(0.3, p, p)
optimal.value <- c(-3, 5, 0)
selection.strength <- diag(3, p, p) + tcrossprod(c(0.1, 0.2, 0.3))
exp.stationary <- optimal.value
var.stationary <- compute_stationary_variance(variance, selection.strength)
root.state <- list(random = TRUE,
stationary.root = TRUE,
value.root = NA,
exp.root = exp.stationary,
var.root = var.stationary)
shifts = list(edges = c(11, 35, 44),
values=cbind(c(4, -10, 3),
c(-5, 5.4, 0),
c(2, -8, 0.3)),
relativeTimes = 0)
X1 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "OU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts)
X1.tips.exp <- extract_simulate_internal(X1, where = "tips", what = "exp")
## Compute expectations with tree matrix
T_tree <- incidence.matrix(tree)
Delta <- shifts.list_to_matrix(tree, shifts)
W <- compute_actualization_matrix_ultrametric(tree, selection.strength)
vec_Y <- kronecker(T_tree, diag(1, p, p)) %*% W %*% as.vector(Delta)
X1.tips.exp.mat <- matrix(vec_Y, p, ntaxa) + optimal.value
expect_that(X1.tips.exp, equals(X1.tips.exp.mat))
## Without simulate
X2 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "OU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts,
simulate_random = FALSE)
X2.tips.exp <- extract_simulate_internal(X2, where = "tips", what = "exp")
expect_that(X2.tips.exp, equals(X2.tips.exp))
})
###############################################################################
test_that("OU - fixed root", {
set.seed(1899)
ntaxa <- 32
tree <- rcoal(ntaxa)
## Simulate Process
p <- 3
variance <- matrix(0.2, p, p) + diag(0.3, p, p)
optimal.value <- c(-3, 5, 0)
selection.strength <- diag(3, p, p) + tcrossprod(c(0.1, 0.2, 0.3))
exp.stationary <- optimal.value
# var.stationary <- compute_stationary_variance(variance, selection.strength)
root.state <- list(random = FALSE,
stationary.root = FALSE,
value.root = exp.stationary,
exp.root = NA,
var.root = NA)
shifts = list(edges = c(11, 35, 44),
values=cbind(c(4, -10, 3),
c(-5, 5.4, 0),
c(2, -8, 0.3)),
relativeTimes = 0)
X1 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "OU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts)
X1.tips.exp <- extract_simulate_internal(X1, where = "tips", what = "exp")
## Compute expectations with tree matrix
T_tree <- incidence.matrix(tree)
Delta <- shifts.list_to_matrix(tree, shifts)
W <- compute_actualization_matrix_ultrametric(tree, selection.strength)
vec_Y <- kronecker(T_tree, diag(1, p, p)) %*% W %*% as.vector(Delta)
X1.tips.exp.mat <- matrix(vec_Y, p, ntaxa) + optimal.value
expect_that(X1.tips.exp, equals(X1.tips.exp.mat))
## Without simulate
X2 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "OU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts,
simulate_random = FALSE)
X2.tips.exp <- extract_simulate_internal(X2, where = "tips", what = "exp")
expect_that(X2.tips.exp, equals(X2.tips.exp))
})
###############################################################################
test_that("Multivariate Scalar (scOU)", {
testthat::skip_on_cran()
testthat::skip_if_not_installed("TreeSim")
set.seed(586)
ntaxa <- 100
tree <- TreeSim::sim.bd.taxa.age(n = ntaxa, numbsim = 1,
lambda = 1, mu = 0,
age = 1, mrca = TRUE)[[1]]
p <- 6
variance <- diag(0.5, p, p)
optimal.value <- c(-3, 5, 0, 7, 3, 0)
selection.strength <- diag(3, p, p)
exp.stationary <- optimal.value
var.stationary <- compute_stationary_variance(variance, selection.strength)
root.state <- list(random = TRUE,
stationary.root = TRUE,
value.root = NA,
exp.root = exp.stationary,
var.root = var.stationary)
shifts = list(edges = c(18, 32),
values=cbind(c(4, -10, 3, 2, 2, 9),
c(-5, 5, 0, -7, 5, -4)),
relativeTimes = 0)
## Forget that it is scalar
Xnot <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "OU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts)
## Use that it is scalar
Xsc <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "scOU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts)
expect_that(Xnot[,,2:3], equals(Xsc[,,2:3]))
## Do not simulate
Xsc2 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "scOU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts,
simulate_random = FALSE)
expect_that(Xnot[,,2:3], equals(Xsc2[,,2:3]))
})
###############################################################################
test_that("Multivariate Scalar (scOU) - Fixed Root", {
testthat::skip_on_cran()
testthat::skip_if_not_installed("TreeSim")
set.seed(586)
ntaxa <- 100
tree <- TreeSim::sim.bd.taxa.age(n = ntaxa, numbsim = 1,
lambda = 1, mu = 0,
age = 1, mrca = TRUE)[[1]]
p <- 6
variance <- diag(0.5, p, p)
optimal.value <- c(-3, 5, 0, 7, 3, 0)
selection.strength <- diag(3, p, p)
exp.stationary <- optimal.value
root.state <- list(random = FALSE,
stationary.root = FALSE,
value.root = exp.stationary,
exp.root = NA,
var.root = NA)
shifts = list(edges = c(18, 32),
values=cbind(c(4, -10, 3, 2, 2, 9),
c(-5, 5, 0, -7, 5, -4)),
relativeTimes = 0)
## Forget that it is scalar
Xnot <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "OU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts)
## Use that it is scalar
Xsc <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "scOU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts)
expect_that(Xnot[,,2:3], equals(Xsc[,,2:3]))
## Do not simulate
Xsc2 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "scOU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts,
simulate_random = FALSE)
expect_that(Xnot[,,2:3], equals(Xsc2[,,2:3]))
})
###############################################################################
test_that("Internal vs simul", {
ntaxa <- 32
tree <- rcoal(ntaxa)
## Simulate Internal
p <- 3
variance <- matrix(0.2, p, p) + diag(0.3, p, p)
optimal.value <- c(-3, 5, 0)
selection.strength <- diag(3, p, p) + tcrossprod(c(0.1, 0.2, 0.3))
exp.stationary <- optimal.value
var.stationary <- compute_stationary_variance(variance, selection.strength)
root.state <- list(random = TRUE,
stationary.root = TRUE,
value.root = NA,
exp.root = exp.stationary,
var.root = var.stationary)
shifts = list(edges = c(11, 35, 44),
values=cbind(c(4, -10, 3),
c(-5, 5.4, 0),
c(2, -8, 0.3)),
relativeTimes = 0)
set.seed(1899)
X1 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "OU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts)
## Simulate External
para <- params_process("OU", p = p, variance = variance,
selection.strength = selection.strength,
optimal.value = optimal.value, random = TRUE,
stationary.root = TRUE, exp.root = exp.stationary,
var.root = var.stationary, edges = shifts$edges,
values = shifts$values)
set.seed(1899)
X2 <- simul_process(para, phylo = tree)
## Comparisons and extractors
X1.tips.exp <- extract_simulate_internal(X1, where = "tips", what = "exp")
X2.tips.exp <- extract(X2, where = "tips", what = "exp")
expect_equal(X1.tips.exp, unname(X2.tips.exp))
X1.nodes.exp <- extract_simulate_internal(X1, where = "nodes", what = "exp")
X2.nodes.exp <- extract(X2, where = "nodes", what = "exp")
expect_equal(X1.nodes.exp, unname(X2.nodes.exp))
X1.tips.states <- extract_simulate_internal(X1, where = "tips", what = "states")
X2.tips.states <- extract(X2, where = "tips", what = "states")
expect_equal(X1.tips.states, unname(X2.tips.states))
X1.nodes.states <- extract_simulate_internal(X1, where = "nodes", what = "states")
X2.nodes.states <- extract(X2, where = "nodes", what = "states")
expect_equal(X1.nodes.states, unname(X2.nodes.states))
X1.tips.optim <- extract_simulate_internal(X1, where = "tips", what = "optim")
X2.tips.optim <- extract(X2, where = "tips", what = "optim")
expect_equal(X1.tips.optim, unname(X2.tips.optim))
X1.nodes.optim <- extract_simulate_internal(X1, where = "nodes", what = "optim")
X2.nodes.optim <- extract(X2, where = "nodes", what = "optim")
expect_equal(X1.nodes.optim, unname(X2.nodes.optim))
})
# test_that("Multivariate Independent (BM)", {
# set.seed(586)
# ntaxa <- 200
# tree <- TreeSim::sim.bd.taxa.age(n = ntaxa, numbsim = 1,
# lambda = 1, mu = 0,
# age = 1, mrca = TRUE)[[1]]
#
# p <- 3
# variance <- diag(0.5, p, p)
# root.state <- list(random = FALSE,
# stationary.root = FALSE,
# value.root = c(1, -6, 2),
# exp.root = NA,
# var.root = NA)
# shifts = list(edges = c(18, 32),
# values=cbind(c(4, -10, 3),
# c(-5, 5, 0)),
# relativeTimes = 0)
# paramsSimu <- list(variance = variance,
# selection.strength = selection.strength,
# shifts = shifts,
# root.state = root.state)
#
# Xnot <- simulate_internal(tree,
# p = p,
# independent = FALSE,
# root.state = root.state,
# process = "BM",
# variance = variance,
# optimal.value = optimal.value,
# selection.strength = selection.strength,
# shifts = shifts)
#
# Xind <- simulate_internal(tree,
# p = p,
# independent = TRUE,
# root.state = root.state,
# process = "OU",
# variance = variance,
# optimal.value = optimal.value,
# selection.strength = selection.strength,
# shifts = shifts)
#
# expect_that(Xnot[,,2], equals(Xind[,,2]))
# })
###############################################################################
test_that("OU/BM", {
set.seed(1899)
ntaxa <- 32
tree <- rcoal(ntaxa)
## Simulate Process
p <- 3
variance <- matrix(0.2, p, p) + diag(0.3, p, p)
optimal.value <- c(-3, 5, 0)
selection.strength <- diag(0:2)
exp.stationary <- optimal.value
# var.stationary <- compute_stationary_variance(variance, selection.strength)
root.state <- list(random = FALSE,
stationary.root = FALSE,
value.root = exp.stationary,
exp.root = NA,
var.root = NA)
shifts = NULL
X1 <- expect_warning(simulate_internal(tree,
p = p,
root.state = root.state,
process = "OUBM",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts),
"All the eigen values of the selection strengh do not have a strictly positive real part.")
X1.tips.exp <- extract_simulate_internal(X1, where = "tips", what = "exp")
## Compute expectations with tree matrix
T_tree <- incidence.matrix(tree)
expect_error(shifts.list_to_matrix(tree, shifts), "the dimension p must be specified when shift is NULL")
Delta <- shifts.list_to_matrix(tree, shifts, p = p)
W <- compute_actualization_matrix_ultrametric(tree, selection.strength)
vec_Y <- kronecker(T_tree, diag(1, p, p)) %*% W %*% as.vector(Delta)
X1.tips.exp.mat <- matrix(vec_Y, p, ntaxa) + optimal.value
expect_that(X1.tips.exp, equals(X1.tips.exp.mat))
## Without simulate
expect_warning(X2 <- simulate_internal(tree,
p = p,
root.state = root.state,
process = "OU",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts,
simulate_random = FALSE),
"strictly positive real part")
X2.tips.exp <- extract_simulate_internal(X2, where = "tips", what = "exp")
expect_that(X2.tips.exp, equals(X2.tips.exp))
## Variances
set.seed(1899)
ntaxa <- 5
p <- 2
tree <- rcoal(ntaxa)
variance <- diag(0.5, p)
optimal.value <- c(-3, 5)
selection.strength <- diag(0:1)
exp.stationary <- optimal.value
root.state <- list(random = FALSE,
stationary.root = FALSE,
value.root = exp.stationary,
exp.root = NA,
var.root = NA)
tree_heigth <- max(node.depth.edgelength(tree))
var_1 <- tree_heigth * variance[1, 1]
var_2 <- variance[2, 2] / (2 * selection.strength[2, 2])
X1.tips.states <- sapply(1:50, function(x) extract_simulate_internal(expect_warning(simulate_internal(tree,
p = p,
root.state = root.state,
process = "OUBM",
variance = variance,
optimal.value = optimal.value,
selection.strength = selection.strength,
shifts = shifts),
"strictly positive real part"),
where = "tips", what = "states"))
X1.tips.states.var <- apply(X1.tips.states, 1, var)
expect_equal(mean(X1.tips.states.var[1 + p * 0:(ntaxa - 1)]), var_1, 0.05)
expect_equal(mean(X1.tips.states.var[2 + p * 0:(ntaxa - 1)]), var_2, 0.05)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.