Introduction

This vignette demonstrates how to use the ShiVa R package to detect evolutionary shifts in both optimal trait values (mean) and evolutionary variance under an Ornstein-Uhlenbeck (OU) model. We illustrate the process using a floral trait dataset from Euphorbiaceae species. It is available at phylolm.

library(ShiVa)
library(phylolm)
#> Warning: package 'phylolm' was built under R version 4.4.3
#> Loading required package: ape
#> Warning: package 'ape' was built under R version 4.4.3

Setup

Load the required packages:

data('flowerTree')
data('flowerSize')

Load example data

We load the phylogenetic tree and trait data. The trait is floral diameter, log-transformed.

Y = flowerSize$log_transformed_size
names(Y) = rownames(flowerSize)
tree = flowerTree
# normalize the tree
tree$edge.length = flowerTree$edge.length/max(node.depth.edgelength(flowerTree))

Run ShiVa

set.seed(111)
result = ShiVa(Y,tree, lambda.type = "lambda.min")
#> ====== Model Selection Round 1 ======
#> Trying lambda2 = 0.003697864 ...
#> [1] "233 steps to converge"
#> Selected lambda1 from CV: 1.127226 
#> [1] "204 steps to converge"
#>   shifts_mean = 6 
#>   shifts_var  = 7,15,16,29,31,37,41,44 
#>   log-likelihood = -14.94914 
#> 
#> ====== Model Selection Round 2 ======
#> Trying lambda2 = 0.005516564 ...
#> [1] "241 steps to converge"
#> Selected lambda1 from CV: 0.5590295 
#> [1] "247 steps to converge"
#>   shifts_mean = 1,6,7 
#>   shifts_var  = 7,15,31,37,41,44 
#>   log-likelihood = -10.27391 
#> 
#> ====== Model Selection Round 3 ======
#> Trying lambda2 = 0.008229747 ...
#> [1] "256 steps to converge"
#> Selected lambda1 from CV: 0.4852793 
#>   shifts_mean = 1,6,7 
#>   shifts_var  = 3,9,15,31,37,41,44 
#>   log-likelihood = -7.947515 
#> 
#> ====== Model Selection Round 4 ======
#> Trying lambda2 = 0.01227734 ...
#> [1] "297 steps to converge"
#> Selected lambda1 from CV: 1.129127 
#> [1] "282 steps to converge"
#>   shifts_mean = 6 
#>   shifts_var  = 7,15,31,37,41,44 
#>   log-likelihood = -16.08287 
#> 
#> ====== Model Selection Round 5 ======
#> Trying lambda2 = 0.01831564 ...
#> Selected lambda1 from CV: 0.7391688 
#> [1] "82 steps to converge"
#>   shifts_mean = 6,7 
#>   shifts_var  = 7,31 
#>   log-likelihood = -18.18387 
#> 
#> ====== Model Selection Round 6 ======
#> Trying lambda2 = 0.02732372 ...
#> [1] "128 steps to converge"
#> Selected lambda1 from CV: 0.7806718 
#> [1] "19 steps to converge"
#>   shifts_mean = 6,7 
#>   shifts_var  = none 
#>   log-likelihood = -21.86493 
#> 
#> ====== Model Selection Round 7 ======
#> Trying lambda2 = 0.0407622 ...
#> [1] "9 steps to converge"
#> Selected lambda1 from CV: 0.8721472 
#> [1] "22 steps to converge"
#>   shifts_mean = 6,7 
#>   shifts_var  = none 
#>   log-likelihood = -21.86493 
#> 
#> ====== Model Selection Round 8 ======
#> Trying lambda2 = 0.06081006 ...
#> [1] "12 steps to converge"
#> Selected lambda1 from CV: 0.8572609 
#> [1] "14 steps to converge"
#>   shifts_mean = 6,7 
#>   shifts_var  = none 
#>   log-likelihood = -21.86493 
#> 
#> ====== Model Selection Round 9 ======
#> Trying lambda2 = 0.09071795 ...
#> [1] "12 steps to converge"
#> Selected lambda1 from CV: 0.8570856 
#> [1] "14 steps to converge"
#>   shifts_mean = 6,7 
#>   shifts_var  = none 
#>   log-likelihood = -21.86493 
#> 
#> ====== Model Selection Round 10 ======
#> Trying lambda2 = 0.1353353 ...
#> [1] "12 steps to converge"
#> Selected lambda1 from CV: 0.8570856 
#> [1] "14 steps to converge"
#>   shifts_mean = 6,7 
#>   shifts_var  = none 
#>   log-likelihood = -21.86493 
#> 
#> 
#> ====== Backward Correction (Top 10 ) ======
#> Correcting model 1 with shifts_mean = 6,7 shifts_var =  ...
#> Correcting model 2 with shifts_mean = 6,7 shifts_var = 7,31 ...
#> Correcting model 3 with shifts_mean = 6 shifts_var = 7,15,31,37,41,44 ...
#> Correcting model 4 with shifts_mean = 1,6,7 shifts_var = 7,15,31,37,41,44 ...
#> Correcting model 5 with shifts_mean = 1,6,7 shifts_var = 3,9,15,31,37,41,44 ...
#> Correcting model 6 with shifts_mean = 6 shifts_var = 7,15,16,29,31,37,41,44 ...
#> ====== Selection Finished. Best BIC = 59.25992 ======

Visualize Detected Shifts

plot(result$best_model,title = "ShiVa")

Summarize Shifts