1 Introduction

This worksheet provides examples of how to estimate stem growth rates. The code chunks are all written in R. The example cleaned dataset (“stems_ABG_clean.csv”) can be found in the workshop Sharepoint folder.

First, load the required packages. If you haven’t installed the packages you should refer to the worksheet from the preparation session.

library(dplyr)
library(ggplot2)
library(lme4)
library(ggeffects)
library(seosawr)
# source("rscripts/seosawr_fn.R")
library(ProdVital)
# source("rscripts/ProdVital_fn.R")

Next, load the data. Either use the data you cleaned in the previous part of the workshop, or load one of the cleaned example datasets,

s <- read.csv("seosaw_data/stems_ABG_clean.csv")
p <- read.csv("seosaw_data/plots_ABG.csv")

2 Growth increment

To calculate diameter growth increment \(G_{Dt}\), simply find the difference in diameter of the stem between the start (\(D_{t}\)) and the end ($D_{t+1}) of the census interval.

\[ G_{Dt} = D_{t+1} - D_{t} \]

You can use the ProdVital::growthAll() function to calculate growth increment for every available consecutive or pairwise census interval for each stem in a dataset. Where a stem measurement in a particular census is missing, the function simply extends the census interval to the next non-missing value.

g_all <- growthAll(s, w = "diam", group = c("plot_id", "stem_id"), 
  census = "census_date", type = "consecutive")

head(g_all)
##          plot_id stem_id   t0   tT int   w0   wT    g
## ABG_10.1  ABG_10       1 2019 2021   2  5.6  6.2  0.6
## ABG_12.1  ABG_12       1 2019 2021   2  4.7 40.6 35.9
## ABG_14.1  ABG_14       1 2019 2021   2  7.5  7.0 -0.5
## ABG_2.1    ABG_2       1 2018 2021   3  6.7  7.6  0.9
## ABG_3.1    ABG_3       1 2018 2021   3 10.6 17.9  7.3
## ABG_4.1    ABG_4       1 2018 2021   3 12.0 11.0 -1.0

The group argument takes a column name, or multiple columns, which define the IDs of individuals between censuses. w is the name of the column which is used to calculate growth. census gives the name of the column where census dates are stored, and type states whether to calculate “consecutive” or “pairwise” growth increments.

The column g gives the growth increment of the stem for that census interval, where t0 is the date of the first census, tT is the date of the final census, int is the census interval length, w0 is the diameter of the stem in the first census, and wT is the diameter in the final census.

You can convert growth increment to a rate by dividing by the census interval (int, \(T\)).

\[ G_{D} = \frac{D_{t+1} - D_{t}}{T} \]

g_all$gr <- g_all$g / g_all$int

You can read the documentation for this function, or any other function in {ProdVital} or {seosawr} in R by typing:

?growthAll

And if you’re really interested, you can read the source code for the function by typing the name of the function without parentheses:

growthAll

3 Average stem growth rates

Here, we present two ways to calculate the average growth rate of a stem across a broader census period with many census intervals. The first way is to simply calculate the mean of all consecutive stem growth rates:

gr_all_mean <- g_all %>% 
  group_by(plot_id, stem_id) %>% 
  summarise(gr_mean = mean(gr))
## `summarise()` has grouped output by 'plot_id'. You can override using the
## `.groups` argument.
head(gr_all_mean)
## # A tibble: 6 × 3
## # Groups:   plot_id [1]
##   plot_id stem_id gr_mean
##   <chr>   <chr>     <dbl>
## 1 ABG_10  1        0.300 
## 2 ABG_10  10       0.0500
## 3 ABG_10  100      0.5   
## 4 ABG_10  101     -0.100 
## 5 ABG_10  102     -0.100 
## 6 ABG_10  103      0

The second, is to fit a linear model to the diameter timeline. This is a rather simplistic method, as stem growth may not always be linear over long time scales, but it is a sufficient starting point. For this, we can use the ProdVital::growthMod() function. The argument full states whether to return a dataframe (TRUE), or a named vector (FALSE).

s_g_mod <- growthMod(s, w = "diam",
  group = c("plot_id", "stem_id"), census = "census_date", full = TRUE)
An example of a linear model fitted through a long stem diameter timelines. The slope of this linear model will provide the annual rate of growth.

Figure 3.1: An example of a linear model fitted through a long stem diameter timelines. The slope of this linear model will provide the annual rate of growth.

