Starting values and bounds for `sigma2_error`

[pbastide]

Hi,

I encountered a small unexpected behavior when trying to a lower value for the sigma2_error parameter.

It seems that bounds or starting values on sigma2_error cannot be set. "Default" bounds are mentioned for sigma2_error/sigma2 in the doc, but apparently they cannot be over-written. If I try to set them, I get (rather uninformative) errors (see example below).

From what I understood reading the code, the lower.bound and upper.bound arguments in phylolm are assumed to be a single value, that can only be used for the parameter to optimize (optpar), and not sigma2_error. Parameter starting.value does accept a list however, but the value on sigma2_error seems to be ignored (from what I understood from l.130 of 'R/phylolm'.)

I made a small local patch to allow for setting bounds on this sigma2_error parameter, I'm happy to do a PR if you want. The only thing is that the bounds are actually on sigma2_error/sigma2 (or sigma2_error/sigma2*2*alpha for an OU).

I hope this helps, and thanks for this great package !

Paul

# Taken from phyolm manual example
set.seed(123456)
tre = rcoal(60)
taxa = sort(tre$tip.label)
b0=0; b1=1;
x <- rTrait(n=1, phy=tre,model="BM",
            parameters=list(ancestral.state=0,sigma2=10))
y <- b0 + b1*x +
  rTrait(n=1,phy=tre,model="kappa",parameters=list(
    ancestral.state=0,sigma2=1,kappa=0.5))
z <- y + rnorm(60,0,1)
dat = data.frame(trait=z[taxa],pred=x[taxa])

# Fitting
fit = phylolm(trait ~ pred, data = dat, phy = tre,
              model = "kappa", measurement_error = TRUE)

fit$optpar; fit$sigma2_error

[1] 0.8276196
[1] 1.062907

# Trying to set a bound on sigma2_error throws an error
fit = phylolm(trait ~ pred, data = dat, phy = tre,
              model = "kappa", measurement_error = TRUE,
              lower.bound = list(sigma2_error = 0.1))

Error in log(lower) : non-numeric argument to mathematical function

# Also when specifying kappa
fit = phylolm(trait ~ pred, data = dat, phy = tre,
              model = "kappa", measurement_error = TRUE,
              lower.bound = list(kappa = 0.5, sigma2_error = 0.1))

Error in log(lower) : non-numeric argument to mathematical function
In addition: Warning message:
In lower > start :
  longer object length is not a multiple of shorter object length

# Trying to set a starting value on sigma2_error throws an error
fit = phylolm(trait ~ pred, data = dat, phy = tre,
              model = "kappa", measurement_error = TRUE,
              starting.value = list(sigma2_error = 0.1))

Error in loglik(prm, y, X) : NA/NaN/Inf in foreign function call (arg 8)

# Does not fail, but starting value on sigma2_error gets ignored (l.130 of R/phylolm)
fit = phylolm(trait ~ pred, data = dat, phy = tre,
              model = "kappa", measurement_error = TRUE,
              starting.value = list(kappa = 0.5, sigma2_error = 0.1))
fit$optpar; fit$sigma2_error

[1] 0.8276196
[1] 1.062907

sessionInfo()

R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

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] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] phylolm_2.6.2  testthat_3.0.0 ape_5.4-1     

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.6         compiler_4.0.3     remotes_2.2.0      prettyunits_1.1.1  tools_4.0.3       
 [6] digest_0.6.27      pkgbuild_1.1.0     pkgload_1.1.0      memoise_1.1.0      nlme_3.1-150      
[11] lattice_0.20-41    rlang_0.4.10       cli_2.2.0          rstudioapi_0.13    parallel_4.0.3    
[16] withr_2.3.0        fs_1.5.0           desc_1.2.0         globals_0.14.0     devtools_2.3.2    
[21] rprojroot_2.0.2    grid_4.0.3         glue_1.4.2         listenv_0.8.0      R6_2.5.0          
[26] processx_3.4.5     future.apply_1.7.0 fansi_0.4.1        parallelly_1.23.0  sessioninfo_1.1.1 
[31] callr_3.5.1        magrittr_2.0.1     usethis_1.6.3      codetools_0.2-18   ps_1.4.0          
[36] ellipsis_0.3.1     assertthat_0.2.1   future_1.21.0      crayon_1.3.4      

[pbastide]

PS: Here is a possible patch commit, that allows for named vectors for bounds and starting values.

[lamho86]

Hi Paul,

It looks like when I add the code, I didn’t let the user to define the bounds for sigma^2_error. Thanks for upgrading this.

Best, Lam

Lam Ho | Assistant Professor FACULTY OF SCIENCE Department of Mathematics and Statistics 902.494.1069

Dalhousie University dal.ca

On Feb 22, 2021, at 9:29 AM, Paul Bastide notifications@github.com wrote:

PS: Here is a possible patch commit https://github.com/pbastide/phylolm/commit/b3787b6bc841d0e5a095519d273073edf67175a2, that allows for named vectors for bounds and starting values.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/lamho86/phylolm/issues/38#issuecomment-783374004, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABW3HPZSYFK7CPJMD3WQDKLTAJL2TANCNFSM4YADCCEA.

[pbastide]

Hi Lam,

Thank you for your answer, and for merging the change. I think it should solve the issue (see test file), so I'm closing this, but do not hesitate to check and warn me if you see any unexpected behavior.

Thanks, Paul