Creating custom models
hatchR also includes basic functionality to generate your own model parameterizations for predicting hatching and emergence phenology using your own data. Importantly, this functionality implements model form #2 of Beacham and Murray (1990), which we chose because of its overall simplicity and negligible loss of accuracy. See Beacham and Murray (1990) and Sparks et al. (2019) for more specific discussion regarding model #2 and the development of the effective value approach.
The model follows the general format of:
\[ Effective Value_i = 1/exp(log_ea - log_e(Temperature_i - b)) \]
Where i is the daily value and a fish hatches or emerges when the cumulative sum reaches one: \[\sum_{i =1}^nEffectiveValue_i = 1\]
fit_model()
The function fit_model() uses data in which average
incubation temperature (°C) and days to phenological event are the
inputs and estimates parameter coefficients for
logea and b using
stats::nls().
Here, we borrow data from Table 8.1 (pg. 183) from Quinn (2018) to generate a custom hatch parameterization for brown trout (Salmo trutta).
You could either create a .csv file with those data and import them
with readr::read_csv() or alternatively, directly input
them as an object in R. We’ll use tibble::tibble() to
create the data.
bt_data <- tibble(
temperature = c(2,5,8,11,14),
days_to_hatch = c(194,87,54,35,28)
)
bt_data
#> # A tibble: 5 × 2
#> temperature days_to_hatch
#> <dbl> <dbl>
#> 1 2 194
#> 2 5 87
#> 3 8 54
#> 4 11 35
#> 5 14 28We can plot our data for a validity check using ggplot2:
We can now use fit_model() to create our custom
parameterization from our data. You must specify a species and a
development type, information which is carried forward in subsequent
functions in hatchR.
bt_fit <- fit_model(temp = bt_data$temperature,
days = bt_data$days_to_hatch,
species = "brown_trout",
development_type = "hatch")The output of fit_model() is a list with several
elements:
bt_fit$model: a model object of class “nls” containing the nonlinear regression modelbt_fit$log_a: a named numeric vector of the estimated coefficient logeabt_fit$b: a named numeric vector of the estimated coefficient bbt_fit$r_squared: pseudo R-squared value (1 - (residual sum of squares / total sum of squares))bt_fit$mse: mean squared error (mean(residuals^2))bt_fit$rmse: root mean squared error (sqrt(mse)bt_fit$expression: a tibble with the species, development type, and the parameterized model expressionbt_fit$pred_plot: a ggplot object showing the observed data and predicted values
bt_fit
#> $model
#> Nonlinear regression model
#> model: y ~ a/(x - b)
#> data: df
#> a b
#> 439.1634 -0.2582
#> residual sum-of-squares: 36.92
#>
#> Number of iterations to convergence: 5
#> Achieved convergence tolerance: 2.812e-07
#>
#> $log_a
#> [1] 6.084872
#>
#> $b
#> [1] -0.258236
#>
#> $r_squared
#> [1] 0.9979987
#>
#> $mse
#> [1] 7.38454
#>
#> $rmse
#> [1] 2.717451
#>
#> $expression
#> # A tibble: 1 × 3
#> species development_type expression
#> <chr> <chr> <chr>
#> 1 brown_trout hatch 1 / exp(6.08487164096323 - log(x + 0.25823604410…
#>
#> $pred_plot
The vast majority of the time, what you will want is the actual
expression with parameter estimates for use in the
model = ... argument of predict_phenology().
This is stored in the expression element of the list.
You can either pass this tibble directly with the $
operator by calling $expression element of the list (e.g.,
model = bt_fit$expression) or set as an object to pass,
like so:
bt_hatch_exp <- bt_fit$expression
bt_hatch_exp
#> # A tibble: 1 × 3
#> species development_type expression
#> <chr> <chr> <chr>
#> 1 brown_trout hatch 1 / exp(6.08487164096323 - log(x + 0.25823604410…predict_phenology() will extract the expression from the
object and use it to predict phenology. We will demonstrate this in the
Predict
phenology: Basic vignette.
Fitting models for other fishes
We demonstrate how the fit_model() function may be used
to create custom parameterizations for species beyond the Salmonids in
the model_table included in the package. We include
parameterizations from three warm-water species to demonstrate the
fit_model() utility for fishes beyond the scope of the
original effective value approach. These parameterizations are for
commonly cultured sportfishes including Smallmouth Bass (Micropterus
dolomieu) (Webster 1948), Channel
Catfish (Ictalurus punctatus) (Small and
Bates 2001) , and Lake Sturgeon (Acipenser fulvescens)
(Smith and King 2005).
We demonstrate the utility of this approach by creating a random thermal regime with an ascending thermograph with a mean temperature of 16 °C, parameterizing models for each species, and demonstrating days to hatch and developmental period for each species with the random thermal regime .
### make temp regime
set.seed(123)
# create random temps and corresponding dates
temps_sim <- sort(rnorm(n =30, mean = 16, sd = 1), decreasing = FALSE)
dates_sim <- seq(from = ymd("2000-07-01"),
to = ymd("2000-07-31"), length.out = 30)
data_sim <- matrix(NA, 30, 2) |> data.frame()
data_sim[,1] <- temps_sim
data_sim[,2] <- dates_sim
# change names so they aren't the same as the vector objects
colnames(data_sim) <- c("temp_sim", "date_sim")Next we’ll parameterize our models for the three different fishes
### smallmouth mod
smallmouth <- matrix(NA, 10, 2) |> data.frame()
colnames(smallmouth) <- c("hours", "temp_F")
smallmouth$hours <- c(52, 54, 70, 78, 90, 98, 150, 167, 238, 234)
smallmouth$temp_F <- c(77, 75, 71, 70, 67, 65, 60, 59, 55, 55)
# change °F to °C and hours to days
smallmouth <- smallmouth |>
mutate(days = ceiling(hours/24),
temp_C = (temp_F -32) * (5/9))
# model object for smallmouth bass
smb_mod <- fit_model(temp = smallmouth$temp_C,
days = smallmouth$days,
species = "smb",
development_type = "hatch")
### catfish mod
catfish <- matrix(NA, 3, 2) |> data.frame()
colnames(catfish) <- c("days", "temp_C")
catfish$days <- c(16,21,26)
catfish$temp_C <- c(22,10,7)
cat_mod <- fit_model(temp = catfish$temp_C,
days = catfish$days,
species = "catfish",
development_type = "hatch")
### lake sturgeon mod
sturgeon <- matrix(NA, 7, 2) |> data.frame()
colnames(sturgeon) <- c("days", "CTU")
sturgeon$days <- c(7,5,6,6,5,11,7)
sturgeon$CTU <- c(58.1, 62.2, 61.1, 57.5, 58.1, 71.4, 54.7)
sturgeon <- sturgeon |>
mutate(temp_C = CTU/days) # change CTUs to average temp and add column
sturgeon_mod <- fit_model(days = sturgeon$days,
temp = sturgeon$temp_C,
species = "sturgeon",
development_type = "hatch")Note the R2 fit from the models below. You can see they generally all preform well and are in line with values from model 2 of Beacham and Murray (1990).
#model fits
smb_mod$r_squared; cat_mod$r_squared; sturgeon_mod$r_squared
#> [1] 0.9868067
#> [1] 0.9433598
#> [1] 0.9217358After we have our fits for each species we can predict phenology
using our data_sim datset we created above.
### predict_phenology
#smallmouth bass
smb_res <- predict_phenology(data = data_sim,
dates = date_sim,
temperature = temp_sim,
spawn.date = "2000-07-01",
model = smb_mod$expression)
# catfish
catfish_res <- predict_phenology(data = data_sim,
dates = date_sim,
temperature = temp_sim,
spawn.date = "2000-07-01",
model = cat_mod$expression)
# sturgeon
# note that 16 C is pretty far out of range of temps for model fit, not best practice
sturgeon_res <- predict_phenology(data = data_sim,
dates = date_sim,
temperature = temp_sim,
spawn.date = "2000-07-01",
model = sturgeon_mod$expression)After we have our predictions we’ll combine our results and do a little data cleaning.
# summary for all species
all_res <- data.frame(matrix(NA, 3, 2))
colnames(all_res) <- c("start", "stop")
all_res$start <- c(catfish_res$dev.period$start,
smb_res$dev.period$start,
sturgeon_res$dev.period$start)
all_res$stop <- c(catfish_res$dev.period$stop,
smb_res$dev.period$stop,
sturgeon_res$dev.period$stop)
all_res <- all_res |>
mutate(days = ceiling(stop-start),
index = c(17,16.5,16)) # index for our horizontal bars
all_res$Species <- c("Channel Catfish", "Smallmouth Bass", "Lake Sturgeon")And then finally we can plot our results.
ggplot() +
geom_point(data = data_sim, aes(x = date_sim, y = temp_sim )) +
geom_line(data = data_sim, aes(x = date_sim, y = temp_sim )) +
geom_rect(data = all_res, aes(xmin = start, xmax = stop, ymax =index-.35, ymin = index-.5, fill = Species)) +
geom_label(data = all_res, aes(x = start + (stop - start) / 1.25, y = (index -0.425), label = days)) +
labs(x = "Date", y = "Temperature (°C)") +
scale_fill_manual(values = c("deepskyblue4", "grey23", "darkolivegreen4")) +
theme_classic() +
theme(legend.position = c(0.75, 0.25))