Supplying formulas to domir programatically #7

MarcRieraDominguez · 2023-11-16T10:26:27Z

MarcRieraDominguez
Nov 16, 2023

Hi Joseph! Thank you for maintaining this awesome package! I would like to consult you on one issue.

I am currently working of Phylogentic Generalized Least Square models (PGLS, phylolm package), and want to run dominance analysis on them (domir::domir()).
The thing is that I need to run multiple variations on the same model, with different response variables, and different explanatory variables. Due to different data availability, the number of species changes between models. Therefore, I end up running 15 models, with around five different combinations of (data frame + phylogenetic tree).

I've found it convenient to use data frames to store formulas, data frames, and phylogenetic trees; and then use rowwise() to feed everything into the modelling function.
However, I am not able to "feed" the formulas into domir::domir(). It seems the issue is related to the formula being stored as a string. To fit models, I just put as.formula() inside the model call, but that doesn't seem to work with domir. See reprex below.
Any help much appreciated!
Incidentally, my knowledge of R is not so good, so perhaps I am not using "programatically" correctly, I am happy to change the title.

library(caper)
#> Loading required package: ape
#> Warning: package 'ape' was built under R version 4.2.3
#> Loading required package: MASS
#> Loading required package: mvtnorm
library(domir)
#> Warning: package 'domir' was built under R version 4.2.3
library(phylolm)
library(MuMIn)
#> Warning: package 'MuMIn' was built under R version 4.2.3
#> Registered S3 methods overwritten by 'MuMIn':
#>   method         from   
#>   nobs.pgls      caper  
#>   nobs.phylolm   phylolm
#>   logLik.phylolm phylolm
library(tidyverse)
#> Warning: package 'ggplot2' was built under R version 4.2.2
#> Warning: package 'tibble' was built under R version 4.2.3
#> Warning: package 'tidyr' was built under R version 4.2.3
#> Warning: package 'purrr' was built under R version 4.2.3
#> Warning: package 'dplyr' was built under R version 4.2.3
#> Warning: package 'stringr' was built under R version 4.2.3

## Create data
data(shorebird)
shorebird.data$M.Mass.std <- (shorebird.data$M.Mass - mean(shorebird.data$M.Mass))/sd(shorebird.data$M.Mass)
str(shorebird.data)
#> 'data.frame':    71 obs. of  7 variables:
#>  $ Species   : chr  "Actophilornis_africanus" "Aphriza_virgata" "Arenaria_interpres" "Arenaria_melanocephala" ...
#>  $ M.Mass    : num  143 186 108 114 151 ...
#>  $ F.Mass    : num  261 216 113 124 164 ...
#>  $ Egg.Mass  : num  8.6 22.4 17.9 17.3 23.5 11.2 10.7 9.6 19.3 13.3 ...
#>  $ Cl.size   : num  4 4 3.5 4 3.99 3.9 3.9 4 3.7 3.9 ...
#>  $ Mat.syst  : Factor w/ 3 levels "MO","PA","PG": 2 1 1 1 1 2 1 1 1 1 ...
#>  $ M.Mass.std: num  -0.1658 0.0666 -0.3557 -0.3255 -0.1238 ...
str(shorebird.tree)
#> List of 4
#>  $ edge       : num [1:132, 1:2] 72 73 74 75 76 77 78 79 79 80 ...
#>  $ edge.length: num [1:132] 18.1 13.4 13.4 2.6 2.3 ...
#>  $ tip.label  : chr [1:71] "Catoptrophorus_semipalmatus" "Tringa_ochropus" "Tringa_stagnatilis" "Tringa_flavipes" ...
#>  $ Nnode      : num 62
#>  - attr(*, "class")= chr "phylo"
#>  - attr(*, "origin")= chr "/Users/dorme/Scripts/CAIC/caic/pkg/data/shorebird.nex"

