In abichat/zazou: Z-Scores As Ornstein-Uhlenbeck

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

library(ggplot2)
library(tictoc)
library(zazou)
theme_set(theme_minimal())
set.seed(42)

data(alcohol)

Abundance

dim(alcohol$X)

z-scores

pval_obs <- test_wilcoxon(alcohol$X, alcohol$Y)$p.value
zsco_obs <- p2z(pval_obs)

Tree

tree <- force_ultrametric(alcohol$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 = 2, tree = tree, alphaOU = 0.1,
                              method = "L-BFGS-B")
toc()
est_lbfgsb
plot(est_lbfgsb)

Shooting (beta)

tic()
est_shooting <- estimate_shifts(zscores = zsco_obs,
                                lambda = 2, tree = tree, alphaOU = 0.1,
                                method = "lasso", constraint_type = "beta")
toc() # 10 times faster
est_shooting # similar results
plot(est_shooting)

Shooting (yhat)

tic()
est_shooting_p <- estimate_shifts(zscores = zsco_obs,
                                  lambda = 2, tree = tree, alphaOU = 0.1,
                                  method = "lasso", constraint_type = "yhat")
toc()
est_shooting_p # similar results
plot(est_shooting_p)

Shooting (none)

tic()
est_shooting_n <- estimate_shifts(zscores = zsco_obs,
                                  lambda = 2, tree = tree, alphaOU = 0.1,
                                  method = "lasso", constraint_type = "none")
toc() # 10 times faster
est_shooting_n # similar results
plot(est_shooting_n)

Comparisons

est_lbfgsb$objective_value
est_shooting$objective_value
est_shooting_p$objective_value
est_shooting_n$objective_value
est_lbfgsb$optim_info
est_shooting$optim_info
est_shooting_p$optim_info
est_shooting_n$optim_info

qplot(est_lbfgsb$zscores_obs, est_lbfgsb$zscores_est) +
  geom_abline() + geom_vline(xintercept = 0) + coord_equal()
qplot(est_shooting$zscores_obs, est_shooting$zscores_est) +
  geom_abline() + geom_vline(xintercept = 0) + coord_equal()
# qplot(est_shooting_p$zscores_obs, est_shooting_p$zscores_est) +
#   geom_abline() + geom_vline(xintercept = 0) + coord_equal()
qplot(-est_shooting$zscores_est, -est_lbfgsb$zscores_est) +
  geom_abline() + scale_x_log10() + scale_y_log10() + coord_equal()
qplot(-est_shooting$shifts_est, -est_lbfgsb$shifts_est) +
  geom_abline() + scale_x_log10() + scale_y_log10() + coord_equal()
# qplot(-est_shooting$shifts_est, -est_shooting_p$shifts_est) +
  # geom_abline() + scale_x_log10() + scale_y_log10() + coord_equal()