lamho86/phylolm

Compatibility between phylolm and visreg

Closed this issue · 2 comments

Hi! Congratulations for the package! And thank you for maintaining it!

I've tried to visualize phylolm objects with the visreg package: the visualization lacks the confidence band around the slope. This lack of features could arise from the fact that model.frame() appears to return a string "lambda", rather than a data frame, when applied to phylolm objects. Could this be fixed? I think that compatibility with visreg would be a useful feature in future releases of phylolm.

See a reproducible example below. Cheers!

library(ape)
library(ade4)
library(phylolm)
library(visreg)

# use dataset from ade4 package: life-history traits for 18 lizards
data(lizards)
# phylogenetic tree
l.tree <- read.tree(text = lizards$hprA)
# traits
l.traits <- lizards$traits


# Fit pgls
mod <- phylolm::phylolm(matur.L ~ age.mat, data = l.traits, phy = l.tree, model = "lambda")
#> Warning in phylolm::phylolm(matur.L ~ age.mat, data = l.traits, phy = l.tree, : the estimation of lambda matches the upper/lower bound for this parameter.
#>                           You may change the bounds using options "upper.bound" and "lower.bound".


summary(mod)
#> 
#> Call:
#> phylolm::phylolm(formula = matur.L ~ age.mat, data = l.traits, 
#>     phy = l.tree, model = "lambda")
#> 
#>    AIC logLik 
#> 153.23 -72.62 
#> 
#> Raw residuals:
#>     Min      1Q  Median      3Q     Max 
#> -18.233 -13.593   1.759   6.276  31.485 
#> 
#> Mean tip height: 33
#> Parameter estimate(s) using ML:
#> lambda : 1e-07
#> sigma2: 5.663128 
#> 
#> Coefficients:
#>             Estimate   StdErr t.value   p.value    
#> (Intercept)  4.17891 11.37494  0.3674 0.7181507    
#> age.mat      5.07030  0.99637  5.0888 0.0001095 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-squared: 0.6181    Adjusted R-squared: 0.5942 
#> 
#> Note: p-values and R-squared are conditional on lambda=1e-07.


# model.frame returns a string: "lambda"
model.frame(mod)
#> [1] "lambda"


# if I understand correctly, it should return something like this
model.frame(lm(matur.L ~ age.mat, data = l.traits))
#>    matur.L age.mat
#> Sa      58      13
#> Sh      42       5
#> Tl     132      19
#> Mc      56      11
#> My      60      10
#> Pb      44       8
#> Ph      39       8
#> Pg      49       8
#> Pa      53       7
#> Pm      49      12
#> Tt      46       7
#> Ts      56      13
#> Ae      61      10
#> Zv      44      11
#> Zo      43      11
#> La      62      15
#> Ls      91      16
#> Lv      84      12


# formula appears to work fine
formula(mod)
#> matur.L ~ age.mat


# a visualization without confidence band
visreg(mod)

sessionInfo()
#> R version 4.2.0 (2022-04-22 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 22000)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
#> [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
#> [5] LC_TIME=Spanish_Spain.utf8    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] visreg_2.7.0  phylolm_2.6.2 ade4_1.7-19   ape_5.6-2    
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.10        parallelly_1.31.1  rstudioapi_0.13    knitr_1.39        
#>  [5] magrittr_2.0.3     MASS_7.3-57        lattice_0.20-45    rlang_1.1.1       
#>  [9] fastmap_1.1.0      stringr_1.5.0      highr_0.9          globals_0.15.0    
#> [13] tools_4.2.0        parallel_4.2.0     grid_4.2.0         nlme_3.1-157      
#> [17] xfun_0.40          cli_3.6.1          withr_2.5.0        htmltools_0.5.6   
#> [21] yaml_2.3.5         digest_0.6.29      lifecycle_1.0.3    codetools_0.2-18  
#> [25] vctrs_0.6.3        fs_1.5.2           glue_1.6.2         evaluate_0.15     
#> [29] rmarkdown_2.14     reprex_2.0.2       future.apply_1.9.0 stringi_1.7.6     
#> [33] compiler_4.2.0     future_1.25.0      listenv_0.8.0

Created on 2023-09-28 with reprex v2.0.2

Hello @MarcRieraDominguez ,

The new version of phylolm should be compatible with visreg. If you note any problems, please let me know.

Best,
Lam

Hello @lamho86
That's great! I've just tested the reprex with phylolm ver. 2.6.5, and visreg now returns a confidence band! Thank you very much!
I'll be closing this.