4 Stems vs. trees

Trees can have multiple stems. Depending on your research question you may want to measure the growth of individual stems, or whole trees. To measure the growth of multi-stemmed trees we can calculate an “effective” diameter (\(D_{e}\)) for the tree. In other words the diameter the tree would have if all the stems were combined into a single stem. We can do this using the seosawr::equivDiamGen() function.

\[ D_{e} = 2 \sqrt{\frac{\pi{}(D\times{}0.5)^2}{\pi{}}} \]

trees <- equivDiamGen(s, plot_id = "plot_id", census_date = "census_date",
  tree_id = "tree_id", diam = "diam")

head(trees)
##   plot_id census_date tree_id equiv_diam
## 1   ABG_1        2018       1       8.30
## 2   ABG_2        2018       1       6.70
## 3   ABG_3        2018       1      10.60
## 4   ABG_4        2018       1      12.00
## 5  ABG_10        2019       1      10.01
## 6  ABG_12        2019       1       4.70

Then we can use these tree-level measures of size to calculate growth rate, using the same growthAll() function we used previously.

g_trees <- growthAll(trees, w = "equiv_diam", group = c("plot_id", "tree_id"), 
  census = "census_date", type = "consecutive")
g_trees$gr <- g_trees$g / g_trees$int

5 Scaling up to plots: basal area growth rates

We can calculate growth rates at the scale of an entire plot by calculating growth by basal area.

First, calculate the total basal area of each plot in each census, using the seosawr::basalAreaGen() function to convert between diameter in cm and basal area in m2.

pba <- s %>%
  mutate(ba = basalAreaGen(., "diam")) %>% 
  group_by(plot_id, census_date) %>% 
  summarise(ba_sum = sum(ba, na.rm = TRUE))
## `summarise()` has grouped output by 'plot_id'. You can override using the
## `.groups` argument.

Then, using growthAll() we can calculate the basal area growth increment.

pbg <- growthAll(pba, w = "ba_sum", group = "plot_id", census = "census_date")

And finally, divide by the census interval and the plot area to get a growth rate in units of m2 ha-1 y-1.

pbg_a <- pbg %>% 
  left_join(., unique(p[,c("plot_id", "plot_area")]), by = "plot_id") %>%
  mutate(gr_ha = g / int / plot_area)

6 More complicated diameter growth models

So far, our estimates of diameter growth have assumed a constant rate of growth. However, diameter growth models can be made more complex by adding predictors to the model; factors which we think modulate the rate of growth.

Commonly, initial stem size is added as a predictor. Diameter growth of trees tends to increase to a maximum early in life and then slowly decreases, approaching zero as the tree matures.

An example from Bowman et al. (2012), showing how growth increment varies with tree age and size.

Figure 6.1: An example from Bowman et al. (2012), showing how growth increment varies with tree age and size.

Proxies for competition among trees within the stand, such as basal area per hectare, are also common.

In this example, we fit a linear diameter growth model which accounts for the effect of tree initial diameter (w0) when predicting growth.

# Fit the model
mod1 <- lm(gr ~ 1 + w0, data = g_all)

# Summarise the model
summary(mod1)
## 
## Call:
## lm(formula = gr ~ 1 + w0, data = g_all)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.439  -0.142  -0.023   0.128  17.918 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.040007   0.015209    2.63   0.0085 **
## w0          -0.001695   0.000993   -1.71   0.0877 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.711 on 6209 degrees of freedom
##   (264 observations deleted due to missingness)
## Multiple R-squared:  0.000469,   Adjusted R-squared:  0.000309 
## F-statistic: 2.92 on 1 and 6209 DF,  p-value: 0.0877
# Extract predictions from the model
mod1_p <- data.frame(w0 = seq(5, 100))
mod1_preds <- predict(mod1, newdata = mod1_p, se.fit = TRUE)
mod1_p$fit <- mod1_preds$fit
mod1_p$se.fit <- mod1_preds$se.fit

Here, the model summary shows us that initial diameter (w0) has a small and non-significant negative effect on stem diameter growth rate.

Or, we may want to include a quadratic term to better represent the compounding decrease in diameter growth with size:

# Fit the model
mod2 <- lm(gr ~ 1 + w0 + I(w0^2), data = g_all)