## Use a data frame to store formula and fit models programatically
# Note that I used rowwise() to fit models
mod <-
  data.frame(response = c("Egg.Mass", "Cl.size")) %>% 
  mutate(formula = paste0(response, " ~ M.Mass.std + Mat.syst"),
         df = list(shorebird.data),
         tree = list(shorebird.tree)) %>% 
  rowwise() %>% 
  mutate(phylolm.mod = list(phylolm::phylolm(as.formula(formula), data = df, phy = tree, model = "lambda")))

# The approach works fine to fit the models
lapply(mod$phylolm.mod, summary)
#> [[1]]
#> 
#> Call:
#> phylolm::phylolm(formula = as.formula(formula), data = df, phy = tree, 
#>     model = "lambda")
#> 
#>    AIC logLik 
#>  417.6 -202.8 
#> 
#> Raw residuals:
#>      Min       1Q   Median       3Q      Max 
#> -19.0780  -2.8232  -0.2333   2.7835  15.4826 
#> 
#> Mean tip height: 81.872
#> Parameter estimate(s) using ML:
#> lambda : 0.8908896
#> sigma2: 0.5773961 
#> 
#> Coefficients:
#>             Estimate   StdErr t.value   p.value    
#> (Intercept) 22.48761  3.38389  6.6455 6.495e-09 ***
#> M.Mass.std  15.29357  0.89905 17.0107 < 2.2e-16 ***
#> Mat.systPA  -2.46079  1.47511 -1.6682   0.09994 .  
#> Mat.systPG   4.21126  2.88749  1.4585   0.14939    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-squared: 0.8229    Adjusted R-squared: 0.815 
#> 
#> Note: p-values and R-squared are conditional on lambda=0.8908896.
#> 
#> [[2]]
#> 
#> Call:
#> phylolm::phylolm(formula = as.formula(formula), data = df, phy = tree, 
#>     model = "lambda")
#> 
#>    AIC logLik 
#>  75.62 -31.81 
#> 
#> Raw residuals:
#>      Min       1Q   Median       3Q      Max 
#> -1.71043 -0.03726  0.14459  0.27334  0.88305 
#> 
#> Mean tip height: 81.872
#> Parameter estimate(s) using ML:
#> lambda : 0.9936611
#> sigma2: 0.00703669 
#> 
#> Coefficients:
#>              Estimate    StdErr t.value   p.value    
#> (Intercept)  3.689908  0.388677  9.4935 5.041e-14 ***
#> M.Mass.std  -0.227910  0.086874 -2.6234   0.01077 *  
#> Mat.systPA  -0.114188  0.123347 -0.9257   0.35790    
#> Mat.systPG  -0.114788  0.234778 -0.4889   0.62650    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-squared: 0.1026    Adjusted R-squared: 0.06244 
#> 
#> Note: p-values and R-squared are conditional on lambda=0.9936611.

# However, I cannot "feed" the formula to domir() inside the data frame
dominance <-
  mod %>% 
  mutate(domir = list(domir::domir(formula,
                                   function(x) {
                                     summary(phylolm::phylolm(as.formula(x), data = df, phy = tree, model = "lambda"))$adj.r.squared})))
#> Error in `mutate()`:
#> ℹ In argument: `domir = list(...)`.
#> ℹ In row 1.
#> Caused by error in `UseMethod()`:
#> ! no applicable method for 'domir' applied to an object of class "character"
#> Backtrace:
#>      ▆
#>   1. ├─mod %>% ...
#>   2. ├─dplyr::mutate(...)
#>   3. ├─dplyr:::mutate.data.frame(...)
#>   4. │ └─dplyr:::mutate_cols(.data, dplyr_quosures(...), by)
#>   5. │   ├─base::withCallingHandlers(...)
#>   6. │   └─dplyr:::mutate_col(dots[[i]], data, mask, new_columns)
#>   7. │     └─mask$eval_all_mutate(quo)
#>   8. │       └─dplyr (local) eval()
#>   9. ├─domir::domir(...)
#>  10. └─base::.handleSimpleError(...)
#>  11.   └─dplyr (local) h(simpleError(msg, call))
#>  12.     └─rlang::abort(message, class = error_class, parent = parent, call = error_call)

