Nothing
R/dm.growth.fit.double.R
In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data
Defines functions
dm.growth.fit.double
Documented in
dm.growth.fit.double
#' Fit bimodal dendrometer growth curves by vegetation season
#'
#' @description
#' Fits bimodal cumulative growth curves to daily dendrometer series using either
#' a double-Gompertz or double-Richards model.
#'
#' Overall growing-season timing is derived from the fitted cumulative-growth
#' curve using \code{growth_fraction}. Pulse-specific timing is derived from the
#' derivative of the fitted curve using a pulse-specific relative rate threshold
#' given by \code{rate_threshold_fraction}.
#'
#' Pulse-specific timing is returned as \code{NA} when the fitted curve does not
#' show a convincing two-pulse pattern according to derivative-based criteria.
#'
#' If \code{fallback_to_single = TRUE} and no convincing two-pulse pattern is
#' found, the function refits a corresponding single growth curve and returns
#' that fit instead. In that case, overall season timing is still returned, but
#' pulse-specific timing remains \code{NA}.
#'
#' For double-Gompertz and double-Richards fits, the total asymptote can
#' optionally be fixed to the observed maximum cumulative growth of the
#' corresponding series and vegetation season. In the double-curve case, this
#' means \code{a + a2} is fixed, while the relative contribution of the first
#' and second pulse is estimated.
#'
#' @param df A data frame or tibble. The first column must be a time stamp
#' \code{Date}, \code{POSIXct}, or parseable date-time string. Remaining
#' selected columns must be numeric dendrometer series.
#' @param TreeNum Either \code{"all"} to use all dendrometer series, a numeric
#' vector selecting dendrometer columns by position, or a character vector of
#' column names.
#' @param method Double-curve fitting method. One of \code{"gompertz"} or
#' \code{"richards"}.
#' @param years Either \code{"all"} to fit all available vegetation seasons, or
#' a character vector of season labels to retain.
#' @param year_mode Either \code{"yearly"} to fit one curve per vegetation
#' season, or \code{"pooled"} to fit one pooled curve across all retained
#' seasons.
#' @param fit_GRO Logical. If \code{TRUE}, processed daily series are converted
#' to cumulative growth using a cumulative maximum within each vegetation
#' season.
#' @param site_mode Vegetation season definition. One of \code{"NH"},
#' \code{"SH"}, or \code{"CS"}.
#' @param custom_veg_season Numeric vector of length two giving the start and
#' end day-of-year for custom vegetation seasons in \code{"CS"} mode.
#' @param growth_fraction Numeric vector of length two giving the lower and
#' upper fractions of final fitted seasonal growth used to define overall
#' growing-season onset and cessation.
#' @param min_unique_growth Minimum number of unique non-missing cumulative
#' growth values required before a fit is attempted.
#' @param rate_threshold_fraction Numeric scalar between 0 and 1. Pulse start
#' and end are defined as the first and last days where fitted growth rate
#' exceeds this fraction of the pulse-specific peak fitted growth rate.
#' @param peak_separation_min Minimum number of days separating the two fitted
#' derivative peaks required to classify a fit as truly two-pulse.
#' @param valley_ratio_max Maximum allowed ratio between the valley rate and the
#' weaker of the two derivative peaks. Smaller values require a deeper valley
#' between pulses.
#' @param min_relative_peak Minimum relative height, expressed as a fraction of
#' the global derivative maximum, for a local derivative peak to be considered.
#' @param fallback_to_single Logical. If \code{TRUE}, and the fitted double curve
#' does not show a convincing two-pulse pattern according to derivative-based
#' criteria, the function refits a corresponding single growth curve
#' \code{"gompertz"} or \code{"richards"} and returns that fit instead.
#' @param fix_a_to_observed_max Logical. If \code{TRUE}, the total asymptote of
#' the double-growth curve is fixed to the observed maximum cumulative growth
#' for that series and vegetation season. For double-Gompertz and
#' double-Richards models, this means \code{a + a2} is fixed, while the
#' relative contribution of the two pulses is estimated. This is most
#' appropriate with \code{year_mode = "yearly"}. For
#' \code{year_mode = "pooled"}, one pooled maximum is used.
#' @param fixed_a_multiplier Numeric multiplier applied to the observed maximum
#' when \code{fix_a_to_observed_max = TRUE}. The default is \code{1}, meaning
#' the total asymptote is fixed exactly to the observed seasonal maximum.
#' Values slightly above 1, for example \code{1.01} or \code{1.05}, allow the
#' fitted asymptote to sit slightly above the observed maximum.
#' @param start_value_double_gompertz_parameters Optional named list of starting
#' values for the double-Gompertz fit. Supported names are \code{a}, \code{k},
#' \code{t0}, \code{a2}, \code{k2}, and \code{t02}. If
#' \code{fix_a_to_observed_max = TRUE}, supplied \code{a} and \code{a2} are
#' used only to initialize the relative pulse contribution.
#' @param start_value_double_richards_parameters Optional named list of starting
#' values for the double-Richards fit. Supported names are \code{a}, \code{k},
#' \code{t0}, \code{v}, \code{a2}, \code{k2}, \code{t02}, and \code{v2}. If
#' \code{fix_a_to_observed_max = TRUE}, supplied \code{a} and \code{a2} are
#' used only to initialize the relative pulse contribution.
#' @param verbose Logical. If \code{TRUE}, prints a short completion message.
#'
#' @return
#' An object of class \code{"dm_growth_fit"} with elements:
#' \describe{
#' \item{call}{The matched function call.}
#' \item{original_daily_data}{Raw daily dendrometer data after resampling to
#' daily maxima and assigning vegetation seasons, but before centering and
#' optional cumulative-growth transformation.}
#' \item{processed_data}{Daily processed data used for fitting.}
#' \item{fitted_data}{Daily fitted values on the full vegetation-season grid.}
#' \item{fit_statistics}{Fit-level statistics and estimated timing.}
#' \item{fit_parameters}{Fit-level model parameters and convergence
#' information.}
#' \item{season_table}{Vegetation seasons retained for fitting.}
#' }
#'
#' @details
#' The second pulse is constrained to occur after the first pulse, which improves
#' numerical stability and reduces label switching between the two components.
#'
#' For double-Gompertz and double-Richards models, the total seasonal asymptote
#' is:
#'
#' \deqn{a_{total} = a + a2}
#'
#' When \code{fix_a_to_observed_max = TRUE}, the function fits:
#'
#' \deqn{a + a2 = max(y_{obs}) \times fixed\_a\_multiplier}
#'
#' where \eqn{y_{obs}} is the observed cumulative growth of the selected
#' dendrometer series and vegetation season.
#'
#' The fitted pulse contributions are still returned as \code{a} and \code{a2},
#' and their sum should equal \code{fixed_a_value}, apart from small numerical
#' rounding.
#'
#' Non-applicable parameters are returned as \code{NA}. For example, \code{b}
#' and \code{b2} are relevant for double-Gompertz, while \code{v} and \code{v2}
#' are relevant for double-Richards.
#'
#' @examples
#' \donttest{
#' # Double-Gompertz with total seasonal asymptote fixed
#' # fit_gomp <- dm.growth.fit.double(
#' # df = dendro_data,
#' # TreeNum = "all",
#' # method = "gompertz",
#' # year_mode = "yearly",
#' # fit_GRO = TRUE,
#' # fix_a_to_observed_max = TRUE
#' # )
#'
#' # Double-Richards with the total asymptote set 2 percent above observed max
#' # fit_rich <- dm.growth.fit.double(
#' # df = dendro_data,
#' # TreeNum = "all",
#' # method = "richards",
#' # year_mode = "yearly",
#' # fit_GRO = TRUE,
#' # fix_a_to_observed_max = TRUE,
#' # fixed_a_multiplier = 1.02
#' # )
#' }
#'
#' @seealso [dm.growth.fit()], [summary.dm_growth_fit()],
#' [print.dm_growth_fit()]
#'
#' @importFrom dplyr %>% select all_of any_of bind_cols bind_rows filter arrange
#' mutate across group_by ungroup distinct left_join
#' @importFrom tibble as_tibble tibble
#' @importFrom lubridate parse_date_time year yday month
#' @importFrom stats lm coef qlogis predict quantile median na.omit
#' @importFrom utils modifyList
#'
#' @export
dm.growth.fit.double <- function(
df,
TreeNum = "all",
method = c("gompertz", "richards"),
years = "all",
year_mode = c("yearly", "pooled"),
fit_GRO = TRUE,
site_mode = c("NH", "SH", "CS"),
custom_veg_season = c(275, 274),
growth_fraction = c(0.1, 0.9),
min_unique_growth = 10,
rate_threshold_fraction = 0.1,
peak_separation_min = 10,
valley_ratio_max = 0.4,
min_relative_peak = 0.05,
fallback_to_single = TRUE,
fix_a_to_observed_max = FALSE,
fixed_a_multiplier = 1,
start_value_double_gompertz_parameters = list(
a = NA_real_, k = NA_real_, t0 = NA_real_,
a2 = NA_real_, k2 = NA_real_, t02 = NA_real_
),
start_value_double_richards_parameters = list(
a = NA_real_, k = NA_real_, t0 = NA_real_, v = 1,
a2 = NA_real_, k2 = NA_real_, t02 = NA_real_, v2 = 1
),
verbose = TRUE
) {
cl <- match.call()
TIME <- any_observed <- doy <- season_day <- season_end <- NULL
season_label <- season_length <- season_start <- NULL
method <- match.arg(method)
year_mode <- match.arg(year_mode)
site_mode <- match.arg(site_mode)
dmgfd_validate_inputs(
df = df,
growth_fraction = growth_fraction,
min_unique_growth = min_unique_growth,
custom_veg_season = custom_veg_season,
rate_threshold_fraction = rate_threshold_fraction,
peak_separation_min = peak_separation_min,
valley_ratio_max = valley_ratio_max,
min_relative_peak = min_relative_peak,
fallback_to_single = fallback_to_single,
fix_a_to_observed_max = fix_a_to_observed_max,
fixed_a_multiplier = fixed_a_multiplier
)
df <- tibble::as_tibble(df)
names(df)[1] <- "TIME"
if (!inherits(df$TIME, "Date") &&
!inherits(df$TIME, "POSIXct") &&
!inherits(df$TIME, "POSIXt")) {
df$TIME <- suppressWarnings(
lubridate::parse_date_time(
df$TIME,
orders = c(
"ymd HMS", "ymd HM", "ymd",
"dmy HMS", "dmy HM", "dmy",
"mdy HMS", "mdy HM", "mdy"
),
quiet = TRUE
)
)
}
if (inherits(df$TIME, "Date")) {
df$TIME <- as.POSIXct(df$TIME)
}
if (any(is.na(df$TIME))) {
stop("Some timestamps in the first column could not be parsed.")
}
obs_years <- sort(unique(lubridate::year(as.Date(df$TIME))))
if (site_mode == "SH" && length(obs_years) < 2) {
stop("For site_mode = 'SH', df must contain observations from at least 2 calendar years.")
}
if (site_mode == "CS" &&
custom_veg_season[1] > custom_veg_season[2] &&
length(obs_years) < 2) {
stop(
"For site_mode = 'CS' with custom_veg_season[1] > custom_veg_season[2], ",
"df must contain observations from at least 2 calendar years."
)
}
all_series <- names(df)[-1]
if (length(all_series) == 0) {
stop("df must contain a time column followed by at least one dendrometer column.")
}
if (is.character(TreeNum) &&
length(TreeNum) == 1 &&
tolower(TreeNum) == "all") {
keep_cols <- all_series
} else if (is.numeric(TreeNum)) {
tree_idx <- as.integer(TreeNum)
if (any(is.na(tree_idx)) ||
any(tree_idx < 1) ||
any(tree_idx > length(all_series))) {
stop(
"'TreeNum' numeric values must be between 1 and ", length(all_series),
". Here, 1 refers to the first dendrometer series after the TIME column."
)
}
keep_cols <- all_series[tree_idx]
} else if (is.character(TreeNum)) {
missing_cols <- setdiff(TreeNum, all_series)
if (length(missing_cols) > 0) {
warning(
"These TreeNum column names were not found and were ignored: ",
paste(missing_cols, collapse = ", ")
)
}
keep_cols <- intersect(TreeNum, all_series)
if (length(keep_cols) == 0) {
stop("None of the requested TreeNum column names were found in df.")
}
} else {
stop("TreeNum must be 'all', a numeric vector, or a character vector of column names.")
}
df <- dplyr::select(df, TIME, dplyr::all_of(keep_cols))
num_cols <- names(df)[vapply(df, is.numeric, logical(1))]
num_cols <- setdiff(num_cols, "TIME")
if (length(num_cols) == 0) {
stop("No numeric dendrometer series found to fit.")
}
df <- dplyr::select(df, TIME, dplyr::all_of(num_cols))
dm_daily <- dendro.resample(
df = df,
by = "D",
value = "max",
method = "aggregate"
)
season_tbl_obs <- dmgfd_build_season_table(
time_vec = dm_daily$TIME,
site_mode = site_mode,
custom_veg_season = custom_veg_season
)
observed_daily <- dplyr::bind_cols(dm_daily, season_tbl_obs) %>%
dplyr::filter(!is.na(season_label)) %>%
dplyr::arrange(TIME)
original_daily_data <- observed_daily %>%
dplyr::select(
TIME, doy, season_label, season_start, season_end, season_day,
dplyr::all_of(num_cols)
)
season_table <- observed_daily %>%
dplyr::distinct(season_label, season_start, season_end) %>%
dplyr::arrange(season_start)
if (!(is.character(years) &&
length(years) == 1 &&
tolower(years) == "all")) {
season_table <- season_table %>%
dplyr::filter(season_label %in% as.character(years))
}
if (nrow(season_table) == 0) {
stop("No seasons remained after applying site_mode/custom_veg_season/years.")
}
full_grids <- lapply(seq_len(nrow(season_table)), function(i) {
ss <- season_table$season_start[i]
se <- season_table$season_end[i]
sl <- season_table$season_label[i]
tt <- seq(ss, se, by = "day")
tibble::tibble(
TIME = tt,
doy = lubridate::yday(tt),
season_label = sl,
season_start = ss,
season_end = se,
season_day = seq_along(tt),
season_length = length(tt)
)
})
full_grid <- dplyr::bind_rows(full_grids)
dat <- full_grid %>%
dplyr::left_join(
observed_daily %>%
dplyr::select(TIME, season_label, dplyr::all_of(num_cols)),
by = c("TIME", "season_label")
) %>%
dplyr::arrange(season_start, TIME)
dat_proc <- dat %>%
dplyr::group_by(season_label) %>%
dplyr::arrange(TIME, .by_group = TRUE) %>%
dplyr::mutate(
dplyr::across(
dplyr::all_of(num_cols),
~ . - dmgfd_first_non_na(.)
)
) %>%
dplyr::ungroup()
if (isTRUE(fit_GRO)) {
dat_proc <- dat_proc %>%
dplyr::group_by(season_label) %>%
dplyr::arrange(TIME, .by_group = TRUE) %>%
dplyr::mutate(
dplyr::across(
dplyr::all_of(num_cols),
dmgfd_cummax_na
)
) %>%
dplyr::ungroup()
}
dat_proc$any_observed <- rowSums(
!is.na(as.data.frame(dat_proc[, num_cols, drop = FALSE]))
) > 0
fitted_dat <- dat_proc %>%
dplyr::select(
TIME, doy, season_label, season_start, season_end,
season_day, season_length, any_observed
)
for (cc in num_cols) {
fitted_dat[[cc]] <- NA_real_
}
param_rows <- list()
row_id <- 1L
if (year_mode == "yearly") {
fit_ids <- unique(dat_proc$season_label)
for (series_name in num_cols) {
for (sid in fit_ids) {
idx <- which(dat_proc$season_label == sid)
x_full <- dat_proc$season_day[idx]
y_full <- dat_proc[[series_name]][idx]
ss <- unique(dat_proc$season_start[idx])[1]
se <- unique(dat_proc$season_end[idx])[1]
full_len <- unique(dat_proc$season_length[idx])[1]
good <- is.finite(y_full)
x_obs <- x_full[good]
y_obs <- y_full[good]
fit_res <- dmgfd_fit_one_curve(
x_obs = x_obs,
y_obs = y_obs,
x_pred = x_full,
season_start_date = ss,
season_length = full_len,
method = method,
growth_fraction = growth_fraction,
min_unique_growth = min_unique_growth,
rate_threshold_fraction = rate_threshold_fraction,
peak_separation_min = peak_separation_min,
valley_ratio_max = valley_ratio_max,
min_relative_peak = min_relative_peak,
fallback_to_single = fallback_to_single,
fix_a_to_observed_max = fix_a_to_observed_max,
fixed_a_multiplier = fixed_a_multiplier,
start_value_double_gompertz_parameters = start_value_double_gompertz_parameters,
start_value_double_richards_parameters = start_value_double_richards_parameters
)
fitted_dat[[series_name]][idx] <- fit_res$pred
param_rows[[row_id]] <- tibble::tibble(
series = series_name,
fit_id = sid,
year_mode = year_mode,
method = fit_res$params$method_used,
site_mode = site_mode,
season_start = ss,
season_end = se,
included_seasons = NA_character_,
season_length = fit_res$params$season_length,
n_obs = fit_res$params$n_obs,
n_days_observed = fit_res$params$n_days_observed,
first_obs_day = fit_res$params$first_obs_day,
last_obs_day = fit_res$params$last_obs_day,
n_missing_days = fit_res$params$n_missing_days,
extrapolated = fit_res$params$extrapolated,
anchor_added = fit_res$params$anchor_added,
growth_start_day = fit_res$params$growth_start_day,
growth_end_day = fit_res$params$growth_end_day,
growth_start_season_day = fit_res$params$growth_start_season_day,
growth_end_season_day = fit_res$params$growth_end_season_day,
growth_start_date = fit_res$params$growth_start_date,
growth_end_date = fit_res$params$growth_end_date,
peak_rate = fit_res$params$peak_rate,
rate_start_day = fit_res$params$rate_start_day,
rate_end_day = fit_res$params$rate_end_day,
rate_start_season_day = fit_res$params$rate_start_season_day,
rate_end_season_day = fit_res$params$rate_end_season_day,
rate_start_date = fit_res$params$rate_start_date,
rate_end_date = fit_res$params$rate_end_date,
fallback_used = fit_res$params$fallback_used,
two_pulse_detected = fit_res$params$two_pulse_detected,
peak1_day = fit_res$params$peak1_day,
peak2_day = fit_res$params$peak2_day,
peak1_season_day = fit_res$params$peak1_season_day,
peak2_season_day = fit_res$params$peak2_season_day,
separator_day = fit_res$params$separator_day,
separator_season_day = fit_res$params$separator_season_day,
separator_date = fit_res$params$separator_date,
valley_rate = fit_res$params$valley_rate,
pulse1_peak_rate = fit_res$params$pulse1_peak_rate,
pulse2_peak_rate = fit_res$params$pulse2_peak_rate,
pulse1_start_day = fit_res$params$pulse1_start_day,
pulse1_end_day = fit_res$params$pulse1_end_day,
pulse1_start_season_day = fit_res$params$pulse1_start_season_day,
pulse1_end_season_day = fit_res$params$pulse1_end_season_day,
pulse1_start_date = fit_res$params$pulse1_start_date,
pulse1_end_date = fit_res$params$pulse1_end_date,
pulse2_start_day = fit_res$params$pulse2_start_day,
pulse2_end_day = fit_res$params$pulse2_end_day,
pulse2_start_season_day = fit_res$params$pulse2_start_season_day,
pulse2_end_season_day = fit_res$params$pulse2_end_season_day,
pulse2_start_date = fit_res$params$pulse2_start_date,
pulse2_end_date = fit_res$params$pulse2_end_date,
converged = fit_res$params$converged,
fixed_a_used = fit_res$params$fixed_a_used,
fixed_a_value = fit_res$params$fixed_a_value,
a = fit_res$params$a,
b = fit_res$params$b,
k = fit_res$params$k,
t0 = fit_res$params$t0,
v = fit_res$params$v,
a2 = fit_res$params$a2,
b2 = fit_res$params$b2,
k2 = fit_res$params$k2,
t02 = fit_res$params$t02,
v2 = fit_res$params$v2,
edf = fit_res$params$edf,
span = fit_res$params$span,
spline_df = fit_res$params$spline_df,
spar = fit_res$params$spar,
model_note = fit_res$params$model_note
)
row_id <- row_id + 1L
}
}
}
if (year_mode == "pooled") {
included_seasons <- paste(unique(dat_proc$season_label), collapse = ", ")
for (series_name in num_cols) {
x_obs_all <- dat_proc$season_day
y_obs_all <- dat_proc[[series_name]]
good <- is.finite(y_obs_all)
x_obs <- x_obs_all[good]
y_obs <- y_obs_all[good]
x_template <- seq_len(max(dat_proc$season_day, na.rm = TRUE))
fit_res <- dmgfd_fit_one_curve(
x_obs = x_obs,
y_obs = y_obs,
x_pred = x_template,
season_start_date = as.Date(NA),
season_length = max(x_template, na.rm = TRUE),
method = method,
growth_fraction = growth_fraction,
min_unique_growth = min_unique_growth,
rate_threshold_fraction = rate_threshold_fraction,
peak_separation_min = peak_separation_min,
valley_ratio_max = valley_ratio_max,
min_relative_peak = min_relative_peak,
fallback_to_single = fallback_to_single,
fix_a_to_observed_max = fix_a_to_observed_max,
fixed_a_multiplier = fixed_a_multiplier,
start_value_double_gompertz_parameters = start_value_double_gompertz_parameters,
start_value_double_richards_parameters = start_value_double_richards_parameters
)
fitted_dat[[series_name]] <- fit_res$pred[match(dat_proc$season_day, x_template)]
param_rows[[row_id]] <- tibble::tibble(
series = series_name,
fit_id = "pooled",
year_mode = year_mode,
method = fit_res$params$method_used,
site_mode = site_mode,
season_start = as.Date(NA),
season_end = as.Date(NA),
included_seasons = included_seasons,
season_length = fit_res$params$season_length,
n_obs = fit_res$params$n_obs,
n_days_observed = fit_res$params$n_days_observed,
first_obs_day = fit_res$params$first_obs_day,
last_obs_day = fit_res$params$last_obs_day,
n_missing_days = fit_res$params$n_missing_days,
extrapolated = fit_res$params$extrapolated,
anchor_added = fit_res$params$anchor_added,
growth_start_day = fit_res$params$growth_start_day,
growth_end_day = fit_res$params$growth_end_day,
growth_start_season_day = fit_res$params$growth_start_season_day,
growth_end_season_day = fit_res$params$growth_end_season_day,
growth_start_date = fit_res$params$growth_start_date,
growth_end_date = fit_res$params$growth_end_date,
peak_rate = fit_res$params$peak_rate,
rate_start_day = fit_res$params$rate_start_day,
rate_end_day = fit_res$params$rate_end_day,
rate_start_season_day = fit_res$params$rate_start_season_day,
rate_end_season_day = fit_res$params$rate_end_season_day,
rate_start_date = fit_res$params$rate_start_date,
rate_end_date = fit_res$params$rate_end_date,
fallback_used = fit_res$params$fallback_used,
two_pulse_detected = fit_res$params$two_pulse_detected,
peak1_day = fit_res$params$peak1_day,
peak2_day = fit_res$params$peak2_day,
peak1_season_day = fit_res$params$peak1_season_day,
peak2_season_day = fit_res$params$peak2_season_day,
separator_day = fit_res$params$separator_day,
separator_season_day = fit_res$params$separator_season_day,
separator_date = fit_res$params$separator_date,
valley_rate = fit_res$params$valley_rate,
pulse1_peak_rate = fit_res$params$pulse1_peak_rate,
pulse2_peak_rate = fit_res$params$pulse2_peak_rate,
pulse1_start_day = fit_res$params$pulse1_start_day,
pulse1_end_day = fit_res$params$pulse1_end_day,
pulse1_start_season_day = fit_res$params$pulse1_start_season_day,
pulse1_end_season_day = fit_res$params$pulse1_end_season_day,
pulse1_start_date = fit_res$params$pulse1_start_date,
pulse1_end_date = fit_res$params$pulse1_end_date,
pulse2_start_day = fit_res$params$pulse2_start_day,
pulse2_end_day = fit_res$params$pulse2_end_day,
pulse2_start_season_day = fit_res$params$pulse2_start_season_day,
pulse2_end_season_day = fit_res$params$pulse2_end_season_day,
pulse2_start_date = fit_res$params$pulse2_start_date,
pulse2_end_date = fit_res$params$pulse2_end_date,
converged = fit_res$params$converged,
fixed_a_used = fit_res$params$fixed_a_used,
fixed_a_value = fit_res$params$fixed_a_value,
a = fit_res$params$a,
b = fit_res$params$b,
k = fit_res$params$k,
t0 = fit_res$params$t0,
v = fit_res$params$v,
a2 = fit_res$params$a2,
b2 = fit_res$params$b2,
k2 = fit_res$params$k2,
t02 = fit_res$params$t02,
v2 = fit_res$params$v2,
edf = fit_res$params$edf,
span = fit_res$params$span,
spline_df = fit_res$params$spline_df,
spar = fit_res$params$spar,
model_note = fit_res$params$model_note
)
row_id <- row_id + 1L
}
}
parameter_table <- dplyr::bind_rows(param_rows)
stats_cols <- c(
"series", "fit_id", "year_mode", "method", "site_mode",
"season_start", "season_end", "included_seasons", "season_length",
"n_obs", "n_days_observed", "first_obs_day", "last_obs_day",
"n_missing_days", "extrapolated", "anchor_added",
"growth_start_day", "growth_end_day",
"growth_start_season_day", "growth_end_season_day",
"growth_start_date", "growth_end_date",
"peak_rate",
"rate_start_day", "rate_end_day",
"rate_start_season_day", "rate_end_season_day",
"rate_start_date", "rate_end_date",
"fallback_used",
"two_pulse_detected",
"peak1_day", "peak2_day",
"peak1_season_day", "peak2_season_day",
"separator_day", "separator_season_day", "separator_date",
"valley_rate",
"pulse1_peak_rate", "pulse2_peak_rate",
"pulse1_start_day", "pulse1_end_day",
"pulse1_start_season_day", "pulse1_end_season_day",
"pulse1_start_date", "pulse1_end_date",
"pulse2_start_day", "pulse2_end_day",
"pulse2_start_season_day", "pulse2_end_season_day",
"pulse2_start_date", "pulse2_end_date"
)
param_cols <- c(
"series", "fit_id", "year_mode", "method", "site_mode",
"fallback_used", "converged",
"fixed_a_used", "fixed_a_value",
"a", "b", "k", "t0", "v",
"a2", "b2", "k2", "t02", "v2",
"edf", "span", "spline_df", "spar", "model_note"
)
fit_statistics <- parameter_table %>%
dplyr::select(dplyr::any_of(stats_cols))
fit_parameters <- parameter_table %>%
dplyr::select(dplyr::any_of(param_cols))
out <- list(
call = cl,
original_daily_data = original_daily_data,
processed_data = dat_proc,
fitted_data = fitted_dat,
fit_statistics = fit_statistics,
fit_parameters = fit_parameters,
season_table = season_table
)
class(out) <- "dm_growth_fit"
if (isTRUE(verbose)) {
message(
"dm.growth.fit.double completed: ",
nrow(fit_parameters), " fit(s) across ",
length(unique(fit_parameters$series)), " series."
)
}
out
}
# helpers ----------------------------------------------------------------------
#' Validate inputs for double dendrometer growth fitting
#'
#' @param df Input data frame.
#' @param growth_fraction Numeric vector of two growth fractions.
#' @param min_unique_growth Minimum number of unique growth values.
#' @param custom_veg_season Custom vegetation-season DOY range.
#' @param rate_threshold_fraction Fraction of peak growth rate.
#' @param peak_separation_min Minimum separation between derivative peaks.
#' @param valley_ratio_max Maximum allowed valley-to-peak ratio.
#' @param min_relative_peak Minimum relative derivative peak height.
#' @param fallback_to_single Logical fallback flag.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param fixed_a_multiplier Multiplier for fixed total asymptote.
#'
#' @return Invisibly returns \code{TRUE} if checks pass.
#'
#' @keywords internal
dmgfd_validate_inputs <- function(df,
growth_fraction,
min_unique_growth,
custom_veg_season,
rate_threshold_fraction,
peak_separation_min,
valley_ratio_max,
min_relative_peak,
fallback_to_single,
fix_a_to_observed_max,
fixed_a_multiplier) {
if (!is.data.frame(df)) {
stop("df must be a data.frame or tibble.")
}
if (ncol(df) < 2) {
stop("df must contain a time column and at least one dendrometer series.")
}
if (!is.numeric(growth_fraction) ||
length(growth_fraction) != 2 ||
any(!is.finite(growth_fraction)) ||
growth_fraction[1] < 0 ||
growth_fraction[2] > 1 ||
growth_fraction[1] >= growth_fraction[2]) {
stop("growth_fraction must be a numeric vector like c(0.1, 0.9).")
}
if (!is.numeric(min_unique_growth) ||
length(min_unique_growth) != 1 ||
!is.finite(min_unique_growth) ||
min_unique_growth < 2) {
stop("min_unique_growth must be a single finite number >= 2.")
}
if (!is.numeric(custom_veg_season) ||
length(custom_veg_season) != 2 ||
any(!is.finite(custom_veg_season))) {
stop("custom_veg_season must be a numeric vector of length 2.")
}
if (!is.numeric(rate_threshold_fraction) ||
length(rate_threshold_fraction) != 1 ||
!is.finite(rate_threshold_fraction) ||
rate_threshold_fraction <= 0 ||
rate_threshold_fraction >= 1) {
stop("rate_threshold_fraction must be a single number between 0 and 1.")
}
if (!is.numeric(peak_separation_min) ||
length(peak_separation_min) != 1 ||
!is.finite(peak_separation_min) ||
peak_separation_min < 1) {
stop("peak_separation_min must be a single number >= 1.")
}
if (!is.numeric(valley_ratio_max) ||
length(valley_ratio_max) != 1 ||
!is.finite(valley_ratio_max) ||
valley_ratio_max <= 0 ||
valley_ratio_max >= 1) {
stop("valley_ratio_max must be a single number between 0 and 1.")
}
if (!is.numeric(min_relative_peak) ||
length(min_relative_peak) != 1 ||
!is.finite(min_relative_peak) ||
min_relative_peak <= 0 ||
min_relative_peak >= 1) {
stop("min_relative_peak must be a single number between 0 and 1.")
}
if (!is.logical(fallback_to_single) ||
length(fallback_to_single) != 1 ||
is.na(fallback_to_single)) {
stop("fallback_to_single must be TRUE or FALSE.")
}
if (!is.logical(fix_a_to_observed_max) ||
length(fix_a_to_observed_max) != 1 ||
is.na(fix_a_to_observed_max)) {
stop("fix_a_to_observed_max must be TRUE or FALSE.")
}
if (!is.numeric(fixed_a_multiplier) ||
length(fixed_a_multiplier) != 1 ||
!is.finite(fixed_a_multiplier) ||
fixed_a_multiplier <= 0) {
stop("fixed_a_multiplier must be a single positive finite number.")
}
invisible(TRUE)
}
#' Check whether a value is a finite numeric scalar
#'
#' @param x Object to test.
#'
#' @return Logical scalar.
#'
#' @keywords internal
dmgfd_is_scalar_finite <- function(x) {
is.numeric(x) && length(x) == 1 && !is.na(x) && is.finite(x)
}
#' Return first non-missing value
#'
#' @param x Numeric vector.
#'
#' @return First non-missing value or \code{NA_real_}.
#'
#' @keywords internal
dmgfd_first_non_na <- function(x) {
idx <- which(!is.na(x))[1]
if (length(idx) == 0 || is.na(idx)) {
return(NA_real_)
}
x[idx]
}
#' Cumulative maximum while preserving missing values
#'
#' @param x Numeric vector.
#'
#' @return Numeric vector.
#'
#' @keywords internal
dmgfd_cummax_na <- function(x) {
out <- x
idx <- which(!is.na(x))
if (length(idx) > 0) {
out[idx] <- cummax(x[idx])
}
out
}
#' Convert year and day-of-year to Date
#'
#' @param year Calendar year.
#' @param doy Day of year.
#'
#' @return Date.
#'
#' @keywords internal
dmgfd_doy_to_date <- function(year, doy) {
as.Date(sprintf("%04d-01-01", year)) + (doy - 1)
}
#' Build vegetation-season table
#'
#' @param time_vec Time vector.
#' @param site_mode One of \code{"NH"}, \code{"SH"}, or \code{"CS"}.
#' @param custom_veg_season Custom vegetation-season DOY range.
#'
#' @return Tibble with season labels and season timing.
#'
#' @keywords internal
dmgfd_build_season_table <- function(time_vec,
site_mode = c("NH", "SH", "CS"),
custom_veg_season = c(275, 274)) {
site_mode <- match.arg(site_mode)
dates <- as.Date(time_vec)
yr <- lubridate::year(dates)
doy <- lubridate::yday(dates)
mo <- lubridate::month(dates)
n <- length(dates)
season_label <- rep(NA_character_, n)
season_start <- as.Date(rep(NA_character_, n))
season_end <- as.Date(rep(NA_character_, n))
season_day <- rep(NA_integer_, n)
if (site_mode == "NH") {
season_start <- as.Date(sprintf("%04d-01-01", yr))
season_end <- as.Date(sprintf("%04d-12-31", yr))
season_label <- as.character(yr)
season_day <- as.integer(dates - season_start) + 1L
}
if (site_mode == "SH") {
start_year <- ifelse(mo >= 7, yr, yr - 1L)
end_year <- start_year + 1L
season_start <- as.Date(sprintf("%04d-07-01", start_year))
season_end <- as.Date(sprintf("%04d-06-30", end_year))
season_label <- sprintf("%04d/%04d", start_year, end_year)
season_day <- as.integer(dates - season_start) + 1L
}
if (site_mode == "CS") {
start_doy <- as.integer(custom_veg_season[1])
end_doy <- as.integer(custom_veg_season[2])
if (start_doy <= end_doy) {
inside <- doy >= start_doy & doy <= end_doy
season_start[inside] <- dmgfd_doy_to_date(yr[inside], start_doy)
season_end[inside] <- dmgfd_doy_to_date(yr[inside], end_doy)
season_label[inside] <- as.character(yr[inside])
season_day[inside] <- as.integer(dates[inside] - season_start[inside]) + 1L
} else {
inside_late <- doy >= start_doy
inside_early <- doy <= end_doy
if (any(inside_late)) {
sy <- yr[inside_late]
ey <- sy + 1L
season_start[inside_late] <- dmgfd_doy_to_date(sy, start_doy)
season_end[inside_late] <- dmgfd_doy_to_date(ey, end_doy)
season_label[inside_late] <- sprintf("%04d/%04d", sy, ey)
season_day[inside_late] <- as.integer(dates[inside_late] - season_start[inside_late]) + 1L
}
if (any(inside_early)) {
ey <- yr[inside_early]
sy <- ey - 1L
season_start[inside_early] <- dmgfd_doy_to_date(sy, start_doy)
season_end[inside_early] <- dmgfd_doy_to_date(ey, end_doy)
season_label[inside_early] <- sprintf("%04d/%04d", sy, ey)
season_day[inside_early] <- as.integer(dates[inside_early] - season_start[inside_early]) + 1L
}
}
}
tibble::tibble(
season_label = season_label,
season_start = season_start,
season_end = season_end,
season_day = season_day
)
}
#' Add zero-growth anchor point
#'
#' @param x Season-day values.
#' @param y Growth values.
#'
#' @return List with \code{x} and \code{y}.
#'
#' @keywords internal
dmgfd_add_anchor_points <- function(x, y) {
if (length(x) == 0) {
return(list(x = x, y = y))
}
if (min(x, na.rm = TRUE) > 1 && !any(x == 1)) {
x <- c(1, x)
y <- c(0, y)
}
ord <- order(x)
list(
x = x[ord],
y = y[ord]
)
}
#' Empty parameter template for double growth fitting
#'
#' @return Named list of model and timing parameters.
#'
#' @keywords internal
dmgfd_empty_model_params <- function() {
list(
converged = FALSE,
fallback_used = FALSE,
fixed_a_used = FALSE,
fixed_a_value = NA_real_,
method_used = NA_character_,
a = NA_real_,
b = NA_real_,
k = NA_real_,
t0 = NA_real_,
v = NA_real_,
a2 = NA_real_,
b2 = NA_real_,
k2 = NA_real_,
t02 = NA_real_,
v2 = NA_real_,
edf = NA_real_,
span = NA_real_,
spline_df = NA_real_,
spar = NA_real_,
model_note = NA_character_,
n_obs = NA_real_,
n_days_observed = NA_real_,
first_obs_day = NA_real_,
last_obs_day = NA_real_,
season_length = NA_real_,
n_missing_days = NA_real_,
extrapolated = NA,
anchor_added = FALSE,
growth_start_day = NA_real_,
growth_end_day = NA_real_,
growth_start_season_day = NA_real_,
growth_end_season_day = NA_real_,
growth_start_date = as.Date(NA),
growth_end_date = as.Date(NA),
peak_rate = NA_real_,
rate_start_day = NA_real_,
rate_end_day = NA_real_,
rate_start_season_day = NA_real_,
rate_end_season_day = NA_real_,
rate_start_date = as.Date(NA),
rate_end_date = as.Date(NA),
two_pulse_detected = FALSE,
peak1_day = NA_real_,
peak2_day = NA_real_,
peak1_season_day = NA_real_,
peak2_season_day = NA_real_,
separator_day = NA_real_,
separator_season_day = NA_real_,
separator_date = as.Date(NA),
valley_rate = NA_real_,
pulse1_peak_rate = NA_real_,
pulse2_peak_rate = NA_real_,
pulse1_start_day = NA_real_,
pulse1_end_day = NA_real_,
pulse1_start_season_day = NA_real_,
pulse1_end_season_day = NA_real_,
pulse1_start_date = as.Date(NA),
pulse1_end_date = as.Date(NA),
pulse2_start_day = NA_real_,
pulse2_end_day = NA_real_,
pulse2_start_season_day = NA_real_,
pulse2_end_season_day = NA_real_,
pulse2_start_date = as.Date(NA),
pulse2_end_date = as.Date(NA)
)
}
#' Guess two pulse locations from growth increments
#'
#' @param x Season-day values.
#' @param y Cumulative growth values.
#'
#' @return Numeric vector of two initial pulse locations.
#'
#' @keywords internal
dmgfd_guess_two_pulses <- function(x, y) {
x <- as.numeric(x)
y <- as.numeric(y)
ord <- order(x)
x <- x[ord]
y <- y[ord]
dy <- c(NA_real_, diff(y))
dy[!is.finite(dy)] <- -Inf
min_sep <- max(14, round(0.15 * diff(range(x, na.rm = TRUE))))
o <- order(dy, decreasing = TRUE)
if (length(o) < 2 || !is.finite(dy[o[1]])) {
qs <- stats::quantile(x, probs = c(0.33, 0.67), na.rm = TRUE, names = FALSE)
return(as.numeric(qs))
}
t1 <- x[o[1]]
t2 <- NA_real_
if (length(o) > 1) {
for (j in o[-1]) {
if (is.finite(dy[j]) && abs(x[j] - t1) >= min_sep) {
t2 <- x[j]
break
}
}
}
if (!is.finite(t2)) {
qs <- stats::quantile(x, probs = c(0.33, 0.67), na.rm = TRUE, names = FALSE)
t1 <- qs[1]
t2 <- qs[2]
}
sort(c(t1, t2))
}
#' Infer single Gompertz starting values
#'
#' @param x Season-day values.
#' @param y Growth values.
#' @param a0 Initial asymptote.
#'
#' @return Named list with \code{a}, \code{b}, and \code{k}.
#'
#' @keywords internal
dmgfd_infer_gompertz_starts <- function(x, y, a0) {
frac <- pmin(pmax(y / a0, 1e-6), 1 - 1e-6)
z <- log(-log(frac))
ok <- is.finite(x) & is.finite(z)
b0 <- 0.5
k0 <- 0.01
if (sum(ok) >= 2) {
lm_fit <- try(stats::lm(z[ok] ~ x[ok]), silent = TRUE)
if (!inherits(lm_fit, "try-error")) {
cf <- stats::coef(lm_fit)
if (length(cf) == 2 && all(is.finite(cf))) {
b_try <- unname(cf[1])
k_try <- -unname(cf[2])
if (is.finite(b_try)) {
b0 <- b_try
}
if (is.finite(k_try) && k_try > 1e-4) {
k0 <- k_try
}
}
}
}
list(
a = a0,
b = b0,
k = k0
)
}
#' Infer logistic/Richards starting values
#'
#' @param x Season-day values.
#' @param y Growth values.
#' @param a0 Initial asymptote.
#'
#' @return Named list with \code{a}, \code{k}, and \code{t0}.
#'
#' @keywords internal
dmgfd_infer_logistic_starts <- function(x, y, a0) {
frac <- pmin(pmax(y / a0, 1e-6), 1 - 1e-6)
z <- stats::qlogis(frac)
ok <- is.finite(x) & is.finite(z)
k0 <- 0.03
t00 <- stats::median(x, na.rm = TRUE)
if (sum(ok) >= 2) {
lm_fit <- try(stats::lm(z[ok] ~ x[ok]), silent = TRUE)
if (!inherits(lm_fit, "try-error")) {
cf <- stats::coef(lm_fit)
if (length(cf) == 2 && all(is.finite(cf))) {
k_try <- abs(unname(cf[2]))
if (is.finite(k_try) && k_try > 1e-4) {
k0 <- k_try
t0_try <- -unname(cf[1]) / k0
if (is.finite(t0_try)) {
t00 <- t0_try
}
}
}
}
}
list(
a = a0,
k = k0,
t0 = t00
)
}
#' Infer double-Gompertz starting values
#'
#' @param x Season-day values.
#' @param y Growth values.
#' @param a0 Initial total asymptote.
#'
#' @return Named list of starting values.
#'
#' @keywords internal
dmgfd_infer_double_gompertz_starts <- function(x, y, a0) {
pulses <- dmgfd_guess_two_pulses(x, y)
base <- dmgfd_infer_gompertz_starts(x, y, a0)
k0 <- if (dmgfd_is_scalar_finite(base$k)) base$k else 0.01
list(
a = a0 * 0.5,
b = k0 * pulses[1],
k = k0,
t0 = pulses[1],
a2 = a0 * 0.5,
b2 = k0 * pulses[2],
k2 = k0,
t02 = pulses[2]
)
}
#' Infer double-Richards starting values
#'
#' @param x Season-day values.
#' @param y Growth values.
#' @param a0 Initial total asymptote.
#'
#' @return Named list of starting values.
#'
#' @keywords internal
dmgfd_infer_double_richards_starts <- function(x, y, a0) {
pulses <- dmgfd_guess_two_pulses(x, y)
base <- dmgfd_infer_logistic_starts(x, y, a0)
k0 <- if (dmgfd_is_scalar_finite(base$k)) base$k else 0.03
list(
a = a0 * 0.5,
k = k0,
t0 = pulses[1],
v = 1,
a2 = a0 * 0.5,
k2 = k0,
t02 = pulses[2],
v2 = 1
)
}
#' Fit double-Gompertz model
#'
#' @param x Observed season-day values.
#' @param y Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param fixed_a_value Optional fixed total asymptote on original scale.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param start_value_double_gompertz_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgfd_fit_model_double_gompertz <- function(
x,
y,
x_pred,
fixed_a_value = NA_real_,
fix_a_to_observed_max = FALSE,
start_value_double_gompertz_parameters = list(
a = NA_real_, k = NA_real_, t0 = NA_real_,
a2 = NA_real_, k2 = NA_real_, t02 = NA_real_
)
) {
params <- dmgfd_empty_model_params()
res <- tryCatch({
cons <- ifelse(
min(y, na.rm = TRUE) < 0,
abs(min(y, na.rm = TRUE)) + 1e-8,
0
)
y_adj <- y + cons
a0 <- max(y_adj, na.rm = TRUE) * 1.05 + 1e-6
start_guess <- dmgfd_infer_double_gompertz_starts(x, y_adj, a0)
start_guess <- utils::modifyList(start_guess, start_value_double_gompertz_parameters)
if (!dmgfd_is_scalar_finite(start_guess$a)) {
start_guess$a <- a0 * 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k)) {
start_guess$k <- 0.01
}
if (!dmgfd_is_scalar_finite(start_guess$t0)) {
start_guess$t0 <- stats::quantile(x, 0.33, na.rm = TRUE, names = FALSE)
}
if (!dmgfd_is_scalar_finite(start_guess$a2)) {
start_guess$a2 <- a0 * 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k2)) {
start_guess$k2 <- start_guess$k
}
if (!dmgfd_is_scalar_finite(start_guess$t02)) {
start_guess$t02 <- stats::quantile(x, 0.67, na.rm = TRUE, names = FALSE)
}
dt_start <- max(5, start_guess$t02 - start_guess$t0)
dat <- data.frame(
x = x,
y = y_adj
)
use_fixed_a <- isTRUE(fix_a_to_observed_max) &&
dmgfd_is_scalar_finite(fixed_a_value) &&
fixed_a_value > 0
if (isTRUE(use_fixed_a)) {
a_fixed_original <- fixed_a_value
a_fixed_fit <- fixed_a_value + cons
p_start <- 0.5
if (dmgfd_is_scalar_finite(start_guess$a) &&
dmgfd_is_scalar_finite(start_guess$a2) &&
(start_guess$a + start_guess$a2) > 0) {
p_start <- start_guess$a / (start_guess$a + start_guess$a2)
}
p_start <- min(max(p_start, 0.05), 0.95)
mod <- minpack.lm::nlsLM(
y ~ a_fixed_fit * p * exp(-exp(-(k * (x - t0)))) +
a_fixed_fit * (1 - p) * exp(-exp(-(k2 * (x - (t0 + dt))))),
data = dat,
start = list(
p = p_start,
k = start_guess$k,
t0 = start_guess$t0,
k2 = start_guess$k2,
dt = dt_start
),
lower = c(
p = 0.001,
k = 1e-6,
t0 = min(x, na.rm = TRUE),
k2 = 1e-6,
dt = 1
),
upper = c(
p = 0.999,
k = 5,
t0 = max(x, na.rm = TRUE),
k2 = 5,
dt = diff(range(x, na.rm = TRUE))
),
control = minpack.lm::nls.lm.control(maxiter = 1000)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- pmax(as.numeric(pred), 0)
pred <- cummax(pred)
cf <- stats::coef(mod)
p_est <- unname(cf["p"])
params$a <- a_fixed_original * p_est
params$a2 <- a_fixed_original * (1 - p_est)
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$b <- params$k * params$t0
params$k2 <- unname(cf["k2"])
params$t02 <- unname(cf["t0"] + cf["dt"])
params$b2 <- params$k2 * params$t02
params$converged <- TRUE
params$fixed_a_used <- TRUE
params$fixed_a_value <- a_fixed_original
params$method_used <- "double_gompertz"
params$model_note <- "Double-Gompertz total asymptote a + a2 fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
mod <- minpack.lm::nlsLM(
y ~ a * exp(-exp(-(k * (x - t0)))) +
a2 * exp(-exp(-(k2 * (x - (t0 + dt))))),
data = dat,
start = list(
a = start_guess$a,
k = start_guess$k,
t0 = start_guess$t0,
a2 = start_guess$a2,
k2 = start_guess$k2,
dt = dt_start
),
lower = c(
a = 0,
k = 1e-6,
t0 = min(x, na.rm = TRUE),
a2 = 0,
k2 = 1e-6,
dt = 1
),
upper = c(
a = Inf,
k = 5,
t0 = max(x, na.rm = TRUE),
a2 = Inf,
k2 = 5,
dt = diff(range(x, na.rm = TRUE))
),
control = minpack.lm::nls.lm.control(maxiter = 1000)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- pmax(as.numeric(pred), 0)
pred <- cummax(pred)
cf <- stats::coef(mod)
params$a <- unname(cf["a"])
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$b <- params$k * params$t0
params$a2 <- unname(cf["a2"])
params$k2 <- unname(cf["k2"])
params$t02 <- unname(cf["t0"] + cf["dt"])
params$b2 <- params$k2 * params$t02
params$converged <- TRUE
params$method_used <- "double_gompertz"
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
params$method_used <- "double_gompertz"
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit double-Richards model
#'
#' @param x Observed season-day values.
#' @param y Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param fixed_a_value Optional fixed total asymptote on original scale.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param start_value_double_richards_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgfd_fit_model_double_richards <- function(
x,
y,
x_pred,
fixed_a_value = NA_real_,
fix_a_to_observed_max = FALSE,
start_value_double_richards_parameters = list(
a = NA_real_, k = NA_real_, t0 = NA_real_, v = 1,
a2 = NA_real_, k2 = NA_real_, t02 = NA_real_, v2 = 1
)
) {
params <- dmgfd_empty_model_params()
res <- tryCatch({
cons <- ifelse(
min(y, na.rm = TRUE) < 0,
abs(min(y, na.rm = TRUE)) + 1e-8,
0
)
y_adj <- y + cons
a0 <- max(y_adj, na.rm = TRUE) * 1.05 + 1e-6
start_guess <- dmgfd_infer_double_richards_starts(x, y_adj, a0)
start_guess <- utils::modifyList(start_guess, start_value_double_richards_parameters)
if (!dmgfd_is_scalar_finite(start_guess$a)) {
start_guess$a <- a0 * 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k)) {
start_guess$k <- 0.03
}
if (!dmgfd_is_scalar_finite(start_guess$t0)) {
start_guess$t0 <- stats::quantile(x, 0.33, na.rm = TRUE, names = FALSE)
}
if (!dmgfd_is_scalar_finite(start_guess$v)) {
start_guess$v <- 1
}
if (!dmgfd_is_scalar_finite(start_guess$a2)) {
start_guess$a2 <- a0 * 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k2)) {
start_guess$k2 <- start_guess$k
}
if (!dmgfd_is_scalar_finite(start_guess$t02)) {
start_guess$t02 <- stats::quantile(x, 0.67, na.rm = TRUE, names = FALSE)
}
if (!dmgfd_is_scalar_finite(start_guess$v2)) {
start_guess$v2 <- 1
}
dt_start <- max(5, start_guess$t02 - start_guess$t0)
dat <- data.frame(
x = x,
y = y_adj
)
use_fixed_a <- isTRUE(fix_a_to_observed_max) &&
dmgfd_is_scalar_finite(fixed_a_value) &&
fixed_a_value > 0
if (isTRUE(use_fixed_a)) {
a_fixed_original <- fixed_a_value
a_fixed_fit <- fixed_a_value + cons
p_start <- 0.5
if (dmgfd_is_scalar_finite(start_guess$a) &&
dmgfd_is_scalar_finite(start_guess$a2) &&
(start_guess$a + start_guess$a2) > 0) {
p_start <- start_guess$a / (start_guess$a + start_guess$a2)
}
p_start <- min(max(p_start, 0.05), 0.95)
mod <- minpack.lm::nlsLM(
y ~ a_fixed_fit * p / ((1 + v * exp(-k * (x - t0)))^(1 / v)) +
a_fixed_fit * (1 - p) / ((1 + v2 * exp(-k2 * (x - (t0 + dt))))^(1 / v2)),
data = dat,
start = list(
p = p_start,
k = start_guess$k,
t0 = start_guess$t0,
v = start_guess$v,
k2 = start_guess$k2,
dt = dt_start,
v2 = start_guess$v2
),
lower = c(
p = 0.001,
k = 1e-6,
t0 = min(x, na.rm = TRUE),
v = 1e-3,
k2 = 1e-6,
dt = 1,
v2 = 1e-3
),
upper = c(
p = 0.999,
k = 5,
t0 = max(x, na.rm = TRUE),
v = 50,
k2 = 5,
dt = diff(range(x, na.rm = TRUE)),
v2 = 50
),
control = minpack.lm::nls.lm.control(maxiter = 1000)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- pmax(as.numeric(pred), 0)
pred <- cummax(pred)
cf <- stats::coef(mod)
p_est <- unname(cf["p"])
params$a <- a_fixed_original * p_est
params$a2 <- a_fixed_original * (1 - p_est)
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$v <- unname(cf["v"])
params$k2 <- unname(cf["k2"])
params$t02 <- unname(cf["t0"] + cf["dt"])
params$v2 <- unname(cf["v2"])
params$converged <- TRUE
params$fixed_a_used <- TRUE
params$fixed_a_value <- a_fixed_original
params$method_used <- "double_richards"
params$model_note <- "Double-Richards total asymptote a + a2 fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
mod <- minpack.lm::nlsLM(
y ~ a / ((1 + v * exp(-k * (x - t0)))^(1 / v)) +
a2 / ((1 + v2 * exp(-k2 * (x - (t0 + dt))))^(1 / v2)),
data = dat,
start = list(
a = start_guess$a,
k = start_guess$k,
t0 = start_guess$t0,
v = start_guess$v,
a2 = start_guess$a2,
k2 = start_guess$k2,
dt = dt_start,
v2 = start_guess$v2
),
lower = c(
a = 0,
k = 1e-6,
t0 = min(x, na.rm = TRUE),
v = 1e-3,
a2 = 0,
k2 = 1e-6,
dt = 1,
v2 = 1e-3
),
upper = c(
a = Inf,
k = 5,
t0 = max(x, na.rm = TRUE),
v = 50,
a2 = Inf,
k2 = 5,
dt = diff(range(x, na.rm = TRUE)),
v2 = 50
),
control = minpack.lm::nls.lm.control(maxiter = 1000)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- pmax(as.numeric(pred), 0)
pred <- cummax(pred)
cf <- stats::coef(mod)
params$a <- unname(cf["a"])
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$v <- unname(cf["v"])
params$a2 <- unname(cf["a2"])
params$k2 <- unname(cf["k2"])
params$t02 <- unname(cf["t0"] + cf["dt"])
params$v2 <- unname(cf["v2"])
params$converged <- TRUE
params$method_used <- "double_richards"
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
params$method_used <- "double_richards"
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit single-Gompertz fallback model
#'
#' @param x Observed season-day values.
#' @param y Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param fixed_a_value Optional fixed asymptote.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param start_value_gompertz_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgfd_fit_model_single_gompertz <- function(
x,
y,
x_pred,
fixed_a_value = NA_real_,
fix_a_to_observed_max = FALSE,
start_value_gompertz_parameters = list(a = NA_real_, b = NA_real_, k = NA_real_)
) {
params <- dmgfd_empty_model_params()
res <- tryCatch({
cons <- ifelse(
min(y, na.rm = TRUE) < 0,
abs(min(y, na.rm = TRUE)) + 1e-8,
0
)
y_adj <- y + cons
a0 <- max(y_adj, na.rm = TRUE) * 1.05 + 1e-6
use_fixed_a <- isTRUE(fix_a_to_observed_max) &&
dmgfd_is_scalar_finite(fixed_a_value) &&
fixed_a_value > 0
if (isTRUE(use_fixed_a)) {
a_fixed_original <- fixed_a_value
a_fixed_fit <- fixed_a_value + cons
start_guess <- dmgfd_infer_gompertz_starts(x, y_adj, a_fixed_fit)
start_guess <- utils::modifyList(start_guess, start_value_gompertz_parameters)
if (!dmgfd_is_scalar_finite(start_guess$b)) {
start_guess$b <- 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k)) {
start_guess$k <- 0.01
}
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a_fixed_fit * exp(-exp(b - k * x)),
data = dat,
start = list(
b = start_guess$b,
k = start_guess$k
),
control = minpack.lm::nls.lm.control(maxiter = 500)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- as.numeric(pred)
cf <- stats::coef(mod)
params$a <- a_fixed_original
params$b <- unname(cf["b"])
params$k <- unname(cf["k"])
params$converged <- TRUE
params$fixed_a_used <- TRUE
params$fixed_a_value <- a_fixed_original
params$method_used <- "single_gompertz"
params$model_note <- "Single-Gompertz fallback asymptote a fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
start_guess <- dmgfd_infer_gompertz_starts(x, y_adj, a0)
start_guess <- utils::modifyList(start_guess, start_value_gompertz_parameters)
if (!dmgfd_is_scalar_finite(start_guess$a)) {
start_guess$a <- a0
}
if (!dmgfd_is_scalar_finite(start_guess$b)) {
start_guess$b <- 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k)) {
start_guess$k <- 0.01
}
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a * exp(-exp(b - k * x)),
data = dat,
start = start_guess,
control = minpack.lm::nls.lm.control(maxiter = 500)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- as.numeric(pred)
cf <- stats::coef(mod)
params$a <- unname(cf["a"])
params$b <- unname(cf["b"])
params$k <- unname(cf["k"])
params$converged <- TRUE
params$method_used <- "single_gompertz"
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
params$method_used <- "single_gompertz"
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit single-Richards fallback model
#'
#' @param x Observed season-day values.
#' @param y Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param fixed_a_value Optional fixed asymptote.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param start_value_richards_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgfd_fit_model_single_richards <- function(
x,
y,
x_pred,
fixed_a_value = NA_real_,
fix_a_to_observed_max = FALSE,
start_value_richards_parameters = list(a = NA_real_, k = NA_real_, t0 = NA_real_, v = 1)
) {
params <- dmgfd_empty_model_params()
res <- tryCatch({
cons <- ifelse(
min(y, na.rm = TRUE) < 0,
abs(min(y, na.rm = TRUE)) + 1e-8,
0
)
y_adj <- y + cons
a0 <- max(y_adj, na.rm = TRUE) * 1.05 + 1e-6
use_fixed_a <- isTRUE(fix_a_to_observed_max) &&
dmgfd_is_scalar_finite(fixed_a_value) &&
fixed_a_value > 0
if (isTRUE(use_fixed_a)) {
a_fixed_original <- fixed_a_value
a_fixed_fit <- fixed_a_value + cons
logi_start <- dmgfd_infer_logistic_starts(x, y_adj, a_fixed_fit)
start_vals <- list(
k = if (dmgfd_is_scalar_finite(start_value_richards_parameters$k)) {
start_value_richards_parameters$k
} else {
logi_start$k
},
t0 = if (dmgfd_is_scalar_finite(start_value_richards_parameters$t0)) {
start_value_richards_parameters$t0
} else {
logi_start$t0
},
v = if (dmgfd_is_scalar_finite(start_value_richards_parameters$v)) {
start_value_richards_parameters$v
} else {
1
}
)
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a_fixed_fit / ((1 + v * exp(-k * (x - t0)))^(1 / v)),
data = dat,
start = start_vals,
control = minpack.lm::nls.lm.control(maxiter = 500)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- as.numeric(pred)
cf <- stats::coef(mod)
params$a <- a_fixed_original
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$v <- unname(cf["v"])
params$converged <- TRUE
params$fixed_a_used <- TRUE
params$fixed_a_value <- a_fixed_original
params$method_used <- "single_richards"
params$model_note <- "Single-Richards fallback asymptote a fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
logi_start <- dmgfd_infer_logistic_starts(x, y_adj, a0)
start_vals <- list(
a = if (dmgfd_is_scalar_finite(start_value_richards_parameters$a)) {
start_value_richards_parameters$a
} else {
logi_start$a
},
k = if (dmgfd_is_scalar_finite(start_value_richards_parameters$k)) {
start_value_richards_parameters$k
} else {
logi_start$k
},
t0 = if (dmgfd_is_scalar_finite(start_value_richards_parameters$t0)) {
start_value_richards_parameters$t0
} else {
logi_start$t0
},
v = if (dmgfd_is_scalar_finite(start_value_richards_parameters$v)) {
start_value_richards_parameters$v
} else {
1
}
)
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a / ((1 + v * exp(-k * (x - t0)))^(1 / v)),
data = dat,
start = start_vals,
control = minpack.lm::nls.lm.control(maxiter = 500)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- as.numeric(pred)
cf <- stats::coef(mod)
params$a <- unname(cf["a"])
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$v <- unname(cf["v"])
params$converged <- TRUE
params$method_used <- "single_richards"
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
params$method_used <- "single_richards"
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Add overall cumulative-growth timing
#'
#' @param params Parameter list.
#' @param pred_full Fitted cumulative growth.
#' @param x_pred Prediction season-day values.
#' @param growth_fraction Growth fractions.
#' @param season_start_date Season start date.
#'
#' @return Updated parameter list.
#'
#' @keywords internal
dmgfd_add_overall_timing <- function(params,
pred_full,
x_pred,
growth_fraction,
season_start_date = as.Date(NA)) {
valid_pred <- is.finite(pred_full) & is.finite(x_pred)
if (sum(valid_pred) == 0) {
return(params)
}
pred_use <- pred_full[valid_pred]
doy_use <- x_pred[valid_pred]
final_val <- tail(stats::na.omit(pred_use), 1)
if (length(final_val) == 1 &&
is.finite(final_val) &&
final_val > 0) {
low_val <- final_val * growth_fraction[1]
high_val <- final_val * growth_fraction[2]
st_pos <- suppressWarnings(min(doy_use[pred_use > low_val], na.rm = TRUE))
en_pos <- suppressWarnings(max(doy_use[pred_use < high_val], na.rm = TRUE))
if (is.finite(st_pos)) {
params$growth_start_season_day <- st_pos
if (!is.na(season_start_date)) {
params$growth_start_date <- season_start_date + (st_pos - 1)
params$growth_start_day <- lubridate::yday(params$growth_start_date)
}
}
if (is.finite(en_pos)) {
params$growth_end_season_day <- en_pos
if (!is.na(season_start_date)) {
params$growth_end_date <- season_start_date + (en_pos - 1)
params$growth_end_day <- lubridate::yday(params$growth_end_date)
}
}
}
params
}
#' Add overall rate-based timing
#'
#' @param params Parameter list.
#' @param pred_full Fitted cumulative growth.
#' @param x_pred Prediction season-day values.
#' @param rate_threshold_fraction Fraction of peak rate.
#' @param season_start_date Season start date.
#'
#' @return Updated parameter list.
#'
#' @keywords internal
dmgfd_add_overall_rate_timing <- function(params,
pred_full,
x_pred,
rate_threshold_fraction = 0.1,
season_start_date = as.Date(NA)) {
rate <- dmgfd_growth_rate(pred_full)
ok <- is.finite(rate) & is.finite(x_pred)
if (sum(ok) == 0) {
return(params)
}
r <- rate[ok]
x <- x_pred[ok]
peak_rate <- suppressWarnings(max(r, na.rm = TRUE))
if (!is.finite(peak_rate) || peak_rate <= 0) {
return(params)
}
thr <- rate_threshold_fraction * peak_rate
above <- which(r >= thr)
if (length(above) == 0) {
return(params)
}
params$peak_rate <- peak_rate
params$rate_start_season_day <- min(x[above], na.rm = TRUE)
params$rate_end_season_day <- max(x[above], na.rm = TRUE)
if (!is.na(season_start_date)) {
params$rate_start_date <- season_start_date + (params$rate_start_season_day - 1)
params$rate_end_date <- season_start_date + (params$rate_end_season_day - 1)
params$rate_start_day <- lubridate::yday(params$rate_start_date)
params$rate_end_day <- lubridate::yday(params$rate_end_date)
}
params
}
#' Fit single-curve fallback
#'
#' @param method Original double method.
#' @param x Observed season-day values.
#' @param y Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param fixed_a_value Optional fixed asymptote.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param start_value_double_gompertz_parameters Double-Gompertz starts.
#' @param start_value_double_richards_parameters Double-Richards starts.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgfd_fit_single_fallback <- function(method,
x,
y,
x_pred,
fixed_a_value = NA_real_,
fix_a_to_observed_max = FALSE,
start_value_double_gompertz_parameters,
start_value_double_richards_parameters) {
params <- dmgfd_empty_model_params()
if (method == "gompertz") {
start_single <- list(
a = if (dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$a) &&
dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$a2)) {
start_value_double_gompertz_parameters$a +
start_value_double_gompertz_parameters$a2
} else if (dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$a)) {
start_value_double_gompertz_parameters$a
} else {
NA_real_
},
b = if (dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$k) &&
dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$t0)) {
start_value_double_gompertz_parameters$k *
start_value_double_gompertz_parameters$t0
} else {
NA_real_
},
k = if (dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$k)) {
start_value_double_gompertz_parameters$k
} else {
NA_real_
}
)
single_res <- dmgfd_fit_model_single_gompertz(
x = x,
y = y,
x_pred = x_pred,
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_gompertz_parameters = start_single
)
params$a <- single_res$params$a
params$b <- single_res$params$b
params$k <- single_res$params$k
params$converged <- single_res$params$converged
params$fixed_a_used <- single_res$params$fixed_a_used
params$fixed_a_value <- single_res$params$fixed_a_value
params$model_note <- single_res$params$model_note
params$method_used <- "single_gompertz"
params$fallback_used <- TRUE
return(list(pred = single_res$pred, params = params))
}
if (method == "richards") {
start_single <- list(
a = if (dmgfd_is_scalar_finite(start_value_double_richards_parameters$a) &&
dmgfd_is_scalar_finite(start_value_double_richards_parameters$a2)) {
start_value_double_richards_parameters$a +
start_value_double_richards_parameters$a2
} else if (dmgfd_is_scalar_finite(start_value_double_richards_parameters$a)) {
start_value_double_richards_parameters$a
} else {
NA_real_
},
k = if (dmgfd_is_scalar_finite(start_value_double_richards_parameters$k)) {
start_value_double_richards_parameters$k
} else {
NA_real_
},
t0 = if (dmgfd_is_scalar_finite(start_value_double_richards_parameters$t0)) {
start_value_double_richards_parameters$t0
} else {
NA_real_
},
v = if (dmgfd_is_scalar_finite(start_value_double_richards_parameters$v)) {
start_value_double_richards_parameters$v
} else {
1
}
)
single_res <- dmgfd_fit_model_single_richards(
x = x,
y = y,
x_pred = x_pred,
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_richards_parameters = start_single
)
params$a <- single_res$params$a
params$k <- single_res$params$k
params$t0 <- single_res$params$t0
params$v <- single_res$params$v
params$converged <- single_res$params$converged
params$fixed_a_used <- single_res$params$fixed_a_used
params$fixed_a_value <- single_res$params$fixed_a_value
params$model_note <- single_res$params$model_note
params$method_used <- "single_richards"
params$fallback_used <- TRUE
return(list(pred = single_res$pred, params = params))
}
stop("Unknown single fallback method: ", method)
}
#' Calculate non-negative fitted growth rate
#'
#' @param pred Fitted cumulative growth.
#'
#' @return Numeric growth-rate vector.
#'
#' @keywords internal
dmgfd_growth_rate <- function(pred) {
pred <- as.numeric(pred)
if (length(pred) < 2) {
return(rep(NA_real_, length(pred)))
}
rate <- c(NA_real_, diff(pred))
rate[!is.finite(rate)] <- NA_real_
rate <- pmax(rate, 0)
rate
}
#' Identify local maxima
#'
#' @param x Numeric vector.
#'
#' @return Integer indices of local maxima.
#'
#' @keywords internal
dmgfd_local_maxima <- function(x) {
x <- as.numeric(x)
n <- length(x)
if (n < 3) {
return(integer(0))
}
which(
is.finite(x[2:(n - 1)]) &
x[2:(n - 1)] >= x[1:(n - 2)] &
x[2:(n - 1)] >= x[3:n]
) + 1L
}
#' Detect two-pulse pattern from fitted growth rate
#'
#' @param rate Fitted growth-rate vector.
#' @param x_pred Prediction season-day values.
#' @param peak_separation_min Minimum separation between peaks.
#' @param valley_ratio_max Maximum allowed valley-to-peak ratio.
#' @param min_relative_peak Minimum relative peak height.
#'
#' @return List describing detected pulse pattern.
#'
#' @keywords internal
dmgfd_detect_two_pulse_pattern <- function(rate,
x_pred,
peak_separation_min = 10,
valley_ratio_max = 0.4,
min_relative_peak = 0.05) {
out <- list(
detected = FALSE,
peak1_day = NA_real_,
peak2_day = NA_real_,
separator_day = NA_real_,
valley_rate = NA_real_,
pulse1_peak_rate = NA_real_,
pulse2_peak_rate = NA_real_
)
ok <- is.finite(rate) & is.finite(x_pred)
if (sum(ok) < 5) {
return(out)
}
r <- as.numeric(rate)
x <- as.numeric(x_pred)
peak_idx <- dmgfd_local_maxima(r)
if (length(peak_idx) < 2) {
return(out)
}
global_peak <- max(r, na.rm = TRUE)
if (!is.finite(global_peak) || global_peak <= 0) {
return(out)
}
peak_idx <- peak_idx[r[peak_idx] >= min_relative_peak * global_peak]
if (length(peak_idx) < 2) {
return(out)
}
best_score <- -Inf
best <- NULL
for (i in seq_len(length(peak_idx) - 1)) {
for (j in (i + 1):length(peak_idx)) {
p1 <- peak_idx[i]
p2 <- peak_idx[j]
if ((x[p2] - x[p1]) < peak_separation_min) {
next
}
valley_idx <- p1:p2
valley_rate <- suppressWarnings(min(r[valley_idx], na.rm = TRUE))
if (!is.finite(valley_rate)) {
next
}
if (valley_rate > valley_ratio_max * min(r[p1], r[p2])) {
next
}
score <- r[p1] + r[p2] - valley_rate
if (score > best_score) {
sep_idx <- valley_idx[which.min(r[valley_idx])]
best_score <- score
best <- list(
p1 = p1,
p2 = p2,
sep = sep_idx,
valley_rate = valley_rate
)
}
}
}
if (is.null(best)) {
return(out)
}
out$detected <- TRUE
out$peak1_day <- x[best$p1]
out$peak2_day <- x[best$p2]
out$separator_day <- x[best$sep]
out$valley_rate <- best$valley_rate
out$pulse1_peak_rate <- r[best$p1]
out$pulse2_peak_rate <- r[best$p2]
out
}
#' Estimate pulse window from fitted growth rate
#'
#' @param rate Fitted growth-rate vector.
#' @param x_pred Prediction season-day values.
#' @param threshold Pulse-specific rate threshold.
#' @param domain_idx Indices belonging to one pulse domain.
#'
#' @return List with start and end season-day.
#'
#' @keywords internal
dmgfd_pulse_window_from_rate <- function(rate,
x_pred,
threshold,
domain_idx) {
out <- list(
start_day = NA_real_,
end_day = NA_real_
)
if (length(domain_idx) == 0) {
return(out)
}
r <- rate[domain_idx]
x <- x_pred[domain_idx]
ok <- is.finite(r) & is.finite(x)
if (sum(ok) == 0) {
return(out)
}
r <- r[ok]
x <- x[ok]
above <- which(r >= threshold)
if (length(above) == 0) {
return(out)
}
out$start_day <- min(x[above], na.rm = TRUE)
out$end_day <- max(x[above], na.rm = TRUE)
out
}
#' Fit one double growth curve
#'
#' @description
#' Internal helper used by [dm.growth.fit.double()] to fit one double curve for
#' one series and one season, or one pooled double curve.
#'
#' @param x_obs Observed season-day values.
#' @param y_obs Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param season_start_date Season start date.
#' @param season_length Season length.
#' @param method Double-growth method.
#' @param growth_fraction Growth fractions used for cumulative timing.
#' @param min_unique_growth Minimum number of unique growth values.
#' @param rate_threshold_fraction Fraction of peak rate.
#' @param peak_separation_min Minimum separation between peaks.
#' @param valley_ratio_max Maximum allowed valley-to-peak ratio.
#' @param min_relative_peak Minimum relative derivative peak height.
#' @param fallback_to_single Logical fallback flag.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param fixed_a_multiplier Multiplier for fixed asymptote.
#' @param start_value_double_gompertz_parameters Double-Gompertz starts.
#' @param start_value_double_richards_parameters Double-Richards starts.
#'
#' @return List with fitted predictions and parameters.
#'
#' @keywords internal
dmgfd_fit_one_curve <- function(x_obs,
y_obs,
x_pred,
season_start_date = as.Date(NA),
season_length = NA_real_,
method,
growth_fraction,
min_unique_growth,
rate_threshold_fraction = 0.1,
peak_separation_min = 10,
valley_ratio_max = 0.4,
min_relative_peak = 0.05,
fallback_to_single = TRUE,
fix_a_to_observed_max = FALSE,
fixed_a_multiplier = 1,
start_value_double_gompertz_parameters,
start_value_double_richards_parameters) {
pred_full <- rep(NA_real_, length(x_pred))
params <- dmgfd_empty_model_params()
params$n_obs <- sum(is.finite(x_obs) & is.finite(y_obs))
params$n_days_observed <- NA_integer_
params$first_obs_day <- NA_real_
params$last_obs_day <- NA_real_
params$season_length <- season_length
params$n_missing_days <- NA_real_
params$extrapolated <- NA
params$anchor_added <- FALSE
params$method_used <- paste0("double_", method)
good <- is.finite(x_obs) & is.finite(y_obs)
if (sum(good) == 0) {
params$model_note <- "No observed data in this vegetation year."
return(list(pred = pred_full, params = params))
}
x_raw <- x_obs[good]
y_raw <- y_obs[good]
params$n_days_observed <- length(x_raw)
params$first_obs_day <- min(x_raw, na.rm = TRUE)
params$last_obs_day <- max(x_raw, na.rm = TRUE)
if (is.finite(season_length)) {
params$n_missing_days <- season_length - length(unique(x_raw))
params$extrapolated <- params$first_obs_day > 1 ||
params$last_obs_day < season_length
}
if (length(unique(y_raw)) < min_unique_growth) {
params$model_note <- paste0(
"Too few unique cumulative-growth values (",
length(unique(y_raw)), " < ", min_unique_growth, ")."
)
return(list(pred = pred_full, params = params))
}
anchored <- dmgfd_add_anchor_points(x_raw, y_raw)
x_fit <- anchored$x
y_fit <- anchored$y
params$anchor_added <- length(x_fit) > length(x_raw)
dyn_range <- max(y_fit, na.rm = TRUE) - min(y_fit, na.rm = TRUE)
if (!is.finite(dyn_range) || dyn_range <= 0) {
params$model_note <- "Observed cumulative-growth range is too small."
return(list(pred = pred_full, params = params))
}
fixed_a_value <- NA_real_
if (isTRUE(fix_a_to_observed_max) &&
method %in% c("gompertz", "richards")) {
fixed_a_value <- suppressWarnings(max(y_raw, na.rm = TRUE)) * fixed_a_multiplier
if (!is.finite(fixed_a_value) || fixed_a_value <= 0) {
fixed_a_value <- NA_real_
}
}
fit_res <- switch(
method,
gompertz = dmgfd_fit_model_double_gompertz(
x = x_fit,
y = y_fit,
x_pred = x_pred,
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_double_gompertz_parameters = start_value_double_gompertz_parameters
),
richards = dmgfd_fit_model_double_richards(
x = x_fit,
y = y_fit,
x_pred = x_pred,
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_double_richards_parameters = start_value_double_richards_parameters
)
)
pred_full <- fit_res$pred
params[names(fit_res$params)] <- fit_res$params
params$method_used <- paste0("double_", method)
params <- dmgfd_add_overall_timing(
params = params,
pred_full = pred_full,
x_pred = x_pred,
growth_fraction = growth_fraction,
season_start_date = season_start_date
)
params <- dmgfd_add_overall_rate_timing(
params = params,
pred_full = pred_full,
x_pred = x_pred,
rate_threshold_fraction = rate_threshold_fraction,
season_start_date = season_start_date
)
rate <- dmgfd_growth_rate(pred_full)
pulse_info <- dmgfd_detect_two_pulse_pattern(
rate = rate,
x_pred = x_pred,
peak_separation_min = peak_separation_min,
valley_ratio_max = valley_ratio_max,
min_relative_peak = min_relative_peak
)
if (isTRUE(pulse_info$detected)) {
params$two_pulse_detected <- TRUE
params$peak1_season_day <- pulse_info$peak1_day
params$peak2_season_day <- pulse_info$peak2_day
params$separator_season_day <- pulse_info$separator_day
params$valley_rate <- pulse_info$valley_rate
params$pulse1_peak_rate <- pulse_info$pulse1_peak_rate
params$pulse2_peak_rate <- pulse_info$pulse2_peak_rate
if (!is.na(season_start_date)) {
peak1_date <- season_start_date + (params$peak1_season_day - 1)
peak2_date <- season_start_date + (params$peak2_season_day - 1)
params$separator_date <- season_start_date + (params$separator_season_day - 1)
params$peak1_day <- lubridate::yday(peak1_date)
params$peak2_day <- lubridate::yday(peak2_date)
params$separator_day <- lubridate::yday(params$separator_date)
}
thr1 <- rate_threshold_fraction * params$pulse1_peak_rate
thr2 <- rate_threshold_fraction * params$pulse2_peak_rate
idx1 <- which(is.finite(x_pred) & x_pred <= params$separator_season_day)
idx2 <- which(is.finite(x_pred) & x_pred > params$separator_season_day)
win1 <- dmgfd_pulse_window_from_rate(
rate = rate,
x_pred = x_pred,
threshold = thr1,
domain_idx = idx1
)
win2 <- dmgfd_pulse_window_from_rate(
rate = rate,
x_pred = x_pred,
threshold = thr2,
domain_idx = idx2
)
if (is.finite(win1$start_day)) {
params$pulse1_start_season_day <- win1$start_day
if (!is.na(season_start_date)) {
params$pulse1_start_date <- season_start_date + (win1$start_day - 1)
params$pulse1_start_day <- lubridate::yday(params$pulse1_start_date)
}
}
if (is.finite(win1$end_day)) {
params$pulse1_end_season_day <- win1$end_day
if (!is.na(season_start_date)) {
params$pulse1_end_date <- season_start_date + (win1$end_day - 1)
params$pulse1_end_day <- lubridate::yday(params$pulse1_end_date)
}
}
if (is.finite(win2$start_day)) {
params$pulse2_start_season_day <- win2$start_day
if (!is.na(season_start_date)) {
params$pulse2_start_date <- season_start_date + (win2$start_day - 1)
params$pulse2_start_day <- lubridate::yday(params$pulse2_start_date)
}
}
if (is.finite(win2$end_day)) {
params$pulse2_end_season_day <- win2$end_day
if (!is.na(season_start_date)) {
params$pulse2_end_date <- season_start_date + (win2$end_day - 1)
params$pulse2_end_day <- lubridate::yday(params$pulse2_end_date)
}
}
return(list(pred = pred_full, params = params))
}
params$two_pulse_detected <- FALSE
if (!isTRUE(fallback_to_single)) {
if (is.na(params$model_note) || !nzchar(params$model_note)) {
params$model_note <- "No convincing two-pulse pattern detected."
} else {
params$model_note <- paste(
params$model_note,
"No convincing two-pulse pattern detected."
)
}
return(list(pred = pred_full, params = params))
}
meta_fields <- c(
"n_obs", "n_days_observed", "first_obs_day", "last_obs_day",
"season_length", "n_missing_days", "extrapolated", "anchor_added"
)
double_note <- params$model_note
fallback_res <- dmgfd_fit_single_fallback(
method = method,
x = x_fit,
y = y_fit,
x_pred = x_pred,
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_double_gompertz_parameters = start_value_double_gompertz_parameters,
start_value_double_richards_parameters = start_value_double_richards_parameters
)
pred_full <- fallback_res$pred
fb_params <- dmgfd_empty_model_params()
fb_params[meta_fields] <- params[meta_fields]
fb_params[names(fallback_res$params)] <- fallback_res$params
fb_params$two_pulse_detected <- FALSE
fb_params$fallback_used <- TRUE
if (is.na(double_note) || !nzchar(double_note)) {
fb_params$model_note <- paste0(
"No convincing two-pulse pattern detected; ",
fb_params$method_used,
" fallback used."
)
} else {
fb_params$model_note <- paste(
double_note,
paste0(
"No convincing two-pulse pattern detected; ",
fb_params$method_used,
" fallback used."
)
)
}
fb_params <- dmgfd_add_overall_timing(
params = fb_params,
pred_full = pred_full,
x_pred = x_pred,
growth_fraction = growth_fraction,
season_start_date = season_start_date
)
fb_params <- dmgfd_add_overall_rate_timing(
params = fb_params,
pred_full = pred_full,
x_pred = x_pred,
rate_threshold_fraction = rate_threshold_fraction,
season_start_date = season_start_date
)
list(
pred = pred_full,
params = fb_params
)
}
Try the dendRoAnalyst package in your browser
Any scripts or data that you put into this service are public.
dendRoAnalyst documentation built on May 20, 2026, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions dm.growth.fit.double
Documented in dm.growth.fit.double
#' Fit bimodal dendrometer growth curves by vegetation season
#'
#' @description
#' Fits bimodal cumulative growth curves to daily dendrometer series using either
#' a double-Gompertz or double-Richards model.
#'
#' Overall growing-season timing is derived from the fitted cumulative-growth
#' curve using \code{growth_fraction}. Pulse-specific timing is derived from the
#' derivative of the fitted curve using a pulse-specific relative rate threshold
#' given by \code{rate_threshold_fraction}.
#'
#' Pulse-specific timing is returned as \code{NA} when the fitted curve does not
#' show a convincing two-pulse pattern according to derivative-based criteria.
#'
#' If \code{fallback_to_single = TRUE} and no convincing two-pulse pattern is
#' found, the function refits a corresponding single growth curve and returns
#' that fit instead. In that case, overall season timing is still returned, but
#' pulse-specific timing remains \code{NA}.
#'
#' For double-Gompertz and double-Richards fits, the total asymptote can
#' optionally be fixed to the observed maximum cumulative growth of the
#' corresponding series and vegetation season. In the double-curve case, this
#' means \code{a + a2} is fixed, while the relative contribution of the first
#' and second pulse is estimated.
#'
#' @param df A data frame or tibble. The first column must be a time stamp
#' \code{Date}, \code{POSIXct}, or parseable date-time string. Remaining
#' selected columns must be numeric dendrometer series.
#' @param TreeNum Either \code{"all"} to use all dendrometer series, a numeric
#' vector selecting dendrometer columns by position, or a character vector of
#' column names.
#' @param method Double-curve fitting method. One of \code{"gompertz"} or
#' \code{"richards"}.
#' @param years Either \code{"all"} to fit all available vegetation seasons, or
#' a character vector of season labels to retain.
#' @param year_mode Either \code{"yearly"} to fit one curve per vegetation
#' season, or \code{"pooled"} to fit one pooled curve across all retained
#' seasons.
#' @param fit_GRO Logical. If \code{TRUE}, processed daily series are converted
#' to cumulative growth using a cumulative maximum within each vegetation
#' season.
#' @param site_mode Vegetation season definition. One of \code{"NH"},
#' \code{"SH"}, or \code{"CS"}.
#' @param custom_veg_season Numeric vector of length two giving the start and
#' end day-of-year for custom vegetation seasons in \code{"CS"} mode.
#' @param growth_fraction Numeric vector of length two giving the lower and
#' upper fractions of final fitted seasonal growth used to define overall
#' growing-season onset and cessation.
#' @param min_unique_growth Minimum number of unique non-missing cumulative
#' growth values required before a fit is attempted.
#' @param rate_threshold_fraction Numeric scalar between 0 and 1. Pulse start
#' and end are defined as the first and last days where fitted growth rate
#' exceeds this fraction of the pulse-specific peak fitted growth rate.
#' @param peak_separation_min Minimum number of days separating the two fitted
#' derivative peaks required to classify a fit as truly two-pulse.
#' @param valley_ratio_max Maximum allowed ratio between the valley rate and the
#' weaker of the two derivative peaks. Smaller values require a deeper valley
#' between pulses.
#' @param min_relative_peak Minimum relative height, expressed as a fraction of
#' the global derivative maximum, for a local derivative peak to be considered.
#' @param fallback_to_single Logical. If \code{TRUE}, and the fitted double curve
#' does not show a convincing two-pulse pattern according to derivative-based
#' criteria, the function refits a corresponding single growth curve
#' \code{"gompertz"} or \code{"richards"} and returns that fit instead.
#' @param fix_a_to_observed_max Logical. If \code{TRUE}, the total asymptote of
#' the double-growth curve is fixed to the observed maximum cumulative growth
#' for that series and vegetation season. For double-Gompertz and
#' double-Richards models, this means \code{a + a2} is fixed, while the
#' relative contribution of the two pulses is estimated. This is most
#' appropriate with \code{year_mode = "yearly"}. For
#' \code{year_mode = "pooled"}, one pooled maximum is used.
#' @param fixed_a_multiplier Numeric multiplier applied to the observed maximum
#' when \code{fix_a_to_observed_max = TRUE}. The default is \code{1}, meaning
#' the total asymptote is fixed exactly to the observed seasonal maximum.
#' Values slightly above 1, for example \code{1.01} or \code{1.05}, allow the
#' fitted asymptote to sit slightly above the observed maximum.
#' @param start_value_double_gompertz_parameters Optional named list of starting
#' values for the double-Gompertz fit. Supported names are \code{a}, \code{k},
#' \code{t0}, \code{a2}, \code{k2}, and \code{t02}. If
#' \code{fix_a_to_observed_max = TRUE}, supplied \code{a} and \code{a2} are
#' used only to initialize the relative pulse contribution.
#' @param start_value_double_richards_parameters Optional named list of starting
#' values for the double-Richards fit. Supported names are \code{a}, \code{k},
#' \code{t0}, \code{v}, \code{a2}, \code{k2}, \code{t02}, and \code{v2}. If
#' \code{fix_a_to_observed_max = TRUE}, supplied \code{a} and \code{a2} are
#' used only to initialize the relative pulse contribution.
#' @param verbose Logical. If \code{TRUE}, prints a short completion message.
#'
#' @return
#' An object of class \code{"dm_growth_fit"} with elements:
#' \describe{
#' \item{call}{The matched function call.}
#' \item{original_daily_data}{Raw daily dendrometer data after resampling to
#' daily maxima and assigning vegetation seasons, but before centering and
#' optional cumulative-growth transformation.}
#' \item{processed_data}{Daily processed data used for fitting.}
#' \item{fitted_data}{Daily fitted values on the full vegetation-season grid.}
#' \item{fit_statistics}{Fit-level statistics and estimated timing.}
#' \item{fit_parameters}{Fit-level model parameters and convergence
#' information.}
#' \item{season_table}{Vegetation seasons retained for fitting.}
#' }
#'
#' @details
#' The second pulse is constrained to occur after the first pulse, which improves
#' numerical stability and reduces label switching between the two components.
#'
#' For double-Gompertz and double-Richards models, the total seasonal asymptote
#' is:
#'
#' \deqn{a_{total} = a + a2}
#'
#' When \code{fix_a_to_observed_max = TRUE}, the function fits:
#'
#' \deqn{a + a2 = max(y_{obs}) \times fixed\_a\_multiplier}
#'
#' where \eqn{y_{obs}} is the observed cumulative growth of the selected
#' dendrometer series and vegetation season.
#'
#' The fitted pulse contributions are still returned as \code{a} and \code{a2},
#' and their sum should equal \code{fixed_a_value}, apart from small numerical
#' rounding.
#'
#' Non-applicable parameters are returned as \code{NA}. For example, \code{b}
#' and \code{b2} are relevant for double-Gompertz, while \code{v} and \code{v2}
#' are relevant for double-Richards.
#'
#' @examples
#' \donttest{
#' # Double-Gompertz with total seasonal asymptote fixed
#' # fit_gomp <- dm.growth.fit.double(
#' # df = dendro_data,
#' # TreeNum = "all",
#' # method = "gompertz",
#' # year_mode = "yearly",
#' # fit_GRO = TRUE,
#' # fix_a_to_observed_max = TRUE
#' # )
#'
#' # Double-Richards with the total asymptote set 2 percent above observed max
#' # fit_rich <- dm.growth.fit.double(
#' # df = dendro_data,
#' # TreeNum = "all",
#' # method = "richards",
#' # year_mode = "yearly",
#' # fit_GRO = TRUE,
#' # fix_a_to_observed_max = TRUE,
#' # fixed_a_multiplier = 1.02
#' # )
#' }
#'
#' @seealso [dm.growth.fit()], [summary.dm_growth_fit()],
#' [print.dm_growth_fit()]
#'
#' @importFrom dplyr %>% select all_of any_of bind_cols bind_rows filter arrange
#' mutate across group_by ungroup distinct left_join
#' @importFrom tibble as_tibble tibble
#' @importFrom lubridate parse_date_time year yday month
#' @importFrom stats lm coef qlogis predict quantile median na.omit
#' @importFrom utils modifyList
#'
#' @export
dm.growth.fit.double <- function(
df,
TreeNum = "all",
method = c("gompertz", "richards"),
years = "all",
year_mode = c("yearly", "pooled"),
fit_GRO = TRUE,
site_mode = c("NH", "SH", "CS"),
custom_veg_season = c(275, 274),
growth_fraction = c(0.1, 0.9),
min_unique_growth = 10,
rate_threshold_fraction = 0.1,
peak_separation_min = 10,
valley_ratio_max = 0.4,
min_relative_peak = 0.05,
fallback_to_single = TRUE,
fix_a_to_observed_max = FALSE,
fixed_a_multiplier = 1,
start_value_double_gompertz_parameters = list(
a = NA_real_, k = NA_real_, t0 = NA_real_,
a2 = NA_real_, k2 = NA_real_, t02 = NA_real_
),
start_value_double_richards_parameters = list(
a = NA_real_, k = NA_real_, t0 = NA_real_, v = 1,
a2 = NA_real_, k2 = NA_real_, t02 = NA_real_, v2 = 1
),
verbose = TRUE
) {
cl <- match.call()
TIME <- any_observed <- doy <- season_day <- season_end <- NULL
season_label <- season_length <- season_start <- NULL
method <- match.arg(method)
year_mode <- match.arg(year_mode)
site_mode <- match.arg(site_mode)
dmgfd_validate_inputs(
df = df,
growth_fraction = growth_fraction,
min_unique_growth = min_unique_growth,
custom_veg_season = custom_veg_season,
rate_threshold_fraction = rate_threshold_fraction,
peak_separation_min = peak_separation_min,
valley_ratio_max = valley_ratio_max,
min_relative_peak = min_relative_peak,
fallback_to_single = fallback_to_single,
fix_a_to_observed_max = fix_a_to_observed_max,
fixed_a_multiplier = fixed_a_multiplier
)
df <- tibble::as_tibble(df)
names(df)[1] <- "TIME"
if (!inherits(df$TIME, "Date") &&
!inherits(df$TIME, "POSIXct") &&
!inherits(df$TIME, "POSIXt")) {
df$TIME <- suppressWarnings(
lubridate::parse_date_time(
df$TIME,
orders = c(
"ymd HMS", "ymd HM", "ymd",
"dmy HMS", "dmy HM", "dmy",
"mdy HMS", "mdy HM", "mdy"
),
quiet = TRUE
)
)
}
if (inherits(df$TIME, "Date")) {
df$TIME <- as.POSIXct(df$TIME)
}
if (any(is.na(df$TIME))) {
stop("Some timestamps in the first column could not be parsed.")
}
obs_years <- sort(unique(lubridate::year(as.Date(df$TIME))))
if (site_mode == "SH" && length(obs_years) < 2) {
stop("For site_mode = 'SH', df must contain observations from at least 2 calendar years.")
}
if (site_mode == "CS" &&
custom_veg_season[1] > custom_veg_season[2] &&
length(obs_years) < 2) {
stop(
"For site_mode = 'CS' with custom_veg_season[1] > custom_veg_season[2], ",
"df must contain observations from at least 2 calendar years."
)
}
all_series <- names(df)[-1]
if (length(all_series) == 0) {
stop("df must contain a time column followed by at least one dendrometer column.")
}
if (is.character(TreeNum) &&
length(TreeNum) == 1 &&
tolower(TreeNum) == "all") {
keep_cols <- all_series
} else if (is.numeric(TreeNum)) {
tree_idx <- as.integer(TreeNum)
if (any(is.na(tree_idx)) ||
any(tree_idx < 1) ||
any(tree_idx > length(all_series))) {
stop(
"'TreeNum' numeric values must be between 1 and ", length(all_series),
". Here, 1 refers to the first dendrometer series after the TIME column."
)
}
keep_cols <- all_series[tree_idx]
} else if (is.character(TreeNum)) {
missing_cols <- setdiff(TreeNum, all_series)
if (length(missing_cols) > 0) {
warning(
"These TreeNum column names were not found and were ignored: ",
paste(missing_cols, collapse = ", ")
)
}
keep_cols <- intersect(TreeNum, all_series)
if (length(keep_cols) == 0) {
stop("None of the requested TreeNum column names were found in df.")
}
} else {
stop("TreeNum must be 'all', a numeric vector, or a character vector of column names.")
}
df <- dplyr::select(df, TIME, dplyr::all_of(keep_cols))
num_cols <- names(df)[vapply(df, is.numeric, logical(1))]
num_cols <- setdiff(num_cols, "TIME")
if (length(num_cols) == 0) {
stop("No numeric dendrometer series found to fit.")
}
df <- dplyr::select(df, TIME, dplyr::all_of(num_cols))
dm_daily <- dendro.resample(
df = df,
by = "D",
value = "max",
method = "aggregate"
)
season_tbl_obs <- dmgfd_build_season_table(
time_vec = dm_daily$TIME,
site_mode = site_mode,
custom_veg_season = custom_veg_season
)
observed_daily <- dplyr::bind_cols(dm_daily, season_tbl_obs) %>%
dplyr::filter(!is.na(season_label)) %>%
dplyr::arrange(TIME)
original_daily_data <- observed_daily %>%
dplyr::select(
TIME, doy, season_label, season_start, season_end, season_day,
dplyr::all_of(num_cols)
)
season_table <- observed_daily %>%
dplyr::distinct(season_label, season_start, season_end) %>%
dplyr::arrange(season_start)
if (!(is.character(years) &&
length(years) == 1 &&
tolower(years) == "all")) {
season_table <- season_table %>%
dplyr::filter(season_label %in% as.character(years))
}
if (nrow(season_table) == 0) {
stop("No seasons remained after applying site_mode/custom_veg_season/years.")
}
full_grids <- lapply(seq_len(nrow(season_table)), function(i) {
ss <- season_table$season_start[i]
se <- season_table$season_end[i]
sl <- season_table$season_label[i]
tt <- seq(ss, se, by = "day")
tibble::tibble(
TIME = tt,
doy = lubridate::yday(tt),
season_label = sl,
season_start = ss,
season_end = se,
season_day = seq_along(tt),
season_length = length(tt)
)
})
full_grid <- dplyr::bind_rows(full_grids)
dat <- full_grid %>%
dplyr::left_join(
observed_daily %>%
dplyr::select(TIME, season_label, dplyr::all_of(num_cols)),
by = c("TIME", "season_label")
) %>%
dplyr::arrange(season_start, TIME)
dat_proc <- dat %>%
dplyr::group_by(season_label) %>%
dplyr::arrange(TIME, .by_group = TRUE) %>%
dplyr::mutate(
dplyr::across(
dplyr::all_of(num_cols),
~ . - dmgfd_first_non_na(.)
)
) %>%
dplyr::ungroup()
if (isTRUE(fit_GRO)) {
dat_proc <- dat_proc %>%
dplyr::group_by(season_label) %>%
dplyr::arrange(TIME, .by_group = TRUE) %>%
dplyr::mutate(
dplyr::across(
dplyr::all_of(num_cols),
dmgfd_cummax_na
)
) %>%
dplyr::ungroup()
}
dat_proc$any_observed <- rowSums(
!is.na(as.data.frame(dat_proc[, num_cols, drop = FALSE]))
) > 0
fitted_dat <- dat_proc %>%
dplyr::select(
TIME, doy, season_label, season_start, season_end,
season_day, season_length, any_observed
)
for (cc in num_cols) {
fitted_dat[[cc]] <- NA_real_
}
param_rows <- list()
row_id <- 1L
if (year_mode == "yearly") {
fit_ids <- unique(dat_proc$season_label)
for (series_name in num_cols) {
for (sid in fit_ids) {
idx <- which(dat_proc$season_label == sid)
x_full <- dat_proc$season_day[idx]
y_full <- dat_proc[[series_name]][idx]
ss <- unique(dat_proc$season_start[idx])[1]
se <- unique(dat_proc$season_end[idx])[1]
full_len <- unique(dat_proc$season_length[idx])[1]
good <- is.finite(y_full)
x_obs <- x_full[good]
y_obs <- y_full[good]
fit_res <- dmgfd_fit_one_curve(
x_obs = x_obs,
y_obs = y_obs,
x_pred = x_full,
season_start_date = ss,
season_length = full_len,
method = method,
growth_fraction = growth_fraction,
min_unique_growth = min_unique_growth,
rate_threshold_fraction = rate_threshold_fraction,
peak_separation_min = peak_separation_min,
valley_ratio_max = valley_ratio_max,
min_relative_peak = min_relative_peak,
fallback_to_single = fallback_to_single,
fix_a_to_observed_max = fix_a_to_observed_max,
fixed_a_multiplier = fixed_a_multiplier,
start_value_double_gompertz_parameters = start_value_double_gompertz_parameters,
start_value_double_richards_parameters = start_value_double_richards_parameters
)
fitted_dat[[series_name]][idx] <- fit_res$pred
param_rows[[row_id]] <- tibble::tibble(
series = series_name,
fit_id = sid,
year_mode = year_mode,
method = fit_res$params$method_used,
site_mode = site_mode,
season_start = ss,
season_end = se,
included_seasons = NA_character_,
season_length = fit_res$params$season_length,
n_obs = fit_res$params$n_obs,
n_days_observed = fit_res$params$n_days_observed,
first_obs_day = fit_res$params$first_obs_day,
last_obs_day = fit_res$params$last_obs_day,
n_missing_days = fit_res$params$n_missing_days,
extrapolated = fit_res$params$extrapolated,
anchor_added = fit_res$params$anchor_added,
growth_start_day = fit_res$params$growth_start_day,
growth_end_day = fit_res$params$growth_end_day,
growth_start_season_day = fit_res$params$growth_start_season_day,
growth_end_season_day = fit_res$params$growth_end_season_day,
growth_start_date = fit_res$params$growth_start_date,
growth_end_date = fit_res$params$growth_end_date,
peak_rate = fit_res$params$peak_rate,
rate_start_day = fit_res$params$rate_start_day,
rate_end_day = fit_res$params$rate_end_day,
rate_start_season_day = fit_res$params$rate_start_season_day,
rate_end_season_day = fit_res$params$rate_end_season_day,
rate_start_date = fit_res$params$rate_start_date,
rate_end_date = fit_res$params$rate_end_date,
fallback_used = fit_res$params$fallback_used,
two_pulse_detected = fit_res$params$two_pulse_detected,
peak1_day = fit_res$params$peak1_day,
peak2_day = fit_res$params$peak2_day,
peak1_season_day = fit_res$params$peak1_season_day,
peak2_season_day = fit_res$params$peak2_season_day,
separator_day = fit_res$params$separator_day,
separator_season_day = fit_res$params$separator_season_day,
separator_date = fit_res$params$separator_date,
valley_rate = fit_res$params$valley_rate,
pulse1_peak_rate = fit_res$params$pulse1_peak_rate,
pulse2_peak_rate = fit_res$params$pulse2_peak_rate,
pulse1_start_day = fit_res$params$pulse1_start_day,
pulse1_end_day = fit_res$params$pulse1_end_day,
pulse1_start_season_day = fit_res$params$pulse1_start_season_day,
pulse1_end_season_day = fit_res$params$pulse1_end_season_day,
pulse1_start_date = fit_res$params$pulse1_start_date,
pulse1_end_date = fit_res$params$pulse1_end_date,
pulse2_start_day = fit_res$params$pulse2_start_day,
pulse2_end_day = fit_res$params$pulse2_end_day,
pulse2_start_season_day = fit_res$params$pulse2_start_season_day,
pulse2_end_season_day = fit_res$params$pulse2_end_season_day,
pulse2_start_date = fit_res$params$pulse2_start_date,
pulse2_end_date = fit_res$params$pulse2_end_date,
converged = fit_res$params$converged,
fixed_a_used = fit_res$params$fixed_a_used,
fixed_a_value = fit_res$params$fixed_a_value,
a = fit_res$params$a,
b = fit_res$params$b,
k = fit_res$params$k,
t0 = fit_res$params$t0,
v = fit_res$params$v,
a2 = fit_res$params$a2,
b2 = fit_res$params$b2,
k2 = fit_res$params$k2,
t02 = fit_res$params$t02,
v2 = fit_res$params$v2,
edf = fit_res$params$edf,
span = fit_res$params$span,
spline_df = fit_res$params$spline_df,
spar = fit_res$params$spar,
model_note = fit_res$params$model_note
)
row_id <- row_id + 1L
}
}
}
if (year_mode == "pooled") {
included_seasons <- paste(unique(dat_proc$season_label), collapse = ", ")
for (series_name in num_cols) {
x_obs_all <- dat_proc$season_day
y_obs_all <- dat_proc[[series_name]]
good <- is.finite(y_obs_all)
x_obs <- x_obs_all[good]
y_obs <- y_obs_all[good]
x_template <- seq_len(max(dat_proc$season_day, na.rm = TRUE))
fit_res <- dmgfd_fit_one_curve(
x_obs = x_obs,
y_obs = y_obs,
x_pred = x_template,
season_start_date = as.Date(NA),
season_length = max(x_template, na.rm = TRUE),
method = method,
growth_fraction = growth_fraction,
min_unique_growth = min_unique_growth,
rate_threshold_fraction = rate_threshold_fraction,
peak_separation_min = peak_separation_min,
valley_ratio_max = valley_ratio_max,
min_relative_peak = min_relative_peak,
fallback_to_single = fallback_to_single,
fix_a_to_observed_max = fix_a_to_observed_max,
fixed_a_multiplier = fixed_a_multiplier,
start_value_double_gompertz_parameters = start_value_double_gompertz_parameters,
start_value_double_richards_parameters = start_value_double_richards_parameters
)
fitted_dat[[series_name]] <- fit_res$pred[match(dat_proc$season_day, x_template)]
param_rows[[row_id]] <- tibble::tibble(
series = series_name,
fit_id = "pooled",
year_mode = year_mode,
method = fit_res$params$method_used,
site_mode = site_mode,
season_start = as.Date(NA),
season_end = as.Date(NA),
included_seasons = included_seasons,
season_length = fit_res$params$season_length,
n_obs = fit_res$params$n_obs,
n_days_observed = fit_res$params$n_days_observed,
first_obs_day = fit_res$params$first_obs_day,
last_obs_day = fit_res$params$last_obs_day,
n_missing_days = fit_res$params$n_missing_days,
extrapolated = fit_res$params$extrapolated,
anchor_added = fit_res$params$anchor_added,
growth_start_day = fit_res$params$growth_start_day,
growth_end_day = fit_res$params$growth_end_day,
growth_start_season_day = fit_res$params$growth_start_season_day,
growth_end_season_day = fit_res$params$growth_end_season_day,
growth_start_date = fit_res$params$growth_start_date,
growth_end_date = fit_res$params$growth_end_date,
peak_rate = fit_res$params$peak_rate,
rate_start_day = fit_res$params$rate_start_day,
rate_end_day = fit_res$params$rate_end_day,
rate_start_season_day = fit_res$params$rate_start_season_day,
rate_end_season_day = fit_res$params$rate_end_season_day,
rate_start_date = fit_res$params$rate_start_date,
rate_end_date = fit_res$params$rate_end_date,
fallback_used = fit_res$params$fallback_used,
two_pulse_detected = fit_res$params$two_pulse_detected,
peak1_day = fit_res$params$peak1_day,
peak2_day = fit_res$params$peak2_day,
peak1_season_day = fit_res$params$peak1_season_day,
peak2_season_day = fit_res$params$peak2_season_day,
separator_day = fit_res$params$separator_day,
separator_season_day = fit_res$params$separator_season_day,
separator_date = fit_res$params$separator_date,
valley_rate = fit_res$params$valley_rate,
pulse1_peak_rate = fit_res$params$pulse1_peak_rate,
pulse2_peak_rate = fit_res$params$pulse2_peak_rate,
pulse1_start_day = fit_res$params$pulse1_start_day,
pulse1_end_day = fit_res$params$pulse1_end_day,
pulse1_start_season_day = fit_res$params$pulse1_start_season_day,
pulse1_end_season_day = fit_res$params$pulse1_end_season_day,
pulse1_start_date = fit_res$params$pulse1_start_date,
pulse1_end_date = fit_res$params$pulse1_end_date,
pulse2_start_day = fit_res$params$pulse2_start_day,
pulse2_end_day = fit_res$params$pulse2_end_day,
pulse2_start_season_day = fit_res$params$pulse2_start_season_day,
pulse2_end_season_day = fit_res$params$pulse2_end_season_day,
pulse2_start_date = fit_res$params$pulse2_start_date,
pulse2_end_date = fit_res$params$pulse2_end_date,
converged = fit_res$params$converged,
fixed_a_used = fit_res$params$fixed_a_used,
fixed_a_value = fit_res$params$fixed_a_value,
a = fit_res$params$a,
b = fit_res$params$b,
k = fit_res$params$k,
t0 = fit_res$params$t0,
v = fit_res$params$v,
a2 = fit_res$params$a2,
b2 = fit_res$params$b2,
k2 = fit_res$params$k2,
t02 = fit_res$params$t02,
v2 = fit_res$params$v2,
edf = fit_res$params$edf,
span = fit_res$params$span,
spline_df = fit_res$params$spline_df,
spar = fit_res$params$spar,
model_note = fit_res$params$model_note
)
row_id <- row_id + 1L
}
}
parameter_table <- dplyr::bind_rows(param_rows)
stats_cols <- c(
"series", "fit_id", "year_mode", "method", "site_mode",
"season_start", "season_end", "included_seasons", "season_length",
"n_obs", "n_days_observed", "first_obs_day", "last_obs_day",
"n_missing_days", "extrapolated", "anchor_added",
"growth_start_day", "growth_end_day",
"growth_start_season_day", "growth_end_season_day",
"growth_start_date", "growth_end_date",
"peak_rate",
"rate_start_day", "rate_end_day",
"rate_start_season_day", "rate_end_season_day",
"rate_start_date", "rate_end_date",
"fallback_used",
"two_pulse_detected",
"peak1_day", "peak2_day",
"peak1_season_day", "peak2_season_day",
"separator_day", "separator_season_day", "separator_date",
"valley_rate",
"pulse1_peak_rate", "pulse2_peak_rate",
"pulse1_start_day", "pulse1_end_day",
"pulse1_start_season_day", "pulse1_end_season_day",
"pulse1_start_date", "pulse1_end_date",
"pulse2_start_day", "pulse2_end_day",
"pulse2_start_season_day", "pulse2_end_season_day",
"pulse2_start_date", "pulse2_end_date"
)
param_cols <- c(
"series", "fit_id", "year_mode", "method", "site_mode",
"fallback_used", "converged",
"fixed_a_used", "fixed_a_value",
"a", "b", "k", "t0", "v",
"a2", "b2", "k2", "t02", "v2",
"edf", "span", "spline_df", "spar", "model_note"
)
fit_statistics <- parameter_table %>%
dplyr::select(dplyr::any_of(stats_cols))
fit_parameters <- parameter_table %>%
dplyr::select(dplyr::any_of(param_cols))
out <- list(
call = cl,
original_daily_data = original_daily_data,
processed_data = dat_proc,
fitted_data = fitted_dat,
fit_statistics = fit_statistics,
fit_parameters = fit_parameters,
season_table = season_table
)
class(out) <- "dm_growth_fit"
if (isTRUE(verbose)) {
message(
"dm.growth.fit.double completed: ",
nrow(fit_parameters), " fit(s) across ",
length(unique(fit_parameters$series)), " series."
)
}
out
}
# helpers ----------------------------------------------------------------------
#' Validate inputs for double dendrometer growth fitting
#'
#' @param df Input data frame.
#' @param growth_fraction Numeric vector of two growth fractions.
#' @param min_unique_growth Minimum number of unique growth values.
#' @param custom_veg_season Custom vegetation-season DOY range.
#' @param rate_threshold_fraction Fraction of peak growth rate.
#' @param peak_separation_min Minimum separation between derivative peaks.
#' @param valley_ratio_max Maximum allowed valley-to-peak ratio.
#' @param min_relative_peak Minimum relative derivative peak height.
#' @param fallback_to_single Logical fallback flag.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param fixed_a_multiplier Multiplier for fixed total asymptote.
#'
#' @return Invisibly returns \code{TRUE} if checks pass.
#'
#' @keywords internal
dmgfd_validate_inputs <- function(df,
growth_fraction,
min_unique_growth,
custom_veg_season,
rate_threshold_fraction,
peak_separation_min,
valley_ratio_max,
min_relative_peak,
fallback_to_single,
fix_a_to_observed_max,
fixed_a_multiplier) {
if (!is.data.frame(df)) {
stop("df must be a data.frame or tibble.")
}
if (ncol(df) < 2) {
stop("df must contain a time column and at least one dendrometer series.")
}
if (!is.numeric(growth_fraction) ||
length(growth_fraction) != 2 ||
any(!is.finite(growth_fraction)) ||
growth_fraction[1] < 0 ||
growth_fraction[2] > 1 ||
growth_fraction[1] >= growth_fraction[2]) {
stop("growth_fraction must be a numeric vector like c(0.1, 0.9).")
}
if (!is.numeric(min_unique_growth) ||
length(min_unique_growth) != 1 ||
!is.finite(min_unique_growth) ||
min_unique_growth < 2) {
stop("min_unique_growth must be a single finite number >= 2.")
}
if (!is.numeric(custom_veg_season) ||
length(custom_veg_season) != 2 ||
any(!is.finite(custom_veg_season))) {
stop("custom_veg_season must be a numeric vector of length 2.")
}
if (!is.numeric(rate_threshold_fraction) ||
length(rate_threshold_fraction) != 1 ||
!is.finite(rate_threshold_fraction) ||
rate_threshold_fraction <= 0 ||
rate_threshold_fraction >= 1) {
stop("rate_threshold_fraction must be a single number between 0 and 1.")
}
if (!is.numeric(peak_separation_min) ||
length(peak_separation_min) != 1 ||
!is.finite(peak_separation_min) ||
peak_separation_min < 1) {
stop("peak_separation_min must be a single number >= 1.")
}
if (!is.numeric(valley_ratio_max) ||
length(valley_ratio_max) != 1 ||
!is.finite(valley_ratio_max) ||
valley_ratio_max <= 0 ||
valley_ratio_max >= 1) {
stop("valley_ratio_max must be a single number between 0 and 1.")
}
if (!is.numeric(min_relative_peak) ||
length(min_relative_peak) != 1 ||
!is.finite(min_relative_peak) ||
min_relative_peak <= 0 ||
min_relative_peak >= 1) {
stop("min_relative_peak must be a single number between 0 and 1.")
}
if (!is.logical(fallback_to_single) ||
length(fallback_to_single) != 1 ||
is.na(fallback_to_single)) {
stop("fallback_to_single must be TRUE or FALSE.")
}
if (!is.logical(fix_a_to_observed_max) ||
length(fix_a_to_observed_max) != 1 ||
is.na(fix_a_to_observed_max)) {
stop("fix_a_to_observed_max must be TRUE or FALSE.")
}
if (!is.numeric(fixed_a_multiplier) ||
length(fixed_a_multiplier) != 1 ||
!is.finite(fixed_a_multiplier) ||
fixed_a_multiplier <= 0) {
stop("fixed_a_multiplier must be a single positive finite number.")
}
invisible(TRUE)
}
#' Check whether a value is a finite numeric scalar
#'
#' @param x Object to test.
#'
#' @return Logical scalar.
#'
#' @keywords internal
dmgfd_is_scalar_finite <- function(x) {
is.numeric(x) && length(x) == 1 && !is.na(x) && is.finite(x)
}
#' Return first non-missing value
#'
#' @param x Numeric vector.
#'
#' @return First non-missing value or \code{NA_real_}.
#'
#' @keywords internal
dmgfd_first_non_na <- function(x) {
idx <- which(!is.na(x))[1]
if (length(idx) == 0 || is.na(idx)) {
return(NA_real_)
}
x[idx]
}
#' Cumulative maximum while preserving missing values
#'
#' @param x Numeric vector.
#'
#' @return Numeric vector.
#'
#' @keywords internal
dmgfd_cummax_na <- function(x) {
out <- x
idx <- which(!is.na(x))
if (length(idx) > 0) {
out[idx] <- cummax(x[idx])
}
out
}
#' Convert year and day-of-year to Date
#'
#' @param year Calendar year.
#' @param doy Day of year.
#'
#' @return Date.
#'
#' @keywords internal
dmgfd_doy_to_date <- function(year, doy) {
as.Date(sprintf("%04d-01-01", year)) + (doy - 1)
}
#' Build vegetation-season table
#'
#' @param time_vec Time vector.
#' @param site_mode One of \code{"NH"}, \code{"SH"}, or \code{"CS"}.
#' @param custom_veg_season Custom vegetation-season DOY range.
#'
#' @return Tibble with season labels and season timing.
#'
#' @keywords internal
dmgfd_build_season_table <- function(time_vec,
site_mode = c("NH", "SH", "CS"),
custom_veg_season = c(275, 274)) {
site_mode <- match.arg(site_mode)
dates <- as.Date(time_vec)
yr <- lubridate::year(dates)
doy <- lubridate::yday(dates)
mo <- lubridate::month(dates)
n <- length(dates)
season_label <- rep(NA_character_, n)
season_start <- as.Date(rep(NA_character_, n))
season_end <- as.Date(rep(NA_character_, n))
season_day <- rep(NA_integer_, n)
if (site_mode == "NH") {
season_start <- as.Date(sprintf("%04d-01-01", yr))
season_end <- as.Date(sprintf("%04d-12-31", yr))
season_label <- as.character(yr)
season_day <- as.integer(dates - season_start) + 1L
}
if (site_mode == "SH") {
start_year <- ifelse(mo >= 7, yr, yr - 1L)
end_year <- start_year + 1L
season_start <- as.Date(sprintf("%04d-07-01", start_year))
season_end <- as.Date(sprintf("%04d-06-30", end_year))
season_label <- sprintf("%04d/%04d", start_year, end_year)
season_day <- as.integer(dates - season_start) + 1L
}
if (site_mode == "CS") {
start_doy <- as.integer(custom_veg_season[1])
end_doy <- as.integer(custom_veg_season[2])
if (start_doy <= end_doy) {
inside <- doy >= start_doy & doy <= end_doy
season_start[inside] <- dmgfd_doy_to_date(yr[inside], start_doy)
season_end[inside] <- dmgfd_doy_to_date(yr[inside], end_doy)
season_label[inside] <- as.character(yr[inside])
season_day[inside] <- as.integer(dates[inside] - season_start[inside]) + 1L
} else {
inside_late <- doy >= start_doy
inside_early <- doy <= end_doy
if (any(inside_late)) {
sy <- yr[inside_late]
ey <- sy + 1L
season_start[inside_late] <- dmgfd_doy_to_date(sy, start_doy)
season_end[inside_late] <- dmgfd_doy_to_date(ey, end_doy)
season_label[inside_late] <- sprintf("%04d/%04d", sy, ey)
season_day[inside_late] <- as.integer(dates[inside_late] - season_start[inside_late]) + 1L
}
if (any(inside_early)) {
ey <- yr[inside_early]
sy <- ey - 1L
season_start[inside_early] <- dmgfd_doy_to_date(sy, start_doy)
season_end[inside_early] <- dmgfd_doy_to_date(ey, end_doy)
season_label[inside_early] <- sprintf("%04d/%04d", sy, ey)
season_day[inside_early] <- as.integer(dates[inside_early] - season_start[inside_early]) + 1L
}
}
}
tibble::tibble(
season_label = season_label,
season_start = season_start,
season_end = season_end,
season_day = season_day
)
}
#' Add zero-growth anchor point
#'
#' @param x Season-day values.
#' @param y Growth values.
#'
#' @return List with \code{x} and \code{y}.
#'
#' @keywords internal
dmgfd_add_anchor_points <- function(x, y) {
if (length(x) == 0) {
return(list(x = x, y = y))
}
if (min(x, na.rm = TRUE) > 1 && !any(x == 1)) {
x <- c(1, x)
y <- c(0, y)
}
ord <- order(x)
list(
x = x[ord],
y = y[ord]
)
}
#' Empty parameter template for double growth fitting
#'
#' @return Named list of model and timing parameters.
#'
#' @keywords internal
dmgfd_empty_model_params <- function() {
list(
converged = FALSE,
fallback_used = FALSE,
fixed_a_used = FALSE,
fixed_a_value = NA_real_,
method_used = NA_character_,
a = NA_real_,
b = NA_real_,
k = NA_real_,
t0 = NA_real_,
v = NA_real_,
a2 = NA_real_,
b2 = NA_real_,
k2 = NA_real_,
t02 = NA_real_,
v2 = NA_real_,
edf = NA_real_,
span = NA_real_,
spline_df = NA_real_,
spar = NA_real_,
model_note = NA_character_,
n_obs = NA_real_,
n_days_observed = NA_real_,
first_obs_day = NA_real_,
last_obs_day = NA_real_,
season_length = NA_real_,
n_missing_days = NA_real_,
extrapolated = NA,
anchor_added = FALSE,
growth_start_day = NA_real_,
growth_end_day = NA_real_,
growth_start_season_day = NA_real_,
growth_end_season_day = NA_real_,
growth_start_date = as.Date(NA),
growth_end_date = as.Date(NA),
peak_rate = NA_real_,
rate_start_day = NA_real_,
rate_end_day = NA_real_,
rate_start_season_day = NA_real_,
rate_end_season_day = NA_real_,
rate_start_date = as.Date(NA),
rate_end_date = as.Date(NA),
two_pulse_detected = FALSE,
peak1_day = NA_real_,
peak2_day = NA_real_,
peak1_season_day = NA_real_,
peak2_season_day = NA_real_,
separator_day = NA_real_,
separator_season_day = NA_real_,
separator_date = as.Date(NA),
valley_rate = NA_real_,
pulse1_peak_rate = NA_real_,
pulse2_peak_rate = NA_real_,
pulse1_start_day = NA_real_,
pulse1_end_day = NA_real_,
pulse1_start_season_day = NA_real_,
pulse1_end_season_day = NA_real_,
pulse1_start_date = as.Date(NA),
pulse1_end_date = as.Date(NA),
pulse2_start_day = NA_real_,
pulse2_end_day = NA_real_,
pulse2_start_season_day = NA_real_,
pulse2_end_season_day = NA_real_,
pulse2_start_date = as.Date(NA),
pulse2_end_date = as.Date(NA)
)
}
#' Guess two pulse locations from growth increments
#'
#' @param x Season-day values.
#' @param y Cumulative growth values.
#'
#' @return Numeric vector of two initial pulse locations.
#'
#' @keywords internal
dmgfd_guess_two_pulses <- function(x, y) {
x <- as.numeric(x)
y <- as.numeric(y)
ord <- order(x)
x <- x[ord]
y <- y[ord]
dy <- c(NA_real_, diff(y))
dy[!is.finite(dy)] <- -Inf
min_sep <- max(14, round(0.15 * diff(range(x, na.rm = TRUE))))
o <- order(dy, decreasing = TRUE)
if (length(o) < 2 || !is.finite(dy[o[1]])) {
qs <- stats::quantile(x, probs = c(0.33, 0.67), na.rm = TRUE, names = FALSE)
return(as.numeric(qs))
}
t1 <- x[o[1]]
t2 <- NA_real_
if (length(o) > 1) {
for (j in o[-1]) {
if (is.finite(dy[j]) && abs(x[j] - t1) >= min_sep) {
t2 <- x[j]
break
}
}
}
if (!is.finite(t2)) {
qs <- stats::quantile(x, probs = c(0.33, 0.67), na.rm = TRUE, names = FALSE)
t1 <- qs[1]
t2 <- qs[2]
}
sort(c(t1, t2))
}
#' Infer single Gompertz starting values
#'
#' @param x Season-day values.
#' @param y Growth values.
#' @param a0 Initial asymptote.
#'
#' @return Named list with \code{a}, \code{b}, and \code{k}.
#'
#' @keywords internal
dmgfd_infer_gompertz_starts <- function(x, y, a0) {
frac <- pmin(pmax(y / a0, 1e-6), 1 - 1e-6)
z <- log(-log(frac))
ok <- is.finite(x) & is.finite(z)
b0 <- 0.5
k0 <- 0.01
if (sum(ok) >= 2) {
lm_fit <- try(stats::lm(z[ok] ~ x[ok]), silent = TRUE)
if (!inherits(lm_fit, "try-error")) {
cf <- stats::coef(lm_fit)
if (length(cf) == 2 && all(is.finite(cf))) {
b_try <- unname(cf[1])
k_try <- -unname(cf[2])
if (is.finite(b_try)) {
b0 <- b_try
}
if (is.finite(k_try) && k_try > 1e-4) {
k0 <- k_try
}
}
}
}
list(
a = a0,
b = b0,
k = k0
)
}
#' Infer logistic/Richards starting values
#'
#' @param x Season-day values.
#' @param y Growth values.
#' @param a0 Initial asymptote.
#'
#' @return Named list with \code{a}, \code{k}, and \code{t0}.
#'
#' @keywords internal
dmgfd_infer_logistic_starts <- function(x, y, a0) {
frac <- pmin(pmax(y / a0, 1e-6), 1 - 1e-6)
z <- stats::qlogis(frac)
ok <- is.finite(x) & is.finite(z)
k0 <- 0.03
t00 <- stats::median(x, na.rm = TRUE)
if (sum(ok) >= 2) {
lm_fit <- try(stats::lm(z[ok] ~ x[ok]), silent = TRUE)
if (!inherits(lm_fit, "try-error")) {
cf <- stats::coef(lm_fit)
if (length(cf) == 2 && all(is.finite(cf))) {
k_try <- abs(unname(cf[2]))
if (is.finite(k_try) && k_try > 1e-4) {
k0 <- k_try
t0_try <- -unname(cf[1]) / k0
if (is.finite(t0_try)) {
t00 <- t0_try
}
}
}
}
}
list(
a = a0,
k = k0,
t0 = t00
)
}
#' Infer double-Gompertz starting values
#'
#' @param x Season-day values.
#' @param y Growth values.
#' @param a0 Initial total asymptote.
#'
#' @return Named list of starting values.
#'
#' @keywords internal
dmgfd_infer_double_gompertz_starts <- function(x, y, a0) {
pulses <- dmgfd_guess_two_pulses(x, y)
base <- dmgfd_infer_gompertz_starts(x, y, a0)
k0 <- if (dmgfd_is_scalar_finite(base$k)) base$k else 0.01
list(
a = a0 * 0.5,
b = k0 * pulses[1],
k = k0,
t0 = pulses[1],
a2 = a0 * 0.5,
b2 = k0 * pulses[2],
k2 = k0,
t02 = pulses[2]
)
}
#' Infer double-Richards starting values
#'
#' @param x Season-day values.
#' @param y Growth values.
#' @param a0 Initial total asymptote.
#'
#' @return Named list of starting values.
#'
#' @keywords internal
dmgfd_infer_double_richards_starts <- function(x, y, a0) {
pulses <- dmgfd_guess_two_pulses(x, y)
base <- dmgfd_infer_logistic_starts(x, y, a0)
k0 <- if (dmgfd_is_scalar_finite(base$k)) base$k else 0.03
list(
a = a0 * 0.5,
k = k0,
t0 = pulses[1],
v = 1,
a2 = a0 * 0.5,
k2 = k0,
t02 = pulses[2],
v2 = 1
)
}
#' Fit double-Gompertz model
#'
#' @param x Observed season-day values.
#' @param y Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param fixed_a_value Optional fixed total asymptote on original scale.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param start_value_double_gompertz_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgfd_fit_model_double_gompertz <- function(
x,
y,
x_pred,
fixed_a_value = NA_real_,
fix_a_to_observed_max = FALSE,
start_value_double_gompertz_parameters = list(
a = NA_real_, k = NA_real_, t0 = NA_real_,
a2 = NA_real_, k2 = NA_real_, t02 = NA_real_
)
) {
params <- dmgfd_empty_model_params()
res <- tryCatch({
cons <- ifelse(
min(y, na.rm = TRUE) < 0,
abs(min(y, na.rm = TRUE)) + 1e-8,
0
)
y_adj <- y + cons
a0 <- max(y_adj, na.rm = TRUE) * 1.05 + 1e-6
start_guess <- dmgfd_infer_double_gompertz_starts(x, y_adj, a0)
start_guess <- utils::modifyList(start_guess, start_value_double_gompertz_parameters)
if (!dmgfd_is_scalar_finite(start_guess$a)) {
start_guess$a <- a0 * 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k)) {
start_guess$k <- 0.01
}
if (!dmgfd_is_scalar_finite(start_guess$t0)) {
start_guess$t0 <- stats::quantile(x, 0.33, na.rm = TRUE, names = FALSE)
}
if (!dmgfd_is_scalar_finite(start_guess$a2)) {
start_guess$a2 <- a0 * 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k2)) {
start_guess$k2 <- start_guess$k
}
if (!dmgfd_is_scalar_finite(start_guess$t02)) {
start_guess$t02 <- stats::quantile(x, 0.67, na.rm = TRUE, names = FALSE)
}
dt_start <- max(5, start_guess$t02 - start_guess$t0)
dat <- data.frame(
x = x,
y = y_adj
)
use_fixed_a <- isTRUE(fix_a_to_observed_max) &&
dmgfd_is_scalar_finite(fixed_a_value) &&
fixed_a_value > 0
if (isTRUE(use_fixed_a)) {
a_fixed_original <- fixed_a_value
a_fixed_fit <- fixed_a_value + cons
p_start <- 0.5
if (dmgfd_is_scalar_finite(start_guess$a) &&
dmgfd_is_scalar_finite(start_guess$a2) &&
(start_guess$a + start_guess$a2) > 0) {
p_start <- start_guess$a / (start_guess$a + start_guess$a2)
}
p_start <- min(max(p_start, 0.05), 0.95)
mod <- minpack.lm::nlsLM(
y ~ a_fixed_fit * p * exp(-exp(-(k * (x - t0)))) +
a_fixed_fit * (1 - p) * exp(-exp(-(k2 * (x - (t0 + dt))))),
data = dat,
start = list(
p = p_start,
k = start_guess$k,
t0 = start_guess$t0,
k2 = start_guess$k2,
dt = dt_start
),
lower = c(
p = 0.001,
k = 1e-6,
t0 = min(x, na.rm = TRUE),
k2 = 1e-6,
dt = 1
),
upper = c(
p = 0.999,
k = 5,
t0 = max(x, na.rm = TRUE),
k2 = 5,
dt = diff(range(x, na.rm = TRUE))
),
control = minpack.lm::nls.lm.control(maxiter = 1000)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- pmax(as.numeric(pred), 0)
pred <- cummax(pred)
cf <- stats::coef(mod)
p_est <- unname(cf["p"])
params$a <- a_fixed_original * p_est
params$a2 <- a_fixed_original * (1 - p_est)
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$b <- params$k * params$t0
params$k2 <- unname(cf["k2"])
params$t02 <- unname(cf["t0"] + cf["dt"])
params$b2 <- params$k2 * params$t02
params$converged <- TRUE
params$fixed_a_used <- TRUE
params$fixed_a_value <- a_fixed_original
params$method_used <- "double_gompertz"
params$model_note <- "Double-Gompertz total asymptote a + a2 fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
mod <- minpack.lm::nlsLM(
y ~ a * exp(-exp(-(k * (x - t0)))) +
a2 * exp(-exp(-(k2 * (x - (t0 + dt))))),
data = dat,
start = list(
a = start_guess$a,
k = start_guess$k,
t0 = start_guess$t0,
a2 = start_guess$a2,
k2 = start_guess$k2,
dt = dt_start
),
lower = c(
a = 0,
k = 1e-6,
t0 = min(x, na.rm = TRUE),
a2 = 0,
k2 = 1e-6,
dt = 1
),
upper = c(
a = Inf,
k = 5,
t0 = max(x, na.rm = TRUE),
a2 = Inf,
k2 = 5,
dt = diff(range(x, na.rm = TRUE))
),
control = minpack.lm::nls.lm.control(maxiter = 1000)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- pmax(as.numeric(pred), 0)
pred <- cummax(pred)
cf <- stats::coef(mod)
params$a <- unname(cf["a"])
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$b <- params$k * params$t0
params$a2 <- unname(cf["a2"])
params$k2 <- unname(cf["k2"])
params$t02 <- unname(cf["t0"] + cf["dt"])
params$b2 <- params$k2 * params$t02
params$converged <- TRUE
params$method_used <- "double_gompertz"
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
params$method_used <- "double_gompertz"
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit double-Richards model
#'
#' @param x Observed season-day values.
#' @param y Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param fixed_a_value Optional fixed total asymptote on original scale.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param start_value_double_richards_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgfd_fit_model_double_richards <- function(
x,
y,
x_pred,
fixed_a_value = NA_real_,
fix_a_to_observed_max = FALSE,
start_value_double_richards_parameters = list(
a = NA_real_, k = NA_real_, t0 = NA_real_, v = 1,
a2 = NA_real_, k2 = NA_real_, t02 = NA_real_, v2 = 1
)
) {
params <- dmgfd_empty_model_params()
res <- tryCatch({
cons <- ifelse(
min(y, na.rm = TRUE) < 0,
abs(min(y, na.rm = TRUE)) + 1e-8,
0
)
y_adj <- y + cons
a0 <- max(y_adj, na.rm = TRUE) * 1.05 + 1e-6
start_guess <- dmgfd_infer_double_richards_starts(x, y_adj, a0)
start_guess <- utils::modifyList(start_guess, start_value_double_richards_parameters)
if (!dmgfd_is_scalar_finite(start_guess$a)) {
start_guess$a <- a0 * 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k)) {
start_guess$k <- 0.03
}
if (!dmgfd_is_scalar_finite(start_guess$t0)) {
start_guess$t0 <- stats::quantile(x, 0.33, na.rm = TRUE, names = FALSE)
}
if (!dmgfd_is_scalar_finite(start_guess$v)) {
start_guess$v <- 1
}
if (!dmgfd_is_scalar_finite(start_guess$a2)) {
start_guess$a2 <- a0 * 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k2)) {
start_guess$k2 <- start_guess$k
}
if (!dmgfd_is_scalar_finite(start_guess$t02)) {
start_guess$t02 <- stats::quantile(x, 0.67, na.rm = TRUE, names = FALSE)
}
if (!dmgfd_is_scalar_finite(start_guess$v2)) {
start_guess$v2 <- 1
}
dt_start <- max(5, start_guess$t02 - start_guess$t0)
dat <- data.frame(
x = x,
y = y_adj
)
use_fixed_a <- isTRUE(fix_a_to_observed_max) &&
dmgfd_is_scalar_finite(fixed_a_value) &&
fixed_a_value > 0
if (isTRUE(use_fixed_a)) {
a_fixed_original <- fixed_a_value
a_fixed_fit <- fixed_a_value + cons
p_start <- 0.5
if (dmgfd_is_scalar_finite(start_guess$a) &&
dmgfd_is_scalar_finite(start_guess$a2) &&
(start_guess$a + start_guess$a2) > 0) {
p_start <- start_guess$a / (start_guess$a + start_guess$a2)
}
p_start <- min(max(p_start, 0.05), 0.95)
mod <- minpack.lm::nlsLM(
y ~ a_fixed_fit * p / ((1 + v * exp(-k * (x - t0)))^(1 / v)) +
a_fixed_fit * (1 - p) / ((1 + v2 * exp(-k2 * (x - (t0 + dt))))^(1 / v2)),
data = dat,
start = list(
p = p_start,
k = start_guess$k,
t0 = start_guess$t0,
v = start_guess$v,
k2 = start_guess$k2,
dt = dt_start,
v2 = start_guess$v2
),
lower = c(
p = 0.001,
k = 1e-6,
t0 = min(x, na.rm = TRUE),
v = 1e-3,
k2 = 1e-6,
dt = 1,
v2 = 1e-3
),
upper = c(
p = 0.999,
k = 5,
t0 = max(x, na.rm = TRUE),
v = 50,
k2 = 5,
dt = diff(range(x, na.rm = TRUE)),
v2 = 50
),
control = minpack.lm::nls.lm.control(maxiter = 1000)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- pmax(as.numeric(pred), 0)
pred <- cummax(pred)
cf <- stats::coef(mod)
p_est <- unname(cf["p"])
params$a <- a_fixed_original * p_est
params$a2 <- a_fixed_original * (1 - p_est)
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$v <- unname(cf["v"])
params$k2 <- unname(cf["k2"])
params$t02 <- unname(cf["t0"] + cf["dt"])
params$v2 <- unname(cf["v2"])
params$converged <- TRUE
params$fixed_a_used <- TRUE
params$fixed_a_value <- a_fixed_original
params$method_used <- "double_richards"
params$model_note <- "Double-Richards total asymptote a + a2 fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
mod <- minpack.lm::nlsLM(
y ~ a / ((1 + v * exp(-k * (x - t0)))^(1 / v)) +
a2 / ((1 + v2 * exp(-k2 * (x - (t0 + dt))))^(1 / v2)),
data = dat,
start = list(
a = start_guess$a,
k = start_guess$k,
t0 = start_guess$t0,
v = start_guess$v,
a2 = start_guess$a2,
k2 = start_guess$k2,
dt = dt_start,
v2 = start_guess$v2
),
lower = c(
a = 0,
k = 1e-6,
t0 = min(x, na.rm = TRUE),
v = 1e-3,
a2 = 0,
k2 = 1e-6,
dt = 1,
v2 = 1e-3
),
upper = c(
a = Inf,
k = 5,
t0 = max(x, na.rm = TRUE),
v = 50,
a2 = Inf,
k2 = 5,
dt = diff(range(x, na.rm = TRUE)),
v2 = 50
),
control = minpack.lm::nls.lm.control(maxiter = 1000)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- pmax(as.numeric(pred), 0)
pred <- cummax(pred)
cf <- stats::coef(mod)
params$a <- unname(cf["a"])
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$v <- unname(cf["v"])
params$a2 <- unname(cf["a2"])
params$k2 <- unname(cf["k2"])
params$t02 <- unname(cf["t0"] + cf["dt"])
params$v2 <- unname(cf["v2"])
params$converged <- TRUE
params$method_used <- "double_richards"
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
params$method_used <- "double_richards"
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit single-Gompertz fallback model
#'
#' @param x Observed season-day values.
#' @param y Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param fixed_a_value Optional fixed asymptote.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param start_value_gompertz_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgfd_fit_model_single_gompertz <- function(
x,
y,
x_pred,
fixed_a_value = NA_real_,
fix_a_to_observed_max = FALSE,
start_value_gompertz_parameters = list(a = NA_real_, b = NA_real_, k = NA_real_)
) {
params <- dmgfd_empty_model_params()
res <- tryCatch({
cons <- ifelse(
min(y, na.rm = TRUE) < 0,
abs(min(y, na.rm = TRUE)) + 1e-8,
0
)
y_adj <- y + cons
a0 <- max(y_adj, na.rm = TRUE) * 1.05 + 1e-6
use_fixed_a <- isTRUE(fix_a_to_observed_max) &&
dmgfd_is_scalar_finite(fixed_a_value) &&
fixed_a_value > 0
if (isTRUE(use_fixed_a)) {
a_fixed_original <- fixed_a_value
a_fixed_fit <- fixed_a_value + cons
start_guess <- dmgfd_infer_gompertz_starts(x, y_adj, a_fixed_fit)
start_guess <- utils::modifyList(start_guess, start_value_gompertz_parameters)
if (!dmgfd_is_scalar_finite(start_guess$b)) {
start_guess$b <- 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k)) {
start_guess$k <- 0.01
}
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a_fixed_fit * exp(-exp(b - k * x)),
data = dat,
start = list(
b = start_guess$b,
k = start_guess$k
),
control = minpack.lm::nls.lm.control(maxiter = 500)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- as.numeric(pred)
cf <- stats::coef(mod)
params$a <- a_fixed_original
params$b <- unname(cf["b"])
params$k <- unname(cf["k"])
params$converged <- TRUE
params$fixed_a_used <- TRUE
params$fixed_a_value <- a_fixed_original
params$method_used <- "single_gompertz"
params$model_note <- "Single-Gompertz fallback asymptote a fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
start_guess <- dmgfd_infer_gompertz_starts(x, y_adj, a0)
start_guess <- utils::modifyList(start_guess, start_value_gompertz_parameters)
if (!dmgfd_is_scalar_finite(start_guess$a)) {
start_guess$a <- a0
}
if (!dmgfd_is_scalar_finite(start_guess$b)) {
start_guess$b <- 0.5
}
if (!dmgfd_is_scalar_finite(start_guess$k)) {
start_guess$k <- 0.01
}
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a * exp(-exp(b - k * x)),
data = dat,
start = start_guess,
control = minpack.lm::nls.lm.control(maxiter = 500)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- as.numeric(pred)
cf <- stats::coef(mod)
params$a <- unname(cf["a"])
params$b <- unname(cf["b"])
params$k <- unname(cf["k"])
params$converged <- TRUE
params$method_used <- "single_gompertz"
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
params$method_used <- "single_gompertz"
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit single-Richards fallback model
#'
#' @param x Observed season-day values.
#' @param y Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param fixed_a_value Optional fixed asymptote.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param start_value_richards_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgfd_fit_model_single_richards <- function(
x,
y,
x_pred,
fixed_a_value = NA_real_,
fix_a_to_observed_max = FALSE,
start_value_richards_parameters = list(a = NA_real_, k = NA_real_, t0 = NA_real_, v = 1)
) {
params <- dmgfd_empty_model_params()
res <- tryCatch({
cons <- ifelse(
min(y, na.rm = TRUE) < 0,
abs(min(y, na.rm = TRUE)) + 1e-8,
0
)
y_adj <- y + cons
a0 <- max(y_adj, na.rm = TRUE) * 1.05 + 1e-6
use_fixed_a <- isTRUE(fix_a_to_observed_max) &&
dmgfd_is_scalar_finite(fixed_a_value) &&
fixed_a_value > 0
if (isTRUE(use_fixed_a)) {
a_fixed_original <- fixed_a_value
a_fixed_fit <- fixed_a_value + cons
logi_start <- dmgfd_infer_logistic_starts(x, y_adj, a_fixed_fit)
start_vals <- list(
k = if (dmgfd_is_scalar_finite(start_value_richards_parameters$k)) {
start_value_richards_parameters$k
} else {
logi_start$k
},
t0 = if (dmgfd_is_scalar_finite(start_value_richards_parameters$t0)) {
start_value_richards_parameters$t0
} else {
logi_start$t0
},
v = if (dmgfd_is_scalar_finite(start_value_richards_parameters$v)) {
start_value_richards_parameters$v
} else {
1
}
)
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a_fixed_fit / ((1 + v * exp(-k * (x - t0)))^(1 / v)),
data = dat,
start = start_vals,
control = minpack.lm::nls.lm.control(maxiter = 500)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- as.numeric(pred)
cf <- stats::coef(mod)
params$a <- a_fixed_original
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$v <- unname(cf["v"])
params$converged <- TRUE
params$fixed_a_used <- TRUE
params$fixed_a_value <- a_fixed_original
params$method_used <- "single_richards"
params$model_note <- "Single-Richards fallback asymptote a fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
logi_start <- dmgfd_infer_logistic_starts(x, y_adj, a0)
start_vals <- list(
a = if (dmgfd_is_scalar_finite(start_value_richards_parameters$a)) {
start_value_richards_parameters$a
} else {
logi_start$a
},
k = if (dmgfd_is_scalar_finite(start_value_richards_parameters$k)) {
start_value_richards_parameters$k
} else {
logi_start$k
},
t0 = if (dmgfd_is_scalar_finite(start_value_richards_parameters$t0)) {
start_value_richards_parameters$t0
} else {
logi_start$t0
},
v = if (dmgfd_is_scalar_finite(start_value_richards_parameters$v)) {
start_value_richards_parameters$v
} else {
1
}
)
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a / ((1 + v * exp(-k * (x - t0)))^(1 / v)),
data = dat,
start = start_vals,
control = minpack.lm::nls.lm.control(maxiter = 500)
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred)) - cons
pred <- as.numeric(pred)
cf <- stats::coef(mod)
params$a <- unname(cf["a"])
params$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$v <- unname(cf["v"])
params$converged <- TRUE
params$method_used <- "single_richards"
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
params$method_used <- "single_richards"
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Add overall cumulative-growth timing
#'
#' @param params Parameter list.
#' @param pred_full Fitted cumulative growth.
#' @param x_pred Prediction season-day values.
#' @param growth_fraction Growth fractions.
#' @param season_start_date Season start date.
#'
#' @return Updated parameter list.
#'
#' @keywords internal
dmgfd_add_overall_timing <- function(params,
pred_full,
x_pred,
growth_fraction,
season_start_date = as.Date(NA)) {
valid_pred <- is.finite(pred_full) & is.finite(x_pred)
if (sum(valid_pred) == 0) {
return(params)
}
pred_use <- pred_full[valid_pred]
doy_use <- x_pred[valid_pred]
final_val <- tail(stats::na.omit(pred_use), 1)
if (length(final_val) == 1 &&
is.finite(final_val) &&
final_val > 0) {
low_val <- final_val * growth_fraction[1]
high_val <- final_val * growth_fraction[2]
st_pos <- suppressWarnings(min(doy_use[pred_use > low_val], na.rm = TRUE))
en_pos <- suppressWarnings(max(doy_use[pred_use < high_val], na.rm = TRUE))
if (is.finite(st_pos)) {
params$growth_start_season_day <- st_pos
if (!is.na(season_start_date)) {
params$growth_start_date <- season_start_date + (st_pos - 1)
params$growth_start_day <- lubridate::yday(params$growth_start_date)
}
}
if (is.finite(en_pos)) {
params$growth_end_season_day <- en_pos
if (!is.na(season_start_date)) {
params$growth_end_date <- season_start_date + (en_pos - 1)
params$growth_end_day <- lubridate::yday(params$growth_end_date)
}
}
}
params
}
#' Add overall rate-based timing
#'
#' @param params Parameter list.
#' @param pred_full Fitted cumulative growth.
#' @param x_pred Prediction season-day values.
#' @param rate_threshold_fraction Fraction of peak rate.
#' @param season_start_date Season start date.
#'
#' @return Updated parameter list.
#'
#' @keywords internal
dmgfd_add_overall_rate_timing <- function(params,
pred_full,
x_pred,
rate_threshold_fraction = 0.1,
season_start_date = as.Date(NA)) {
rate <- dmgfd_growth_rate(pred_full)
ok <- is.finite(rate) & is.finite(x_pred)
if (sum(ok) == 0) {
return(params)
}
r <- rate[ok]
x <- x_pred[ok]
peak_rate <- suppressWarnings(max(r, na.rm = TRUE))
if (!is.finite(peak_rate) || peak_rate <= 0) {
return(params)
}
thr <- rate_threshold_fraction * peak_rate
above <- which(r >= thr)
if (length(above) == 0) {
return(params)
}
params$peak_rate <- peak_rate
params$rate_start_season_day <- min(x[above], na.rm = TRUE)
params$rate_end_season_day <- max(x[above], na.rm = TRUE)
if (!is.na(season_start_date)) {
params$rate_start_date <- season_start_date + (params$rate_start_season_day - 1)
params$rate_end_date <- season_start_date + (params$rate_end_season_day - 1)
params$rate_start_day <- lubridate::yday(params$rate_start_date)
params$rate_end_day <- lubridate::yday(params$rate_end_date)
}
params
}
#' Fit single-curve fallback
#'
#' @param method Original double method.
#' @param x Observed season-day values.
#' @param y Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param fixed_a_value Optional fixed asymptote.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param start_value_double_gompertz_parameters Double-Gompertz starts.
#' @param start_value_double_richards_parameters Double-Richards starts.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgfd_fit_single_fallback <- function(method,
x,
y,
x_pred,
fixed_a_value = NA_real_,
fix_a_to_observed_max = FALSE,
start_value_double_gompertz_parameters,
start_value_double_richards_parameters) {
params <- dmgfd_empty_model_params()
if (method == "gompertz") {
start_single <- list(
a = if (dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$a) &&
dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$a2)) {
start_value_double_gompertz_parameters$a +
start_value_double_gompertz_parameters$a2
} else if (dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$a)) {
start_value_double_gompertz_parameters$a
} else {
NA_real_
},
b = if (dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$k) &&
dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$t0)) {
start_value_double_gompertz_parameters$k *
start_value_double_gompertz_parameters$t0
} else {
NA_real_
},
k = if (dmgfd_is_scalar_finite(start_value_double_gompertz_parameters$k)) {
start_value_double_gompertz_parameters$k
} else {
NA_real_
}
)
single_res <- dmgfd_fit_model_single_gompertz(
x = x,
y = y,
x_pred = x_pred,
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_gompertz_parameters = start_single
)
params$a <- single_res$params$a
params$b <- single_res$params$b
params$k <- single_res$params$k
params$converged <- single_res$params$converged
params$fixed_a_used <- single_res$params$fixed_a_used
params$fixed_a_value <- single_res$params$fixed_a_value
params$model_note <- single_res$params$model_note
params$method_used <- "single_gompertz"
params$fallback_used <- TRUE
return(list(pred = single_res$pred, params = params))
}
if (method == "richards") {
start_single <- list(
a = if (dmgfd_is_scalar_finite(start_value_double_richards_parameters$a) &&
dmgfd_is_scalar_finite(start_value_double_richards_parameters$a2)) {
start_value_double_richards_parameters$a +
start_value_double_richards_parameters$a2
} else if (dmgfd_is_scalar_finite(start_value_double_richards_parameters$a)) {
start_value_double_richards_parameters$a
} else {
NA_real_
},
k = if (dmgfd_is_scalar_finite(start_value_double_richards_parameters$k)) {
start_value_double_richards_parameters$k
} else {
NA_real_
},
t0 = if (dmgfd_is_scalar_finite(start_value_double_richards_parameters$t0)) {
start_value_double_richards_parameters$t0
} else {
NA_real_
},
v = if (dmgfd_is_scalar_finite(start_value_double_richards_parameters$v)) {
start_value_double_richards_parameters$v
} else {
1
}
)
single_res <- dmgfd_fit_model_single_richards(
x = x,
y = y,
x_pred = x_pred,
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_richards_parameters = start_single
)
params$a <- single_res$params$a
params$k <- single_res$params$k
params$t0 <- single_res$params$t0
params$v <- single_res$params$v
params$converged <- single_res$params$converged
params$fixed_a_used <- single_res$params$fixed_a_used
params$fixed_a_value <- single_res$params$fixed_a_value
params$model_note <- single_res$params$model_note
params$method_used <- "single_richards"
params$fallback_used <- TRUE
return(list(pred = single_res$pred, params = params))
}
stop("Unknown single fallback method: ", method)
}
#' Calculate non-negative fitted growth rate
#'
#' @param pred Fitted cumulative growth.
#'
#' @return Numeric growth-rate vector.
#'
#' @keywords internal
dmgfd_growth_rate <- function(pred) {
pred <- as.numeric(pred)
if (length(pred) < 2) {
return(rep(NA_real_, length(pred)))
}
rate <- c(NA_real_, diff(pred))
rate[!is.finite(rate)] <- NA_real_
rate <- pmax(rate, 0)
rate
}
#' Identify local maxima
#'
#' @param x Numeric vector.
#'
#' @return Integer indices of local maxima.
#'
#' @keywords internal
dmgfd_local_maxima <- function(x) {
x <- as.numeric(x)
n <- length(x)
if (n < 3) {
return(integer(0))
}
which(
is.finite(x[2:(n - 1)]) &
x[2:(n - 1)] >= x[1:(n - 2)] &
x[2:(n - 1)] >= x[3:n]
) + 1L
}
#' Detect two-pulse pattern from fitted growth rate
#'
#' @param rate Fitted growth-rate vector.
#' @param x_pred Prediction season-day values.
#' @param peak_separation_min Minimum separation between peaks.
#' @param valley_ratio_max Maximum allowed valley-to-peak ratio.
#' @param min_relative_peak Minimum relative peak height.
#'
#' @return List describing detected pulse pattern.
#'
#' @keywords internal
dmgfd_detect_two_pulse_pattern <- function(rate,
x_pred,
peak_separation_min = 10,
valley_ratio_max = 0.4,
min_relative_peak = 0.05) {
out <- list(
detected = FALSE,
peak1_day = NA_real_,
peak2_day = NA_real_,
separator_day = NA_real_,
valley_rate = NA_real_,
pulse1_peak_rate = NA_real_,
pulse2_peak_rate = NA_real_
)
ok <- is.finite(rate) & is.finite(x_pred)
if (sum(ok) < 5) {
return(out)
}
r <- as.numeric(rate)
x <- as.numeric(x_pred)
peak_idx <- dmgfd_local_maxima(r)
if (length(peak_idx) < 2) {
return(out)
}
global_peak <- max(r, na.rm = TRUE)
if (!is.finite(global_peak) || global_peak <= 0) {
return(out)
}
peak_idx <- peak_idx[r[peak_idx] >= min_relative_peak * global_peak]
if (length(peak_idx) < 2) {
return(out)
}
best_score <- -Inf
best <- NULL
for (i in seq_len(length(peak_idx) - 1)) {
for (j in (i + 1):length(peak_idx)) {
p1 <- peak_idx[i]
p2 <- peak_idx[j]
if ((x[p2] - x[p1]) < peak_separation_min) {
next
}
valley_idx <- p1:p2
valley_rate <- suppressWarnings(min(r[valley_idx], na.rm = TRUE))
if (!is.finite(valley_rate)) {
next
}
if (valley_rate > valley_ratio_max * min(r[p1], r[p2])) {
next
}
score <- r[p1] + r[p2] - valley_rate
if (score > best_score) {
sep_idx <- valley_idx[which.min(r[valley_idx])]
best_score <- score
best <- list(
p1 = p1,
p2 = p2,
sep = sep_idx,
valley_rate = valley_rate
)
}
}
}
if (is.null(best)) {
return(out)
}
out$detected <- TRUE
out$peak1_day <- x[best$p1]
out$peak2_day <- x[best$p2]
out$separator_day <- x[best$sep]
out$valley_rate <- best$valley_rate
out$pulse1_peak_rate <- r[best$p1]
out$pulse2_peak_rate <- r[best$p2]
out
}
#' Estimate pulse window from fitted growth rate
#'
#' @param rate Fitted growth-rate vector.
#' @param x_pred Prediction season-day values.
#' @param threshold Pulse-specific rate threshold.
#' @param domain_idx Indices belonging to one pulse domain.
#'
#' @return List with start and end season-day.
#'
#' @keywords internal
dmgfd_pulse_window_from_rate <- function(rate,
x_pred,
threshold,
domain_idx) {
out <- list(
start_day = NA_real_,
end_day = NA_real_
)
if (length(domain_idx) == 0) {
return(out)
}
r <- rate[domain_idx]
x <- x_pred[domain_idx]
ok <- is.finite(r) & is.finite(x)
if (sum(ok) == 0) {
return(out)
}
r <- r[ok]
x <- x[ok]
above <- which(r >= threshold)
if (length(above) == 0) {
return(out)
}
out$start_day <- min(x[above], na.rm = TRUE)
out$end_day <- max(x[above], na.rm = TRUE)
out
}
#' Fit one double growth curve
#'
#' @description
#' Internal helper used by [dm.growth.fit.double()] to fit one double curve for
#' one series and one season, or one pooled double curve.
#'
#' @param x_obs Observed season-day values.
#' @param y_obs Observed cumulative growth values.
#' @param x_pred Prediction season-day values.
#' @param season_start_date Season start date.
#' @param season_length Season length.
#' @param method Double-growth method.
#' @param growth_fraction Growth fractions used for cumulative timing.
#' @param min_unique_growth Minimum number of unique growth values.
#' @param rate_threshold_fraction Fraction of peak rate.
#' @param peak_separation_min Minimum separation between peaks.
#' @param valley_ratio_max Maximum allowed valley-to-peak ratio.
#' @param min_relative_peak Minimum relative derivative peak height.
#' @param fallback_to_single Logical fallback flag.
#' @param fix_a_to_observed_max Logical fixed-asymptote flag.
#' @param fixed_a_multiplier Multiplier for fixed asymptote.
#' @param start_value_double_gompertz_parameters Double-Gompertz starts.
#' @param start_value_double_richards_parameters Double-Richards starts.
#'
#' @return List with fitted predictions and parameters.
#'
#' @keywords internal
dmgfd_fit_one_curve <- function(x_obs,
y_obs,
x_pred,
season_start_date = as.Date(NA),
season_length = NA_real_,
method,
growth_fraction,
min_unique_growth,
rate_threshold_fraction = 0.1,
peak_separation_min = 10,
valley_ratio_max = 0.4,
min_relative_peak = 0.05,
fallback_to_single = TRUE,
fix_a_to_observed_max = FALSE,
fixed_a_multiplier = 1,
start_value_double_gompertz_parameters,
start_value_double_richards_parameters) {
pred_full <- rep(NA_real_, length(x_pred))
params <- dmgfd_empty_model_params()
params$n_obs <- sum(is.finite(x_obs) & is.finite(y_obs))
params$n_days_observed <- NA_integer_
params$first_obs_day <- NA_real_
params$last_obs_day <- NA_real_
params$season_length <- season_length
params$n_missing_days <- NA_real_
params$extrapolated <- NA
params$anchor_added <- FALSE
params$method_used <- paste0("double_", method)
good <- is.finite(x_obs) & is.finite(y_obs)
if (sum(good) == 0) {
params$model_note <- "No observed data in this vegetation year."
return(list(pred = pred_full, params = params))
}
x_raw <- x_obs[good]
y_raw <- y_obs[good]
params$n_days_observed <- length(x_raw)
params$first_obs_day <- min(x_raw, na.rm = TRUE)
params$last_obs_day <- max(x_raw, na.rm = TRUE)
if (is.finite(season_length)) {
params$n_missing_days <- season_length - length(unique(x_raw))
params$extrapolated <- params$first_obs_day > 1 ||
params$last_obs_day < season_length
}
if (length(unique(y_raw)) < min_unique_growth) {
params$model_note <- paste0(
"Too few unique cumulative-growth values (",
length(unique(y_raw)), " < ", min_unique_growth, ")."
)
return(list(pred = pred_full, params = params))
}
anchored <- dmgfd_add_anchor_points(x_raw, y_raw)
x_fit <- anchored$x
y_fit <- anchored$y
params$anchor_added <- length(x_fit) > length(x_raw)
dyn_range <- max(y_fit, na.rm = TRUE) - min(y_fit, na.rm = TRUE)
if (!is.finite(dyn_range) || dyn_range <= 0) {
params$model_note <- "Observed cumulative-growth range is too small."
return(list(pred = pred_full, params = params))
}
fixed_a_value <- NA_real_
if (isTRUE(fix_a_to_observed_max) &&
method %in% c("gompertz", "richards")) {
fixed_a_value <- suppressWarnings(max(y_raw, na.rm = TRUE)) * fixed_a_multiplier
if (!is.finite(fixed_a_value) || fixed_a_value <= 0) {
fixed_a_value <- NA_real_
}
}
fit_res <- switch(
method,
gompertz = dmgfd_fit_model_double_gompertz(
x = x_fit,
y = y_fit,
x_pred = x_pred,
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_double_gompertz_parameters = start_value_double_gompertz_parameters
),
richards = dmgfd_fit_model_double_richards(
x = x_fit,
y = y_fit,
x_pred = x_pred,
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_double_richards_parameters = start_value_double_richards_parameters
)
)
pred_full <- fit_res$pred
params[names(fit_res$params)] <- fit_res$params
params$method_used <- paste0("double_", method)
params <- dmgfd_add_overall_timing(
params = params,
pred_full = pred_full,
x_pred = x_pred,
growth_fraction = growth_fraction,
season_start_date = season_start_date
)
params <- dmgfd_add_overall_rate_timing(
params = params,
pred_full = pred_full,
x_pred = x_pred,
rate_threshold_fraction = rate_threshold_fraction,
season_start_date = season_start_date
)
rate <- dmgfd_growth_rate(pred_full)
pulse_info <- dmgfd_detect_two_pulse_pattern(
rate = rate,
x_pred = x_pred,
peak_separation_min = peak_separation_min,
valley_ratio_max = valley_ratio_max,
min_relative_peak = min_relative_peak
)
if (isTRUE(pulse_info$detected)) {
params$two_pulse_detected <- TRUE
params$peak1_season_day <- pulse_info$peak1_day
params$peak2_season_day <- pulse_info$peak2_day
params$separator_season_day <- pulse_info$separator_day
params$valley_rate <- pulse_info$valley_rate
params$pulse1_peak_rate <- pulse_info$pulse1_peak_rate
params$pulse2_peak_rate <- pulse_info$pulse2_peak_rate
if (!is.na(season_start_date)) {
peak1_date <- season_start_date + (params$peak1_season_day - 1)
peak2_date <- season_start_date + (params$peak2_season_day - 1)
params$separator_date <- season_start_date + (params$separator_season_day - 1)
params$peak1_day <- lubridate::yday(peak1_date)
params$peak2_day <- lubridate::yday(peak2_date)
params$separator_day <- lubridate::yday(params$separator_date)
}
thr1 <- rate_threshold_fraction * params$pulse1_peak_rate
thr2 <- rate_threshold_fraction * params$pulse2_peak_rate
idx1 <- which(is.finite(x_pred) & x_pred <= params$separator_season_day)
idx2 <- which(is.finite(x_pred) & x_pred > params$separator_season_day)
win1 <- dmgfd_pulse_window_from_rate(
rate = rate,
x_pred = x_pred,
threshold = thr1,
domain_idx = idx1
)
win2 <- dmgfd_pulse_window_from_rate(
rate = rate,
x_pred = x_pred,
threshold = thr2,
domain_idx = idx2
)
if (is.finite(win1$start_day)) {
params$pulse1_start_season_day <- win1$start_day
if (!is.na(season_start_date)) {
params$pulse1_start_date <- season_start_date + (win1$start_day - 1)
params$pulse1_start_day <- lubridate::yday(params$pulse1_start_date)
}
}
if (is.finite(win1$end_day)) {
params$pulse1_end_season_day <- win1$end_day
if (!is.na(season_start_date)) {
params$pulse1_end_date <- season_start_date + (win1$end_day - 1)
params$pulse1_end_day <- lubridate::yday(params$pulse1_end_date)
}
}
if (is.finite(win2$start_day)) {
params$pulse2_start_season_day <- win2$start_day
if (!is.na(season_start_date)) {
params$pulse2_start_date <- season_start_date + (win2$start_day - 1)
params$pulse2_start_day <- lubridate::yday(params$pulse2_start_date)
}
}
if (is.finite(win2$end_day)) {
params$pulse2_end_season_day <- win2$end_day
if (!is.na(season_start_date)) {
params$pulse2_end_date <- season_start_date + (win2$end_day - 1)
params$pulse2_end_day <- lubridate::yday(params$pulse2_end_date)
}
}
return(list(pred = pred_full, params = params))
}
params$two_pulse_detected <- FALSE
if (!isTRUE(fallback_to_single)) {
if (is.na(params$model_note) || !nzchar(params$model_note)) {
params$model_note <- "No convincing two-pulse pattern detected."
} else {
params$model_note <- paste(
params$model_note,
"No convincing two-pulse pattern detected."
)
}
return(list(pred = pred_full, params = params))
}
meta_fields <- c(
"n_obs", "n_days_observed", "first_obs_day", "last_obs_day",
"season_length", "n_missing_days", "extrapolated", "anchor_added"
)
double_note <- params$model_note
fallback_res <- dmgfd_fit_single_fallback(
method = method,
x = x_fit,
y = y_fit,
x_pred = x_pred,
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_double_gompertz_parameters = start_value_double_gompertz_parameters,
start_value_double_richards_parameters = start_value_double_richards_parameters
)
pred_full <- fallback_res$pred
fb_params <- dmgfd_empty_model_params()
fb_params[meta_fields] <- params[meta_fields]
fb_params[names(fallback_res$params)] <- fallback_res$params
fb_params$two_pulse_detected <- FALSE
fb_params$fallback_used <- TRUE
if (is.na(double_note) || !nzchar(double_note)) {
fb_params$model_note <- paste0(
"No convincing two-pulse pattern detected; ",
fb_params$method_used,
" fallback used."
)
} else {
fb_params$model_note <- paste(
double_note,
paste0(
"No convincing two-pulse pattern detected; ",
fb_params$method_used,
" fallback used."
)
)
}
fb_params <- dmgfd_add_overall_timing(
params = fb_params,
pred_full = pred_full,
x_pred = x_pred,
growth_fraction = growth_fraction,
season_start_date = season_start_date
)
fb_params <- dmgfd_add_overall_rate_timing(
params = fb_params,
pred_full = pred_full,
x_pred = x_pred,
rate_threshold_fraction = rate_threshold_fraction,
season_start_date = season_start_date
)
list(
pred = pred_full,
params = fb_params
)
}
Try the dendRoAnalyst package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.