Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Error with model-averaged phylolm object #387

Closed
MarcRieraDominguez opened this issue Oct 9, 2023 · 1 comment
Closed

Error with model-averaged phylolm object #387

MarcRieraDominguez opened this issue Oct 9, 2023 · 1 comment
Labels
bug 🐛 Something isn't working

Comments

@MarcRieraDominguez
Copy link

Hi! Congratulations on the great package! I've com across an error when trying to plot predictions from model-averaged phylolm objects. While ggpredict will plot a phylolm model just fine, it won't plot a model-averaged one. For completion, I show that the error message also appears with model.avg(fit = FALSE).
Cant his be fixed?

Note that the order in which you install the packages is potentially important: this reprex won't work if you load phylolm after MuMIn, I think due to differences in the logLik method. I'll notify the developpers of phylolm of such behaviour.

library(caper)
#> Loading required package: ape
#> Loading required package: MASS
#> Loading required package: mvtnorm
library(ggeffects)
#> Warning: package 'ggeffects' 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

## Create data
data(shorebird)
summary(shorebird.data)
#>    Species              M.Mass           F.Mass         Egg.Mass    
#>  Length:71          Min.   : 20.30   Min.   : 22.2   Min.   : 5.80  
#>  Class :character   1st Qu.: 54.45   1st Qu.: 62.0   1st Qu.: 9.85  
#>  Mode  :character   Median :108.30   Median :117.0   Median :16.50  
#>                     Mean   :173.95   Mean   :197.3   Mean   :22.45  
#>                     3rd Qu.:181.25   3rd Qu.:243.3   3rd Qu.:29.20  
#>                     Max.   :740.30   Max.   :788.0   Max.   :76.00  
#>     Cl.size      Mat.syst
#>  Min.   :1.700   MO:54   
#>  1st Qu.:3.700   PA:15   
#>  Median :3.900   PG: 2   
#>  Mean   :3.658           
#>  3rd Qu.:4.000           
#>  Max.   :4.100

summary(shorebird.tree)
#> 
#> Phylogenetic tree: shorebird.tree 
#> 
#>   Number of tips: 71 
#>   Number of nodes: 62 
#>   Branch lengths:
#>     mean: 10.12843 
#>     variance: 84.77508 
#>     distribution summary:
#>    Min. 1st Qu.  Median 3rd Qu.    Max. 
#>   0.435   3.483   7.400  13.762  52.100 
#>   No root edge.
#>   First ten tip labels: Catoptrophorus_semipalmatus 
#>                         Tringa_ochropus
#>                         Tringa_stagnatilis
#>                         Tringa_flavipes
#>                         Tringa_nebularia
#>                         Tringa_totanus
#>                         Tringa_erythropus
#>                         Tringa_melanoleuca
#>                         Tringa_glareola
#>                         Steganopus_tricolor
#>   No node labels.


## Fit model
mod <- phylolm::phylolm(Egg.Mass ~ M.Mass + F.Mass, data = shorebird.data, phy = shorebird.tree, model = "lambda")


## Model average
options(na.action = "na.fail")
mod.d <- MuMIn::dredge(mod, rank = "AICc")
#> Fixed term is "(Intercept)"
mod.avg.fit <- MuMIn::model.avg(mod.d, revised.var = TRUE, fit = TRUE)
mod.avg.no.fit <- MuMIn::model.avg(mod.d, revised.var = TRUE, fit = FALSE)


## Compare coefficients
cbind(coef(mod), confint(mod))
#>                               2.5 %     97.5 %
#> (Intercept)  6.41712259 -0.93856939 13.7728146
#> M.Mass      -0.01096584 -0.07074487  0.0488132
#> F.Mass       0.08434847  0.03085706  0.1378399

cbind(coef(mod.avg.fit, full = TRUE), confint(mod.avg.fit, full = TRUE))
#>                                2.5 %      97.5 %
#> (Intercept)  6.437285624 -0.79029863 13.66486988
#> F.Mass       0.076339943  0.04345065  0.10922924
#> M.Mass      -0.001736286 -0.03770974  0.03423716
cbind(coef(mod.avg.no.fit, full = TRUE), confint(mod.avg.no.fit, full = TRUE))
#>                                2.5 %      97.5 %
#> (Intercept)  6.437285624 -0.79029863 13.66486988
#> F.Mass       0.076339943  0.04345065  0.10922924
#> M.Mass      -0.001736286 -0.03770974  0.03423716


plot( ggeffects::ggpredict(mod, terms = c("M.Mass")) )

plot( ggeffects::ggpredict(mod.avg.fit, terms = c("M.Mass")) )
#> Warning: Could not access model information.
#> Warning: Could not access model information.
#> Error in !is.null(model_info) && model_info$is_trial: invalid 'y' type in 'x && y'


plot( ggeffects::ggpredict(mod.avg.no.fit, terms = c("M.Mass")) )
#> Warning: Can't calculate covariance matrix. Please use 'fit = TRUE' in
#>   'model.avg()'.
#> Warning: Can't retrieve data. Please use `fit = TRUE` in `model.avg()`.

