系统发育线性模型（PGLS）：如何可视化模型预测的不确定性？

Question

我正在使用模型平均系统发育广义最小二乘模型（PGLS），我想可视化模型的预测，并且预测存在不确定性。我使用

phylolm::phylolm()

来表示速度，使用

MuMIn

来表示模型平均过程。我尝试使用

averaging

为

ggeffects::ggpredict()

对象建模时遇到错误。此外，似乎

predict()

的

phylolm::phylolm()

方法不喜欢

se.fit = TRUE

选项。谁能推荐一种方法来绘制模型平均 PGLS 的预测趋势及其不确定性？

请参阅下面的可重现示例。

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

^{由 reprex 包于 2023 年 10 月 20 日创建（v2.0.1）}

Answer 1

从 2.6.5 版本开始，phylolm 的

predict()

方法现在包含

se.fit = TRUE

选项。这允许可视化模型平均预测的不确定性。