# Summarise the model
summary(mod2)
## 
## Call:
## lm(formula = gr ~ 1 + w0 + I(w0^2), data = g_all)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.451  -0.155  -0.023   0.127  17.928 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.68e-02   2.36e-02    0.71     0.48
## w0           1.39e-03   2.59e-03    0.54     0.59
## I(w0^2)     -6.34e-05   4.92e-05   -1.29     0.20
## 
## Residual standard error: 0.711 on 6208 degrees of freedom
##   (264 observations deleted due to missingness)
## Multiple R-squared:  0.000737,   Adjusted R-squared:  0.000415 
## F-statistic: 2.29 on 2 and 6208 DF,  p-value: 0.101
# Extract predictions from the model
mod2_p <- data.frame(w0 = seq(5, 100))
mod2_preds <- predict(mod2, newdata = mod2_p, se.fit = TRUE)
mod2_p$fit <- mod2_preds$fit
mod2_p$se.fit <- mod2_preds$se.fit

To illustrate the differences between the fits of these models, we can then plot the predicted values from these models to visualise how they differ.

ggplot() + 
  geom_ribbon(data = mod2_p, 
    aes(x = w0, ymin = fit - se.fit, ymax = fit + se.fit),
    alpha = 0.2, fill = "red") + 
  geom_line(data = mod2_p, 
    aes(x = w0, y = fit), colour = "red") + 
  geom_ribbon(data = mod1_p, 
    aes(x = w0, ymin = fit - se.fit, ymax = fit + se.fit),
    alpha = 0.2, fill = "blue") + 
  geom_line(data = mod1_p, 
    aes(x = w0, y = fit), colour = "blue")

Using this framework, hopefully you get a sense of how models of diameter increment can be extended to much more complex models including many predictors.

7 Aggregating to demographic and taxonomic groups

We often have populations of stems within plots, with variation among individuals in their growth rate due to various factors, some measured, others unmeasured. You may want to estimate the growth rate of the average stem, with some measure of uncertainty, from a given stem size class, species, or some other grouping that makes sense for your research question.

We can estimate the sample mean and uncertainty of a distribution of growth rates. In this example, we will explore variation in stem growth rates within different diameter size classes, across all plots.

First, use cut() to define diameter size classes (dsc), and remove measurements below the minimum diameter threshold (5 cm).

s_dsc <- s %>% 
  mutate(dsc = cut(diam, breaks = c(5, 10, 25, 50, 200), 
      include.lowest = TRUE)) %>% 
  filter(!is.na(dsc))

Then, calculate the consecutive growth rates for each stem. Notice that we include dsc as a grouping factor, so it appears in the output dataframe.

g_dsc <- growthAll(s_dsc, w = "diam", group = c("plot_id", "dsc", "stem_id"), 
  census = "census_date")
g_dsc$gr <- g_dsc$g / g_dsc$int

Calculate the sample means (gr_mean) and 95% confidence interval (gr_95ci) for each diameter size class.

g_group_summ <- g_dsc %>% 
  group_by(dsc) %>% 
  summarise(
    gr_mean = mean(gr),
    gr_sd = sd(gr),
    n = n()) %>% 
  mutate(
    gr_sem = gr_sd / sqrt(n),
    gr_ci95 = gr_sem * 1.96)

And finally visualise the difference in distributions and sample means among diameter size classes.

ggplot() + 
  geom_histogram(data = g_dsc, 
    aes(x = gr, y = after_stat(density)), 
    colour = "black", fill = "grey", binwidth = 0.2) + 
  geom_vline(data = g_group_summ, aes(xintercept = gr_mean), 
    linetype = 2, colour = "red") + 
  geom_point(data = g_group_summ, aes(x = gr_mean, y = 0.5)) + 
  geom_errorbarh(data = g_group_summ, 
    aes(xmin = gr_mean - gr_ci95, xmax = gr_mean + gr_ci95, y = 0.5),
    height = 0.1) + 
  facet_wrap(~dsc, scales = "fixed", ncol = 1) + 
  xlim(-2, 2)
## Warning: Removed 33 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 8 rows containing missing values (`geom_bar()`).

An alternative way to estimate population rates of growth and their variance is to fit a linear mixed effects model, with the response of growth rate. By including random effects for individual or plot, we can account for non-independence of measurements when we have very unequal sampling across these groups.

mod3 <- lmer(gr ~ 1 + (1|plot_id), data = g_all)

