context("Function merge rotation") ############################################################################### test_that("Rotations", { testthat::skip_on_cran() data(monkeys) # Fit data res <- PhyloEM(Y_data = monkeys$dat, phylo = monkeys$phy, process = "scOU", nbr_alpha = 2, random.root = FALSE, K_max = 10, method.selection = "BGHmlraw", parallel_alpha = F) # rotate data rot <- matrix(c(cos(pi/4), -sin(pi/4), sin(pi/4), cos(pi/4)), nrow= 2, ncol = 2) Yrot <- t(rot) %*% monkeys$dat # Fit rot res_rot <- PhyloEM(Y_data = Yrot, phylo = monkeys$phy, process = "scOU", nbr_alpha = 2, random.root = FALSE, K_max = 10, method.selection = "BGHmlraw", parallel_alpha = F) # Find rotation expect_equal(find_rotation(res, res_rot), rot) # rotate params expect_equivalent(res$alpha_max$params_estim$`3`, rotate_params(res_rot$alpha_max$params_estim$`3`, rot)) expect_equal(log_likelihood(params_process(res, K = 3), phylo = monkeys$phy, Y_data = monkeys$dat), log_likelihood(rotate_params(params_process(res_rot, K = 3), rot), phylo = monkeys$phy, Y_data = monkeys$dat)) # another rotation rot2 <- matrix(c(cos(pi/3), -sin(pi/3), sin(pi/3), cos(pi/3)), nrow= 2, ncol = 2) Yrot2 <- t(rot2) %*% monkeys$dat # Fit rot res_rot2 <- PhyloEM(Y_data = Yrot2, phylo = monkeys$phy, process = "scOU", nbr_alpha = 2, random.root = FALSE, K_max = 10, method.selection = "BGHmlraw", parallel_alpha = F) # merge solution res_merge <- merge_rotations(res, res_rot, res_rot2) expect_equal(res_merge$alpha_max$results_summary$log_likelihood, apply(sapply(list(res, res_rot, res_rot2), function(x) x$alpha_max$results_summary[["log_likelihood"]]), 1, max)) expect_equal(log_likelihood(params_process(res_merge, K = 3), phylo = monkeys$phy, Y_data = monkeys$dat), res_merge$alpha_max$results_summary$log_likelihood[4]) ############################################################################### ## Not rotations ## transform data norot <- matrix(c(1, 0.5, 0.1, 2), 2) Ynorot <- t(norot) %*% monkeys$dat # Fit rot expect_warning(res_norot <- PhyloEM(Y_data = Ynorot, phylo = monkeys$phy, process = "scOU", nbr_alpha = 2, random.root = FALSE, K_max = 10, method.selection = "BGHmlraw", parallel_alpha = F), "The solution fund by lasso had non independent vectors. Had to modify this solution.") # Rotation ? expect_error(find_rotation(res, res_norot), "The datasets are not linked by a rotation.") ## Other completely unrelated data Ynorot2 <- matrix(rnorm(length(monkeys$dat)), dim(monkeys$dat)) colnames(Ynorot2) <- colnames(Ynorot) # Fit rot res_norot2 <- PhyloEM(Y_data = Ynorot2, phylo = monkeys$phy, process = "scOU", nbr_alpha = 2, random.root = FALSE, K_max = 10, method.selection = "BGHmlraw", parallel_alpha = F) # Rotation ? expect_error(find_rotation(res, res_norot2), "The datasets are not linearly mapped.") ############################################################################### ## Missing data dat <- monkeys$dat dat[, c(12, 23, 32)] <- NA # Fit data res <- PhyloEM(Y_data = dat, phylo = monkeys$phy, process = "scOU", nbr_alpha = 2, random.root = FALSE, K_max = 10, method.selection = "BGHmlraw", parallel_alpha = F) # rotate data rot <- matrix(c(cos(pi/4), -sin(pi/4), sin(pi/4), cos(pi/4)), nrow= 2, ncol = 2) Yrot <- t(rot) %*% dat # Fit rot res_rot <- PhyloEM(Y_data = Yrot, phylo = monkeys$phy, process = "scOU", nbr_alpha = 2, random.root = FALSE, K_max = 10, method.selection = "BGHmlraw", parallel_alpha = F) # Find rotation expect_equal(find_rotation(res, res_rot), rot) # rotate params expect_equivalent(res$alpha_max$params_estim$`3`, rotate_params(res_rot$alpha_max$params_estim$`3`, rot), 1e-3) expect_equal(log_likelihood(params_process(res, K = 3), phylo = monkeys$phy, Y_data = monkeys$dat), log_likelihood(rotate_params(params_process(res_rot, K = 3), rot), phylo = monkeys$phy, Y_data = monkeys$dat), tolerance = 1e-4) # another rotation rot2 <- matrix(c(cos(pi/3), -sin(pi/3), sin(pi/3), cos(pi/3)), nrow= 2, ncol = 2) Yrot2 <- t(rot2) %*% dat # Fit rot res_rot2 <- PhyloEM(Y_data = Yrot2, phylo = monkeys$phy, process = "scOU", nbr_alpha = 2, random.root = FALSE, K_max = 10, method.selection = "BGHmlraw", parallel_alpha = F) # merge solution res_merge <- merge_rotations(res, res_rot, res_rot2) expect_equal(res_merge$alpha_max$results_summary$log_likelihood, apply(sapply(list(res, res_rot, res_rot2), function(x) x$alpha_max$results_summary[["log_likelihood"]]), 1, max)) expect_equal(log_likelihood(params_process(res_merge, K = 3), phylo = monkeys$phy, Y_data = dat), res_merge$alpha_max$results_summary$log_likelihood[4]) ############################################################################### ## Missing data - Errors dat <- monkeys$dat dat[, c(12, 23, 32)] <- NA dat[1, 13] <- NA # Fit data res <- PhyloEM(Y_data = dat, phylo = monkeys$phy, process = "scOU", nbr_alpha = 2, random.root = FALSE, K_max = 10, method.selection = "BGHmlraw", parallel_alpha = F) # rotate data rot <- matrix(c(cos(pi/4), -sin(pi/4), sin(pi/4), cos(pi/4)), nrow= 2, ncol = 2) Yrot <- t(rot) %*% dat # Fit rot res_rot <- PhyloEM(Y_data = Yrot, phylo = monkeys$phy, process = "scOU", nbr_alpha = 2, random.root = FALSE, K_max = 10, method.selection = "BGHmlraw", parallel_alpha = F) # Find rotation expect_error(find_rotation(res, res_rot), "Rotations can only be applied to datasets that have entire species missing") expect_error(find_rotation(res_rot, res), "Rotations can only be applied to datasets that have entire species missing") # all missing plus another dat[2, 13] <- NA dat[, 14] <- NA res <- PhyloEM(Y_data = dat, phylo = monkeys$phy, process = "scOU", nbr_alpha = 2, random.root = FALSE, K_max = 10, method.selection = "BGHmlraw", parallel_alpha = F) expect_error(find_rotation(res_rot, res), "The two datasets used in the analyses do not have the same missing data.") })