06-Guyafor.Rmd

```{r guyafor, include=FALSE, eval=T}
rm(list = ls()) ; invisible(gc()) ; set.seed(42)
library(knitr)
library(tidyverse)
library(bayesplot)
library(cmdstanr)
library(ggtree)
library(phylosignal)
library(adephylo)
library(ape)
library(phylobase)
theme_set(bayesplot::theme_default())
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)
```

# Guyafor

In this chapter, I repeated the model fit for the whole Guyafor network.

## Data

I used only recruited trees in the censuses with at least 4 measurements of diameter at breast height (DBH, cm).
I used only species with at least 4 trees following previous requirements (Tab. \@ref(tab:guyafortab) & Fig. \@ref(fig:mdataguyaforfig)).

```{r guyaforlim}
n_ind_species <- 6
n_years <- 10
n_census <- 4
```

```{r guyaforfull, eval=F}
trees <- vroom::vroom("data/Selection_Guyafor_Gmax_4_4.csv", locale=locale(decimal_mark=',')) %>% 
  mutate(species = paste(Genus, Species)) %>% 
  filter(!grepl("Indet", species)) %>% 
  filter(BotaSource == "Bota") %>% 
  group_by(idTree) %>% 
  arrange(CensusYear) %>% 
  mutate(FirstDead = first(CensusYear[CodeAlive == 0])) %>% 
  mutate(FirstDead = ifelse(is.na(FirstDead), max(CensusYear)+1, FirstDead)) %>% 
  filter(CensusYear < FirstDead) %>% 
  ungroup() %>% 
  mutate(DBH = CircCorr/pi) %>%
  group_by(Plot) %>% 
  mutate(StartYear = min(CensusYear)) %>% 
  group_by(idTree) %>% 
  arrange(CensusYear) %>% 
  filter(first(CensusYear) > StartYear) %>% 
  filter(first(DBH) < 15) %>%
  filter(last(DBH) > first(DBH)) %>% 
  group_by(idTree) %>% 
  filter(max(CensusYear) - min(CensusYear) >= n_years) %>% 
  filter(n() >= n_census) %>% 
  group_by(species) %>% 
  filter(length(unique(idTree)) >= n_ind_species)
rm(n_census, n_ind_species)
vroom::vroom_write(trees, "save/guyafor/trees_guyafor.tsv")
```

```{r  guyafortab}
options(knitr.kable.NA = '')
vroom::vroom("save/guyafor/trees_guyafor.tsv") %>% 
  group_by(Family, Genus, species, idTree) %>% 
  summarise(census = n(), year_start = min(CensusYear), year_end = max(CensusYear),
            dbh_start = min(DBH), dbh_end = max(DBH)) %>% 
  ungroup() %>% 
  summarise(n_families = length(unique(Family)),
            n_genera = length(unique(Genus)),
            n_species = length(unique(species)),
            n_individuals = length(unique(idTree)),
            n_observations = nrow(vroom::vroom("save/guyafor/trees_guyafor.tsv")),
            min_census = min(census),
            med_census = median(census),
            max_census = max(census),
            min_year0 = min(year_start),
            med_year0 = median(year_start),
            max_year0 = max(year_start),
            min_yearmax = min(year_end),
            med_yearmax = median(year_end),
            max_yearmax = max(year_end),
            min_dbh0 = min(dbh_start),
            med_dbh0 = median(dbh_start),
            max_dbh0 = max(dbh_start),
            min_dbhmax = min(dbh_end),
            med_dbhmax = median(dbh_end),
            max_dbhmax = max(dbh_end)) %>% 
  reshape2::melt() %>% 
  separate(variable, c("measure", "variable")) %>% 
  reshape2::dcast(variable ~ measure) %>% 
  dplyr::select(variable, n, med, min, max) %>% 
  mutate(variable = factor(variable, levels = c("families", "genera", "species", "individuals", "observations",
                                                "census", "year0", "yearmax", "dbh0", "dbhmax"))) %>% 
  arrange(variable) %>% 
  kable(col.names = c("", "N", "Median", "Minimum", "Maximum"), format.args = list(big.mark = " "), digits = 0,
        caption = "Metrics on inventory data used to fit the full model including sample size (N), memdian, minimum and maximum values for families, genera, species, individuals, observations, cenusus, recruitment year (year0), last censused year (yearmax), recruitment diameter (dbh0) and last censused diameter (dbhmax).")    
```


