04-EcoEvoAnalyses.Rmd

```{r analyses, include=FALSE, eval=T}
rm(list = ls()) ; invisible(gc()) ; set.seed(42)
library(knitr)
library(tidyverse)
library(bayesplot)
library(ggtree)
library(phylosignal)
library(adephylo)
library(ape)
library(phylobase)
theme_set(bayesplot::theme_default())
opts_chunk$set(
  echo = F, message = F, warning = F, fig.height = 6, fig.width = 8,
  cache = T, cache.lazy = F, eval=T)
```

```{r env, eval=F}
rm(list = ls()) ; invisible(gc()) ; set.seed(42)
library(parallel)
library(doSNOW)
library(foreach)
library(tidyverse)
data <- vroom::vroom("save/full/mdatafull.tsv") %>% 
  dplyr::select(idTree, Plot, Xutm, Yutm) %>% 
  unique()
cl <- makeCluster(4, outfile = "/dev/null")
registerDoSNOW(cl)
pb <- txtProgressBar(max = nrow(data), style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress = progress)
NC <- foreach(ind=1:nrow(data),
              .packages = "dplyr",
              .options.snow = opts) %dopar% {
                con <- DBI::dbConnect(RSQLite::SQLite(), dbname = "data/guyafor.sql")
                NC <- tbl(con, "inventory") %>% 
                  filter(Plot == local(data$Plot[ind])) %>% 
                  filter(idTree != local(data$idTree[ind])) %>% 
                  mutate(dij = sqrt((local(data$Xutm[ind]) - Xutm)^2+(local(data$Yutm[ind]) - Yutm)^2)) %>% 
                  filter(dij < 20) %>% 
                  mutate(DBH = CircCorr/pi) %>% 
                  collect() %>% 
                  group_by(CensusYear) %>% 
                  summarise(NCI = sum(DBH*DBH*exp(-0.25*dij))) %>% 
                  ungroup() %>% 
                  summarise(idTree = local(data$idTree[ind]),
                            NCI = mean(NCI))
                DBI::dbDisconnect(con)
                return(NC)
              } 
close(pb)
stopCluster(cl)
NC <- bind_rows(NC)
XY <- NC %>% 
  left_join(data) %>% 
  sf::st_as_sf(coords = c("Xutm", "Yutm"),
             crs = '+proj=utm +zone=22 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0')
NC$twi <- raster::extract(raster::raster("data/TWI_1m.tif"), XY)
write_tsv(NC, file = "save/ecoevo/env_full.tsv")
```

```{r phylo, eval=F}
library(V.PhyloMaker)
splist <- vroom::vroom("save/full/mdatafull.tsv") %>% 
  dplyr::select(Family, Genus, species) %>% 
  unique() %>% 
  mutate(genus = Genus, family = Family) %>% 
  dplyr::select(species, genus, family)
tree <- phylo.maker(sp.list = splist, tree = GBOTB.extended, nodes = nodes.info.1, scenarios = "S3")
ape::write.tree(tree$scenario.3, "save/ecoevo/phylogeny_full.tree")
```

```{r parsan, eval=F}
load("save/growth.Rdata")  
pars <- as.data.frame(fit_full, pars = c('gmax')) %>% 
  reshape2::melt(variable.name = "parameter") %>% 
  group_by(parameter) %>% 
  summarise(value = median(value)) %>% 
  separate(parameter, c("parameter", "ind"), convert = T) %>% 
  reshape2::dcast(ind ~ parameter)
rm(fit_full) ; invisible(gc())
mdata <- vroom::vroom("save/full/mdatafull.tsv")
growth <- pars %>% 
  left_join(dplyr::select(mdata, ind, idTree, Family, Genus, species)) %>% 
  left_join(vroom::vroom(file = "save/ecoevo/env_full.tsv")) %>% 
  unique()
vroom::vroom_write(growth, file = "save/ecoevo/growthfull.tsv") 
```