Using summary(), we can then extract the population level estimate of the annual growth rate.

mod3_summ <- summary(mod2)
mod3_summ$coefficients[1]
## [1] 0.01679
mod3_summ$coefficients[2]
## [1] 0.001394

We can also then add fixed effects to investigate how these factors affect rates of growth, and how growth rates vary among groups. In this example, we add the effects of initial diameter w0, wood density wood_density_mean, and their interaction (with *).

# Add wood density from the original stems table
g_all_wd <- g_all %>% 
  left_join(., s[,c("plot_id", "stem_id", "census_date", "wood_density_mean")], 
    by = c("plot_id", "stem_id", "t0" = "census_date")) 

# Fit a model with the fixed effect of wood density
mod4 <- lmer(gr ~ w0*wood_density_mean + (1|plot_id), data = g_all_wd)

# Summarise the model
summary(mod4)
## Linear mixed model fit by REML ['lmerMod']
## Formula: gr ~ w0 * wood_density_mean + (1 | plot_id)
##    Data: g_all_wd
## 
## REML criterion at convergence: 13391
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -17.400  -0.212  -0.021   0.190  25.160 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plot_id  (Intercept) 0.00335  0.0579  
##  Residual             0.50231  0.7087  
## Number of obs: 6211, groups:  plot_id, 14
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)            0.2966     0.1228    2.42
## w0                    -0.0170     0.0101   -1.69
## wood_density_mean     -0.3389     0.1699   -2.00
## w0:wood_density_mean   0.0209     0.0143    1.46
## 
## Correlation of Fixed Effects:
##             (Intr) w0     wd_dn_
## w0          -0.787              
## wd_dnsty_mn -0.982  0.803       
## w0:wd_dnst_  0.770 -0.995 -0.798
# Extract and plot predicted values from the model
plot(ggpredict(mod4, 
    terms = c("w0", "wood_density_mean[0.25,0.5,0.75,1]"), type = "fe"))

8 Forecasting growth

Now that we have reliable sample estimates of growth rate, we can use those estimates to forecast growth into the future. A plausible research scenario would be to estimate how many years it would take a hypothetical stem from a given species to reach harvestable size (e.g. 30 cm).

First, join species information to each stem in g_all, filter to species of interest, and filter out any negative diameter increments.

g_sp <- g_all %>% 
  left_join(., 
    s[,c("plot_id", "stem_id", "census_date", "species_name_clean")],
    by = c("plot_id", "stem_id", "t0" = "census_date"))

sp_dom <- c("Baikiaea plurijuga", "Julbernardia paniculata", 
  "Brachystegia spiciformis", "Burkea africana")


g_sp_fil <- g_sp %>% 
  filter(
    species_name_clean %in% sp_dom,
    gr > 0)

Then, construct a mixed effects model with the fixed effect of species to predict growth rates, and the random effect of plot to account for non-independence of measurements:

mod5 <- lmer(gr ~ species_name_clean + (1|plot_id), data = g_sp_fil)

summary(mod5)
## Linear mixed model fit by REML ['lmerMod']
## Formula: gr ~ species_name_clean + (1 | plot_id)
##    Data: g_sp_fil
## 
## REML criterion at convergence: 4075
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -0.991 -0.281 -0.169  0.020 20.656 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plot_id  (Intercept) 0.028    0.167   
##  Residual             0.684    0.827   
## Number of obs: 1645, groups:  plot_id, 14
## 
## Fixed effects:
##                                            Estimate Std. Error t value
## (Intercept)                                   0.647      0.092    7.03
## species_name_cleanBrachystegia spiciformis   -0.245      0.106   -2.31
## species_name_cleanBurkea africana            -0.424      0.107   -3.97
## species_name_cleanJulbernardia paniculata    -0.385      0.106   -3.63
## 
## Correlation of Fixed Effects:
##             (Intr) sp__Bs sp__Ba
## spcs_nm_cBs -0.677              
## spcs_nm_cBa -0.747  0.736       
## spcs_nm_cJp -0.796  0.660  0.776
mod5_pred <- ggpredict(mod5, terms = "species_name_clean", type = "fe")

