03-ModelFit.Rmd

```{r modelfull, include=FALSE, eval=T}
rm(list = ls()) ; invisible(gc()) ; set.seed(42)
library(knitr)
library(tidyverse)
library(bayesplot)
library(cmdstanr)
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)
```

# Model fit

In this chapter, I fitted the selected model.

## Data

I focused on trees at 20 meters from any plot edges for neighbourhood effect.
I used only recruited trees in the censuses with at least 10 measurements of diameter at breast height (DBH, cm).
I used only species with at least 10 trees following previous requirements (Tab. \@ref(tab:datafulltab) & Fig. \@ref(fig:mdatafullfig)).

```{r datafulllim}
dist_edge <- 20
n_census <- 10
n_ind_species <- 10
```

```{r datafull, eval=F}
guyafor <- DBI::dbConnect(RSQLite::SQLite(), dbname = "data/guyafor.sql")
trees <- tbl(guyafor, "inventory") %>% 
  filter(Forest == "Paracou") %>% 
  filter(Xfield > local(dist_edge), Xfield < 250-local(dist_edge), 
         Yfield > local(dist_edge), Yfield < 250-local(dist_edge)) %>% 
  mutate(species = paste(Genre, Espece)) %>% 
  collect() %>% 
  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(n() > n_census) %>% 
  group_by(species) %>% 
  filter(length(unique(idTree)) > n_ind_species)
DBI::dbDisconnect(guyafor) ; rm(guyafor, dist_edge, n_census, n_ind_species)
vroom::vroom_write(trees, "save/full/trees_full.tsv")
```

```{r datafulltab}
options(knitr.kable.NA = '')
vroom::vroom("save/full/trees_full.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/full/trees_full.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 mdatafull, eval=F}
vroom::vroom("save/full/trees_full.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/full/mdatafull.tsv")
```

```{r mdatafullfig, fig.cap="Tree diameter trajectories in reduced data. Color represent individuals."}
dplyr::select(vroom::vroom("save/full/mdatafull.tsv"), sp, ind) %>% 
  unique() %>% 
  filter(sp %in% sample(unique(.$sp), 9)) %>% 
  group_by(sp) %>% 
  sample_n(10, replace = T) %>%
  unique() %>% 
  left_join(vroom::vroom("save/full/mdatafull.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 growthfit, eval=F}
mdata <- vroom::vroom("save/full/mdatafull.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/full/growth")
fit_full <- rstan::read_stan_csv(fit$output_files())
save(fit_full, file = "save/full/growth.Rdata") 
```

## Fit

The model correctly converged ($\hat R < 1.1$) for the majority of $Gmax_i$.
All $\sigma$ have a small posterior but difficulties to converge.
The correlation between $Dopt$ and $Ks$ is acceptable but marked.
$Gmax_i$ posteriors have logical uncertainty but are varying widely among individuals.

```{r growthfullsampling, eval=F}
mdata <- vroom::vroom("save/full/mdatafull.tsv")
load("save/full/growth.Rdata")
```

```{r growthfullrhat}
# r <- rhat(fit_full)
# data.frame(parameter = names(r), rhat = r) %>% 
#   vroom::vroom_write("save/full/rhatfull.tsv")
vroom::vroom("save/full/rhatfull.tsv") %>% 
  separate(parameter, c("parameter", "indsp")) %>% 
  filter(parameter == "gmax") %>% 
  # arrange(desc(rhat)) %>% filter(rhat > 1.1)
  select(rhat) %>% 
  unlist() %>% 
  mcmc_rhat() +
  ggtitle("7,961 Gmax_i", "64 with rhat > 1.1 (0.8%)")
```

```{r growthfulltrace}
# g <- mcmc_trace(as.array(fit_full, pars = c("gmax[1]", "dopt[1]", "ks[1]", "sigma", "sigmaD", "sigmaK")))
# ggsave(plot = g, filename = "save/full/figs/growthfulltrace.png", bg = "white")
include_graphics("save/full/figs/growthfulltrace.png")
```

