library(tidyverse)
library(ape)
library(tidybayes)
library(brms)
library(ggtree)
library(sjPlot)

Mixed model

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

Reduced models

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

phylogenetic effect

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

Palatability coefficients

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)