I am working with model-averaged phylogenetic generalized least square models (PGLS), and I want to visualize model's predictions, with uncertainty around the prediction. I use phylolm::phylolm() for its speed, and MuMIn for the model-averaging procedure. I get an error with I try to model the averaging object with ggeffects::ggpredict(). Moreover, it appears that the predict() method of phylolm::phylolm() does not fave a se.fit = TRUE option. Can anyone recommend a way to plot the predicted trends from model-averaged PGLS, together with their uncertainty?
See below a reproducible example.
library(caper)
#> Warning: package 'caper' was built under R version 4.1.3
#> Loading required package: ape
#> Warning: package 'ape' was built under R version 4.1.3
#> Loading required package: MASS
#> Loading required package: mvtnorm
#> Warning: package 'mvtnorm' was built under R version 4.1.1
library(ggeffects)
library(phylolm)
#> Warning: package 'phylolm' was built under R version 4.1.3
library(MuMIn)
#> Warning: package 'MuMIn' was built under R version 4.1.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
## 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)
# Plot model-averaged object
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'
# Predict() with se.fit
summary(predict(mod, se.fit = TRUE))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 8.067 11.045 15.189 21.152 25.001 65.624
summary(predict(mod.avg.fit, se.fit = TRUE))
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 8.097 11.075 15.195 21.197 24.707 65.444
sessionInfo()
#> R version 4.1.0 (2021-05-18)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19045)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=Spanish_Spain.1252 LC_CTYPE=Spanish_Spain.1252
#> [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C
#> [5] LC_TIME=Spanish_Spain.1252
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] MuMIn_1.46.0 phylolm_2.6.2 ggeffects_1.3.1 caper_1.0.1
#> [5] mvtnorm_1.1-3 MASS_7.3-54 ape_5.6-2
#>
#> loaded via a namespace (and not attached):
#> [1] Rcpp_1.0.8.3 compiler_4.1.0 highr_0.9 tools_4.1.0
#> [5] digest_0.6.29 evaluate_0.15 lifecycle_1.0.3 nlme_3.1-152
#> [9] lattice_0.20-44 rlang_1.1.0 Matrix_1.3-3 reprex_2.0.1
#> [13] cli_3.6.0 rstudioapi_0.13 yaml_2.3.5 parallel_4.1.0
#> [17] xfun_0.30 fastmap_1.1.0 withr_2.5.0 stringr_1.5.0
#> [21] knitr_1.38 fs_1.5.2 vctrs_0.6.1 globals_0.15.0
#> [25] stats4_4.1.0 grid_4.1.0 glue_1.6.2 listenv_0.8.0
#> [29] future.apply_1.9.0 parallelly_1.31.1 rmarkdown_2.13 magrittr_2.0.3
#> [33] codetools_0.2-18 htmltools_0.5.2 insight_0.19.1 future_1.25.0
#> [37] stringi_1.7.6
Created on 2023-10-20 by the reprex package (v2.0.1)
The predict() method of phylolm now includes a se.fit = TRUE option, from version 2.6.5 onwards. That allows to visualize uncertainty around model-averaged predictions.