tests/testthat/test-utilities-simulate.R

In PhylogeneticEM: Automatic Shift Detection using a Phylogenetic EM

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)
})