This worksheet provides research question scenarios which utilise estimates of stem growth rates, rates of productivity, biomass change, and vital rates. These scenarios are designed to strengthen understanding of the code presented in the previous part of the workshop, and to stimulate your thinking on how you could use these techniques in your own research.
The code chunks are all written in R.
First, load the required packages. If you haven’t installed the packages you
should do so with install.packages("<pkg>"), where <pkg> is the name of the
package.
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. For example:
s <- read.csv("seosaw_data/stems_ABG_clean.csv")
p <- read.csv("seosaw_data/plots_ABG.csv")Some tree species grow fast and die young, while others grow slowly and live long. In this scenario, you should aim to calculate average diameter growth rates for the most common species in your dataset, and compare across species.
Which species grow the fastest? Which grow the slowest? Does this match your understanding of those species’ functional traits?
Are there certain genera which grow particularly quickly, or slowly?
How do growth rates vary with stem size among common species? How could you visualise this?
Do the fast growing species also have high rates of recruitment, or high rates of mortality?
What traits do you think should correlate with a fast-growing, short-lived tree?
To find the n (5) most common species or genera in your dataset:
s %>%
group_by(species_name_clean) %>%
tally() %>%
slice_max(n, n = 5)The genus of each individual is found in column genus_clean.
To visualise variation in growth rates among common species, by diameter size class, you could do something like this:
# Find the 3 most common species
common_sp <- s %>%
group_by(species_name_clean) %>%
tally() %>%
slice_max(n, n = 3) %>%
pull(species_name_clean)
# Filter stem data to common species, add diameter size class bins
s_dsc <- s %>%
filter(species_name_clean %in% common_sp) %>%
mutate(dsc = cut(diam, breaks = c(5, 10, 25, 50, 200),
include.lowest = TRUE)) %>%
filter(!is.na(dsc))
# Calculate growth for each individual, with groups of species and size class
g_dsc <- growthAll(s_dsc, w = "diam",
group = c("plot_id", "species_name_clean", "dsc", "stem_id"),
census = "census_date")
g_dsc$gr <- g_dsc$g / g_dsc$int
# Summarise to mean growth and 95% CI
g_group_summ <- g_dsc %>%
group_by(species_name_clean, 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)
# Plot points by diameter size class and species, with error bars
ggplot(g_group_summ, aes(x = dsc, y = gr_mean, colour = species_name_clean)) +
geom_point(position = position_dodge(width = 0.5)) +
geom_line(aes(group = species_name_clean),
position = position_dodge(width = 0.5)) +
geom_errorbar(aes(ymin = gr_mean - gr_ci95, ymax = gr_mean + gr_ci95),
width = 0.2, position = position_dodge(width = 0.5))To calculate recruitment rates for each species within each plot:
# 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_m <- s %>%
filter(
plot_id %in% s_multi_plot_id,
species_name_clean %in% common_sp)
# Split by plot and common species
s_split <- split(s_m, list(s_m$plot_id, s_m$species_name_clean), drop = TRUE)
s_split_fil <- s_split[unlist(lapply(s_split, nrow)) > 50]
vital_all <- do.call(rbind, lapply(s_split_fil, 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
vK_census <- lapply(comb_list, function(y) {
# Estimate vital rates according to Kohyama
vitalKohyama(x, t0 = y[1], tT = y[2],
group = "stem_id", census = "census_date", plot_area = 1)
})
# Create dataframes of metrics
vK <- data.frame()[seq_along(vK_census), ]
for (y in names(vK_census[[1]])) {
vK[[y]] <- unlist(lapply(vK_census, "[[", y))
}
# Add plot ID and species name
vK$plot_id <- unique(x$plot_id)
vK$species_name_clean <- unique(x$species_name_clean)
# Return
vK
}))
# Plot growth rate vs. annual mortality rate.
g_dsc %>%
left_join(., vital_all[,c("plot_id", "species_name_clean", "Ma")],
by = c("plot_id", "species_name_clean")) %>%
ggplot(., aes(x = gr, y = Ma)) +
geom_point()In this scenario you should aim to calculate diameter growth rates for a subset of species from your plots which provide commercially valuable timber. Then, find out how many years it would take for the average stem of each species to reach harvestable size (e.g. 35 cm DBH). Finally, determine which plots are experiencing the highest rates of basal area growth, to prioritise plots for selected timber removal.
Which species are growing the fastest?
Which plots are growing the fastest?
Do multi-stemmed trees grow faster than single stemmed trees?
What other factors other than growth rates should be taken into account when determining sustainable harvest rates?
To find out how many years it would take for a species to grow from 5 cm (minimum diameter threshold) to harvestable size (35 cm):
# Choose a single species
s_sp <- s %>%
filter(species_name_clean == "Baikiaea plurijuga")
# Calculate growth rates
g_sp <- growthAll(s_sp, w = "diam", group = c("plot_id", "stem_id"),
census = "census_date", type = "consecutive")
g_sp$gr <- g_sp$g / g_sp$int
# Estimate mean growth rate
g_sp_summ <- g_sp %>%
filter(!is.na(gr)) %>%
summarise(
gr_mean = mean(gr),
gr_sd = sd(gr),
n = n()) %>%
mutate(
gr_sem = gr_sd / sqrt(n),
gr_ci95 = gr_sem * 1.96)
# How many years to grow 30 cm (5-35 cm)?
30 / g_sp_summ$gr_meanTo find which plots are growing the fastest, by basal area:
# Find basal area of plots per census
pba <- s %>%
group_by(plot_id, census_date) %>%
summarise(ba_sum = sum(ba, na.rm = TRUE))
# Find consecutive growth intervals of plots, and sort descending order
growthAll(pba, w = "ba_sum", group = "plot_id", census = "census_date") %>%
left_join(., unique(p[,c("plot_id", "plot_area")]), by = "plot_id") %>%
mutate(gr_ha = g / int / plot_area) %>%
arrange(desc(gr_ha))To compare single-stemmed and multi-stemmed trees:
# Find number of stems per tree
tr_n <- s %>%
group_by(plot_id, tree_id, census_date) %>%
tally(name = "n_stems")
# Calculate effective diameter for each tree
tr <- equivDiamGen(s, plot_id = "plot_id", census_date = "census_date",
tree_id = "tree_id", diam = "diam") %>%
left_join(., tr_n, by = c("plot_id", "tree_id", "census_date"))
# Calculate tree diameter growth rates
tr_g <- growthAll(tr, w = "equiv_diam", group = c("plot_id", "tree_id", "n_stems"),
census = "census_date", type = "consecutive") %>%
mutate(
gr = g / int,
multi_stem = ifelse(n_stems > 1, "1", "0"))
# Run mixed model with binary multi-stemmed fixed effect
mod <- lmer(gr ~ multi_stem + (1|plot_id), data = tr_g)
# Summarise model
summary(mod)
# Plot model predictions
plot(ggpredict(mod))Competition among individuals has been shown to be a key determinant of stem growth rates. In this scenario, you should aim to calculate diameter growth rates for all individuals in each plot in your dataset, then compare their growth with plot basal area per hectare, which is a commonly used proxy for competition.
Is there a relationship between diameter growth rates and plot basal area?, Why? Why not?
How does the effect of plot basal area on growth rates vary between small and large stems?
Do plots with high basal area experience higher or lower rates of biomass productivity?
How else might you quantify tree-tree competition within a plot?
To compare plot basal area and diameter growth rates:
# Find plot basal area per census
pba <- s %>%
group_by(plot_id, census_date) %>%
summarise(ba_sum = sum(ba, na.rm = TRUE))
# Calculate consecutive stem growth rates
sg <- growthAll(s, w = "diam", group = c("plot_id", "stem_id"),
census = "census_date", type = "consecutive") %>%
left_join(., pba, by = c("plot_id", "t0" = "census_date"))
# Linear model of plot basal area vs. stem growth rates
mod1 <- lmer(g ~ ba_sum + (1|plot_id), data = sg)
summary(mod1)
plot(ggpredict(mod1))To extend the model to consider small and large trees:
# Add a binary variable for small and large trees
sg$tree_size <- ifelse(sg$w0 < 20, "sm", "lg")
mod2 <- lmer(g ~ ba_sum + ba_sum:tree_size + (1|plot_id), data = sg)
summary(mod2)
plot(ggpredict(mod2, terms = c("ba_sum[5:12]", "tree_size")))
# Add a continuous variable for tree size
mod3 <- lmer(g ~ ba_sum + ba_sum:w0 + (1|plot_id), data = sg)
summary(mod3)
plot(ggpredict(mod3, terms = c("ba_sum[5:12]", "w0"))) To calculate biomass productivity and compare to plot basal area:
# 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 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",
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
}))
# Find plot basal area per census
pba <- s %>%
group_by(plot_id, census_date) %>%
summarise(ba_sum = sum(ba, na.rm = TRUE)) %>%
left_join(., prod_all, by = c("plot_id", "census_date" = "t0")) %>%
filter(!is.na(int))
mod1 <- lm(P_ann ~ ba_sum, data = pba)
summary(mod1)
plot(ggpredict(mod1))