library(tidyverse)
library(ape)
library(tidybayes)
library(brms)
library(ggtree)
library(sjPlot)
speciesdat <- read_csv("completedataR.csv")
## Parsed with column specification:
## cols(
## Family = col_character(),
## Genus = col_character(),
## Species = col_character(),
## TAG = col_double(),
## DBH = col_double(),
## Colony = col_double(),
## Foraged = col_double(),
## Distance = col_double(),
## Succession = col_character()
## )
tre <- read.nexus("treenew.nexus") #Output from Phylomatic
inv.phylo <- MCMCglmm::inverseA(tre, nodes = "TIPS", scale = FALSE)
A <- solve(inv.phylo$Ainv)
rownames(A) <- rownames(inv.phylo$Ainv)
priorscomplete <- get_prior(Foraged ~ DBH + Distance + (1 | Colony) + (1 | Species) + (1 | TAG), data = speciesdat, family = bernoulli())
## Warning: Rows containing NAs were excluded from the model.
complete <- brm(Foraged ~ DBH + Distance + (1 | Colony) + (1 | Species) + (1 | TAG), data = speciesdat, family = bernoulli, cov_ranef = list(phylo = A), save_all_pars = TRUE, cores = 4, control = list(adapt_delta = 0.9), prior = priorscomplete, iter = 4000)
## Warning: Rows containing NAs were excluded from the model.
## Compiling the C++ model
## Start sampling
## Warning: There were 4 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
summary(complete)
## Warning: There were 4 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help.
## See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Family: bernoulli
## Links: mu = logit
## Formula: Foraged ~ DBH + Distance + (1 | Colony) + (1 | Species) + (1 | TAG)
## Data: speciesdat (Number of observations: 820)
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup samples = 8000
##
## Group-Level Effects:
## ~Colony (Number of levels: 8)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 0.62 0.46 0.05 1.74 1023 1.00
##
## ~Species (Number of levels: 224)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 1.45 0.93 0.15 3.77 435 1.01
##
## ~TAG (Number of levels: 724)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 3.30 1.84 0.67 7.88 251 1.01
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept -4.40 1.92 -9.25 -2.05 298 1.01
## DBH 0.07 0.03 0.03 0.15 312 1.01
## Distance -0.07 0.04 -0.17 -0.02 348 1.01
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(complete)
r2 <- list()
(r2[["complete"]] <- bayes_R2(complete))
## Estimate Est.Error Q2.5 Q97.5
## R2 0.5149481 0.1474828 0.1918668 0.7728898
tab_model(complete, show.intercept = FALSE, show.r2 = TRUE)
 | Foraged | ||
