knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", 
  fig.align = "center",
  eval = TRUE
)
library(zazou)
library(tictoc)
data(chlamydiae)

z-scores

pval_obs <- test_kruskalwallis(chlamydiae$X, chlamydiae$Y)$p.value
zsco_obs <- p2z(pval_obs)

Tree

tree <- force_ultrametric(chlamydiae$tree)
N_branch <- length(tree$edge.length)
plot_shifts(tree, NA, obs_scores = zsco_obs)

Shift estimation

L-BFGS-B

tic()
est_lbfgsb <- estimate_shifts(zscores = zsco_obs,
                              lambda = 0.1, tree = tree, alphaOU = 1,
                              method = "L-BFGS-B")
toc()
est_lbfgsb
plot(est_lbfgsb)

Shooting

set.seed(42)
tic()
est_shooting <- estimate_shifts(zscores = zsco_obs,
                                lambda = 0.1, tree = tree, alphaOU = 1,
                                method = "lasso")
toc()
est_shooting
plot(est_shooting)

Comparisons

est_lbfgsb$objective_value
est_shooting$objective_value
est_lbfgsb$optim_info
est_shooting$optim_info
library(ggplot2)
theme_set(theme_minimal())
qplot(est_lbfgsb$zscores_obs, est_lbfgsb$zscores_est) +
  geom_abline() + geom_vline(xintercept = 0)
qplot(est_shooting$zscores_obs, est_shooting$zscores_est) +
  geom_abline() + geom_vline(xintercept = 0)
qplot(-est_shooting$zscores_est, -est_lbfgsb$zscores_est) +
  geom_abline() + scale_x_log10() + scale_y_log10()


abichat/zazou documentation built on Sept. 8, 2021, 6:53 a.m.