```{r ecoevo, eval=F}
rm(list = ls()) ; invisible(gc())
# pars
niter_in <- 500
niter_out <- 1000
cores <- 10

# libs
library(brms)
library(tidyverse)

# data
load("save/full/growth.Rdata")  
mdata <- vroom::vroom("save/full/mdatafull.tsv")
env <- vroom::vroom(file = "save/ecoevo/env_full.tsv")

# code
iters <- sample(1:4000, niter_in)
fits <- list()
for(i in 1:length(iters)){
  iter <- iters[i]
  data <- as.data.frame(fit_full, pars = c('gmax'))[iter,] %>% 
    reshape2::melt(variable.name = "parameter") %>% 
    separate(parameter, c("parameter", "ind"), convert = T) %>% 
    reshape2::dcast(ind ~ parameter) %>% 
    left_join(dplyr::select(mdata, ind, idTree, Family, Genus, species)) %>% 
    left_join(env) %>% 
    unique() %>% 
    dplyr::select(Family, Genus, species, NCI, twi, gmax) %>% 
    na.omit() %>% 
    mutate(twi = ifelse(twi <= 0, 10^-6, twi))
  fit <- brm(log(gmax) ~ log(NCI) + log(twi+1) + (1  | Family/Genus/species), 
             data, chain = 1, 
             threads = cores, 
             iter = niter_out, 
             warmup = niter_out - 1,
             backend = "cmdstanr")
  fits[[i]] <- fit
}
allfits <- combine_models(mlist = fits)

# save
save(allfits, file = "save/ecoevo/ecoevo.Rdata")
```

```{r dataan}
growth <- vroom::vroom("save/ecoevo/growthfull.tsv")
phylo <- ape::read.tree("save/ecoevo/phylogeny_full.tree")
p4d <- phylo4d(phylo,
               data.frame(species = gsub("_", " ", phylo$tip.label)) %>% 
                 left_join(group_by(growth, species) %>% 
                             summarise(gmax  = median(gmax), loggmax = log(gmax))
                           ) %>% 
                 dplyr::select(-species))   
```

# Eco-evolutionary analyses

In this chapter, I investigated effects of phylogeny and ecological processes on individual growth, using phylogeny, topography and neighbourhood indices.

## Methods

### Environmental descriptors

I used the mean neighbourhood crowding index [$NCI$; @Uriarte2004] over the last 30 years,
an indirect measurement of access to light and forest gap dynamics.
The mean neighbourhood crowding index $NCI_i$ from tree individual $i$ was calculated as follows:

$$NCI_i=\overline{\sum_{j|\delta_{i,j}<20m}DBH^2_{j,t}.e^{-\frac14\delta_{i,j}}}$$

with $DBH_{j,t}$ the diameter of the neighbouring tree $j$ in year $t$ and $\delta_{i,j}$ its distance to the individual tree $i$.
$NCI_i$ is computed for all neighbours at a distance $\delta_{i,j}$ inferior to the maximum neighbouring distance of 20 meters.
The power of neighbours $DBH_{j,t}$ effect was set to 2 to represent a surface.
The decrease of neighbours' diameter effect with distance was set to -0.25 to represent trees at 20 meters of the focal trees having 1% of the effect of the same tree at 0 meters.
$NCI_i$ is computed as the mean of yearly $NCI_{i,t}$ over the last 30 years denoted by the overline.

I used the topographic wetness index ($TWI$) as proxies of the distribution of soil water and nutrients in Paracou.
Waterlogging and topography have been highlighted as crucial for forest dynamics [@Ferry2010],
species-habitat relationships [@Engelbrecht2007], and phenotypic variation [@Schmitt2020].
Topography, through the dissolution of iron oxides, litter- and tree-fall transfers and waterlogging,
shapes soil nutrient distribution in tropical forests [@John2007; @Ferry2010].
TWI was derived from a 1-m-resolution digital elevation model using SAGA-GIS [@Conrad2015]
based on a LiDAR campaign of the whole Paracou field station done in 2015.

### Analyses

To study the effect of phylogeny and environment,
we investigated the effects of family, genus, species & topography $TWI$ and neighbourhood $NCI$ indices
on individual growth potential $Gmax$ with the following linear mixed model:

$$log(Gmax) ~ \sim \mathcal N (\alpha_{species} + \beta_{NCI} \times log(NCI),\sigma) \\ |~  \alpha_{species} \sim \mathcal N(\alpha_{genus},\sigma_{species})\\ |~  \alpha_{genus} \sim \mathcal N(\alpha_{family},\sigma_{genus})\\ |~  \alpha_{family} \sim \mathcal N(\alpha,\sigma_{family})$$

We reported the corresponding marginal and conditional $R^2$ [@Nakagawa2013].
We further plotted individual growth potential $Gmax$ along environmental gradients and across the phylogeny.
We used `phylosingal` to test for phylogenetic signal [@keck_phylosignal:_2016],
before computing phylogenetic correlogram and local indicator of phylogenetic association.

We used R version 3.6 for all statistical analyses (<www.r-project.org>). 

## Results

### General