mod5_pred
## # Predicted values of gr
## 
## species_name_clean       | Predicted |       95% CI
## ---------------------------------------------------
## Baikiaea plurijuga       |      0.65 | [0.47, 0.83]
## Brachystegia spiciformis |      0.40 | [0.24, 0.56]
## Burkea africana          |      0.22 | [0.08, 0.36]
## Julbernardia paniculata  |      0.26 | [0.13, 0.39]
## 
## Adjusted for:
## * plot_id = 0 (population-level)
plot(mod5_pred)

Then we can calculate the number of years to a harvestable size of 30 cm:

mod5_pred <- as.data.frame(ggpredict(mod5, 
    "species_name_clean", type = "fe"))

# Mean growth rate over 30 years
mod5_pred$t30 <- 30 / mod5_pred$predicted

# Confidence intervals
mod5_pred$t30_lo <- 30 / mod5_pred$conf.low
mod5_pred$t30_hi <- 30 / mod5_pred$conf.high

… and finally plot the results:

ggplot(mod5_pred, aes(x = x, y = t30)) + 
  geom_bar(stat = "identity") + 
  geom_errorbar(aes(ymin = t30_lo, ymax = t30_hi), width = 0.2) + 
  labs(x = NULL, y = "Years to 30 cm")

9 Estimating biomass

In the next part of the session, we will move on from diameter growth rates towards estimates of productivity based on woody biomass. So, we need to estimate the biomass of each stem.

We will use one of the biomass allometries from (Chave et al. 2014), which uses wood density, diameter, and plot location to estimate above-ground woody biomass for each stem. We already have wood density estimates in the SEOSAW dataset for each individual, but if you ever need to calculate wood density again, you can use the seosawr::woodDensity() function.

Then, we can use the seosawr::chave14() function, passing the stems table and plots table as arguments, to estimate biomass (Mg).

s$agb_chave <- chave14(s, p, diam = "diam", 
  wood_density_mean = "wood_density_mean",
  longitude_of_centre = "longitude_of_centre", 
  latitude_of_centre = "latitude_of_centre")

10 Productivity, loss, biomass change

There is a close link between tree growth and higher-level estimations of tree community dynamics. Woody Net Primary Productivity (NPP) can be calculated as the biomass growth of individuals over a given interval. Loss can be calculated as the biomass removed from the living pool through mortality. Productivity minus loss gives us a rate of biomass change over the census interval, also known as Net Biome Production (NBP).

The ProdVital::prodKohyama() function uses methods from (Kohyama, Kohyama, and Sheil 2019) to estimate annual rates of productivity (\(P_a\)), loss (\(L_a\)), and biomass change (\(\Delta{}B\)) while accounting for the productivity not captured within the census interval from stems which died and/or recruited within the interval. These annual estimates are in units of Mg ha^-1$ y-1.

\[ B_w = \left( B_T - B_0 \right) \left( {B_T / B_0}^{\frac{1}{T}} \right) \]

\[ P_a = B_w \left( ( {B_T / B_0}^{\frac{1}{T}} ) ( 1 - {B_{s0}B_T}^{\frac{1}{T}} ) \right) \]

\[ L_a = B_w \left( 1 - {B_{s0}B_0}^{\frac{1}{T}} \right) \]

\[ \Delta{}B = P_a - L_a \]

The prodKohyama() function only measures productivity within a single census interval for a single plot, but it’s fairly easy to extend the function to measure for all consecutive intervals for all plots. First, for a single interval in a single plot.

s_ABG_5 <- s %>% 
  filter(
    plot_id == "ABG_5",
    census_date %in% c(2019, 2021))

pK_ABG_5 <- prodKohyama(s_ABG_5, 2019, 2021, "agb_chave", "stem_id", "census_date")

And for all censuses in all plots:

# Remove plots with only one census
s_multi_plot_id <- s %>% 
  dplyr::select(plot_id, census_date) %>% 
  distinct() %>% 
  group_by(plot_id) %>%
  tally() %>% 
  filter(n > 1) %>% 
  pull(plot_id) 

s_multi <- s %>% 
  filter(plot_id %in% s_multi_plot_id)

s_split <- split(s_multi, s_multi$plot_id)