```{r mdataguyafor, eval=F}
vroom::vroom("save/guyafor/trees_guyafor.tsv") %>% 
  group_by(idTree) %>% 
  mutate(Year = CensusYear - min(CensusYear)) %>%
  ungroup() %>% 
  mutate(ind = as.numeric(as.factor(as.character(idTree)))) %>% 
  mutate(sp = as.numeric(as.factor(species))) %>% 
  vroom::vroom_write("save/guyafor/mdataguyafor.tsv")
```

```{r mdataguyaforfig, fig.cap="Tree diameter trajectories in reduced data. Color represent individuals."}
dplyr::select(vroom::vroom("save/guyafor/mdataguyafor.tsv"), sp, ind) %>% 
  unique() %>% 
  filter(sp %in% sample(unique(.$sp), 9)) %>% 
  group_by(sp) %>% 
  sample_n(5, replace = T) %>%
  unique() %>% 
  left_join(vroom::vroom("save/guyafor/mdataguyafor.tsv")) %>% 
  ggplot(aes(Year, DBH, group = as.factor(ind))) +
  geom_point(col = "grey") +
  geom_smooth(se = F, aes(col = as.factor(ind))) +
  xlim(0,NA) +
  facet_wrap(~ species, scales = "free") +
  scale_color_discrete(guide = "none")  
```

## Model

I used a Gompertz model [@Herault2011], were the diameter of individual $i$ at year $t$ is the sum of annual growth from $t0$ to $t$:

$$  DBH_{t,i,s} \sim \mathcal N  (10 + Gmax_i \times \sum _{y=1|DBH_{t=0}} ^{y=t} exp(-\frac12.[\frac{log(\frac{DBH_{t,i}}{100.Dopt_i})}{Ks_i}]^2)), \sigma) \\| Dopt_i \sim \mathcal N(Dopt_s,\sigma_D), Ks_i \sim \mathcal N(Ks_s,\sigma_K) $$

The annual growth rate for individual $i$ at year $y$ with a diameter of $DBH_{y,i}$ is defined following a Gompertz model [@Gompertz1825] already identified as the best model for growth-trajectories in Paracou [@Herault2011],
where $Gmax_i$ is the fixed maximum growth potential of every individual,
$Dopt_i$ is the optimal diameter at which the individual reaches its maximum growth potential,
and $Ks_i$ is the kurtosis defining the width of the bell-shaped growth-trajectory [see figure 1 in @Herault2011].
$Dopt_i$ and $Ks_i$ are random effects centered on species parameters $Dopt_s$ and $Ks_s$ with associated variances $\sigma_D$ and $\sigma_K$.

```{r growthguyaforfit, eval=F}
mdata <- vroom::vroom("save/mdataguyafor.tsv")
growth <- cmdstan_model("model/growth.stan")
fit <- growth$sample(
  data = list(
    N = nrow(filter(mdata, Year > 0)),
    I = max(mdata$ind),
    S = max(mdata$sp),
    Y = max(mdata$Year),
    year = filter(mdata, Year > 0)$Year,
    dbh = filter(mdata, Year > 0)$DBH,
    dbh0 = filter(mdata, Year == 0)$DBH,
    dmax = summarise(group_by(mdata, sp), dmax = max(DBH))$dmax,
    ind = filter(mdata, Year > 0)$ind,
    indsp = arrange(unique(mdata[c("ind", "sp")]), ind)$sp
  ),
  chains = 4, 
  parallel_chains = 4,
  refresh = 10,
  save_warmup = F,
  max_treedepth = 12
)
fit$save_output_files(dir = "save/growthguyafor")
```


```{r, eval=F, echo=T}
fit <- read_cmdstan_csv("save/growthguyafor/growth-202210280057-1-361cf6.csv")
draws <- drop(fit$post_warmup_draws) %>% 
  as_data_frame() %>% 
  mutate(iteration = 1:n()) %>% 
  gather(parameter, value, -iteration)
vroom::vroom_write(draws, file = 'save/growthguyafor/chain1.tsv')
rm(fit, draws)
gc()
```