#> Warning: Can't retrieve data. Please use `fit = TRUE` in `model.avg()`.
#> Warning: Can't calculate covariance matrix. Please use 'fit = TRUE' in
#>   'model.avg()'.
#> Error in !is.null(model_info) && model_info$is_trial: invalid 'y' type in 'x && y'


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] MuMIn_1.47.5      phylolm_2.6.2     ggeffects_1.3.1.5 caper_1.0.1      
#> [5] mvtnorm_1.1-3     MASS_7.3-57       ape_5.6-2        
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.0   sjlabelled_1.2.0   xfun_0.40          listenv_0.8.0     
#>  [5] haven_2.5.0        lattice_0.20-45    snakecase_0.11.0   generics_0.1.2    
#>  [9] colorspace_2.0-3   vctrs_0.6.3        htmltools_0.5.6    stats4_4.2.0      
#> [13] yaml_2.3.5         utf8_1.2.2         rlang_1.1.1        pillar_1.9.0      
#> [17] glue_1.6.2         withr_2.5.0        lifecycle_1.0.3    stringr_1.5.0     
#> [21] munsell_0.5.0      gtable_0.3.0       future_1.25.0      codetools_0.2-18  
#> [25] evaluate_0.15      labeling_0.4.2     knitr_1.39         forcats_0.5.1     
#> [29] fastmap_1.1.0      datawizard_0.7.0   parallel_4.2.0     fansi_1.0.3       
#> [33] highr_0.9          Rcpp_1.0.10        scales_1.2.1       farver_2.1.0      
#> [37] parallelly_1.31.1  fs_1.5.2           ggplot2_3.4.0      hms_1.1.1         
#> [41] digest_0.6.29      stringi_1.7.6      dplyr_1.1.2        insight_0.19.5    
#> [45] grid_4.2.0         cli_3.6.1          tools_4.2.0        magrittr_2.0.3    
#> [49] tibble_3.2.1       future.apply_1.9.0 pkgconfig_2.0.3    ellipsis_0.3.2    
#> [53] Matrix_1.4-1       reprex_2.0.2       rmarkdown_2.14     rstudioapi_0.13   
#> [57] R6_2.5.1           globals_0.15.0     nlme_3.1-157       compiler_4.2.0

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

@strengejacke
Copy link
Owner

Thanks! This issue will be fixed once this PR is merged:
easystats/insight#975

library(caper)
#> Loading required package: ape
#> Loading required package: MASS
#> Loading required package: mvtnorm
library(ggeffects)
library(phylolm)
library(MuMIn)
#> Registered S3 methods overwritten by 'MuMIn':
#>   method         from   
#>   nobs.pgls      caper  
#>   nobs.phylolm   phylolm
#>   logLik.phylolm phylolm
## Create data
data(shorebird)
## Fit model
mod <- phylolm::phylolm(Egg.Mass ~ M.Mass + F.Mass, data = shorebird.data, phy = shorebird.tree, model = "lambda")

## Model average
options(na.action = "na.fail")
mod.d <- MuMIn::dredge(mod, rank = "AICc")
#> Fixed term is "(Intercept)"
mod.avg.fit <- MuMIn::model.avg(mod.d, revised.var = TRUE, fit = TRUE)
mod.avg.no.fit <- MuMIn::model.avg(mod.d, revised.var = TRUE, fit = FALSE)

ggeffects::ggpredict(mod.avg.fit, terms = c("M.Mass"))
#> # Predicted values of Egg.Mass
#> 
#> M.Mass | Predicted |       95% CI
#> ---------------------------------
#>      0 |     19.38 | 10.69, 28.08
#>    100 |     19.21 | 12.05, 26.38
#>    200 |     19.04 | 11.76, 26.31
#>    300 |     18.86 |  9.90, 27.83
#>    350 |     18.78 |  8.59, 28.96
#>    450 |     18.60 |  5.56, 31.64
#>    550 |     18.43 |  2.24, 34.62
#>    750 |     18.08 | -4.81, 40.98
#> 
#> Adjusted for:
#> * F.Mass = 169.60
#> 
#> Not all rows are shown in the output. Use `print(..., n = Inf)` to show
#>   all rows.
ggeffects::ggpredict(mod.avg.no.fit, terms = c("M.Mass"))
#> Warning: Can't access model information. Please use 'fit = TRUE' in
#>   'model.avg()'.
#> Warning: Can't access model information. Please use 'fit = TRUE' in
#>   'model.avg()'.
#> Warning: Can't retrieve data. Please use `fit = TRUE` in `model.avg()`.
#> Warning: Can't retrieve data. Please use `fit = TRUE` in `model.avg()`.
#> Warning: Can't access model information. Please use 'fit = TRUE' in
#>   'model.avg()'.
#> Error in names(datlist) <- names(focal_terms): 'names' attribute [1] must be the same length as the vector [0]

Created on 2024-11-26 with reprex v2.1.1

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug 🐛 Something isn't working
Projects
None yet
Development

No branches or pull requests

2 participants