logLik method and compatibility with MuMIn functions

[MarcRieraDominguez]

Hi! It looks like the logLik method implemented by phylolm reduces compatibility with MuMIn::model.avg(). If phylolm is loaded after MuMIn, it overwrittes the logLik method, and the model.avg(fit = TRUE) function won't work. Issues with the logLikmethod have been broached before: https://github.com/lamho86/phylolm/issues/56; https://github.com/lamho86/phylolm/pull/55, perhaps changes have not been pushed yet into the latest version?

Thank you for developping and maintaining the package!

1st MuMIn, 2nd phlyolm - problematic

library(caper)
#> Loading required package: ape
#> Loading required package: MASS
#> Loading required package: mvtnorm
library(MuMIn)
#> Warning: package 'MuMIn' was built under R version 4.2.3
#> Registered S3 method overwritten by 'MuMIn':
#>   method    from 
#>   nobs.pgls caper
library(phylolm)
#> Registered S3 methods overwritten by 'phylolm':
#>   method         from 
#>   logLik.phylolm MuMIn
#>   nobs.phylolm   MuMIn

## 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)
#> Error in vapply(logLiks, attr, 0, "df"): values must be length 1,
#>  but FUN(X[[1]]) result is length 0
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 -1.07181850 13.9060637
#> M.Mass      -0.01096584 -0.07182777  0.0498961
#> F.Mass       0.08434847  0.02988806  0.1388089

cbind(coef(mod.avg.fit, full = TRUE), confint(mod.avg.fit, full = TRUE))
#> Error in coef(mod.avg.fit, full = TRUE): object 'mod.avg.fit' not found
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

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] phylolm_2.6.4 MuMIn_1.47.5  caper_1.0.1   mvtnorm_1.1-3 MASS_7.3-57  
#> [6] ape_5.6-2    
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.10         compiler_4.2.0      highr_0.9          
#>  [4] tools_4.2.0         digest_0.6.29       evaluate_0.15      
#>  [7] lifecycle_1.0.3     nlme_3.1-157        lattice_0.20-45    
#> [10] rlang_1.1.1         reprex_2.0.2        Matrix_1.4-1       
#> [13] cli_3.6.1           rstudioapi_0.13     yaml_2.3.5         
#> [16] parallel_4.2.0      xfun_0.40           fastmap_1.1.0      
#> [19] withr_2.5.0         stringr_1.5.0       knitr_1.39         
#> [22] fs_1.5.2            vctrs_0.6.3         globals_0.16.2     
#> [25] stats4_4.2.0        grid_4.2.0          glue_1.6.2         
#> [28] listenv_0.9.0       parallelly_1.36.0   future.apply_1.11.0
#> [31] rmarkdown_2.14      magrittr_2.0.3      codetools_0.2-18   
#> [34] htmltools_0.5.6     future_1.33.0       stringi_1.7.6

^{Created on 2023-10-09 with reprex v2.0.2}

1st phylolm, 2nd MuMIn - works fine

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

## 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 -1.07181850 13.9060637
#> M.Mass      -0.01096584 -0.07182777  0.0498961
#> F.Mass       0.08434847  0.02988806  0.1388089

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

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.4 caper_1.0.1   mvtnorm_1.1-3 MASS_7.3-57  
#> [6] ape_5.6-2    
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.10         compiler_4.2.0      highr_0.9          
#>  [4] tools_4.2.0         digest_0.6.29       evaluate_0.15      
#>  [7] lifecycle_1.0.3     nlme_3.1-157        lattice_0.20-45    
#> [10] rlang_1.1.1         Matrix_1.4-1        reprex_2.0.2       
#> [13] cli_3.6.1           rstudioapi_0.13     yaml_2.3.5         
#> [16] parallel_4.2.0      xfun_0.40           fastmap_1.1.0      
#> [19] withr_2.5.0         stringr_1.5.0       knitr_1.39         
#> [22] fs_1.5.2            vctrs_0.6.3         globals_0.16.2     
#> [25] stats4_4.2.0        grid_4.2.0          glue_1.6.2         
#> [28] listenv_0.9.0       future.apply_1.11.0 parallelly_1.36.0  
#> [31] rmarkdown_2.14      magrittr_2.0.3      codetools_0.2-18   
#> [34] htmltools_0.5.6     future_1.33.0       stringi_1.7.6

^{Created on 2023-10-09 with reprex v2.0.2}

[Ax3man]

This happens because both phylolm and MuMIn provide their own versions of logLik.phylolm and logLik.phyloglm. They aren't compatible because the former is a list object, while the latter uses a numeric vector. Looking at stats::logLik.default, it seems that MuMIn follows the convention.

One solution is to update the logLik (and associated) functions in phylolm, which should improve it's compatibility with other packages in general. PR #66 implements this.

[lamho86]

Thanks @Ax3man for the fix. Although this would be a breaking change as mentioned in #66 but it would make phylolm more compatible with other packages.