Generalized linear models
We also compare results among software for fitting phylogenetic generalized linear models (PGLM).
Poisson-distributed response
First, we specifically explore a Poisson-distributed PGLM, comparing
phylosem against phylolm::phyloglm (which uses
Generalized Estimating Equations) and phyr::pglmm_compare
(which uses maximum likelihood).
# Settings
Ntree = 100
sd_x = 0.3
sd_y = 0.3
b0_x = 1
b0_y = 0
b_xy = 1
# Simulate tree
set.seed(1)
tree = ape::rtree(n=Ntree)
# Simulate data
x = b0_x + sd_x * phylolm::rTrait(n = 1, phy=tree)
ybar = b0_y + b_xy*x
y_normal = ybar + sd_y * phylolm::rTrait(n = 1, phy=tree)
y_pois = rpois( n=Ntree, lambda=exp(y_normal) )
# Construct, re-order, and reduce data
Data = data.frame(x=x,y=y_pois)
# Compare using phylolm::phyloglm
pglm = phylolm::phyloglm(y ~ 1 + x, data=Data, phy=tree, method="poisson_GEE" )
knitr::kable(summary(pglm)$coefficients, digits=3)| Estimate | StdErr | z.value | p.value | |
|---|---|---|---|---|
| (Intercept) | -1.098 | 0.633 | -1.736 | 0.083 |
| x | 1.314 | 0.247 | 5.320 | 0.000 |
#
pglmm = phyr::pglmm_compare(
y ~ 1 + x,
family = "poisson",
data = Data,
phy = tree )
knitr::kable(summary(pglmm), digits=3)
#> Generalized linear mixed model for poisson data fit by restricted maximum likelihood
#>
#> Call:y ~ 1 + x
#>
#> logLik AIC BIC
#> -173.4 354.7 360.6
#>
#> Phylogenetic random effects variance (s2):
#> Variance Std.Dev
#> s2 0.05511 0.2348
#>
#> Fixed effects:
#> Value Std.Error Zscore Pvalue
#> (Intercept) -0.57009 0.30469 -1.8710 0.06134 .
#> x 1.18137 0.19807 5.9645 2.454e-09 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#
pgsem = phylosem( sem = "x -> y, p",
data = Data,
family = c("fixed","poisson"),
tree = tree,
control = phylosem_control(quiet = TRUE) )
knitr::kable(summary(pgsem)$coefficients, digits=3)| Path | VarName | Estimate | StdErr | t.value | p.value |
|---|---|---|---|---|---|
| NA | Intercept_x | 1.087 | 0.183 | 5.945 | 0.000 |
| NA | Intercept_y | -0.581 | 0.305 | 1.904 | 0.057 |
| x -> y | p | 1.190 | 0.199 | 5.971 | 0.000 |
| x <-> x | V[x] | 0.315 | 0.022 | 14.072 | 0.000 |
| y <-> y | V[y] | 0.232 | 0.054 | 4.310 | 0.000 |
We also compare results against brms (which fits a
Bayesian hierarchical model), although we load results from compiled run
of brms to avoid users having to install STAN to run
vignettes for phylosem:
# Comare using Bayesian implementation in brms
library(brms)
Amat <- ape::vcv.phylo(tree)
Data$tips <- rownames(Data)
mcmc <- brm(
y ~ 1 + x + (1 | gr(tips, cov = A)),
data = Data, data2 = list(A = Amat),
family = 'poisson',
cores = 4
)
knitr::kable(fixef(mcmc), digits = 3)
# Plot them together
library(ggplot2)
pdat <- rbind.data.frame(
coef(summary(pglm))[, 1:2],
data.frame(Estimate = pglmm$B, StdErr = pglmm$B.se),
setNames(as.data.frame(fixef(mcmc))[1:2], c('Estimate', 'StdErr')),
setNames(summary(pgsem)$coefficients[2:3, 3:4], c('Estimate', 'StdErr'))
)
pdat$Param <- rep(c('Intercept', 'Slope'), 4)
pdat$Method <- rep( c('phylolm::phyloglm', 'phyr::pglmm_compare',
'brms::brm', 'phylosem::phylosem'), each = 2)
figure = ggplot(pdat, aes(
x = Estimate, xmin = Estimate - StdErr,
xmax = Estimate + StdErr, y = Param, color = Method
)) +
geom_pointrange(position = position_dodge(width = 0.6)) +
theme_classic() +
theme(panel.grid.major.x = element_line(), panel.grid.minor.x = element_line())#> Warning: package 'ggplot2' was built under R version 4.3.3
In this instance (and in others we have explored), results from
phylolm::phyloglm are generally different while those from
phylosem, phyr::pglmm_compare, and
brms are close but not quite identical.
Binomial regression
We also compare results for a Bernoulli-distributed response using
PGLM. We again compare phylosem against
phyr::pglmm_compare, and do not explore threshold models
which we expect to give different results due differences in assumptions
about how latent variables affect measurements.
# Settings
Ntree = 100
sd_x = 0.3
sd_y = 0.3
b0_x = 1
b0_y = 0
b_xy = 1
# Simulate tree
set.seed(1)
tree = ape::rtree(n=Ntree)
# Simulate data
x = b0_x + sd_x * phylolm::rTrait(n = 1, phy=tree)
ybar = b0_y + b_xy*x
y_normal = ybar + sd_y * phylolm::rTrait(n = 1, phy=tree)
y_binom = rbinom( n=Ntree, size=1, prob=plogis(y_normal) )
# Construct, re-order, and reduce data
Data = data.frame(x=x,y=y_binom)
#
pglmm = phyr::pglmm_compare(
y ~ 1 + x,
family = "binomial",
data = Data,
phy = tree )
knitr::kable(summary(pglmm), digits=3)
#> Generalized linear mixed model for binomial data fit by restricted maximum likelihood
#>
#> Call:y ~ 1 + x
#>
#> logLik AIC BIC
#> -63.74 135.47 141.32
#>
#> Phylogenetic random effects variance (s2):
#> Variance Std.Dev
#> s2 0.1076 0.328
#>
#> Fixed effects:
#> Value Std.Error Zscore Pvalue
#> (Intercept) 0.23179 0.60507 0.3831 0.7017
#> x 0.44548 0.45708 0.9746 0.3297
#
pgsem = phylosem( sem = "x -> y, p",
data = Data,
family = c("fixed","binomial"),
tree = tree,
control = phylosem_control(quiet = TRUE) )
knitr::kable(summary(pgsem)$coefficients, digits=3)| Path | VarName | Estimate | StdErr | t.value | p.value |
|---|---|---|---|---|---|
| NA | Intercept_x | 1.087 | 0.183 | 5.945 | 0.000 |
| NA | Intercept_y | 0.204 | 0.589 | 0.346 | 0.730 |
| x -> y | p | 0.458 | 0.468 | 0.977 | 0.328 |
| x <-> x | V[x] | -0.315 | 0.022 | 14.072 | 0.000 |
| y <-> y | V[y] | 0.290 | 0.284 | 1.020 | 0.308 |
In this instance, phylosem and
phyr::pglmm_compare give similar estimates and standard
errors for the slope term.
Summary of PGLM results
Based on these two comparisons, we conclude that phylosem provides an interface for maximum-likelihood estimate of phylogenetic generalized linear models (PGLM), and extends this class to include mixed data (i.e., a combination of different measurement types), missing data, and non-recursive structural linkages. However, we also encourage further cross-testing of different software for fitting phylogenetic generalized linear models.
