To install this package:
devtools::install_github("daijiang/phyr")
# or install from binary file (may not be the latest version)
# macOS
install.packages("https://raw.githubusercontent.com/daijiang/phyr/master/phyr_1.0.3.tgz", repos = NULL)
# Windows
install.packages("https://raw.githubusercontent.com/daijiang/phyr/master/phyr_0.1.6.zip", repos = NULL)The phyr package has three groups of functions:
psv,
psr, pse, etc. and beta: pcd),
which were included in the picante package originally. They
were updated with c++ to improve speed.cor_phylo),
which was included in the ape package originally. It has
new syntax, much improved performance (c++), and bootstrapping
option.pglmm),
which was originally included in the pez package. It has
new model formula syntax that allows straightforward model set up, a
faster version of maximum likelihood implementation via c++, and a
Bayesian model fitting framework based on INLA.
brms for
Stan).pglmm.compare), which was
originally from ape::binaryPGLMM() but has more
features.pglmm()pglmm use similar syntax as lme4::lmer to
specify random terms: add __ (two underscores) at the end
of grouping variable (e.g. sp) to specify both phylogenetic
and non-phylogenetic random terms; use (1|sp__@site) to
specify nested term (i.e. species phylogenetic matrix V_sp
nested within the diagonal of site matrix I_site) to test
phylogenetic overdispersion or underdispersion. This should be the most
commonly used one and is equal to
kronecker(I_site, V_sp).
We can also use a second phylogeny for bipartite questions. For
example, (1|parasite@host__) will be converted to
kronecker(V_host, I_parasite);
(1|parasite__@host__) will be converted to
kronecker(V_host, V_parasite).
For details about model formula, see documentation
?phyr::pglmm. More application examples can be found in Ives 2018 Chapter 4.
library(phyr)library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
comm = comm_a
comm$site = row.names(comm)
dat = tidyr::gather(comm, key = "sp", value = "freq", -site) %>%
left_join(envi, by = "site") %>%
left_join(traits, by = "sp")
dat$pa = as.numeric(dat$freq > 0)
head(dat)
## site sp freq sand shade precip tmin sla veg.height
## 1 s3293 Acer_rubrum 0 80.75 20.9 1.902397 0.1288019 294 170.5
## 2 s3294 Acer_rubrum 3 83.36 45.1 1.902397 0.1288019 294 170.5
## 3 s3295 Acer_rubrum 8 88.83 58.9 1.922669 -0.1061756 294 170.5
## 4 s3296 Acer_rubrum 0 91.24 19.7 1.922669 -0.1061756 294 170.5
## 5 s3297 Acer_rubrum 0 90.04 56.6 1.922669 -0.1061756 294 170.5
## 6 s3299 Acer_rubrum 15 81.87 87.0 1.899665 0.1736423 294 170.5
## disp.mode pa
## 1 Wind 0
## 2 Wind 1
## 3 Wind 1
## 4 Wind 0
## 5 Wind 0
## 6 Wind 1
# phy-LMM
test1 = phyr::pglmm(freq ~ 1 + shade + (1|sp__) + (1|site) + (1|sp__@site),
data = dat, family = "gaussian", REML = FALSE,
cov_ranef = list(sp = phylotree))
## Warning: Drop species from the phylogeny that are not in the variable sp
test1
## Linear mixed model fit by maximum likelihood
##
## Call:freq ~ 1 + shade
##
## logLik AIC BIC
## -463.3 940.6 956.5
##
## Random effects:
## Variance Std.Dev
## 1|sp 7.345e-01 0.8570105
## 1|sp__ 1.800e-04 0.0134157
## 1|site 1.035e-07 0.0003217
## 1|sp__@site 2.138e-05 0.0046238
## residual 3.261e+00 1.8058430
##
## Fixed effects:
## Value Std.Error Zscore Pvalue
## (Intercept) -0.1911039 0.3920853 -0.4874 0.625972
## shade 0.0226917 0.0067263 3.3736 0.000742 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# phy-GLMM
test2 = phyr::pglmm(pa ~ 1 + shade + (1|sp__) + (1|site) + (1|sp__@site),
data = dat, family = "binomial", REML = FALSE,
cov_ranef = list(sp = phylotree))
## Warning: Drop species from the phylogeny that are not in the variable sp
test2
## Generalized linear mixed model for binomial data fit by maximum likelihood
##
## Call:pa ~ 1 + shade
##
##
## Random effects:
## Variance Std.Dev
## 1|sp 1.786e-06 0.001336
## 1|sp__ 4.441e-01 0.666389
## 1|site 4.496e-06 0.002120
## 1|sp__@site 8.689e-06 0.002948
##
## Fixed effects:
## Value Std.Error Zscore Pvalue
## (Intercept) -2.0835724 0.5744500 -3.6271 0.0002867 ***
## shade 0.0165916 0.0087165 1.9035 0.0569784 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# bipartite
tree_site = ape::rtree(n = n_distinct(dat$site), tip.label = sort(unique(dat$site)))
z_bipartite = phyr::pglmm(freq ~ 1 + shade + (1|sp__) + (1|site__) +
(1|sp__@site) + (1|sp@site__) + (1|sp__@site__),
data = dat, family = "gaussian",REML = TRUE,
cov_ranef = list(sp = phylotree, site = tree_site))
## Warning: Drop species from the phylogeny that are not in the variable sp
z_bipartite
## Linear mixed model fit by restricted maximum likelihood
##
## Call:freq ~ 1 + shade
##
## logLik AIC BIC
## -466.0 952.1 974.8
##
## Random effects:
## Variance Std.Dev
## 1|sp 1.648e-02 0.128377
## 1|sp__ 1.173e+00 1.082923
## 1|site 2.792e-02 0.167098
## 1|site__ 8.659e-03 0.093052
## 1|sp__@site 1.965e+00 1.401671
## 1|sp@site__ 7.968e-02 0.282273
## 1|sp__@site__ 8.041e-05 0.008967
## residual 9.625e-01 0.981064
##
## Fixed effects:
## Value Std.Error Zscore Pvalue
## (Intercept) -0.127328 0.815075 -0.1562 0.8759
## shade 0.019393 0.011889 1.6311 0.1029Licensed under the GPL-3 license.
Contributions are welcome. You can provide comments and feedback or ask questions by filing an issue on Github here or making pull requests.
Please note that the ‘phyr’ project is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.