khabbazian / l1ou Goto Github PK

Detection of evolutionary shifts in Ornstein-Uhlenbeck models

License: GNU General Public License v2.0

R 91.48% C++ 5.24% C 3.28%

phylogeny ornstein-uhlenbeck study-trait-evolution evolutionary-shift

l1ou's Introduction

Tools for detecting past changes in the expected mean trait values and studying trait evolution from comparative data

The l1ou package provides functions to study trait evolution from comparative data and detect past changes in the expected mean trait values, as well as convergent evolution. It uses the Ornstein-Uhlenbeck process along a phylogenetic tree, which can model a changing adaptive landscape over time and over lineages.

Shifts can be detected from multiple traits, assuming that all traits shifted along the same lineages. Estimation is very fast thanks to lasso techniques, and the user can choose from multiple information criteria for model selection, including a phylognetic BIC. Citation:

M. Khabbazian, R. Kriebel, K. Rohe, and Cécile Ané. Fast and accurate detection of evolutionary shifts in Ornstein-Uhlenbeck models. Methods in Ecology and Evolution, 7(7):811–824. doi:10.1111/2041-210X.12534

l1ou Reference manual

Install using the devtools package

From within R:

install.packages("devtools")
library(devtools)
install_github("glmgen/genlasso")
install_github("khabbazian/l1ou")

Windows users will first need to install Rtools.

Install without the devtools package

To resolve dependencies, first install the following packages from CRAN, then the knitr package. From within R:

install.packages(c("igraph", "phylolm", "lars", "grplasso", "magic", "Rcpp"))
install.packages("knitr")

Download genlasso version 1.3 R package from CRAN archive (link). From within R:

 install.packages("genlasso_1.3.tar.gz", repos=NULL, type="source")

Now in the shell, with asterisks to be replaced with the correct version number:

git clone https://github.com/khabbazian/l1ou.git 
R CMD build l1ou 
R -e 'install.packages("l1ou_*.**.tar.gz")'

Version notes

major changes are indicated below.

v1.43 (2022-08-05): compatibility updates
v1.42 (2019-02-10): bug fix
v1.41 (2017-06-18): small bug fixes
v1.40 (2017-03-27):
- intercept correctly handled after noise-whitening (results may change for variables with a mean far from 0)
- bug fix in the function calculating the square-root (and inverse) of the phylogenetic covariance matrix.
v1.25: the penalty term in the AICc score now considers the intercept as a free variable. The change only affects the final value of the AICc score.
v1.23: "estimate_convergent_regimes" function accepts multiple traits.
v1.22:
- the scores returned by "estimate_shift_configuration” function are now for the non-normalized, original data.
- "estimate_shift_configuration" function also accepts multiple traits with missing values.

l1ou's People

Contributors

Stargazers

Watchers

Forkers

cecileane yuqing19118 kfuku52 lzhangss c-g-l keffy luisdva

l1ou's Issues

Error with argument alpha.upper=0

Hi,

I was trying to fit a straight BM model with l1ou package implementing alpha=0 in the estimate_shift_configuration function. To do that, I used the argument alpha.upper=0, but I got the following error message:

eModel <- estimate_shift_configuration(tree = lizard$tree, 
                                        lizard$Y, criterion="AICc",  
                                        alpha.upper = 0, 
                                        nCores = 3)