---|---|---|---|
Predictors | Odds Ratios | CI (50%) | CI (95%) |
DBH | 1.06 | 1.05 – 1.09 | 1.03 – 1.17 |
Distance | 0.94 | 0.92 – 0.96 | 0.84 – 0.98 |
Random Effects | |||
σ2 | 0.07 | ||
Ï„00 | 0.05 | ||
ICC | 0.60 | ||
N Colony | 8 | ||
N Species | 224 | ||
N TAG | 724 | ||
Observations | 820 | ||
Marginal R2 / Conditional R2 | 0.096 / 0.526 |
Compute R2 for reduced models to estimate how big individual effects are
# colony
form <- formula(Foraged ~ (1 | Colony))
priors <- get_prior(form, data = speciesdat, family = bernoulli())
reduced <- brm(form, data = speciesdat, family = bernoulli, cov_ranef = list(phylo = A), cores = 4, control = list(adapt_delta = 0.9), prior = priors, iter = 4000)
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
r2[["colony"]] <- bayes_R2(reduced)
# species effect
form <- formula(Foraged ~ (1 | Species) )
priors <- get_prior(form, data = speciesdat, family = bernoulli())
## Warning: Rows containing NAs were excluded from the model.
reduced <- brm(form, data = speciesdat, family = bernoulli, cov_ranef = list(phylo = A), cores = 4, control = list(adapt_delta = 0.9), prior = priors, iter = 4000)
## Warning: Rows containing NAs were excluded from the model.
r2[["species"]] <- bayes_R2(reduced)
# tag effect
form <- formula(Foraged ~ (1 | TAG) )
priors <- get_prior(form, data = speciesdat, family = bernoulli())
reduced <- brm(form, data = speciesdat, family = bernoulli, cov_ranef = list(phylo = A), cores = 4, control = list(adapt_delta = 0.9), prior = priors, iter = 4000)
## Warning: There were 46 divergent transitions after warmup. Increasing adapt_delta above 0.9 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
r2[["tag"]] <- bayes_R2(reduced)
# DBH effect
form <- formula(Foraged ~ DBH)
priors <- get_prior(form, data = speciesdat, family = bernoulli())
## Warning: Rows containing NAs were excluded from the model.
reduced <- brm(form, data = speciesdat, family = bernoulli, cov_ranef = list(phylo = A), cores = 4, control = list(adapt_delta = 0.9), prior = priors, iter = 4000)
## Warning: Rows containing NAs were excluded from the model.
r2[["dbh"]] <- bayes_R2(reduced)
# distance effect
form <- formula(Foraged ~ Distance)
priors <- get_prior(form, data = speciesdat, family = bernoulli())
reduced <- brm(form, data = speciesdat, family = bernoulli, cov_ranef = list(phylo = A), cores = 4, control = list(adapt_delta = 0.9), prior = priors, iter = 4000)
r2[["distance"]] <- bayes_R2(reduced)
r2
## $complete
## Estimate Est.Error Q2.5 Q97.5
## R2 0.5085471 0.1432578 0.2190053 0.7706338
##
## $colony
## Estimate Est.Error Q2.5 Q97.5
## R2 0.007515573 0.006968773 2.787144e-05 0.02507872
##
## $species
## Estimate Est.Error Q2.5 Q97.5
## R2 0.07357428 0.03192131 0.01580775 0.139037
##
## $tag
## Estimate Est.Error Q2.5 Q97.5
## R2 0.2170842 0.1192893 0.008944221 0.4542796
##
## $dbh
## Estimate Est.Error Q2.5 Q97.5
## R2 0.07468829 0.02105632 0.03746641 0.1186829
##
## $distance
## Estimate Est.Error Q2.5 Q97.5
## R2 0.01014808 0.006682576 0.0008406406 0.02690465
We are are also testing the hypothesis that there is no phylogenetic effect (ignoring all other factors), following Nakagawa and Schielzeth’s (2012) Table 2 for distribution-specific variance.
phylo_only <- brm(Foraged ~ (1 | Species), data = speciesdat, family = bernoulli, cov_ranef = list(phylo = A), cores = 4, control = list(adapt_delta = 0.9))
## Warning: Rows containing NAs were excluded from the model.
## Compiling the C++ model
## Start sampling
bayes_R2(phylo_only)
## Estimate Est.Error Q2.5 Q97.5
## R2 0.07320662 0.03105433 0.01903532 0.1368142
hyp <- paste("sd_Species__Intercept /", "(sd_Species__Intercept + 3.141593^(2/3)) = 0")
(hyp <- hypothesis(phylo_only, hyp, class = NULL))
## Hypothesis Tests for class :
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (sd_Species__Inte... = 0 0.28 0.06 0.16 0.38 NA
## Post.Prob Star
## 1 NA *
## ---
## '*': The expected value under the hypothesis lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
plot(hyp)
We find that the complete model has the best fit to the data
UsingWe can look a the coefficients of individual species, which correspond to how likely an they are to be foraged, aftera accounting for all the otehr factors, like the phylogenetic effect.
palatcoef <- complete %>% spread_draws(b_Intercept, r_Species[Species ]) %>% group_by(Species) %>%
median_qi(Species_mean = r_Species) %>% mutate(Species = gsub("_"," ", Species))
## Warning: Expected 1 pieces. Additional pieces discarded in 1792000 rows [1,
## 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...].
ggplot(palatcoef, aes(Species_mean)) +geom_histogram() + theme_minimal() + xlab("Palatability coefficients") + ylab("Number of tree species")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#ggsave("figures/palatbility.pdf", height = 5, width = 5)
plotdat <- as.data.frame(speciesdat %>% group_by(Species) %>% summarize(Foraged = as.factor(max(Foraged))) %>% ungroup() %>% mutate(Species = gsub("_", " ", Species)) %>% select(Species, Foraged) %>% na.omit() %>% left_join(palatcoef %>% select(Species, Species_mean)))
## Joining, by = "Species"
row.names(plotdat) <- plotdat$Species
t2 <- tre
t2$tip.label <- gsub("_", " ", t2$tip.label)
p1 <- ggtree(t2, layout = 'fan') %<+% plotdat + geom_tiplab2(aes(label=paste0('italic(\'', label, '\')'), color=Foraged), size = 2, offset=0.5, parse = TRUE) + scale_color_manual(values = c("black", "red"))
p2 <- gheatmap(p1, plotdat[,"Species_mean", drop=FALSE], offset = -1, width=0.05, colnames=FALSE, colnames_position = "bottom") + scale_fill_gradient2(low = "black", mid = "white", high = "red", midpoint = 0) + theme(legend.position = "bottom")
## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.
print(p2)
#ggsave("figures/tree.pdf", height = 8, width = 8)