# Changing the type to formula inside the function won't work
dominance <-
  mod %>% 
  mutate(domir = list(domir::domir(formula,
                                   function(x) {
                                     fml <- as.formula(formula)
                                     summary(phylolm::phylolm(fml, data = df, phy = tree, model = "lambda"))$adj.r.squared})))
#> Error in `mutate()`:
#> ℹ In argument: `domir = list(...)`.
#> ℹ In row 1.
#> Caused by error in `UseMethod()`:
#> ! no applicable method for 'domir' applied to an object of class "character"
#> Backtrace:
#>      ▆
#>   1. ├─mod %>% ...
#>   2. ├─dplyr::mutate(...)
#>   3. ├─dplyr:::mutate.data.frame(...)
#>   4. │ └─dplyr:::mutate_cols(.data, dplyr_quosures(...), by)
#>   5. │   ├─base::withCallingHandlers(...)
#>   6. │   └─dplyr:::mutate_col(dots[[i]], data, mask, new_columns)
#>   7. │     └─mask$eval_all_mutate(quo)
#>   8. │       └─dplyr (local) eval()
#>   9. ├─domir::domir(...)
#>  10. └─base::.handleSimpleError(...)
#>  11.   └─dplyr (local) h(simpleError(msg, call))
#>  12.     └─rlang::abort(message, class = error_class, parent = parent, call = error_call)



# For completion, domir() works on a phylolm model
domir::domir(Egg.Mass ~ M.Mass.std + Mat.syst,
             function(x) {
               summary(phylolm::phylolm(x, data = shorebird.data, phy = shorebird.tree, model = "lambda"))$adj.r.squared})
#> Overall Value:      0.8149597 
#> 
#> General Dominance Values:
#>            General Dominance Standardized Ranks
#> M.Mass.std         0.7982642   0.97951371     1
#> Mat.syst           0.0166955   0.02048629     2
#> 
#> Conditional Dominance Values:
#>            Subset Size: 1 Subset Size: 2
#> M.Mass.std      0.8051429    0.791385499
#> Mat.syst        0.0235742    0.009816789
#> 
#> Complete Dominance Designations:
#>                    Dmnated?M.Mass.std Dmnated?Mat.syst
#> Dmnates?M.Mass.std                 NA             TRUE
#> Dmnates?Mat.syst                FALSE               NA

# Converting to formula in the modelling function won't work
domir::domir("Egg.Mass ~ M.Mass.std + Mat.syst",
             function(x) {
               summary(phylolm::phylolm(as.formula(x), data = shorebird.data, phy = shorebird.tree, model = "lambda"))$adj.r.squared})
#> Error in UseMethod("domir"): no applicable method for 'domir' applied to an object of class "character"

# Changing the type to formula inside the function won't work
domir::domir("Egg.Mass ~ M.Mass.std + Mat.syst",
             function(x) {
               fml <- as.formula(x)
               summary(phylolm::phylolm(fml, data = shorebird.data, phy = shorebird.tree, model = "lambda"))$adj.r.squared})