```{r, eval=F, echo=T}
fit <- read_cmdstan_csv("save/growthguyafor/growth-202210280057-2-361cf6.csv")
draws <- drop(fit$post_warmup_draws) %>% 
  as_data_frame() %>% 
  mutate(iteration = 1:n()) %>% 
  gather(parameter, value, -iteration)
vroom::vroom_write(draws, file = 'save/growthguyafor/chain2.tsv')
rm(fit, draws)
gc()
```

```{r, eval=F, echo=T}
fit <- read_cmdstan_csv("save/growthguyafor/growth-202210280057-3-361cf6.csv")
draws <- drop(fit$post_warmup_draws) %>% 
  as_data_frame() %>% 
  mutate(iteration = 1:n()) %>% 
  gather(parameter, value, -iteration)
vroom::vroom_write(draws, file = 'save/growthguyafor/chain3.tsv')
rm(fit, draws)
gc()
```

```{r, eval=F, echo=T}
fit <- read_cmdstan_csv("save/growthguyafor/growth-202210280057-4-361cf6.csv")
draws <- drop(fit$post_warmup_draws) %>% 
  as_data_frame() %>% 
  mutate(iteration = 1:n()) %>% 
  gather(parameter, value, -iteration)
vroom::vroom_write(draws, file = 'save/growthguyafor/chain4.tsv')
rm(fit, draws)
gc()
```

```{bash, eval=F, echo=T}
grep gmax chain1.tsv > gmax1.tsv
grep gmax chain2.tsv > gmax2.tsv
grep gmax chain3.tsv > gmax3.tsv
grep gmax chain4.tsv > gmax4.tsv
```

```{r, eval=F, echo=T}
bind_rows(
  vroom::vroom('save/growthguyafor/gmax1.tsv', col_names = c("iteration", "parameter", "value")) %>%
  mutate(chain = 1),
  vroom::vroom('save/growthguyafor/gmax2.tsv', col_names = c("iteration", "parameter", "value")) %>%
  mutate(chain = 2),
  vroom::vroom('save/growthguyafor/gmax3.tsv', col_names = c("iteration", "parameter", "value")) %>%
  mutate(chain = 3),
  vroom::vroom('save/growthguyafor/gmax4.tsv', col_names = c("iteration", "parameter", "value")) %>%
  mutate(chain = 4)
) %>% 
  separate(parameter, c("gmax", "ind")) %>% 
  vroom::vroom_write(file = 'save/growthguyafor/gmax.tsv')
```

```{r, eval=F, echo=T}
library(csv2sql)
gmax <- vroom::vroom(file = 'save/growthguyafor/gmax.tsv')
write_csv(gmax, file = 'save/growthguyafor/gmax.csv')
csv_to_sqlite(csv_name = 'save/growthguyafor/gmax.csv', 
              db_name = 'save/growthguyafor/gmax.sql', 
              table_name = "gmax")
unlink('save/growthguyafor/gmax.csv')
```

```{r draws, eval=F}
rm(list = ls()) ; gc()
library(foreach)
I <- max(vroom::vroom("save/mdataguyafor.tsv")$ind)
cl <- parallel::makeCluster(10, outfile = "")
doSNOW::registerDoSNOW(cl)
pb <- utils::txtProgressBar(max = I, style = 3)
progress <- function(n) utils::setTxtProgressBar(pb, n)
opts <- list(progress = progress)
gmax <- foreach(i=1:I, .options.snow = opts) %dopar% {
  suppressMessages(library(dplyr))
  gmaxdb <- DBI::dbConnect(RSQLite::SQLite(), 
                           dbname = "save/growthguyafor/gmax.sql")
  gmax <- tbl(gmaxdb, "gmax") %>% 
    filter(ind == i) %>% 
    collect() %>% 
    summarise(q5 = quantile(value, 0.05), q25 = quantile(value, 0.25), 
              median= quantile(value, 0.5), 
              q75 = quantile(value, 0.75), q95 = quantile(value, 0.95))
  DBI::dbDisconnect(gmaxdb) ; rm(gmaxdb)
  gmax
}
parallel::stopCluster(cl) ; rm(cl)
save(gmax, file = "gmax.Rdata")
names(gmax) <- 1:I 
gmax_df <- bind_rows(gmax, .id = "ind")
write_tsv(gmax_df, file = "save/growthguyafor/gmax_sum.tsv")
```