prod_all <- do.call(rbind, lapply(s_split, function(x) { 
  # Find all census dates
  census_date_all <- sort(unique(x$census_date))

  # Create all pairwise combinations of censuses
  comb_list <- combn(census_date_all, 2, simplify = FALSE)

  # Alternatively, every consecutive interval
  # tT <- c(census_date_all, NA_real_)
  # t0 <- c(NA_real_, census_date_all)
  # comb_df <- data.frame(t0, tT)
  # comb_df <- comb_df[!is.na(comb_df$t0) & !is.na(comb_df$tT),]
  # comb_list <- apply(comb_df, 1, "[", simplify = FALSE)
  
  # For each of the pairwise combinations of censuses
  pK_census <- lapply(comb_list, function(y) { 
    # Estimate productivity according to Kohyama 
    prodKohyama(x, t0 = y[1], tT = y[2], w = "agb_chave", 
      group = "stem_id", census = "census_date")
  })

  # Create dataframes of metrics
  pK <- data.frame()[seq_along(pK_census), ]
  for (y in names(pK_census[[1]])) {
    pK[[y]] <- unlist(lapply(pK_census, "[[", y))
  }

  # Add plot ID
  pK$plot_id <- unique(x$plot_id)

  # Return
  pK
}))
paged_table(head(prod_all))

As with all the other functions, you can check the description of the values returned by the function by typing ?prodKohyama in R.

Probably the most important values returned are P_ann, the annual rate of productivity (\(P_a\)), L_ann, the annual rate of loss (\(L_a\)), and dB_ann, the annual rate of change (\(\Delta{}B\)) (i.e. P_ann - L_ann).

11 Recruitment and mortality rates

Simulation models of community dynamics frequently combine rates of growth with rates of recruitment and mortality. The ProdVital::vitalKohyama() function uses methods from (Kohyama, Kohyama, and Sheil 2018) to estimate rates of recruitment an mortality, either as per capita (individual-based, \(m\), \(l\)) or areal (area-based, \(M_a\), \(R_a\)) rates. Areal rates are in units of N stems ha-1 y-1, while instantaneous rates can be interpreted as the proportional likelihood of a given stem recruiting or dying per year.

\[ M_a = (N_0 / A) \left( 1 - ( N_{S_{T}} / N_0 )^{\frac{1}{T}} \right) \]

\[ R_a = M_a ( N_T - N_{S_{T}} ) / ( N_0 - N_{S_{T}} ) \]

As with prodKohyama(), vitalKohyama() only estimates recruitment and mortality within a single census interval for a single plot by default.

vK_ABG_5 <- vitalKohyama(s_ABG_5, 2019, 2021, group = "stem_id", 
  census = "census_date", plot_area = 1)

unlist(vK_ABG_5)
##       int        t0        tT        N0        NT       NST        nT         g 
## 2.000e+00 2.019e+03 2.021e+03 6.560e+02 6.470e+02 6.160e+02 3.100e+01 2.455e-02 
##         m         r         l        ma        ra       raf       ras       raz 
## 3.146e-02 2.455e-02 9.931e-01 3.097e-02 2.647e-02 2.425e-02 2.485e-02 2.336e-02 
##         M         R        Ma        Ra       Ras 
## 2.064e+01 1.599e+01 2.031e+01 1.574e+01 1.625e+01

vitalKohyama() can be scaled up to measure the vital rates of multiple census combinations in exactly the same way as with the prodKohyama() function above.

12 Finishing off

Remember to save the growth rates estimates to a file so they can be used in other parts of the workshop.

write.csv(g_all, "seosaw_data/ABG_gr.csv", row.names = FALSE)

References

Chave, J., M. Réjou-Méchain, A. Búrquez, E. Chidumayo, M. S. Colgan, W. B. C. Delitti, A. Duque, et al. 2014. “Improved Allometric Models to Estimate the Aboveground Biomass of Tropical Trees.” Global Change Biology 20 (10): 3177–90. https://doi.org/http://dx.doi.org/10.1111/gcb.12629.
Kohyama, Takashi S., Tetsuo I. Kohyama, and Douglas Sheil. 2018. “Definition and Estimation of Vital Rates from Repeated Censuses: Choices, Comparisons and Bias Corrections Focusing on Trees.” Edited by Mark Rees. Methods in Ecology and Evolution 9 (4): 809–21. https://doi.org/10.1111/2041-210x.12929.
———. 2019. “Estimating Net Biomass Production and Loss from Repeated Measurements of Trees in Forests and Woodlands: Formulae, Biases and Recommendations.” Forest Ecology and Management 433 (February): 729–40. https://doi.org/10.1016/j.foreco.2018.11.010.