Case study 2: von Bertalanffy growth with lizard size data
Source:vignettes/von-bertalanffy.Rmd
von-bertalanffy.Rmd
# remotes::install_github("traitecoevo/hmde")
# install.packages(c("dplyr", "ggplot2"))
library(hmde)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(ggplot2)Our second demo introduces size-dependent growth based on the von Bertalanffy function where is the asymptotic maximum size and controls the growth rate. We have implemented the analytic solution which is independent of age at the starting size and instead uses the first size as the initial condition. The key behaviour of the von Bertalanffy model is a high growth rate at small sizes that declines linearly as the size approaches . This manifests as growth slowing as a creature matures with a hard finite limit on the eventual size. We restrict and to be positive, and to be larger than the observed sizes. As a result the growth rate is non-negative.
In the following code we plot an example of the growth function and the solution to get a feel for the behaviour.
#Analytic solution in function form
solution <- function(t, pars = list(y_0, beta, S_max)){
return(
pars$S_max + (y_0 - pars$S_max)*exp(-t * pars$beta)
)
}
#Parameters
beta <- 0.35 #Growth rate
y_0 <- 1 #Starting size
S_max <- 20 #Asymptotic max size
time <- c(0,30)
pars_list <- list(y_0 = y_0,
beta = beta,
S_max = S_max)
y_final <- solution(time[2], pars_list)
#Plot of growth function
ggplot() +
xlim(y_0, y_final) +
ylim(0, beta*(S_max-y_0)*1.1) +
labs(x = "Y(t)", y = "f", title = "von Berralanffy growth") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title=element_text(size=18,face="bold")) +
geom_function(fun=hmde_model_des("vb_single_ind"),
args=list(pars = list(S_max, beta)),
colour="green4", linewidth=1,
xlim=c(y_0, y_final))
#Size over time
ggplot() +
geom_function(fun=solution,
args=list(pars = pars_list),
colour="green4", linewidth=1,
xlim=c(time)) +
xlim(time) +
ylim(0, y_final*1.05) +
labs(x = "Time", y = "Y(t)", title = "von Bertalanffy growth") +
theme_classic() +
theme(axis.text=element_text(size=16),
axis.title=element_text(size=18,face="bold"))
The von Bertalanffy model is commonly used in fishery management , but has also been used in reptile studies such as .
Lizard size data
Our data is sourced from which measured mass and snout-vent-length (SVL) of delicate skinks – – under experimental conditions to examine the effect of temperature on development. We are going to use the SVL metric for size.
We took a simple random sample without replacement of 50 individuals with at least 5 observations each. The von Bertalanffy model can be fit to shorter observation lengths, but fewer than 3 observations is not advised as there are two growth parameters per individual.
Implementation
The workflow for the second example is the same as the first, with the change in model name and data object.
lizard_vb_fit <- hmde_model("vb_multi_ind") |>
hmde_assign_data(data = Lizard_Size_Data) |>
hmde_run(chains = 4, cores = 1, iter = 2000)
#>
#> SAMPLING FOR MODEL 'vb_multi_ind' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000171 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.71 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 24.03 seconds (Warm-up)
#> Chain 1: 5.988 seconds (Sampling)
#> Chain 1: 30.018 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'vb_multi_ind' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 0.000132 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.32 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 7.591 seconds (Warm-up)
#> Chain 2: 7.402 seconds (Sampling)
#> Chain 2: 14.993 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'vb_multi_ind' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 0.000131 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.31 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3:
#> Chain 3:
#> Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 3:
#> Chain 3: Elapsed Time: 23.282 seconds (Warm-up)
#> Chain 3: 3.779 seconds (Sampling)
#> Chain 3: 27.061 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'vb_multi_ind' NOW (CHAIN 4).
#> Chain 4:
#> Chain 4: Gradient evaluation took 0.000151 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.51 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4:
#> Chain 4:
#> Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 4:
#> Chain 4: Elapsed Time: 20.8 seconds (Warm-up)
#> Chain 4: 3.843 seconds (Sampling)
#> Chain 4: 24.643 seconds (Total)
#> Chain 4:
#> Warning: There were 120 divergent transitions after warmup. See
#> https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See
#> https://mc-stan.org/misc/warnings.html#bfmi-low
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess
lizard_vb_estimates <- hmde_extract_estimates(model = "vb_multi_ind",
fit = lizard_vb_fit,
input_measurement_data = Lizard_Size_Data)As before, we can compare the observed sizes over time to those predicted by the model.
measurement_data_transformed <- lizard_vb_estimates$measurement_data %>%
group_by(ind_id) %>%
mutate(
delta_y_obs = y_obs - lag(y_obs),
obs_interval = time - lag(time),
obs_growth_rate = delta_y_obs/obs_interval,
delta_y_est = y_hat - lag(y_hat),
est_growth_rate = delta_y_est/obs_interval
) %>%
ungroup()
#Distributions of estimated growth and size
hist(measurement_data_transformed$y_hat,
main = "Estimated size distribution",
xlab = "Size (cm)")
hist(measurement_data_transformed$delta_y_est,
main = "Estimated growth increments",
xlab = "Growth increment (cm)")
hist(measurement_data_transformed$est_growth_rate,
main = "Estimated annualised growth rate distribution",
xlab = "Growth rate (cm/yr)")
#Quantitative R^2
cor(measurement_data_transformed$y_obs, measurement_data_transformed$y_hat)^2
#> [1] 0.7470669
#Plots of size over time for a sample of 5 individuals
sample_ids <- sample(1:nrow(lizard_vb_estimates$individual_data), size=5)
plot_data <- measurement_data_transformed %>%
filter(ind_id %in% sample_ids)
ggplot(data=plot_data, aes(group = ind_id)) +
geom_point(aes(x = time, y=y_obs, colour = as.factor(ind_id)),
shape = 1) +
geom_line(aes(x = time, y=y_obs, colour = as.factor(ind_id)),
linetype = "dashed") +
geom_point(aes(x = time, y=y_hat, colour = as.factor(ind_id)),
shape = 2) +
geom_line(aes(x = time, y=y_hat, colour = as.factor(ind_id)),
linetype = "solid") +
labs(x="Time (days)", y="Size (mm)", colour="Ind. ID") +
theme_classic()
We have two parameters at the individual level and are interested in both their separate distributions, and if we see evidence of a relationship between them. We can also use the individual parameter estimates and estimated sizes to plot the growth function pieces.
#1-dimensional parameter distributions
hist(lizard_vb_estimates$individual_data$ind_max_size_mean,
main = "Individual max size parameters",
xlab = "Smax estimate")
hist(lizard_vb_estimates$individual_data$ind_growth_rate_mean,
main = "Individual growth rate parameters",
xlab = "beta estimate")
#2-dimensional parameter distribution
ggplot(data = lizard_vb_estimates$individual_data,
aes(x = ind_max_size_mean, y = ind_growth_rate_mean)) +
geom_point(shape = 16, size = 1, colour = "green4") +
xlab("Individual max sizes (mm)") +
ylab("Individual growth rate parameters") +
theme_classic()
#Correlation of parameters
cor(lizard_vb_estimates$individual_data$ind_max_size_mean,
lizard_vb_estimates$individual_data$ind_growth_rate_mean)
#> [1] 0.4733599
#Plot function pieces over estimated sizes.
hmde_plot_de_pieces(model = "vb_multi_ind",
individual_data = lizard_vb_estimates$individual_data,
measurement_data = lizard_vb_estimates$measurement_data)
At the hyper-parameter level for the whole population we have centre and spread parameters for the log-normal distributions of and . As before, we can look at these as species-level features.
#Max size
lizard_vb_estimates$population_data$mean[1] #Raw value
#> [1] 3.198664
print(paste0("95% CI for mean log max size: (",
lizard_vb_estimates$population_data$CI_lower[1], " , ",
lizard_vb_estimates$population_data$CI_upper[1], ")")) #Raw CI
#> [1] "95% CI for mean log max size: (3.17906559168809 , 3.21938603555537)"
exp(lizard_vb_estimates$population_data$mean[1]) #In mm units
#> [1] 24.49979
print(paste0("95% CI for mean max size in mm: (",
exp(lizard_vb_estimates$population_data$CI_lower[1]), " , ",
exp(lizard_vb_estimates$population_data$CI_upper[1]), ")"))
#> [1] "95% CI for mean max size in mm: (24.0242945602394 , 25.0127585216695)"
#Standard deviation of underlying normal distribution
lizard_vb_estimates$population_data$mean[2]
#> [1] 0.03385361
print(paste0("95% CI for log max size standard deviation: (",
lizard_vb_estimates$population_data$CI_lower[2], " , ",
lizard_vb_estimates$population_data$CI_upper[2], ")")) #Raw CI
#> [1] "95% CI for log max size standard deviation: (0.0193396389437248 , 0.049790035571581)"
#Beta
lizard_vb_estimates$population_data$mean[3] #Raw value
#> [1] -4.000668
print(paste0("95% CI for mean log growth par: (",
lizard_vb_estimates$population_data$CI_lower[3], " , ",
lizard_vb_estimates$population_data$CI_upper[3], ")")) #Raw CI
#> [1] "95% CI for mean log growth par: (-4.19076174579693 , -3.80345960817351)"
exp(lizard_vb_estimates$population_data$mean[3]) #In cm/yr units
#> [1] 0.0183034
print(paste0("95% CI for mean growth par mm/yr: (",
exp(lizard_vb_estimates$population_data$CI_lower[3]), " , ",
exp(lizard_vb_estimates$population_data$CI_upper[3]), ")"))
#> [1] "95% CI for mean growth par mm/yr: (0.0151347516474566 , 0.0222935114733902)"
#Standard deviation of underlying normal distribution
lizard_vb_estimates$population_data$mean[4]
#> [1] 0.1122711
print(paste0("95% CI for log growth par standard deviation: (",
lizard_vb_estimates$population_data$CI_lower[4], " , ",
lizard_vb_estimates$population_data$CI_upper[4], ")")) #Raw CI
#> [1] "95% CI for log growth par standard deviation: (0.038047243081422 , 0.259142006998395)"