#> Error in UseMethod("domir"): no applicable method for 'domir' applied to an object of class "character"

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] forcats_0.5.1   stringr_1.5.0   dplyr_1.1.2     purrr_1.0.1    
#>  [5] readr_2.1.2     tidyr_1.3.0     tibble_3.2.1    ggplot2_3.4.0  
#>  [9] tidyverse_1.3.1 MuMIn_1.47.5    phylolm_2.6.5   domir_1.1.0    
#> [13] caper_1.0.1     mvtnorm_1.1-3   MASS_7.3-57     ape_5.7-1      
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.11         lubridate_1.8.0     lattice_0.20-45    
#>  [4] listenv_0.9.0       assertthat_0.2.1    digest_0.6.33      
#>  [7] utf8_1.2.2          parallelly_1.36.0   cellranger_1.1.0   
#> [10] R6_2.5.1            backports_1.4.1     reprex_2.0.2       
#> [13] stats4_4.2.0        evaluate_0.15       httr_1.4.3         
#> [16] highr_0.9           pillar_1.9.0        rlang_1.1.1        
#> [19] readxl_1.4.0        rstudioapi_0.13     Matrix_1.4-1       
#> [22] rmarkdown_2.14      munsell_0.5.0       broom_0.8.0        
#> [25] modelr_0.1.8        compiler_4.2.0      xfun_0.40          
#> [28] pkgconfig_2.0.3     globals_0.16.2      htmltools_0.5.6    
#> [31] tidyselect_1.2.0    codetools_0.2-18    fansi_1.0.3        
#> [34] future_1.33.0       crayon_1.5.1        tzdb_0.3.0         
#> [37] dbplyr_2.1.1        withr_2.5.0         grid_4.2.0         
#> [40] jsonlite_1.8.0      nlme_3.1-157        gtable_0.3.0       
#> [43] lifecycle_1.0.3     DBI_1.1.2           magrittr_2.0.3     
#> [46] scales_1.2.1        future.apply_1.11.0 cli_3.6.1          
#> [49] stringi_1.7.6       fs_1.5.2            xml2_1.3.3         
#> [52] ellipsis_0.3.2      generics_0.1.2      vctrs_0.6.3        
#> [55] tools_4.2.0         glue_1.6.2          hms_1.1.1          
#> [58] parallel_4.2.0      fastmap_1.1.0       yaml_2.3.5         
#> [61] colorspace_2.0-3    rvest_1.0.2         knitr_1.39         
#> [64] haven_2.5.0

^{Created on 2023-11-16 with reprex v2.0.2}

Answered by jluchman

Nov 18, 2023

Hi again Marc,

I've shifted this issue over to be a discussion (issues are more for bug reports or improvement requests).

domir requires a formula (or formula_list) classed object is passed to it in the .obj argument and it does not have a way to work with a string, even if it looks like a formula. That is, the first argument/.obj to domir must be a formula and even if it is parsed as a formula within the function/.fct argument.

The reason for why this is the case gets technical but it suffices to say that if any formula-like string is passed to domir it just needs to be pulled through as.formula first before it is submitted.

For example, something like:

dominance <-
  mod %>% 
  mutate(d…

View full answer

jluchman · 2023-11-18T15:09:32Z

jluchman
Nov 18, 2023
Maintainer

Hi again Marc,

I've shifted this issue over to be a discussion (issues are more for bug reports or improvement requests).

domir requires a formula (or formula_list) classed object is passed to it in the .obj argument and it does not have a way to work with a string, even if it looks like a formula. That is, the first argument/.obj to domir must be a formula and even if it is parsed as a formula within the function/.fct argument.

The reason for why this is the case gets technical but it suffices to say that if any formula-like string is passed to domir it just needs to be pulled through as.formula first before it is submitted.

For example, something like:

dominance <-
  mod %>% 
  mutate(domir = list(domir::domir(as.formula(formula),
                                   function(x) {
                                     fml <- as.formula(formula)
                                     summary(phylolm::phylolm(fml, data = df, phy = tree, model = "lambda"))$adj.r.squared})))

May avoid the error you've been seeing.

joe

1 reply

MarcRieraDominguez Nov 28, 2023
Author

Hi Joe! Apologies for the delayed reply! The solution you suggetsed fix the issue, thank you very much! I'll be closing the dicussion.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Supplying formulas to domir programatically #7

{{title}}

Replies: 1 comment 1 reply

{{title}}

{{title}}

Select a reply

Supplying formulas to domir programatically #7

MarcRieraDominguez Nov 16, 2023

Replies: 1 comment · 1 reply

jluchman Nov 18, 2023 Maintainer

MarcRieraDominguez Nov 28, 2023 Author

MarcRieraDominguez
Nov 16, 2023

Replies: 1 comment 1 reply

jluchman
Nov 18, 2023
Maintainer

MarcRieraDominguez Nov 28, 2023
Author