```{r growthfullpairs}
# g <- mcmc_pairs(as.array(fit_full, pars = c("gmax[1]", "dopt[1]", "ks[1]", "sigma", "sigmaD", "sigmaK")))
# ggsave(plot = g, filename = "save/full/figs/growthfullpairs.png", bg = "white")
include_graphics("save/full/figs/growthfullpairs.png")
```

```{r growthfullgmax}
# g <- mcmc_intervals(as.array(fit_full, pars = paste0("gmax[", 1:100, "]")))
# g <- g + xlab(expression(gmax[i])) +
#   theme(axis.text.y = element_text(size = 5))
# ggsave(plot = g, filename = "save/full/figs/growthfullgmax.png", bg = "white")
include_graphics("save/full/figs/growthfullgmax.png")
```

```{r growthfulldoptdmax, eval=F}
g <- mcmc_intervals_data(as.array(fit_full, pars = "dopt_s")) %>%
  separate(parameter, c("parameter", "X1", "sp"), convert = T) %>%
  left_join(mdata %>%
              group_by(sp) %>%
              summarise(dmax = max(DBH))) %>%
  ggplot(aes(sp, (m*100)/dmax)) +
  geom_point() +
  coord_flip() +
  scale_y_sqrt() +
  ylab(expression(frac(Dopt,Dmax)*100)) + xlab("Species")
ggsave(plot = g, filename = "growthfulldoptdmax.png", bg = "white")
```

```{r predsfull, eval=F}
preds <- as.data.frame(fit_full, pars = "mu")
preds <- preds %>% 
  reshape2::melt(variable.name = "parameter") %>% 
  group_by(parameter) %>% 
  summarise(ll = quantile(value, 0.05), 
            l = quantile(value, 0.25), 
            m = median(value), 
            h = quantile(value, 0.75),
            hh = quantile(value, 0.95))
mdata <- vroom::vroom("save/full/mdatafull.tsv")
allpreds <- filter(mdata, Year > 0) %>% 
  bind_cols(preds)
vroom::vroom_write(allpreds, "save/full/predsfull.tsv")
```

```{r predsfullfig}
preds <- vroom::vroom("save/full/predsfull.tsv") %>% 
  filter(m > 10)
rmse <- sqrt(mean((preds$m - preds$DBH)^2))
g <- preds %>% 
  ggplot(aes(DBH, m)) +
  geom_point(alpha = 0.25) +
  geom_errorbar(aes(ymin = ll, ymax = hh), alpha = 0.25) +
  geom_abline(col = "red", linetype = "dashed") +
  scale_y_log10() + scale_x_log10() +
  xlab("Observed diameter (Diameter at breast height, cm)") +
  ylab("Predicted diameter (Diameter at breast height, cm)") +
  ggtitle(paste("RMSE =", round(rmse, 2), "cm"))
ggsave(plot = g, filename = "figs/fig1.png", dpi = 300, bg = "white")
g
```

```{r predfullfig}
dplyr::select(vroom::vroom("save/full/predsfull.tsv"), sp, ind) %>% 
  unique() %>% 
  filter(sp %in% sample(unique(.$sp), 9)) %>% 
  group_by(sp) %>% 
  sample_n(10, replace = T) %>%
  unique() %>% 
  left_join(vroom::vroom("save/full/predsfull.tsv")) %>% 
  ggplot(aes(x = Year, group = as.factor(ind))) +
  geom_ribbon(aes(ymin = ll, ymax = hh, fill = as.factor(ind)), 
              alpha = 0.5, col = NA) +
  geom_line(aes(y = m, col = as.factor(ind))) +
  geom_point(aes(y = DBH, col = as.factor(ind))) +
  xlim(0, NA) +
  facet_wrap(~ species) +
  scale_color_discrete(guide = "none") +
  scale_fill_discrete(guide = "none") +
  ylab("Diameter at breast height (cm)") +
  theme(strip.text = element_text(face = "italic"))
```