## Fit

**Unfortunately, gmax upper limit was set to 5 instead of 10 (overconstrained).**
**But the model with more stringent data had a better behavior with gmax limited to 5.**

```{r growthguyaforgmaxtrace}
gmax <- DBI::dbConnect(RSQLite::SQLite(), dbname = "save/guyafor/gmax.sql")
t <- tbl(gmax, "gmax") %>% 
  filter(ind %in% 1:4) %>% 
  collect()
DBI::dbDisconnect(gmax) ; rm(gmax)
t %>% 
  ggplot(aes(iteration, value, col = as.factor(chain))) +
  geom_line() +
  facet_wrap(~ ind, scales = "free_y") + 
  viridis::scale_color_viridis(discrete = T, guide = "none") +
  ylab(expression(g[max]))
```

```{r growthguyaforgmax}
gmax <- DBI::dbConnect(RSQLite::SQLite(), dbname = "save/guyafor/gmax.sql")
t <- tbl(gmax, "gmax") %>% 
  filter(ind %in% 1:100) %>% 
  collect() %>% 
  group_by(ind) %>% 
  summarise(ll = quantile(value, 0.05), l = quantile(value, 0.25), 
            m = quantile(value, 0.5), 
            h = quantile(value, 0.75), hh = quantile(value, 0.95))
DBI::dbDisconnect(gmax) ; rm(gmax)
ggplot(t, aes(x = ind, xend = ind)) +
  geom_point(aes(y = m), shape = 21, size = 3, alpha = 0.5) +
  geom_segment(aes(y = ll, yend = hh),
               size = 1, show.legend = F, alpha = 0.5) +
  geom_segment(aes(y = l, yend = h), size = 2, alpha = 0.5) +
  coord_flip() +
  xlab("First 100 individuals") +
  ylab("gmax") +
  scale_y_log10()
```

## Evolutionary analyses

```{r parsanc, eval=F}
draws <- vroom::vroom('save/guyafor/gmax_sum.tsv')
inds <- vroom::vroom("save/guyafor/mdataguyafor.tsv") %>% 
              group_by(Forest, Plot, SubPlot, TreeFieldNum, idTree, Xfield, Yfield, Xutm, Yutm, Lat, Lon, Family, Genus, Species, ind) %>% 
                summarise(FirstCensus = min(CensusYear), LastCensus = max(CensusYear), NCensuses = n(),
                          FirstDBH = min(DBH), LastDBH = max(DBH))
gmax <- left_join(inds, draws) %>% 
  select(-ind)
vroom::vroom_write(gmax, 'save/guyafor/growthguyafor.tsv') 
```

```{r phyloc, eval=F}
library(V.PhyloMaker)
splist <- vroom::vroom("save/mdataguyafor.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/guyafor/phylogeny_guyafor.tree")
```

```{r dataanc}
growth <- vroom::vroom("save/guyafor/growthguyafor.tsv") %>% 
  mutate(species = paste(Genus, Species))
phylo <- ape::read.tree("save/guyafor/phylogeny_guyafor.tree")
p4d <- phylo4d(phylo,
               data.frame(species = gsub("_", " ", phylo$tip.label)) %>% 
                 left_join(
                   growth %>% 
                     group_by(species) %>% 
                     summarise(gmax = median(median), loggmax = log(gmax))
                 ) %>% 
                 dplyr::select(-species))   
```

```{r gmaxphyloc, fig.cap="Distribution of species growth potential (Gmax, cm/yr) in the phylogeny."}
g <- fortify(phylo) %>% 
   mutate(species = gsub("_", " ", label)) %>%
   mutate(label = species) %>% 
   left_join(group_by(growth, species) %>% summarise(gmax = median(median))) %>% 
   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")
ggsave(g, filename = "save/figs/guyafor/phylogenyguyafor.png", dpi = 500, width = 25, height = 25, units = "cm", bg = 'white')
include_graphics("save/figs/guyafor/phylogenyguyafor.png")
```

## Discussion

* Within-chain parallelization available with [gompertzsumpar](https://github.com/sylvainschmitt/treegrowth/blob/main/model/gompertzsumpar.stan) model.
* Next steps:
    * Stop?
    * Use?
    * Run again? 
        * Cluster? 
        * Who? 
        * When?
   * What for? Data paper?