ShiVa Example

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 = 10.83296 
#> 
#> ====== 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 = 19.75432 
#> 
#> ====== 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 = 24.40707 
#> 
#> ====== 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 = 8.136406 
#> 
#> ====== 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 = 3.934407 
#> 
#> ====== 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 = -3.42772 
#> 
#> ====== 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 = -3.42772 
#> 
#> ====== 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 = -3.42772 
#> 
#> ====== 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 = -3.42772 
#> 
#> ====== 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 = -3.42772 
#> 
#> 
#> ====== Backward Correction (Top 10 ) ======
#> Correcting model 1 with shifts_mean = 1,6,7 shifts_var = 3,9,15,31,37,41,44 ...
#> Correcting model 2 with shifts_mean = 6,7 shifts_var = 7,31 ...
#> Correcting model 3 with shifts_mean = 1,6,7 shifts_var = 7,15,31,37,41,44 ...
#> Correcting model 4 with shifts_mean = 6,7 shifts_var =  ...
#> Correcting model 5 with shifts_mean = 6 shifts_var = 7,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 = 7.875518 ======

Visualize Detected Shifts

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

Summarize Shifts

print(summary(result$best_model))
#> $alpha
#> [1] 1.1e-07
#> 
#> $sigma2
#> [1] 0.0067
#> 
#> $loglik
#> [1] 20.2
#> 
#> $BIC
#> [1] 7.88
#> 
#> $mean_shifts
#>   edge node      size
#> 1    6   29 0.2668755
#> 2    7   30 4.1463832
#> 
#> $var_shifts
#>   edge node      size
#> 1    9   31 0.1624866
#> 2   15   35 0.4933164
#> 3   37   46 0.1293234
#> 4   44   49 0.3449323
#> 
#> attr(,"class")
#> [1] "summary.ShiVaShiftModel"