Error in estimate_shift_configuration(tree = lizard$tree, lizard$Y, criterion = "AICc",  : 
  alpha.upper must be strictly positive.

Follow to read this error message, I fit the alpha.upper=0.001 but I got an similar error:

eModel <- estimate_shift_configuration(tree = lizard$tree, 
                                        lizard$Y, criterion="AICc",  
                                        alpha.upper = 0.001,
                                        nCores = 3)

Starting first LASSO (alpha=0) to find a list of candidate configurations.
Error in fit_OU_model(tree, Y, result$shift.configuration, opt) : 
  model score is NA in fit_OU_model function! 
		 This should not happen. Please set quietly to false to see the reason.

Thanks,

Jaiber

problem when runing estimate_convergent_regimes with pBIC option

Hi,
I am running the following piece of code

lasso1TreeB1_pBIC<-estimate_shift_configuration(AdjTree, Y[,1], criterion="pBIC")
lasso1TreeB1_pBIC$nShifts

[1] 4

lasso2TreeB1_pBIC <- estimate_convergent_regimes(lasso1TreeB1_pBIC, criterion="pBIC")

Error: $ operator is invalid for atomic vectors

However, when I use AICc or BIC instead of pBIC, it seems to work and $nShifts gives:

[1] 4

I have checked the manual for the function but not sure what is not working. Is there any problem with pBIC when using estimate_convergent_regimes? Something I am overlooking? Thank you!

Error of phylolm_CR in convergence estimation

Hi, I got an error in estimate_convergent_regimes(). This happens only when alpha parameters are optimized like estimate_convergent_regimes(..., fixed.alpha=TRUE). I would appreciate any comments or suggestions.

Error message:

Error in phylolm_CR(Y ~ preds - 1, phy = tr, model = "OUfixedRoot", sc = shift.configuration, : The starting value is not within the bounds of the parameter.
Traceback:

1. estimate_convergent_regimes(fit_ind, criterion = "AICc", method = "backward", 
 .     fixed.alpha = FALSE)
2. estimate_convergent_regimes_surface(model, opt = opt)
3. cmp_model_score_CR(tr, Y, regimes, model$alpha, opt = opt)
4. cmp_AICc_CR(tree, Y, conv.regimes = regimes, alpha = alpha, opt = opt)
5. phylolm_interface_CR(tree, matrix(Y[, i]), conv.regimes, alpha = alpha[[i]], 
 .     opt = opt)
6. phylolm_CR(Y ~ preds - 1, phy = tr, model = "OUfixedRoot", sc = shift.configuration, 
 .     cr = conv.regimes, starting.value = alpha, upper.bound = opt$alpha.upper.bound, 
 .     lower.bound = alpha/100)
7. stop("The starting value is not within the bounds of the parameter.")

Code to reproduce the error:

tree_text="(((((((Mus_musculus_ENSMUSG00000089661:12.89838925,Mus_musculus_ENSMUSG00000095538:12.89838925)n2:7.989028152,Rattus_norvegicus_ENSRNOG00000001499:20.8874174)n4:61.25338149,Oryctolagus_cuniculus_ENSOCUG00000022977:82.14079889)n6:7.68238853,((Homo_sapiens_ENSG00000268975:13.90337806,Homo_sapiens_ENSG00000261857:13.90337806)n9:15.53816876,Macaca_mulatta_ENSMMUG00000013632:29.44154682)n11:60.3816406)n12:6.63920175,(Sus_scrofa_ENSSSCG00000002997:61.96598852,Bos_taurus_ENSBTAG00000009078:61.96598852)n15:34.49640065)n16:255.2960941,Xenopus_tropicalis_ENSXETG00000002588:351.7584832)n18:723.0614697,(((((((Rattus_norvegicus_ENSRNOG00000028721:20.8874174,Mus_musculus_ENSMUSG00000027416:20.8874174)n21:61.25338149,Oryctolagus_cuniculus_ENSOCUG00000001424:82.14079889)n23:7.68238853,(Macaca_mulatta_ENSMMUG00000014494:29.44154682,Homo_sapiens_ENSG00000125879:29.44154682)n26:60.3816406)n27:6.63920175,((Sus_scrofa_ENSSSCG00000007082:61.96598852,Bos_taurus_ENSBTAG00000000516:61.96598852)n30:15.78896789,Canis_lupus_ENSCAFG00000030197:77.75495641)n32:18.70743276)n33:62.13519841,Monodelphis_domestica_ENSMODG00000005398:158.5975876)n35:153.3063338,Gallus_gallus_ENSGALG00000008732:311.9039214)n37:39.85456185,Xenopus_tropicalis_ENSXETG00000003548:351.7584832)n39:723.0614697)n40;"
tree = read.tree(text=tree_text)

trait_text="brain,heart,kidney,liver,ovary,testis
Mus_musculus_ENSMUSG00000089661,2.1769406578477,0.202036051693925,0.63488272196565,0.638510496834611,1.31250675784555,1.45913097497683
Mus_musculus_ENSMUSG00000095538,0.103377188952908,0.0392406207305018,0.341134950053078,0.267010753389883,1.20089284030784,0.0804824077175943
Rattus_norvegicus_ENSRNOG00000001499,1.21464073013422,0.367395910873769,0.52235098372207,0.286406906103611,2.01061005264727,2.54677844995076
Oryctolagus_cuniculus_ENSOCUG00000022977,1.50193008898029,0.419327680366938,0.37886749408802,0.495140302749664,6.93602946184919,2.18609877961886
Homo_sapiens_ENSG00000268975,0,0.266988677902243,0.0531581474038556,0.194942717189735,0,0.0411595822452399
Homo_sapiens_ENSG00000261857,0.738193035867579,1.12596946521021,0.571242331193359,0.647248115527252,1.91034890419106,1.41385417335533
Macaca_mulatta_ENSMMUG00000013632,1.58742742928768,1.18179700137321,1.51345681691155,0.969525697852313,1.82109626650635,2.23218724717014
Sus_scrofa_ENSSSCG00000002997,1.14184323490317,1.06198735908003,0.707240258247191,0.736592284442087,0.65577709045319,0.814151034319946
Bos_taurus_ENSBTAG00000009078,4.09085164303143,2.40058424595315,2.9874525565945,2.85228369236132,3.12337006791316,2.90337478427231
Xenopus_tropicalis_ENSXETG00000002588,0.862939940626679,0.360229223073727,0.327912837312163,0.912571974051781,0.16464914584949,0.570870724658469
Rattus_norvegicus_ENSRNOG00000028721,0.793865212543779,0.126495504481227,0.189564135427474,0.0205186446505775,0.404758033920202,0.622541958944324
Mus_musculus_ENSMUSG00000027416,0.0539642463275866,0.0133687235135755,0.0119717717558467,0.026844934044235,1.89059101437825,0
Oryctolagus_cuniculus_ENSOCUG00000001424,0,0.443929837825391,1.09544972717812,0,0,0.550209062958829
Macaca_mulatta_ENSMMUG00000014494,0.0960701959183422,0.0365536296958946,0.0460032798433054,0,0.0805454117819435,3.56218496133299
Homo_sapiens_ENSG00000125879,0,0,0.00201377894184635,0,0.0816434192670939,0.0899198236644605
Sus_scrofa_ENSSSCG00000007082,1.02235798411815,0.745458667732567,0.754893805042313,0.632903316261303,0.0960255936389955,0.452407969797526
Bos_taurus_ENSBTAG00000000516,3.4163923290383,0.723236705462928,0,0.226080404205593,2.80622903364123,3.21201931204057
Canis_lupus_ENSCAFG00000030197,0.586061750953316,0.202319615318614,0.112868933178583,0.031461327458019,0.00608259828675725,0
Monodelphis_domestica_ENSMODG00000005398,2.00670644992751,0.345354276333881,0.228500546652594,0.310639667267954,2.17171017530506,2.98230624273514
Gallus_gallus_ENSGALG00000008732,1.99335564029219,0.448788658290314,0.140094253253896,0.123957568649032,3.10671747036773,0.383263265231511
Xenopus_tropicalis_ENSXETG00000003548,3.55950533293704,0.774471483114139,0.561099010697649,0.0882024773033921,0.200575784207915,0.963046405029511"
trait_matrix = read.table(text=trait_text, sep=",")

adj_data = adjust_data(tree=tree, Y=trait_matrix, normalize = FALSE)
fit_ind = estimate_shift_configuration(tree=adj_data$tree, Y=adj_data$Y, criterion="AICc", root.model="OUrandomRoot", nCores=1, rescale=TRUE)

# no error, fixed.alpha=TRUE
fit_conv1 = estimate_convergent_regimes(fit_ind, criterion="AICc", method="backward", fixed.alpha=TRUE)

# error, fixed.alpha=FALSE
fit_conv2 = estimate_convergent_regimes(fit_ind, criterion="AICc", method="backward", fixed.alpha=FALSE)

Many thanks.
Kenji

issue with candidate edges

I've been running a few experiments with l1ou and it seems like analysis of multiple traits takes much longer when candidate edges are specified (ie, a total set of edges that is less than all edges in the tree). Is this expected?

I'm currently running a particularly extreme test with 8 traits, 200 tips -- the model is using about 70gb of memory on the first LASSO (and > 20 hours of computing so far). Running the same model without candidate edges completes in less time.

Thanks,
Jake Berv

Error in running estimate_shift_configuration - the condition has length > 1

I have a error when running estimate_shift_configuration ():
Error in if (class(Y) != "matrix") { : the condition has length > 1
Even with the function example eModel <- estimate_shift_configuration(lizard$tree, lizard$Y) . Possibly because I'm in new R version 4.2.0?

Empty profile when using multiple cores

Hello!

When running the estimate_shift_configuration function with multiple cores the profiles element in the output, where the different configurations and their scores should be, comes back empty. This doesn't happen when using a single core. Below is an example of how I'm getting this issue with example data:

library(l1ou)

data(lizard.tree, lizard.traits)
lizard <- adjust_data(lizard.tree, lizard.traits[,1])

# the new tree is normalized: each tip at distance 1 from the root.
# new Y: matrix of size 100 x 1 

eModel <- estimate_shift_configuration(lizard$tree, lizard$Y)

# Starting first LASSO (alpha=0) to find a list of candidate configurations.
# Starting second LASSO (alpha= 0.61 ) for another list of candidates.

str(eModel$profile)

# List of 2
#  $ scores        : num [1:2399] 16.8 18.3 18.7 19.3 19.9 ...
#  $ configurations:List of 2399
#   ..$ : num [1:8] 14 32 55 74 77 98 118 164
#   ..$ : num [1:7] 14 32 55 77 98 118 164
#   ..$ : num [1:6] 14 32 77 98 118 164
#   ..$ : num [1:9] 14 32 37 55 74 77 98 118 164
 
# [etc...]

Now multithreaded, produces an empty profile

eModel_mc <- estimate_shift_configuration(lizard$tree, lizard$Y, nCores = 4)

str(eModel_mc$profile)

# List of 2
#  $ scores        : num(0)
#  $ configurations: list()

I'm running R v. 4.22 and l1ou v. 1.43.

Thanks!

matrix class warning

Hi,

I assume this is relatively minor, but I am receiving the following warning message when I run functions that detect the class of input data:

Warning message:
In if (class(Y) != "matrix") { :
the condition has length > 1 and only the first element will be used

I believe this is related to how R v4 now has the class of matrices as a vector including "matrix" and "array".

CPU overuse

I just realized that estimate_shift_configuration() tries to use almost all CPU resources no matter what nCores is. I guess this is not expected behavior. For example, when this code (nCores=1) is executed...

library(l1ou)

tree = rtree(100)
tree = ape::chronoMPL(tree)

trait = matrix(runif(1000, 0, 3), 100, 10)
rownames(trait) = tree$tip.label
colnames(trait) = paste0('trait', 1:ncol(trait))

adj_data = l1ou::adjust_data(tree=tree, Y=trait, normalize=FALSE)

fit_ind = l1ou::estimate_shift_configuration(tree=adj_data$tree, Y=adj_data$Y, nCores=1, rescale=FALSE)

top says the R session is using almost all CPU resources (6236/6400%).

  PID USER      PR  NI    VIRT    RES    SHR S  %CPU %MEM     TIME+ COMMAND                                                                                              
 3710 kfuku     20   0 2836076 481560  10460 R  6236  0.1  71:33.71 R                                                                                                    
 4899 root      20   0 5221332   3.0g  66304 S   3.6  0.6 213:35.75 cmd                                                                                                  
10244 root      20   0   51560   7836   2520 S   1.3  0.0   0:00.04 sample_ilo                                                                                           
10227 kfuku     20   0  173380   3532   1704 R   1.0  0.0   0:00.06 top

Here is the result from nCores=8. It initiated with the single process with ~6400/6400 %CPU...

  PID USER      PR  NI    VIRT    RES    SHR S  %CPU %MEM     TIME+ COMMAND                                                                                              
 3710 kfuku     20   0 2838736 490788  10460 R  6332  0.1 115:28.19 R                                                                                                    
17267 root      20   0   89352   7308   1556 S   1.7  0.0 196:08.84 zabbix_agentd                                                                                        
11635 kfuku     20   0  173380   3524   1700 R   0.7  0.0   0:00.12 top

Each of 8 processes started using more than one core without limit.

  PID USER      PR  NI    VIRT    RES    SHR S  %CPU %MEM     TIME+ COMMAND                                                                                              
11666 kfuku     20   0 2838772 481720   2392 R  1100  0.1   1:44.52 R                                                                                                    
11664 kfuku     20   0 2838772 481784   2456 R 843.0  0.1   0:53.74 R                                                                                                    
11662 kfuku     20   0 2838772 481704   2376 R 802.0  0.1   0:49.33 R                                                                                                    
11663 kfuku     20   0 2838772 481720   2392 R 782.7  0.1   0:52.42 R                                                                                                    
11921 kfuku     20   0 2838772 481708   2380 R 750.5  0.1   0:45.36 R                                                                                                    
11730 kfuku     20   0 2838772 481724   2396 R 740.1  0.1   0:48.67 R                                                                                                    
11665 kfuku     20   0 2838772 481716   2388 R 700.0  0.1   0:47.18 R                                                                                                    
11736 kfuku     20   0 2838772 481720   2392 R 645.9  0.1   0:41.93 R                                                                                                    
12180 kfuku     20   0  173380   3528   1700 R   1.0  0.0   0:00.06 top

My session info:

> sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /lustre7/home/lustre4/kfuku/.pyenv/versions/miniconda3-4.3.30/lib/R/lib/libRblas.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] l1ou_1.42      genlasso_1.3   Matrix_1.2-17  MASS_7.3-51.4  magic_1.5-9   
 [6] abind_1.4-5    grplasso_0.4-6 lars_1.2       phylolm_2.6    ape_5.3       
[11] igraph_1.2.4.1

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1         lattice_0.20-38    codetools_0.2-16   listenv_0.7.0     
 [5] future_1.12.0      digest_0.6.18      grid_3.5.1         nlme_3.1-139      
 [9] magrittr_1.5       future.apply_1.2.0 tools_3.5.1        compiler_3.5.1    
[13] pkgconfig_2.0.2    globals_0.12.4

I would appreciate it if you could fix this problem.

Error: object 'elist.min' not found when estimating convergent regimes using "pBIC" and "pBICess"

Hi,

I am running the following piece of code

lasso1TreeB1_pBIC<-estimate_shift_configuration(AdjTree, Y[,1], criterion="pBIC")
lasso2TreeB1_pBIC <- estimate_convergent_regimes(lasso1TreeB1_pBIC, criterion="pBIC")

I get the following error

Error: object 'elist.min' not found

This seems to be also the case when I choose the criterion "pBICess" for estimating shifts. Is there any problem with my code?

Also, this is probably a silly question, but is the criterion for estimating shifts and convergence supposed to be the same? Or can it be, say pBIC for estimating shifts and BIC/AICc/ for estimating convergence (I illustrate it with an example bellow):

lasso1TreeB1_pBIC<-estimate_shift_configuration(AdjTree, Y[,1], criterion="pBIC")
lasso2TreeB1_pBIC <- estimate_convergent_regimes(lasso1TreeB1_pBIC, criterion="AICc")

Does this make any sense? Thank you!

Gonzalo

plot.l1ou: improve documentation about the trait rescaling

The scale on the plot shows the values of the standardized variables, but this is undocumented and is therefore quite confusing for users (as I've heard).

That is, the scale shows the values of: (trait - mean)/standard deviation. This is done in the code here. The values printed on the left and right side of the axis are the min and max of the standardized trait, which are: (min and max of the original values - mean)/sd. The place where the bars start corresponds to the mean of the original trait values (and corresponds to 0 for the standardized values).

This is all useful for visually highlighting the species with low / average / high values, but should be documented.

github installation issues

user request received by email: an error was received when attempting to install the package with install_github(). The error message said "cannot load shared object, invalid ELF header". Multiple people got this error: on an ubuntu machine and on a Mac machine.

Error message in adjust_tree()

Should line 49 of the function read?
Y <- Y[order(order(tree$tip.label)), ]

estimate_shift_configuration error

I'm trying to use the program to estimate optimal regimens over a large phylogeny (tips=1804). However, regardless of the maximum number of shift allowed, the program gives me the following error:

Starting first LASSO (alpha=0) to find a list of candidate configurations.
Error in if (min.score > res$score) { : argument has zero length

The tree is ultrametric and binary.

measurement error

Is there a way to include standard error for input traits in an l1ou analysis?

estimate_convergent_regimes error

the function estimate_convergent_regimes appears to be broken - running it returns the following error:
Error in sqrt_OU_covariance(tr, alpha = alpha, root.model = "OUfixedRoot", :
unused argument (normalize.tree.height = TRUE)

It can be traced to the function find_convergent_regimes, which calls sqrt_OU_covariance but with the argument normalize.tree.height=TRUE, which is missing from that function's definition.

Error in select_best_solution(tree, Y, sol, opt = opt)

Hi, I got a following error when I tried estimate_shift_configuration() on my dataset.

Warning message in matrix(ifelse(abs(x[!is.na(grpIdx)]) > 0, 1, 0), ncol = nVariables):
“data length [25] is not a sub-multiple or multiple of the number of columns [6]”Warning message in matrix(ifelse(abs(x[!is.na(grpIdx)]) > 0, 1, 0), ncol = nVariables):
...
Error in select_best_solution(tree, Y, sol, opt = opt): object 'best.shift.configuration' not found
Traceback:

1. estimate_shift_configuration(tree = adj_data$tree, Y = adj_data$Y, 
 .     criterion = "AICc", root.model = "OUrandomRoot", nCores = 1, 
 .     rescale = TRUE)
2. estimate_shift_configuration_known_alpha_multivariate(tree, Y, 
 .     est.alpha = TRUE, opt = l1ou.options)
3. select_best_solution(tree, Y, sol, opt = opt)

My session info is:

R version 3.3.2 (2016-10-31)
Platform: x86_64-apple-darwin11.4.2 (64-bit)
Running under: macOS Sierra 10.12.3

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] l1ou_1.40      genlasso_1.3   Matrix_1.2-7.1 MASS_7.3-45    magic_1.5-6   
 [6] abind_1.4-5    grplasso_0.4-5 lars_1.2       phylolm_2.5    ape_4.0       
[11] igraph_1.0.1  

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.9     magrittr_1.5    uuid_0.1-2      lattice_0.20-34
 [5] R6_2.2.0        stringr_1.1.0   tools_3.3.2     grid_3.3.2     
 [9] nlme_3.1-128    digest_0.6.10   crayon_1.3.2    IRdisplay_0.4.4
[13] repr_0.10       IRkernel_0.7.1  evaluate_0.10   pbdZMQ_0.2-4   
[17] stringi_1.1.2   jsonlite_1.1

You can reproduce the error using this code.

tree_text = '((Canis_lupus_ENSCAFG00000007369:77.75495641,(Bos_taurus_ENSBTAG00000031373:61.96598852,Sus_scrofa_ENSSSCG00000028270:61.96598852)n2:15.78896789)n4:18.70743276,Oryctolagus_cuniculus_ENSOCUG00000010555:96.46238917)n6;'
tree = read.tree(text=tree_text)

trait_text="brain,heart,kidney,liver,ovary,testis
Sus_scrofa_ENSSSCG00000028270,0,0.00613622199739828,0,0,0,0.086648021213729
Bos_taurus_ENSBTAG00000031373,0.0037485584053104,0.0025221135879922,0,0,0,0.0987999096017663
Canis_lupus_ENSCAFG00000007369,0,0.0611830531983022,0.0671139398094019,0.0692757761941058,0,0
Oryctolagus_cuniculus_ENSOCUG00000010555,0,4.48056104371082e-17,0,0,0,1.5554235315619"
trait_matrix = read.table(text=trait_text, sep=",")

adj_data = adjust_data(tree=tree, Y=trait_matrix, normalize = FALSE)
fit_ind = estimate_shift_configuration(tree=adj_data$tree, Y=adj_data$Y, criterion="AICc", root.model="OUrandomRoot", nCores=1, rescale=TRUE)

I would appreciate it if you would give me some advice.

Kenji

sqrt_OU_covariance: tree reordering

on line 38 of sqrt_OU_covariance (in my fork), the command re-orders the edges in the tree:

tree <- reorder(tree, "prun")

It would be safer to use a postorder instead of a 'pruningwise' order (in most cases it would not matter, but in some cases it could)
When the user needs to use the function over and over on the same tree, but different alpha values, it would be faster if this re-ordering could be avoided: add an option check.order=T or F, and turn off the re-ordering of edges if check.order=F. The function would then be a lot faster when re-used on the same tree, if that tree is known to have its edges in post-order traversal already.
Add a comment: the function assumes that reorder(...) does not change the order of the nodes, just the order of edges, so that column i in each matrix still corresponds to internal node i, etc. (and the last column corresponds to the root edge instead of an internal node).

Then on line 41 the branch lengths are transformed, (only if alpha>0 though):

tre <- transf.branch.lengths(tree, model=root.model, parameters=list(alpha=alpha))$tree

Add a comment: this is the step that requires an ultrametric tree. If the tree is not ultrametric, the function returns a wrong result with no warning (see example below).
For now, add to the documentation that the tree needs not be ultrametric for the BM model, i.e. with default value alpha=0.
There should be a way to allow for non-ultrametric trees, for the BM at the very least, and for the OU model in general too: see related trick in phylolm. Otherwise, add code to check that the tree is ultrametric if alpha>0, but also an option to turn off this check (for same purpose as above).

Example of wrong result:

tre1 <- read.tree(text="((Strix_aluco:4.2,Asio_otus:4.2):3.1,Athene_noctua:7.3);")
tre1$edge.length[3] <- 8 # was 4.2 originally: tre1 no longer ultrametric
res1 <- sqrt_OU_covariance(tre1, alpha=50); # large alpha: almost like iid model
res1$sqrtSigma %*% t(res1$sqrtSigma) # should be close to identity matrix, but far from it!!

Error in solve.default(comp$XX): system is computationally singular

Sorry to keep bothering you. I found another error which happens to some of my relatively large datasets (typically >150 leaves). Like the problem in #29, the error can be circumvented by specifying fixed.alpha=TRUE in estimate_convergent_regimes(), although the error message is different. I would appreciate any comments or suggestion.

Input data:
traits.txt
tree.txt

Error message:

Error in solve.default(comp$XX): system is computationally singular: reciprocal condition number = 4.72733e-22
Traceback:

1. estimate_convergent_regimes(fit_ind, criterion = "AICc", method = "backward", 
 .     fixed.alpha = FALSE, nCores = 1)
2. estimate_convergent_regimes_surface(model, opt = opt)
3. cmp_model_score_CR(tr, Y, regimes, model$alpha, opt = opt)
4. cmp_AICc_CR(tree, Y, conv.regimes = regimes, alpha = alpha, opt = opt)
5. phylolm_interface_CR(tree, matrix(Y[, i]), conv.regimes, alpha = alpha[[i]], 
 .     opt = opt)
6. phylolm_CR(Y ~ preds - 1, phy = tr, model = "OUfixedRoot", sc = shift.configuration, 
 .     cr = conv.regimes, starting.value = alpha, upper.bound = opt$alpha.upper.bound, 
 .     lower.bound = alpha/100)
7. optim(logstart, fn = minus2llh, method = "L-BFGS-B", lower = loglower, 
 .     upper = logupper, ...)
8. (function (par) 
 . fn(par, ...))(-22.855602561775)
9. fn(par, ...)
10. loglik(prm, y, X = preds)
11. solve(comp$XX)
12. solve.default(comp$XX)

Code to reproduce:

library(l1ou)
tree = read.tree("tree.txt")
trait_matrix = read.table("traits.txt", sep="\t")
adj_data = adjust_data(tree=tree, Y=trait_matrix, normalize = FALSE)
fit_ind = estimate_shift_configuration(tree=adj_data$tree, Y=adj_data$Y, criterion="AICc", root.model="OUrandomRoot", nCores=4, rescale=TRUE)

# no error
fit_conv = estimate_convergent_regimes(fit_ind, criterion="AICc", method="backward", fixed.alpha=TRUE, nCores=4)

# error
fit_conv = estimate_convergent_regimes(fit_ind, criterion="AICc", method="backward", fixed.alpha=FALSE, nCores=4)

Many thanks.

Kenji

3.4.2 compilation error

Hi,

I'm having trouble installing this package in R version 3.4.2. I've tried both methods in the readme, but both result in compilation failures.
I'd love to be able to use this software, so let me know what other information you need from me!

Best,
Will

Installation of l1ou

Hi, I have problem in installing l10u. I'm window user, please help me in using this package.
Thank you!

criterion arg in estimate_convergent_regimes

This is really a two-problem issue:
In the past I found that with some versions of the package pBIC worked, and with others AICc worked, but both didn't always work. I found that if I let l1ou decide what was default it would work. With the current version, when I did that I received the error message:

system.time(backward <- estimate_convergent_regimes(model=forward))
Error in match.arg(criterion) : 'arg' must be of length 1
Calls: system.time ... estimate_convergent_regimes -> estimate_convergent_regimes_surface -> match.arg

When I tried re-running it explicitly setting criterion to pBIC, I got this error message

system.time(backward <- estimate_convergent_regimes(model=forward, criterion="pBIC"))
Error in cmp_model_score_CR(tr, Y, model$shift.configuration, regimes, :
undefined criterion for convergent evolution.
Calls: system.time ... estimate_convergent_regimes_surface -> cmp_model_score_CR

It appears to be running when I explicitly set it to AICc.

The second part of this problem is that the estimate_convergent_regimes function suggests that only AICc and pBIC are options, but the documentation refers back to configuration_ic, which has 5 options. It should be made clear in the definition of the function that those other options either are or not options (i.e. the part where the function is defined and the "default" for criterion is set to c("AICc", "pBIC") needs to be updated if those other 3 options are indeed possible.

Lastly, and slightly unrelatedly, there isn't much documentation given for the new argument "method", and it doesn't say what the default is. This should be made clear. (as should it be made clear what the default criterion is for all functions that use it--since it does choose something whether provided with one or not).

Thanks!

manual unavailable

There is a broken link to the manual from the readme file: l1ou Reference manual. Could this manual be posted some other place? (just passing on user request received via email)

in estimate_shift_configuration: error in if (zmin < gamhat)

The error shows this, after calling estimate_shift_configuration:

Starting first LASSO (alpha=0) to find a list of candidate configurations.
Starting second LASSO (alpha= 0.97 ) for another list of candidates.
Error in if (zmin < gamhat) { : missing value where TRUE/FALSE needed

Here is code to reproduce the error. (I tried hard to find a relatively fast & reproducible example).

library(l1ou)
data(lizard.tree)
Y = c(0.1944112326082775843439,0.2615068389462029685433,-0.07637984205668357784447,-0.1727580692981197652003,0.2934480262963777841279,-0.05445719914907685976768,0.02650650042204652154232,0.3928677951246398736274,-0.03066280657201297943359,0.1400230714268974885339,0.3111385068720688984456,0.2181363661125211295122,0.1454875513324073998955,0.1160210239105412266536,0.04071687564916981472152,0.4185011496745297820965,-0.3291055094818656767686,0.2883096640636917951106,0.3121786434458098113964,0.2462950059644864742037,0.2304067511391049150049,-0.04037490293265803137368,-0.04700945743145562571996,0.2378876919219531649308,-0.08614324242048232438407,0.07222972184190346034427,-0.170195694921206452932,0.1650936975015285246293,0.2192309259413985633724,0.2581594313692500697321,0.1321251831067828952371,0.1526127388876840074161,0.1015299368103808885788,0.2359338510881836270539,0.30855097082924637375,0.311676065613906649876,0.2173648621135936820359,0.3282197952720358746781,-0.8509565110655578079601,-0.001985024323349571240271,-0.2233235872028019319835,-0.1795363601903273709226,-0.2962850975932650454681,-0.2308048749146208167282,-0.2068763128374506432561,-0.04057174583271036527599,-0.04749797662037748280373,0.1586246256821226829903,-0.1344704141733906732625,-0.07948665511398328442638,-0.2462935522959820433542,-0.06303738479084433632416,-0.01857153931405388855302,0.1663626253977315561094,-0.7498727809457496062961,-0.8179861788370244024549,-0.7003474280515201710884,-0.7119943617081601061614,-0.4730439692236153570448,-0.5015238315496822751882,-0.5400803197131716082424,-0.6970199862183987793429,-0.2385252570273939110024,-0.3685431580852820410144,-0.5290574059715392740699,-0.5523618635942365573399,-0.5121314496197453269843,-1.374570027092596236074,-1.365040272048220426626,-1.683909521118079899438,-0.5656336055030394271981,-0.6861931343335359034796,-0.7320270959066771387924,0.3507437055324111874199,0.282431420417687040203,0.113453673033593854802,0.05762758896325523988446,0.1054005028098378704549,0.1513609292613641299496,0.1667395796168281429939,0.2258102347045503055512,-0.208204209901056869203,-0.2599131879421122115481,0.2178061517235754529498,0.1508268013584183608877,0.4308035251951539135185,0.5016602589937741996096,0.5106318642667988516592,0.3951991244587529927834,0.3629133822297955780378,0.1809739468049765820368,0.18504628339645168289,0.2994298963219166331839,0.23710754696231783889,-0.562054641655930997679,-0.1166817102120809163113,0.1353760968930994013082,0.06631886308094243898115,0.3852605661204347442528,-0.9656860100773794197693)
names(Y) = lizard.tree$tip.label
suppressWarnings(estimate_shift_configuration(lizard.tree, Y, max.nShifts=10))

Super strange: rounding the values down (sometimes up) eliminates the error. The code below

fit = estimate_shift_configuration(lizard.tree, round(Y,15), max.nShifts=10) # no error
fit
library(phytools) # for the plotting function below: to show tree & data
plotTree.wBars(lizard.tree, Y)

produces no error, with this (edited) output:

number of shifts: 7
pBIC score: 5.02116
edge indices of the shift configuration (column names) and the corresponding shift values:
           77        32      55       118        98        74        14
s.v -1.245648 -1.899369 -2.2836 -3.451839 -1.867629 -1.802358 -1.602436

convergent regimes and edge indices of the shift configuration
regime 1-> 14
regime 2-> 74
regime 3-> 98
regime 4-> 118
regime 5-> 55
regime 6-> 32
regime 7-> 77

adaptation rate (alpha)      0.96764282
variance (sigma2)            0.07336021
stationary variance (gamma)  0.03790666
logLik                      48.62335357

error related with R version

Dear developers of l1ou package

I have use your package l1ou before with previous data, but now I am getting this error: Error in cmp_sqrt_OU_covariance(my.edge.list, length(tre$tip.label), tre$root.edge) : internal error!@sqrt_OU_covariance.cpp.

I am almost certain that the error has to be related with R version, because after I got the error, I tried to run the analysis in a older R version and it work. But now week after, I am getting the error again, so I think that is some kind of bug error between the package and the R version. I would really appreciatte any help that you can provided me, because this analysis is a very important part of my PhD.

thanks in advance

Error in starting.value$alpha : $ operator is invalid for atomic vectors

Hi,

I get an error when estimating the convergent regimes with a tree. The function estimate_shift_configuration seems to work out fine (I get a tree with the regimes on it)

>AdjData<-adjust_data(MyTree, MyDataSet)

>lasso1Tree<-estimate_shift_configuration(AdjData$tree, MyDataSet[,5:6], criterion="pBIC", lars.alg = "lasso")
>lasso1Tree$nShifts
[1] 3
lasso2Tree<- estimate_convergent_regimes(lasso1Tree, criterion="pBIC")
Error in starting.value$alpha : $ operator is invalid for atomic vectors

Not sure what is not working. Seems like there are indeed regimes found with my tree and data set but for some reason the model produced does not seem to work when analyzing for convergence. I tried tweaking the options offered by estimate_convergent_regimes but method needs to be = "backward" since there is more than one trait, and changing fixed.alpha from TRUE to FALSE doesn't change the error message. Any pointers to what could be happening?

Thank you very much!

Fit Brownian motion

Hello, I am writing to ask if it is possible to fit a Brownian model? I have tried setting the value of alpha.upper to 0 or a small value (eg. 0.001), as well as the parameters alpha.lower and alpha.starting. However, for any combination of these parameters, the program gives me the following error:

Error: $ operator is invalid for atomic vectors
Error in fit_OU_model(tree, Y, shift.configuration, opt) :
model score is NA in fit_OU_model function!
This should not happen. Please set quietly to false to see the reason.
Also: Warning message:
In my_phylolm_interface(tr, y, s.c, opt) :
phylolm returned error with a shift configuration
of size 0. You may
want to change alpha.upper/alpha.lower!

Greetings