Most of the variation in growth potential is individual, then explained by genus before species and family. 
The taxonomic structure explains almost a third of the observed variation in individual growth potential.
Finally, the neighbourhood ($NCI$) has a marked negative significant effect (Fig. \@ref(fig:gmaxncib)) which explains 10% of the observed variation.

```{r}
load("save/ecoevo/ecoevo.Rdata")
```

```{r, eval=FALSE}
growth %>% 
  group_by(Family, species) %>% 
  summarise(`Median Gmax` = median(gmax), `CV[log] Gmax` = sqrt(exp(sd(log(gmax))^2)-1), 
            `Median NCI` = median(NCI), `Median TWI` = median(twi, na.rm = TRUE),
            `N trees`= n()) %>% 
  left_join(vroom::vroom("save/full/mdatafull.tsv") %>% 
              mutate(species = paste(Genus, Species)) %>% 
              group_by(species) %>% 
              summarise(`Median DBH` = median(DBH), `DBH 0.95` = quantile(DBH, 0.95))) %>% 
  write_tsv("~/Téléchargements/growth.tsv")
```


```{r lmm}
allfits %>%
  sjPlot::tab_model(show.icc = F)
```

```{r, eval=F}
library(tidyverse)
load("save/ecoevo/ecoevo.Rdata")
res <- residuals(allfits) %>% 
  as_tibble() %>% 
  mutate(ind = 1:n()) %>% 
  left_join(vroom::vroom("save/full/mdatafull.tsv") %>% 
              group_by(ind) %>% 
              summarise(Xutm = mean(Xutm), Yutm = mean(Yutm))) %>% 
  group_by(Xutm, Yutm) %>% 
  sample_n(1)
cor <- ncf::correlog(res$Xutm, res$Yutm, res$Estimate,
                     resamp = 0, increment = 100,
                     latlon = TRUE)
data.frame(N = cor$n, distance = cor$mean.of.class, I = cor$correlation) %>% 
  ggplot(aes(distance, I)) +
  geom_point(aes(size = N)) +
  theme_bw() +
  geom_smooth(method = "lm") +
  ylim(-0.02, 0.02) +
  xlab("Distance (m)") + ylab("Moran's I") +
  ggtitle("Model 2 residuals", "log(Gmax) ~ log(NCI) + log(TWI+1) + (1  | Family/Genus/Species)")
```

```{r}
growth %>% 
  left_join(vroom::vroom("save/full/mdatafull.tsv") %>% 
              dplyr::select(idTree, Plot, SubPlot) %>% 
              unique()) %>% 
  mutate(Quadrat = paste0(Plot, "-", SubPlot)) %>% 
  lme4::lmer(log(gmax) ~ log(NCI) + log(twi + 1) + (1 | Family/Genus/species) + (1 | Quadrat), data = .) %>% 
  sjPlot::tab_model()
```


```{r}
# brms::parnames(allfits)
bayesplot::mcmc_dens(allfits, pars = c("b_logNCI", "b_logtwiP1"))
```

```{r}
(g.var <-data.frame(level = c("Family", "Genus", "Species", "Individual"), 
           variance = c(0.03, 0.12, 0.04, 0.72)) %>% 
  mutate(pct = variance / sum(variance)*100) %>% 
  mutate(level = factor(level, levels = c("Family", "Genus", "Species", "Individual"))) %>% 
  mutate(pct_text = paste0(round(pct), "%")) %>% 
  mutate(pct_text = gsub("^0%", "", pct_text)) %>% 
  ggplot(aes(x = 1, y = variance, fill = level)) +
  geom_bar(stat = "identity", position = "stack") +
  geom_text(aes(y = variance, label = pct_text), col = "white", 
            position = position_stack(vjust = 0.5), size = 3) +
  scale_fill_manual(expression(sigma^2), 
                    values = ggpubfigs::friendly_pal(4, name = "nickel_five")) +
  theme(axis.title = element_blank(), axis.line = element_blank(), axis.text = element_blank(), 
        axis.ticks = element_blank(), legend.position = "bottom")+
  coord_flip() +
  scale_y_reverse() +
  ggtitle("B."))   
```

```{r gmaxncib, fig.cap="Relation between Neighbourhood crowding index (NCI) and individual growth potential (Gmax, cm/yr)."}
(g.nci <-ggplot(growth, aes(NCI, gmax, col = species, group = NA)) +
  geom_bin2d(bins = 100) +
  scale_fill_continuous(type = "viridis") +
  geom_smooth(method = MASS::rlm, aes(group = NA), se = T, col = "red") +
  viridis::scale_fill_viridis(guide = "none") +
  scale_x_log10() + scale_y_log10() +
  xlab("Neighbourhood crowding index (NCI)") +
  ylab(expression(Gmax[i]~cm.year^-1)) +
  ggtitle("A.", expression(beta==-0.51~R^2~Marginal==0.019)))
```

