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

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

library(zazou)
library(ape)
set.seed(1234)

Tree

m <- 17
M <- 2 * (m - 1)
tree <- stree(m, type = "left")
edlength <- rep(1, M)
edlength[1:M %% 2 == 1] <- (m - 1):1
tree$edge.length <- edlength
mat_incid <- incidence_matrix(tree)
plot(tree)

All shifts are negative

true_shifts <- rep(0, M)
true_shifts[9 * M / 16] <- -3
true_zscores <- mat_incid %*% true_shifts
true_zscores <- true_zscores[, 1]
plot_shifts(tree, true_shifts, true_scores = true_zscores)

Constraint on beta

estimation1 <- estimate_shifts(zscores = true_zscores, 
                               tree = tree, alphaOU = 1, lambda = 0.1, 
                               method = "lasso")
estimation1
plot(estimation1)

All z-scores are negative but shifts could be positive

true_shifts[13 * M / 16] <- 1
true_zscores <- mat_incid %*% true_shifts
true_zscores <- true_zscores[, 1]
plot_shifts(tree, true_shifts, true_scores = true_zscores)

Constraint on beta

estimation2 <- estimate_shifts(zscores = true_zscores, 
                               tree = tree, alphaOU = 1, lambda = 0.1,
                               method = "lasso")
estimation2
plot(estimation2)

Z-score are well recovered but the solution is absolutely not sparse. We should prioritize ancestral branches before

Constraint on yhat

estimation3 <- estimate_shifts(zscores = true_zscores, 
                               tree = tree, alphaOU = 1, lambda = 0.1,  
                               method = "lasso",
                               constraint_type = "yhat")
estimation3
plot(estimation3)

All z-scores could be positive

true_shifts[2 * M / 16 + 1] <- 1
true_zscores <- mat_incid %*% true_shifts
true_zscores <- true_zscores[, 1]
plot_shifts(tree, true_shifts, true_scores = true_zscores)

Constraint on yhat

estimation4 <- estimate_shifts(zscores = true_zscores, 
                               tree = tree, alphaOU = 1, lambda = 0.1,  
                               method = "lasso",
                               constraint_type = "yhat")
estimation4
plot(estimation4)

No constraint

estimation5 <- estimate_shifts(zscores = true_zscores, 
                               tree = tree, alphaOU = 1, lambda = 0.1, 
                               method = "lasso",
                               constraint_type = "none")
estimation5
plot(estimation5)