### Phylogenetic

Species growth potential in the phylogeny (Fig. \@ref(fig:gmaxphylob)) is significantly structured (Tab. \@ref(tab:gmaxphylosignalb)),
with a short distance significant positive association and a long distant significant negative association (Fig. \@ref(fig:gmaxcorrelogb)).
Local indicator of phylogenetic association highlight the conservation of species growth potential at the genus level (Fig. \@ref(fig:lipab)),
as illustrated for instance with fast growing *Cecropia* opposed to slow growing *Eschweilera*.
However, a few species have different growth potential in the same genus, 
such as slow growing *Drypetes variabilis* opposed to fast growing *Drypetes fanshawei*.


```{r gmaxphylo, fig.cap="Distribution of species growth potential (Gmax, cm/yr) in the phylogeny."}
(g.phylo <- fortify(phylo) %>% 
  mutate(species = gsub("_", " ", label)) %>%
  mutate(label = species) %>% 
  left_join(group_by(growth, species) %>% summarise(gmax = median(gmax))) %>% 
  ggtree(aes(color = gmax), layout="circular") + 
  geom_tiplab2(size = 2) +
  theme_tree(legend.position='right', legend.text = element_text(face = "italic")) +
  scale_alpha_manual("taxon", values = c(0.2, 1)) +
  scale_size_manual("taxon", values = c(1, 2)) +
  viridis::scale_color_viridis(expression(Median(Gmax[i])~cm.year^-1), trans = "log", labels=scales::comma) +
  theme(legend.position = "bottom") + xlim(NA, 200) +
  ggtitle("C.", expression(Pagel~lambda==0.78~p<0.001)))
```

```{r gmaxphylosignalb}
phyloSignal(p4d = p4d, method = "all") %>% 
  lapply(as.data.frame) %>% 
  lapply(rownames_to_column, "parameter") %>% 
  bind_rows(.id = "type") %>% 
  reshape2::melt(c("type", "parameter")) %>% 
  reshape2::dcast(parameter + variable ~ type) %>% 
  mutate(value = paste0(round(stat, 4), " (p=", round(pvalue, 4), ")")) %>% 
  reshape2::dcast(parameter ~ variable) %>% 
  kable(caption = "Phylogenetic signal of species growth potential (Gmax, cm/yr) with different methods.")
```

```{r gmaxcorrelog, fig.cap="phylogenetic correlogram of species growth potential (Gmax, cm/yr)."}
# crlg <- phyloCorrelogram(p4d, trait = "loggmax", ci.bs	= 100)
# png(file = "save/figs/gmaxcorrelog_full.png", width = 1000, height = 1000, units = "px")
# plot(crlg, main="log(gmax) phylogenetic correlogram")
# dev.off()
include_graphics("save/figs/gmaxcorrelog_full.png") 
```

```{r lipa, fig.cap="Local phylogenetic associations of of species growth potential (Gmax, cm/yr) among taxa in the phylogeny."}
fortify(phylo) %>% 
  mutate(species = gsub("_", " ", label)) %>%
  mutate(label = species) %>% 
  left_join(lipaMoran(p4d)$p.value %>% 
              as.data.frame() %>% 
              rownames_to_column("species") %>% 
              mutate(species = gsub("_", " ", species)) %>% 
              dplyr::rename(pval = loggmax)) %>% 
  ggtree(aes(color = pval < 0.05), layout="circular") + 
  geom_tiplab2(size = 2) +
  theme_tree(legend.position='right', legend.text = element_text(face = "italic")) +
  scale_alpha_manual("taxon", values = c(0.2, 1)) +
  scale_size_manual("taxon", values = c(1, 2)) +
  scale_color_manual("Local\nIndicator\nof\nPhylogenetic\nAssociation\np<0.05", 
                     values =  c("grey", "red", "grey")) +
  xlim(NA, 200) 
```

```{r figall, fig.width=10, fig.height=7}
gridExtra::grid.arrange(
  grobs = list(g.phylo,  g.nci, g.var),
  layout_matrix = rbind(c(2, 1),
                        c(3, 1)),
  heights = c(3, 1), widths = c(1.5, 2))
# ggsave(plot = g, filename = "figs/fig2.png", dpi = 300, bg = "white", width = 10, height = 7)
```