Nothing
R/dm.growth.fit.R
In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data
Defines functions
dm.growth.fit
Documented in
dm.growth.fit
#' Fit dendrometer growth curves by vegetation season
#'
#' @description
#' Fits cumulative growth curves to daily dendrometer series using one of
#' several supported methods: generalized additive model \code{"gam"},
#' Gompertz \code{"gompertz"}, logistic \code{"logistic"}, Richards
#' \code{"richards"}, local regression \code{"loess"}, or smoothing spline
#' \code{"spline"}.
#'
#' The function:
#' \itemize{
#' \item resamples dendrometer data to daily maxima using [dendro.resample()],
#' \item assigns each daily observation to a vegetation season,
#' \item optionally converts daily stem-size values to cumulative growth
#' with \code{fit_GRO = TRUE},
#' \item fits one curve per series and season with
#' \code{year_mode = "yearly"} or one pooled curve across retained seasons
#' with \code{year_mode = "pooled"},
#' \item estimates growing-season timing from the fitted cumulative-growth
#' curve,
#' \item additionally identifies active-growth timing from the fitted
#' growth-rate curve.
#' }
#'
#' For Gompertz, logistic, and Richards fits, the asymptote parameter
#' \code{a}, representing maximum seasonal growth, can optionally be fixed to
#' the observed maximum cumulative growth for that series and vegetation season.
#' This is controlled by \code{fix_a_to_observed_max}.
#'
#' @param df A data frame or tibble. The first column must be a timestamp
#' \code{Date}, \code{POSIXct}, or parseable date-time string. All selected
#' remaining 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
#' with series names.
#' @param method Fitting method. One of \code{"gam"}, \code{"gompertz"},
#' \code{"logistic"}, \code{"richards"}, \code{"loess"}, or \code{"spline"}.
#' @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.
#' Cross-year seasons are supported, for example \code{c(275, 274)}.
#' @param growth_fraction Numeric vector of length two giving the lower and
#' upper fractions of final fitted seasonal growth used to define cumulative
#' growth onset and cessation, for example \code{c(0.1, 0.9)}.
#' @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. Active-growth
#' dates are derived from the fitted growth-rate curve as the first and last
#' days where fitted growth rate exceeds this fraction of the peak fitted
#' growth rate.
#' @param fix_a_to_observed_max Logical. If \code{TRUE} and
#' \code{method = "gompertz"}, \code{"logistic"}, or \code{"richards"}, the
#' asymptote parameter \code{a}, representing maximum seasonal growth, is fixed
#' to the observed maximum cumulative growth of the corresponding series and
#' vegetation season. 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
#' \code{a} 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 just above the observed maximum.
#' @param start_value_gompertz_parameters Optional named list with starting
#' values for Gompertz fits. Supported names are \code{a}, \code{b}, and
#' \code{k}. If \code{fix_a_to_observed_max = TRUE}, supplied \code{a} is
#' ignored.
#' @param start_value_richards_parameters Optional named list with starting
#' values for Richards and logistic fits. Supported names are \code{a},
#' \code{k}, \code{t0}, and \code{v}. If
#' \code{fix_a_to_observed_max = TRUE}, supplied \code{a} is ignored.
#' @param loess_span Span used when \code{method = "loess"}.
#' @param spline_df Degrees of freedom used when \code{method = "spline"}.
#' @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}{Table with fit-level statistics and estimated timing.}
#' \item{fit_parameters}{Table with fit-level model parameters and convergence
#' information.}
#' \item{season_table}{Vegetation seasons retained for fitting.}
#' }
#'
#' @details
#' The first column of \code{df} is treated as the time column and renamed to
#' \code{"TIME"} internally. If it is not already a date-time object, the
#' function attempts to parse it using [lubridate::parse_date_time()].
#'
#' \code{growth_start_*} and \code{growth_end_*} are based on cumulative fitted
#' growth. \code{rate_start_*} and \code{rate_end_*} are based on the fitted
#' growth-rate curve.
#'
#' For Gompertz, logistic, and Richards models, \code{a} is the upper asymptote.
#' When \code{fix_a_to_observed_max = TRUE}, this parameter is not estimated by
#' nonlinear least squares. Instead, it is fixed to:
#'
#' \deqn{a = max(y_{obs}) \times fixed\_a\_multiplier}
#'
#' where \eqn{y_{obs}} is the observed cumulative growth for that series and
#' vegetation season.
#'
#' @examples
#' \donttest{
#' # fit <- dm.growth.fit(
#' # df = dendro_data,
#' # TreeNum = "all",
#' # method = "gompertz",
#' # year_mode = "yearly",
#' # fit_GRO = TRUE,
#' # fix_a_to_observed_max = TRUE
#' # )
#'
#' # fit2 <- dm.growth.fit(
#' # 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 [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 summarise
#' @importFrom tibble as_tibble tibble
#' @importFrom lubridate parse_date_time year yday month
#' @importFrom stats lm coef qlogis predict smooth.spline loess na.omit
#' @importFrom utils modifyList head
#' @export
dm.growth.fit <- function(df,
TreeNum = "all",
method = c("gam", "gompertz", "logistic", "richards", "loess", "spline"),
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,
fix_a_to_observed_max = FALSE,
fixed_a_multiplier = 1,
start_value_gompertz_parameters = list(a = NA_real_, b = NA_real_, k = NA_real_),
start_value_richards_parameters = list(a = NA_real_, k = NA_real_, t0 = NA_real_, v = 1),
loess_span = 0.2,
spline_df = 10,
verbose = TRUE) {
cl <- match.call()
TIME <- doy <- season_label <- season_start <- season_end <- season_day <- season_length <- any_observed <- NULL
method <- match.arg(method)
year_mode <- match.arg(year_mode)
site_mode <- match.arg(site_mode)
dmgf_validate_fit_growth_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,
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 two 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 two 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 <- dmgf_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),
~ . - dmgf_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),
dmgf_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]
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 <- dmgf_fit_one_curve(
x_obs = x_obs,
y_obs = y_obs,
x_pred = x_full,
season_start_date = ss,
season_length = full_len,
season_id = NULL,
method = method,
loess_span = loess_span,
spline_df = spline_df,
growth_fraction = growth_fraction,
min_unique_growth = min_unique_growth,
rate_threshold_fraction = rate_threshold_fraction,
fix_a_to_observed_max = fix_a_to_observed_max,
fixed_a_multiplier = fixed_a_multiplier,
start_value_gompertz_parameters = start_value_gompertz_parameters,
start_value_richards_parameters = start_value_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 = method,
site_mode = site_mode,
season_start = ss,
season_end = unique(dat_proc$season_end[idx])[1],
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,
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,
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_full <- dat_proc$season_day
y_full <- dat_proc[[series_name]]
sid_full <- dat_proc$season_label
good <- is.finite(y_full)
x_obs <- x_full[good]
y_obs <- y_full[good]
sid_obs <- sid_full[good]
fit_res <- dmgf_fit_one_curve(
x_obs = x_obs,
y_obs = y_obs,
x_pred = x_full,
season_start_date = as.Date(NA),
season_length = NA_real_,
season_id = sid_obs,
method = method,
loess_span = loess_span,
spline_df = spline_df,
growth_fraction = growth_fraction,
min_unique_growth = min_unique_growth,
rate_threshold_fraction = rate_threshold_fraction,
fix_a_to_observed_max = fix_a_to_observed_max,
fixed_a_multiplier = fixed_a_multiplier,
start_value_gompertz_parameters = start_value_gompertz_parameters,
start_value_richards_parameters = start_value_richards_parameters
)
fitted_dat[[series_name]] <- fit_res$pred
param_rows[[row_id]] <- tibble::tibble(
series = series_name,
fit_id = "pooled",
year_mode = year_mode,
method = method,
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,
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,
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"
)
param_cols <- c(
"series", "fit_id", "year_mode", "method", "site_mode",
"converged", "fixed_a_used", "fixed_a_value",
"a", "b", "k", "t0", "v",
"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 completed: ",
nrow(fit_parameters), " fit(s) across ",
length(unique(fit_parameters$series)), " series."
)
}
out
}
# helpers ----------------------------------------------------------------------
#' Validate inputs for 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 vector.
#' @param rate_threshold_fraction Fraction of peak rate used for rate window.
#' @param fixed_a_multiplier Multiplier used when fixing asymptote \code{a}.
#'
#' @return Invisibly returns \code{TRUE} if checks pass.
#'
#' @keywords internal
dmgf_validate_fit_growth_inputs <- function(df,
growth_fraction,
min_unique_growth,
custom_veg_season,
rate_threshold_fraction,
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(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
dmgf_is_scalar_finite <- function(x) {
is.numeric(x) && length(x) == 1 && !is.na(x) && is.finite(x)
}
#' Return the first non-missing value
#'
#' @param x Numeric vector.
#'
#' @return First non-missing value, or \code{NA_real_}.
#'
#' @keywords internal
dmgf_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
dmgf_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
dmgf_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
dmgf_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] <- dmgf_doy_to_date(yr[inside], start_doy)
season_end[inside] <- dmgf_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] <- dmgf_doy_to_date(sy, start_doy)
season_end[inside_late] <- dmgf_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] <- dmgf_doy_to_date(sy, start_doy)
season_end[inside_early] <- dmgf_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 points
#'
#' @param x Season-day values.
#' @param y Growth values.
#' @param season_id Optional season identifier for pooled fits.
#'
#' @return List with \code{x} and \code{y}.
#'
#' @keywords internal
dmgf_add_anchor_points <- function(x, y, season_id = NULL) {
if (length(x) == 0) {
return(list(x = x, y = y))
}
if (is.null(season_id)) {
if (min(x, na.rm = TRUE) > 1 && !any(x == 1)) {
x <- c(1, x)
y <- c(0, y)
}
ord <- order(x)
return(list(
x = x[ord],
y = y[ord]
))
}
season_id <- as.character(season_id)
xs <- list()
ys <- list()
uids <- unique(season_id)
for (i in seq_along(uids)) {
sid <- uids[i]
idx <- season_id == sid
x_i <- x[idx]
y_i <- y[idx]
if (length(x_i) == 0) {
next
}
if (min(x_i, na.rm = TRUE) > 1 && !any(x_i == 1)) {
x_i <- c(1, x_i)
y_i <- c(0, y_i)
}
ord <- order(x_i)
xs[[i]] <- x_i[ord]
ys[[i]] <- y_i[ord]
}
list(
x = unlist(xs, use.names = FALSE),
y = unlist(ys, use.names = FALSE)
)
}
#' Infer logistic 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
dmgf_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 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
dmgf_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
)
}
#' Empty model-parameter template
#'
#' @return Named list of empty model parameters.
#'
#' @keywords internal
dmgf_empty_model_params <- function() {
list(
converged = FALSE,
fixed_a_used = FALSE,
fixed_a_value = NA_real_,
a = NA_real_,
b = NA_real_,
k = NA_real_,
t0 = NA_real_,
v = NA_real_,
edf = NA_real_,
span = NA_real_,
spline_df = NA_real_,
spar = NA_real_,
model_note = NA_character_
)
}
#' Dispatch growth-curve model fitting
#'
#' @param method Fitting method.
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param loess_span Span for loess.
#' @param spline_df Degrees of freedom for spline.
#' @param extrapolated Logical flag for extrapolation.
#' @param fixed_a_value Optional fixed asymptote.
#' @param fix_a_to_observed_max Logical.
#' @param start_value_gompertz_parameters Starting values for Gompertz.
#' @param start_value_richards_parameters Starting values for logistic/Richards.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model <- function(method,
x,
y,
x_pred,
loess_span = 0.2,
spline_df = 10,
extrapolated = FALSE,
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_),
start_value_richards_parameters = list(a = NA_real_, k = NA_real_, t0 = NA_real_, v = 1)) {
switch(
method,
gam = dmgf_fit_model_gam(
x = x,
y = y,
x_pred = x_pred,
extrapolated = extrapolated
),
gompertz = dmgf_fit_model_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_value_gompertz_parameters
),
logistic = dmgf_fit_model_logistic(
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_value_richards_parameters
),
richards = dmgf_fit_model_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_value_richards_parameters
),
loess = dmgf_fit_model_loess(
x = x,
y = y,
x_pred = x_pred,
loess_span = loess_span,
extrapolated = extrapolated
),
spline = dmgf_fit_model_spline(
x = x,
y = y,
x_pred = x_pred,
spline_df = spline_df,
extrapolated = extrapolated
),
stop("Unknown method: ", method)
)
}
#' Fit GAM growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param extrapolated Logical flag for extrapolation.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_gam <- function(x, y, x_pred, extrapolated = FALSE) {
params <- dmgf_empty_model_params()
res <- tryCatch({
if (!requireNamespace("mgcv", quietly = TRUE)) {
stop("Package 'mgcv' is required for method = 'gam'. Please install it.")
}
x <- as.numeric(x)
y <- as.numeric(y)
x_pred <- as.numeric(x_pred)
dat <- data.frame(x = x, y = y)
ux <- sort(unique(dat$x))
if (length(ux) < 4) {
mod <- stats::lm(y ~ x, data = dat)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred))
pred <- as.numeric(pred)
pred <- pmax(pred, 0)
pred <- cummax(pred)
params$converged <- TRUE
params$model_note <- "lm fallback used because too few unique x values for GAM."
return(list(pred = pred, params = params))
}
k_gam <- min(10L, length(ux) - 1L)
k_gam <- max(4L, k_gam)
mod <- mgcv::gam(
y ~ s(x, bs = "cs", k = k_gam),
data = dat,
method = "REML",
select = TRUE
)
pred <- stats::predict(
mod,
newdata = data.frame(x = x_pred),
type = "response"
)
pred <- as.numeric(pred)
pred <- pmax(pred, 0)
pred <- cummax(pred)
st <- summary(mod)$s.table
if (!is.null(st)) {
st <- as.matrix(st)
if ("edf" %in% colnames(st)) {
params$edf <- unname(st[1, "edf"])
}
}
params$converged <- TRUE
if (isTRUE(extrapolated)) {
params$model_note <- "GAM predictions outside observed DOY range may be unstable."
}
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit Gompertz growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param fixed_a_value Optional fixed asymptote value.
#' @param fix_a_to_observed_max Logical.
#' @param start_value_gompertz_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_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 <- dmgf_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) &&
dmgf_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 <- dmgf_infer_gompertz_starts(x, y_adj, a_fixed_fit)
start_guess <- utils::modifyList(start_guess, start_value_gompertz_parameters)
if (!dmgf_is_scalar_finite(start_guess$b)) {
start_guess$b <- 0.5
}
if (!dmgf_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$model_note <- "Gompertz asymptote a fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
start_guess <- dmgf_infer_gompertz_starts(x, y_adj, a0)
start_guess <- utils::modifyList(start_guess, start_value_gompertz_parameters)
if (!dmgf_is_scalar_finite(start_guess$a)) {
start_guess$a <- a0
}
if (!dmgf_is_scalar_finite(start_guess$b)) {
start_guess$b <- 0.5
}
if (!dmgf_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
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit logistic growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param fixed_a_value Optional fixed asymptote value.
#' @param fix_a_to_observed_max Logical.
#' @param start_value_richards_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_logistic <- 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 <- dmgf_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) &&
dmgf_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 <- dmgf_infer_logistic_starts(x, y_adj, a_fixed_fit)
if (dmgf_is_scalar_finite(start_value_richards_parameters$k)) {
start_guess$k <- start_value_richards_parameters$k
}
if (dmgf_is_scalar_finite(start_value_richards_parameters$t0)) {
start_guess$t0 <- start_value_richards_parameters$t0
}
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a_fixed_fit / (1 + exp(-(k * (x - t0)))),
data = dat,
start = list(
k = start_guess$k,
t0 = start_guess$t0
),
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$converged <- TRUE
params$fixed_a_used <- TRUE
params$fixed_a_value <- a_fixed_original
params$model_note <- "Logistic asymptote a fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
start_guess <- dmgf_infer_logistic_starts(x, y_adj, a0)
if (dmgf_is_scalar_finite(start_value_richards_parameters$a)) {
start_guess$a <- start_value_richards_parameters$a
}
if (dmgf_is_scalar_finite(start_value_richards_parameters$k)) {
start_guess$k <- start_value_richards_parameters$k
}
if (dmgf_is_scalar_finite(start_value_richards_parameters$t0)) {
start_guess$t0 <- start_value_richards_parameters$t0
}
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a / (1 + exp(-(k * (x - t0)))),
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$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$converged <- TRUE
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit Richards growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param fixed_a_value Optional fixed asymptote value.
#' @param fix_a_to_observed_max Logical.
#' @param start_value_richards_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_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 <- dmgf_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) &&
dmgf_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 <- dmgf_infer_logistic_starts(x, y_adj, a_fixed_fit)
start_vals <- list(
k = if (dmgf_is_scalar_finite(start_value_richards_parameters$k)) {
start_value_richards_parameters$k
} else {
logi_start$k
},
t0 = if (dmgf_is_scalar_finite(start_value_richards_parameters$t0)) {
start_value_richards_parameters$t0
} else {
logi_start$t0
},
v = if (dmgf_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$model_note <- "Richards asymptote a fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
logi_start <- dmgf_infer_logistic_starts(x, y_adj, a0)
start_vals <- list(
a = if (dmgf_is_scalar_finite(start_value_richards_parameters$a)) {
start_value_richards_parameters$a
} else {
logi_start$a
},
k = if (dmgf_is_scalar_finite(start_value_richards_parameters$k)) {
start_value_richards_parameters$k
} else {
logi_start$k
},
t0 = if (dmgf_is_scalar_finite(start_value_richards_parameters$t0)) {
start_value_richards_parameters$t0
} else {
logi_start$t0
},
v = if (dmgf_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
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit loess growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param loess_span Loess span.
#' @param extrapolated Logical flag for extrapolation.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_loess <- function(x,
y,
x_pred,
loess_span = 0.2,
extrapolated = FALSE) {
params <- dmgf_empty_model_params()
res <- tryCatch({
dat <- data.frame(x = x, y = y)
mod <- stats::loess(
y ~ x,
data = dat,
span = loess_span,
degree = 2
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred))
pred <- as.numeric(pred)
params$span <- loess_span
params$edf <- if (!is.null(mod$enp)) mod$enp else NA_real_
params$converged <- TRUE
if (isTRUE(extrapolated)) {
params$model_note <- "loess extrapolation outside observed DOY range may be unreliable."
}
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit smoothing-spline growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param spline_df Spline degrees of freedom.
#' @param extrapolated Logical flag for extrapolation.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_spline <- function(x,
y,
x_pred,
spline_df = 10,
extrapolated = FALSE) {
params <- dmgf_empty_model_params()
res <- tryCatch({
mod <- stats::smooth.spline(x, y, df = spline_df)
pred <- stats::predict(mod, x_pred)$y
pred <- as.numeric(pred)
params$spline_df <- mod$df
params$spar <- mod$spar
params$converged <- TRUE
if (isTRUE(extrapolated)) {
params$model_note <- "Spline predictions outside observed DOY range may be uncertain."
}
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Reduce prediction curve to unique x values
#'
#' @param x_pred Prediction x values.
#' @param pred Predicted y values.
#'
#' @return List with unique \code{x} and mean \code{pred}.
#'
#' @keywords internal
dmgf_unique_prediction_curve <- function(x_pred, pred) {
ok <- is.finite(x_pred) & is.finite(pred)
if (sum(ok) == 0) {
return(list(x = numeric(0), pred = numeric(0)))
}
xp <- as.numeric(x_pred[ok])
pp <- as.numeric(pred[ok])
ux <- sort(unique(xp))
mp <- vapply(
ux,
function(xx) mean(pp[xp == xx]),
numeric(1)
)
list(
x = ux,
pred = mp
)
}
#' Calculate non-negative fitted growth rate
#'
#' @param pred Fitted cumulative growth.
#'
#' @return Numeric growth-rate vector.
#'
#' @keywords internal
dmgf_growth_rate <- function(pred) {
pred <- as.numeric(pred)
if (length(pred) < 2) {
return(rep(NA_real_, length(pred)))
}
pred_mono <- cummax(pmax(pred, 0))
rate <- c(NA_real_, diff(pred_mono))
rate[!is.finite(rate)] <- NA_real_
rate <- pmax(rate, 0)
rate
}
#' Estimate active-growth window from growth rate
#'
#' @param rate Fitted growth-rate vector.
#' @param x_pred Prediction x values.
#' @param threshold_fraction Fraction of peak growth rate.
#'
#' @return List with peak rate, start day, and end day.
#'
#' @keywords internal
dmgf_rate_window <- function(rate,
x_pred,
threshold_fraction = 0.1) {
out <- list(
peak_rate = NA_real_,
start_day = NA_real_,
end_day = NA_real_
)
ok <- is.finite(rate) & is.finite(x_pred)
if (sum(ok) == 0) {
return(out)
}
r <- rate[ok]
x <- x_pred[ok]
peak_rate <- suppressWarnings(max(r, na.rm = TRUE))
if (!is.finite(peak_rate) || peak_rate <= 0) {
return(out)
}
thr <- threshold_fraction * peak_rate
above <- which(r >= thr)
if (length(above) == 0) {
return(out)
}
out$peak_rate <- peak_rate
out$start_day <- min(x[above], na.rm = TRUE)
out$end_day <- max(x[above], na.rm = TRUE)
out
}
#' Fit one growth curve
#'
#' @description
#' Internal helper used by [dm.growth.fit()] to fit one curve for one series and
#' one season, or one pooled curve.
#'
#' @param x_obs Observed x values.
#' @param y_obs Observed growth values.
#' @param x_pred Prediction x values.
#' @param season_start_date Season start date.
#' @param season_length Season length.
#' @param season_id Optional season ID for pooled fits.
#' @param method Fitting method.
#' @param loess_span Loess span.
#' @param spline_df Spline degrees of freedom.
#' @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 used for rate timing.
#' @param fix_a_to_observed_max Logical.
#' @param fixed_a_multiplier Multiplier for fixed asymptote.
#' @param start_value_gompertz_parameters Starting values for Gompertz.
#' @param start_value_richards_parameters Starting values for logistic/Richards.
#'
#' @return List with fitted predictions and parameters.
#'
#' @keywords internal
dmgf_fit_one_curve <- function(x_obs,
y_obs,
x_pred,
season_start_date = as.Date(NA),
season_length = NA_real_,
season_id = NULL,
method,
loess_span,
spline_df,
growth_fraction,
min_unique_growth,
rate_threshold_fraction = 0.1,
fix_a_to_observed_max = FALSE,
fixed_a_multiplier = 1,
start_value_gompertz_parameters,
start_value_richards_parameters) {
pred_full <- rep(NA_real_, length(x_pred))
params <- list(
converged = FALSE,
fixed_a_used = FALSE,
fixed_a_value = NA_real_,
n_obs = sum(is.finite(x_obs) & is.finite(y_obs)),
n_days_observed = NA_integer_,
first_obs_day = NA_real_,
last_obs_day = NA_real_,
season_length = season_length,
n_missing_days = NA_real_,
extrapolated = NA,
anchor_added = FALSE,
a = NA_real_,
b = NA_real_,
k = NA_real_,
t0 = NA_real_,
v = NA_real_,
edf = NA_real_,
span = NA_real_,
spline_df = NA_real_,
spar = NA_real_,
model_note = NA_character_,
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)
)
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]
sid_raw <- if (is.null(season_id)) NULL else season_id[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 <- dmgf_add_anchor_points(
x = x_raw,
y = y_raw,
season_id = sid_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", "logistic", "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 <- dmgf_fit_model(
method = method,
x = x_fit,
y = y_fit,
x_pred = x_pred,
loess_span = loess_span,
spline_df = spline_df,
extrapolated = isTRUE(params$extrapolated),
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_gompertz_parameters = start_value_gompertz_parameters,
start_value_richards_parameters = start_value_richards_parameters
)
pred_full <- as.numeric(fit_res$pred)
params[names(fit_res$params)] <- fit_res$params
timing_curve <- dmgf_unique_prediction_curve(
x_pred = x_pred,
pred = pred_full
)
timing_x <- timing_curve$x
timing_pred <- timing_curve$pred
if (length(timing_pred) > 0) {
timing_pred_mono <- cummax(pmax(timing_pred, 0))
final_val <- tail(stats::na.omit(timing_pred_mono), 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(timing_x[timing_pred_mono > low_val], na.rm = TRUE)
)
en_pos <- suppressWarnings(
max(timing_x[timing_pred_mono < 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)
}
}
}
rate <- dmgf_growth_rate(timing_pred)
rate_win <- dmgf_rate_window(
rate = rate,
x_pred = timing_x,
threshold_fraction = rate_threshold_fraction
)
params$peak_rate <- rate_win$peak_rate
if (is.finite(rate_win$start_day)) {
params$rate_start_season_day <- rate_win$start_day
if (!is.na(season_start_date)) {
params$rate_start_date <- season_start_date + (rate_win$start_day - 1)
params$rate_start_day <- lubridate::yday(params$rate_start_date)
}
}
if (is.finite(rate_win$end_day)) {
params$rate_end_season_day <- rate_win$end_day
if (!is.na(season_start_date)) {
params$rate_end_date <- season_start_date + (rate_win$end_day - 1)
params$rate_end_day <- lubridate::yday(params$rate_end_date)
}
}
}
list(
pred = pred_full,
params = params
)
}
# S3 methods -------------------------------------------------------------------
#' Print a dm_growth_fit object
#'
#' @description
#' Prints a compact overview of an object returned by [dm.growth.fit()] or
#' [dm.growth.fit.double()].
#'
#' @param x An object of class \code{"dm_growth_fit"}.
#' @param ... Further arguments passed to or from other methods.
#'
#' @return The input object, invisibly.
#'
#' @method print dm_growth_fit
#' @export
print.dm_growth_fit <- function(x, ...) {
fit_parameters <- x$fit_parameters
fit_statistics <- x$fit_statistics
if (is.null(fit_parameters) || nrow(fit_parameters) == 0) {
cat("<dm_growth_fit>\n")
cat("No fit results available.\n")
return(invisible(x))
}
methods_used <- if ("method" %in% names(fit_parameters)) {
unique(stats::na.omit(fit_parameters$method))
} else {
character(0)
}
n_series <- if ("series" %in% names(fit_parameters)) {
length(unique(stats::na.omit(fit_parameters$series)))
} else {
NA_integer_
}
n_fits <- nrow(fit_parameters)
converged_count <- if ("converged" %in% names(fit_parameters)) {
sum(fit_parameters$converged %in% TRUE, na.rm = TRUE)
} else {
NA_integer_
}
fixed_a_count <- if ("fixed_a_used" %in% names(fit_parameters)) {
sum(fit_parameters$fixed_a_used %in% TRUE, na.rm = TRUE)
} else {
0L
}
fallback_count <- if ("fallback_used" %in% names(fit_parameters)) {
sum(fit_parameters$fallback_used %in% TRUE, na.rm = TRUE)
} else {
0L
}
two_pulse_count <- if (!is.null(fit_statistics) &&
"two_pulse_detected" %in% names(fit_statistics)) {
sum(fit_statistics$two_pulse_detected %in% TRUE, na.rm = TRUE)
} else {
0L
}
cat("Call:\n")
print(x$call)
cat("\n")
cat("<dm_growth_fit>\n")
cat("Series:", n_series, "\n")
cat("Fits:", n_fits, "\n")
cat("Methods used:", paste(methods_used, collapse = ", "), "\n")
cat("Converged fits:", converged_count, "\n")
cat("Fixed-a fits:", fixed_a_count, "\n")
if ("fallback_used" %in% names(fit_parameters)) {
cat("Fallback to single curve:", fallback_count, "\n")
}
if (!is.null(fit_statistics) &&
"two_pulse_detected" %in% names(fit_statistics)) {
cat("True two-pulse fits:", two_pulse_count, "\n")
}
invisible(x)
}
#' Summarize a dm_growth_fit object
#'
#' @description
#' Summarizes an object returned by [dm.growth.fit()] or
#' [dm.growth.fit.double()].
#'
#' @param object An object of class \code{"dm_growth_fit"}.
#' @param ... Further arguments passed to or from other methods.
#'
#' @return An object of class \code{"summary.dm_growth_fit"}.
#'
#' @method summary dm_growth_fit
#' @export
summary.dm_growth_fit <- function(object, ...) {
fit_parameters <- object$fit_parameters
fit_statistics <- object$fit_statistics
season_table <- object$season_table
if (is.null(fit_parameters)) {
fit_parameters <- tibble::tibble()
}
if (is.null(fit_statistics)) {
fit_statistics <- tibble::tibble()
}
if (is.null(season_table)) {
season_table <- tibble::tibble()
}
methods_used <- if ("method" %in% names(fit_parameters)) {
unique(stats::na.omit(fit_parameters$method))
} else {
character(0)
}
fixed_a_count <- if ("fixed_a_used" %in% names(fit_parameters)) {
sum(fit_parameters$fixed_a_used %in% TRUE, na.rm = TRUE)
} else {
0L
}
fallback_count <- if ("fallback_used" %in% names(fit_parameters)) {
sum(fit_parameters$fallback_used %in% TRUE, na.rm = TRUE)
} else {
0L
}
nonfallback_count <- if ("fallback_used" %in% names(fit_parameters)) {
sum(fit_parameters$fallback_used %in% FALSE, na.rm = TRUE)
} else {
NA_integer_
}
two_pulse_count <- if ("two_pulse_detected" %in% names(fit_statistics)) {
sum(fit_statistics$two_pulse_detected %in% TRUE, na.rm = TRUE)
} else {
0L
}
not_two_pulse_count <- if ("two_pulse_detected" %in% names(fit_statistics)) {
sum(fit_statistics$two_pulse_detected %in% FALSE, na.rm = TRUE)
} else {
NA_integer_
}
out <- list(
call = object$call,
n_series = if ("series" %in% names(fit_parameters)) {
length(unique(stats::na.omit(fit_parameters$series)))
} else {
0L
},
n_fits = nrow(fit_parameters),
methods_used = methods_used,
site_mode = if ("site_mode" %in% names(fit_parameters)) {
unique(stats::na.omit(fit_parameters$site_mode))
} else {
character(0)
},
converged = if ("converged" %in% names(fit_parameters)) {
sum(fit_parameters$converged %in% TRUE, na.rm = TRUE)
} else {
0L
},
failed = if ("converged" %in% names(fit_parameters)) {
sum(fit_parameters$converged %in% FALSE, na.rm = TRUE)
} else {
0L
},
fixed_a_used = fixed_a_count,
fallback_used = fallback_count,
no_fallback = nonfallback_count,
true_two_pulse_fits = two_pulse_count,
not_two_pulse_fits = not_two_pulse_count,
season_table = season_table,
statistics_head = utils::head(fit_statistics, 10),
parameter_head = utils::head(fit_parameters, 10),
fit_statistics = fit_statistics,
fit_parameters = fit_parameters
)
class(out) <- "summary.dm_growth_fit"
out
}
#' Print a summary.dm_growth_fit object
#'
#' @description
#' Prints a formatted summary of a \code{"summary.dm_growth_fit"} object.
#'
#' @param x An object of class \code{"summary.dm_growth_fit"}.
#' @param ... Further arguments passed to or from other methods.
#'
#' @return The input object, invisibly.
#'
#' @method print summary.dm_growth_fit
#' @export
print.summary.dm_growth_fit <- function(x, ...) {
cat("Call:\n")
print(x$call)
cat("\n")
cat("dm.growth.fit summary\n")
cat("---------------------\n")
cat("Number of fitted series:", x$n_series, "\n")
cat("Number of fits:", x$n_fits, "\n")
cat("Methods used:", paste(x$methods_used, collapse = ", "), "\n")
cat("Site mode:", paste(x$site_mode, collapse = ", "), "\n")
cat("Converged fits:", x$converged, "\n")
cat("Failed fits:", x$failed, "\n")
cat("Fixed-a fits:", x$fixed_a_used, "\n")
if (!is.null(x$fallback_used)) {
cat("Fits using fallback to single curve:", x$fallback_used, "\n")
}
if (!is.null(x$true_two_pulse_fits)) {
cat("True two-pulse fits:", x$true_two_pulse_fits, "\n")
}
cat("\nSeason table:\n")
print(x$season_table)
cat("\nFit statistics preview:\n")
print(x$statistics_head)
cat("\nParameter overview:\n")
print(x$parameter_head)
invisible(x)
}
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
Documented in dm.growth.fit
#' Fit dendrometer growth curves by vegetation season
#'
#' @description
#' Fits cumulative growth curves to daily dendrometer series using one of
#' several supported methods: generalized additive model \code{"gam"},
#' Gompertz \code{"gompertz"}, logistic \code{"logistic"}, Richards
#' \code{"richards"}, local regression \code{"loess"}, or smoothing spline
#' \code{"spline"}.
#'
#' The function:
#' \itemize{
#' \item resamples dendrometer data to daily maxima using [dendro.resample()],
#' \item assigns each daily observation to a vegetation season,
#' \item optionally converts daily stem-size values to cumulative growth
#' with \code{fit_GRO = TRUE},
#' \item fits one curve per series and season with
#' \code{year_mode = "yearly"} or one pooled curve across retained seasons
#' with \code{year_mode = "pooled"},
#' \item estimates growing-season timing from the fitted cumulative-growth
#' curve,
#' \item additionally identifies active-growth timing from the fitted
#' growth-rate curve.
#' }
#'
#' For Gompertz, logistic, and Richards fits, the asymptote parameter
#' \code{a}, representing maximum seasonal growth, can optionally be fixed to
#' the observed maximum cumulative growth for that series and vegetation season.
#' This is controlled by \code{fix_a_to_observed_max}.
#'
#' @param df A data frame or tibble. The first column must be a timestamp
#' \code{Date}, \code{POSIXct}, or parseable date-time string. All selected
#' remaining 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
#' with series names.
#' @param method Fitting method. One of \code{"gam"}, \code{"gompertz"},
#' \code{"logistic"}, \code{"richards"}, \code{"loess"}, or \code{"spline"}.
#' @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.
#' Cross-year seasons are supported, for example \code{c(275, 274)}.
#' @param growth_fraction Numeric vector of length two giving the lower and
#' upper fractions of final fitted seasonal growth used to define cumulative
#' growth onset and cessation, for example \code{c(0.1, 0.9)}.
#' @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. Active-growth
#' dates are derived from the fitted growth-rate curve as the first and last
#' days where fitted growth rate exceeds this fraction of the peak fitted
#' growth rate.
#' @param fix_a_to_observed_max Logical. If \code{TRUE} and
#' \code{method = "gompertz"}, \code{"logistic"}, or \code{"richards"}, the
#' asymptote parameter \code{a}, representing maximum seasonal growth, is fixed
#' to the observed maximum cumulative growth of the corresponding series and
#' vegetation season. 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
#' \code{a} 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 just above the observed maximum.
#' @param start_value_gompertz_parameters Optional named list with starting
#' values for Gompertz fits. Supported names are \code{a}, \code{b}, and
#' \code{k}. If \code{fix_a_to_observed_max = TRUE}, supplied \code{a} is
#' ignored.
#' @param start_value_richards_parameters Optional named list with starting
#' values for Richards and logistic fits. Supported names are \code{a},
#' \code{k}, \code{t0}, and \code{v}. If
#' \code{fix_a_to_observed_max = TRUE}, supplied \code{a} is ignored.
#' @param loess_span Span used when \code{method = "loess"}.
#' @param spline_df Degrees of freedom used when \code{method = "spline"}.
#' @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}{Table with fit-level statistics and estimated timing.}
#' \item{fit_parameters}{Table with fit-level model parameters and convergence
#' information.}
#' \item{season_table}{Vegetation seasons retained for fitting.}
#' }
#'
#' @details
#' The first column of \code{df} is treated as the time column and renamed to
#' \code{"TIME"} internally. If it is not already a date-time object, the
#' function attempts to parse it using [lubridate::parse_date_time()].
#'
#' \code{growth_start_*} and \code{growth_end_*} are based on cumulative fitted
#' growth. \code{rate_start_*} and \code{rate_end_*} are based on the fitted
#' growth-rate curve.
#'
#' For Gompertz, logistic, and Richards models, \code{a} is the upper asymptote.
#' When \code{fix_a_to_observed_max = TRUE}, this parameter is not estimated by
#' nonlinear least squares. Instead, it is fixed to:
#'
#' \deqn{a = max(y_{obs}) \times fixed\_a\_multiplier}
#'
#' where \eqn{y_{obs}} is the observed cumulative growth for that series and
#' vegetation season.
#'
#' @examples
#' \donttest{
#' # fit <- dm.growth.fit(
#' # df = dendro_data,
#' # TreeNum = "all",
#' # method = "gompertz",
#' # year_mode = "yearly",
#' # fit_GRO = TRUE,
#' # fix_a_to_observed_max = TRUE
#' # )
#'
#' # fit2 <- dm.growth.fit(
#' # 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 [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 summarise
#' @importFrom tibble as_tibble tibble
#' @importFrom lubridate parse_date_time year yday month
#' @importFrom stats lm coef qlogis predict smooth.spline loess na.omit
#' @importFrom utils modifyList head
#' @export
dm.growth.fit <- function(df,
TreeNum = "all",
method = c("gam", "gompertz", "logistic", "richards", "loess", "spline"),
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,
fix_a_to_observed_max = FALSE,
fixed_a_multiplier = 1,
start_value_gompertz_parameters = list(a = NA_real_, b = NA_real_, k = NA_real_),
start_value_richards_parameters = list(a = NA_real_, k = NA_real_, t0 = NA_real_, v = 1),
loess_span = 0.2,
spline_df = 10,
verbose = TRUE) {
cl <- match.call()
TIME <- doy <- season_label <- season_start <- season_end <- season_day <- season_length <- any_observed <- NULL
method <- match.arg(method)
year_mode <- match.arg(year_mode)
site_mode <- match.arg(site_mode)
dmgf_validate_fit_growth_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,
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 two 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 two 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 <- dmgf_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),
~ . - dmgf_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),
dmgf_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]
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 <- dmgf_fit_one_curve(
x_obs = x_obs,
y_obs = y_obs,
x_pred = x_full,
season_start_date = ss,
season_length = full_len,
season_id = NULL,
method = method,
loess_span = loess_span,
spline_df = spline_df,
growth_fraction = growth_fraction,
min_unique_growth = min_unique_growth,
rate_threshold_fraction = rate_threshold_fraction,
fix_a_to_observed_max = fix_a_to_observed_max,
fixed_a_multiplier = fixed_a_multiplier,
start_value_gompertz_parameters = start_value_gompertz_parameters,
start_value_richards_parameters = start_value_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 = method,
site_mode = site_mode,
season_start = ss,
season_end = unique(dat_proc$season_end[idx])[1],
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,
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,
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_full <- dat_proc$season_day
y_full <- dat_proc[[series_name]]
sid_full <- dat_proc$season_label
good <- is.finite(y_full)
x_obs <- x_full[good]
y_obs <- y_full[good]
sid_obs <- sid_full[good]
fit_res <- dmgf_fit_one_curve(
x_obs = x_obs,
y_obs = y_obs,
x_pred = x_full,
season_start_date = as.Date(NA),
season_length = NA_real_,
season_id = sid_obs,
method = method,
loess_span = loess_span,
spline_df = spline_df,
growth_fraction = growth_fraction,
min_unique_growth = min_unique_growth,
rate_threshold_fraction = rate_threshold_fraction,
fix_a_to_observed_max = fix_a_to_observed_max,
fixed_a_multiplier = fixed_a_multiplier,
start_value_gompertz_parameters = start_value_gompertz_parameters,
start_value_richards_parameters = start_value_richards_parameters
)
fitted_dat[[series_name]] <- fit_res$pred
param_rows[[row_id]] <- tibble::tibble(
series = series_name,
fit_id = "pooled",
year_mode = year_mode,
method = method,
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,
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,
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"
)
param_cols <- c(
"series", "fit_id", "year_mode", "method", "site_mode",
"converged", "fixed_a_used", "fixed_a_value",
"a", "b", "k", "t0", "v",
"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 completed: ",
nrow(fit_parameters), " fit(s) across ",
length(unique(fit_parameters$series)), " series."
)
}
out
}
# helpers ----------------------------------------------------------------------
#' Validate inputs for 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 vector.
#' @param rate_threshold_fraction Fraction of peak rate used for rate window.
#' @param fixed_a_multiplier Multiplier used when fixing asymptote \code{a}.
#'
#' @return Invisibly returns \code{TRUE} if checks pass.
#'
#' @keywords internal
dmgf_validate_fit_growth_inputs <- function(df,
growth_fraction,
min_unique_growth,
custom_veg_season,
rate_threshold_fraction,
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(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
dmgf_is_scalar_finite <- function(x) {
is.numeric(x) && length(x) == 1 && !is.na(x) && is.finite(x)
}
#' Return the first non-missing value
#'
#' @param x Numeric vector.
#'
#' @return First non-missing value, or \code{NA_real_}.
#'
#' @keywords internal
dmgf_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
dmgf_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
dmgf_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
dmgf_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] <- dmgf_doy_to_date(yr[inside], start_doy)
season_end[inside] <- dmgf_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] <- dmgf_doy_to_date(sy, start_doy)
season_end[inside_late] <- dmgf_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] <- dmgf_doy_to_date(sy, start_doy)
season_end[inside_early] <- dmgf_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 points
#'
#' @param x Season-day values.
#' @param y Growth values.
#' @param season_id Optional season identifier for pooled fits.
#'
#' @return List with \code{x} and \code{y}.
#'
#' @keywords internal
dmgf_add_anchor_points <- function(x, y, season_id = NULL) {
if (length(x) == 0) {
return(list(x = x, y = y))
}
if (is.null(season_id)) {
if (min(x, na.rm = TRUE) > 1 && !any(x == 1)) {
x <- c(1, x)
y <- c(0, y)
}
ord <- order(x)
return(list(
x = x[ord],
y = y[ord]
))
}
season_id <- as.character(season_id)
xs <- list()
ys <- list()
uids <- unique(season_id)
for (i in seq_along(uids)) {
sid <- uids[i]
idx <- season_id == sid
x_i <- x[idx]
y_i <- y[idx]
if (length(x_i) == 0) {
next
}
if (min(x_i, na.rm = TRUE) > 1 && !any(x_i == 1)) {
x_i <- c(1, x_i)
y_i <- c(0, y_i)
}
ord <- order(x_i)
xs[[i]] <- x_i[ord]
ys[[i]] <- y_i[ord]
}
list(
x = unlist(xs, use.names = FALSE),
y = unlist(ys, use.names = FALSE)
)
}
#' Infer logistic 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
dmgf_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 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
dmgf_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
)
}
#' Empty model-parameter template
#'
#' @return Named list of empty model parameters.
#'
#' @keywords internal
dmgf_empty_model_params <- function() {
list(
converged = FALSE,
fixed_a_used = FALSE,
fixed_a_value = NA_real_,
a = NA_real_,
b = NA_real_,
k = NA_real_,
t0 = NA_real_,
v = NA_real_,
edf = NA_real_,
span = NA_real_,
spline_df = NA_real_,
spar = NA_real_,
model_note = NA_character_
)
}
#' Dispatch growth-curve model fitting
#'
#' @param method Fitting method.
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param loess_span Span for loess.
#' @param spline_df Degrees of freedom for spline.
#' @param extrapolated Logical flag for extrapolation.
#' @param fixed_a_value Optional fixed asymptote.
#' @param fix_a_to_observed_max Logical.
#' @param start_value_gompertz_parameters Starting values for Gompertz.
#' @param start_value_richards_parameters Starting values for logistic/Richards.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model <- function(method,
x,
y,
x_pred,
loess_span = 0.2,
spline_df = 10,
extrapolated = FALSE,
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_),
start_value_richards_parameters = list(a = NA_real_, k = NA_real_, t0 = NA_real_, v = 1)) {
switch(
method,
gam = dmgf_fit_model_gam(
x = x,
y = y,
x_pred = x_pred,
extrapolated = extrapolated
),
gompertz = dmgf_fit_model_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_value_gompertz_parameters
),
logistic = dmgf_fit_model_logistic(
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_value_richards_parameters
),
richards = dmgf_fit_model_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_value_richards_parameters
),
loess = dmgf_fit_model_loess(
x = x,
y = y,
x_pred = x_pred,
loess_span = loess_span,
extrapolated = extrapolated
),
spline = dmgf_fit_model_spline(
x = x,
y = y,
x_pred = x_pred,
spline_df = spline_df,
extrapolated = extrapolated
),
stop("Unknown method: ", method)
)
}
#' Fit GAM growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param extrapolated Logical flag for extrapolation.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_gam <- function(x, y, x_pred, extrapolated = FALSE) {
params <- dmgf_empty_model_params()
res <- tryCatch({
if (!requireNamespace("mgcv", quietly = TRUE)) {
stop("Package 'mgcv' is required for method = 'gam'. Please install it.")
}
x <- as.numeric(x)
y <- as.numeric(y)
x_pred <- as.numeric(x_pred)
dat <- data.frame(x = x, y = y)
ux <- sort(unique(dat$x))
if (length(ux) < 4) {
mod <- stats::lm(y ~ x, data = dat)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred))
pred <- as.numeric(pred)
pred <- pmax(pred, 0)
pred <- cummax(pred)
params$converged <- TRUE
params$model_note <- "lm fallback used because too few unique x values for GAM."
return(list(pred = pred, params = params))
}
k_gam <- min(10L, length(ux) - 1L)
k_gam <- max(4L, k_gam)
mod <- mgcv::gam(
y ~ s(x, bs = "cs", k = k_gam),
data = dat,
method = "REML",
select = TRUE
)
pred <- stats::predict(
mod,
newdata = data.frame(x = x_pred),
type = "response"
)
pred <- as.numeric(pred)
pred <- pmax(pred, 0)
pred <- cummax(pred)
st <- summary(mod)$s.table
if (!is.null(st)) {
st <- as.matrix(st)
if ("edf" %in% colnames(st)) {
params$edf <- unname(st[1, "edf"])
}
}
params$converged <- TRUE
if (isTRUE(extrapolated)) {
params$model_note <- "GAM predictions outside observed DOY range may be unstable."
}
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit Gompertz growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param fixed_a_value Optional fixed asymptote value.
#' @param fix_a_to_observed_max Logical.
#' @param start_value_gompertz_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_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 <- dmgf_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) &&
dmgf_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 <- dmgf_infer_gompertz_starts(x, y_adj, a_fixed_fit)
start_guess <- utils::modifyList(start_guess, start_value_gompertz_parameters)
if (!dmgf_is_scalar_finite(start_guess$b)) {
start_guess$b <- 0.5
}
if (!dmgf_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$model_note <- "Gompertz asymptote a fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
start_guess <- dmgf_infer_gompertz_starts(x, y_adj, a0)
start_guess <- utils::modifyList(start_guess, start_value_gompertz_parameters)
if (!dmgf_is_scalar_finite(start_guess$a)) {
start_guess$a <- a0
}
if (!dmgf_is_scalar_finite(start_guess$b)) {
start_guess$b <- 0.5
}
if (!dmgf_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
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit logistic growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param fixed_a_value Optional fixed asymptote value.
#' @param fix_a_to_observed_max Logical.
#' @param start_value_richards_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_logistic <- 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 <- dmgf_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) &&
dmgf_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 <- dmgf_infer_logistic_starts(x, y_adj, a_fixed_fit)
if (dmgf_is_scalar_finite(start_value_richards_parameters$k)) {
start_guess$k <- start_value_richards_parameters$k
}
if (dmgf_is_scalar_finite(start_value_richards_parameters$t0)) {
start_guess$t0 <- start_value_richards_parameters$t0
}
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a_fixed_fit / (1 + exp(-(k * (x - t0)))),
data = dat,
start = list(
k = start_guess$k,
t0 = start_guess$t0
),
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$converged <- TRUE
params$fixed_a_used <- TRUE
params$fixed_a_value <- a_fixed_original
params$model_note <- "Logistic asymptote a fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
start_guess <- dmgf_infer_logistic_starts(x, y_adj, a0)
if (dmgf_is_scalar_finite(start_value_richards_parameters$a)) {
start_guess$a <- start_value_richards_parameters$a
}
if (dmgf_is_scalar_finite(start_value_richards_parameters$k)) {
start_guess$k <- start_value_richards_parameters$k
}
if (dmgf_is_scalar_finite(start_value_richards_parameters$t0)) {
start_guess$t0 <- start_value_richards_parameters$t0
}
dat <- data.frame(x = x, y = y_adj)
mod <- minpack.lm::nlsLM(
y ~ a / (1 + exp(-(k * (x - t0)))),
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$k <- unname(cf["k"])
params$t0 <- unname(cf["t0"])
params$converged <- TRUE
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit Richards growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param fixed_a_value Optional fixed asymptote value.
#' @param fix_a_to_observed_max Logical.
#' @param start_value_richards_parameters Starting values.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_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 <- dmgf_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) &&
dmgf_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 <- dmgf_infer_logistic_starts(x, y_adj, a_fixed_fit)
start_vals <- list(
k = if (dmgf_is_scalar_finite(start_value_richards_parameters$k)) {
start_value_richards_parameters$k
} else {
logi_start$k
},
t0 = if (dmgf_is_scalar_finite(start_value_richards_parameters$t0)) {
start_value_richards_parameters$t0
} else {
logi_start$t0
},
v = if (dmgf_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$model_note <- "Richards asymptote a fixed to observed seasonal maximum."
return(list(pred = pred, params = params))
}
logi_start <- dmgf_infer_logistic_starts(x, y_adj, a0)
start_vals <- list(
a = if (dmgf_is_scalar_finite(start_value_richards_parameters$a)) {
start_value_richards_parameters$a
} else {
logi_start$a
},
k = if (dmgf_is_scalar_finite(start_value_richards_parameters$k)) {
start_value_richards_parameters$k
} else {
logi_start$k
},
t0 = if (dmgf_is_scalar_finite(start_value_richards_parameters$t0)) {
start_value_richards_parameters$t0
} else {
logi_start$t0
},
v = if (dmgf_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
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit loess growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param loess_span Loess span.
#' @param extrapolated Logical flag for extrapolation.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_loess <- function(x,
y,
x_pred,
loess_span = 0.2,
extrapolated = FALSE) {
params <- dmgf_empty_model_params()
res <- tryCatch({
dat <- data.frame(x = x, y = y)
mod <- stats::loess(
y ~ x,
data = dat,
span = loess_span,
degree = 2
)
pred <- stats::predict(mod, newdata = data.frame(x = x_pred))
pred <- as.numeric(pred)
params$span <- loess_span
params$edf <- if (!is.null(mod$enp)) mod$enp else NA_real_
params$converged <- TRUE
if (isTRUE(extrapolated)) {
params$model_note <- "loess extrapolation outside observed DOY range may be unreliable."
}
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Fit smoothing-spline growth curve
#'
#' @param x Observed x values.
#' @param y Observed y values.
#' @param x_pred Prediction x values.
#' @param spline_df Spline degrees of freedom.
#' @param extrapolated Logical flag for extrapolation.
#'
#' @return List with \code{pred} and \code{params}.
#'
#' @keywords internal
dmgf_fit_model_spline <- function(x,
y,
x_pred,
spline_df = 10,
extrapolated = FALSE) {
params <- dmgf_empty_model_params()
res <- tryCatch({
mod <- stats::smooth.spline(x, y, df = spline_df)
pred <- stats::predict(mod, x_pred)$y
pred <- as.numeric(pred)
params$spline_df <- mod$df
params$spar <- mod$spar
params$converged <- TRUE
if (isTRUE(extrapolated)) {
params$model_note <- "Spline predictions outside observed DOY range may be uncertain."
}
list(pred = pred, params = params)
}, error = function(e) {
params$model_note <- conditionMessage(e)
list(pred = rep(NA_real_, length(x_pred)), params = params)
})
res
}
#' Reduce prediction curve to unique x values
#'
#' @param x_pred Prediction x values.
#' @param pred Predicted y values.
#'
#' @return List with unique \code{x} and mean \code{pred}.
#'
#' @keywords internal
dmgf_unique_prediction_curve <- function(x_pred, pred) {
ok <- is.finite(x_pred) & is.finite(pred)
if (sum(ok) == 0) {
return(list(x = numeric(0), pred = numeric(0)))
}
xp <- as.numeric(x_pred[ok])
pp <- as.numeric(pred[ok])
ux <- sort(unique(xp))
mp <- vapply(
ux,
function(xx) mean(pp[xp == xx]),
numeric(1)
)
list(
x = ux,
pred = mp
)
}
#' Calculate non-negative fitted growth rate
#'
#' @param pred Fitted cumulative growth.
#'
#' @return Numeric growth-rate vector.
#'
#' @keywords internal
dmgf_growth_rate <- function(pred) {
pred <- as.numeric(pred)
if (length(pred) < 2) {
return(rep(NA_real_, length(pred)))
}
pred_mono <- cummax(pmax(pred, 0))
rate <- c(NA_real_, diff(pred_mono))
rate[!is.finite(rate)] <- NA_real_
rate <- pmax(rate, 0)
rate
}
#' Estimate active-growth window from growth rate
#'
#' @param rate Fitted growth-rate vector.
#' @param x_pred Prediction x values.
#' @param threshold_fraction Fraction of peak growth rate.
#'
#' @return List with peak rate, start day, and end day.
#'
#' @keywords internal
dmgf_rate_window <- function(rate,
x_pred,
threshold_fraction = 0.1) {
out <- list(
peak_rate = NA_real_,
start_day = NA_real_,
end_day = NA_real_
)
ok <- is.finite(rate) & is.finite(x_pred)
if (sum(ok) == 0) {
return(out)
}
r <- rate[ok]
x <- x_pred[ok]
peak_rate <- suppressWarnings(max(r, na.rm = TRUE))
if (!is.finite(peak_rate) || peak_rate <= 0) {
return(out)
}
thr <- threshold_fraction * peak_rate
above <- which(r >= thr)
if (length(above) == 0) {
return(out)
}
out$peak_rate <- peak_rate
out$start_day <- min(x[above], na.rm = TRUE)
out$end_day <- max(x[above], na.rm = TRUE)
out
}
#' Fit one growth curve
#'
#' @description
#' Internal helper used by [dm.growth.fit()] to fit one curve for one series and
#' one season, or one pooled curve.
#'
#' @param x_obs Observed x values.
#' @param y_obs Observed growth values.
#' @param x_pred Prediction x values.
#' @param season_start_date Season start date.
#' @param season_length Season length.
#' @param season_id Optional season ID for pooled fits.
#' @param method Fitting method.
#' @param loess_span Loess span.
#' @param spline_df Spline degrees of freedom.
#' @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 used for rate timing.
#' @param fix_a_to_observed_max Logical.
#' @param fixed_a_multiplier Multiplier for fixed asymptote.
#' @param start_value_gompertz_parameters Starting values for Gompertz.
#' @param start_value_richards_parameters Starting values for logistic/Richards.
#'
#' @return List with fitted predictions and parameters.
#'
#' @keywords internal
dmgf_fit_one_curve <- function(x_obs,
y_obs,
x_pred,
season_start_date = as.Date(NA),
season_length = NA_real_,
season_id = NULL,
method,
loess_span,
spline_df,
growth_fraction,
min_unique_growth,
rate_threshold_fraction = 0.1,
fix_a_to_observed_max = FALSE,
fixed_a_multiplier = 1,
start_value_gompertz_parameters,
start_value_richards_parameters) {
pred_full <- rep(NA_real_, length(x_pred))
params <- list(
converged = FALSE,
fixed_a_used = FALSE,
fixed_a_value = NA_real_,
n_obs = sum(is.finite(x_obs) & is.finite(y_obs)),
n_days_observed = NA_integer_,
first_obs_day = NA_real_,
last_obs_day = NA_real_,
season_length = season_length,
n_missing_days = NA_real_,
extrapolated = NA,
anchor_added = FALSE,
a = NA_real_,
b = NA_real_,
k = NA_real_,
t0 = NA_real_,
v = NA_real_,
edf = NA_real_,
span = NA_real_,
spline_df = NA_real_,
spar = NA_real_,
model_note = NA_character_,
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)
)
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]
sid_raw <- if (is.null(season_id)) NULL else season_id[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 <- dmgf_add_anchor_points(
x = x_raw,
y = y_raw,
season_id = sid_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", "logistic", "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 <- dmgf_fit_model(
method = method,
x = x_fit,
y = y_fit,
x_pred = x_pred,
loess_span = loess_span,
spline_df = spline_df,
extrapolated = isTRUE(params$extrapolated),
fixed_a_value = fixed_a_value,
fix_a_to_observed_max = fix_a_to_observed_max,
start_value_gompertz_parameters = start_value_gompertz_parameters,
start_value_richards_parameters = start_value_richards_parameters
)
pred_full <- as.numeric(fit_res$pred)
params[names(fit_res$params)] <- fit_res$params
timing_curve <- dmgf_unique_prediction_curve(
x_pred = x_pred,
pred = pred_full
)
timing_x <- timing_curve$x
timing_pred <- timing_curve$pred
if (length(timing_pred) > 0) {
timing_pred_mono <- cummax(pmax(timing_pred, 0))
final_val <- tail(stats::na.omit(timing_pred_mono), 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(timing_x[timing_pred_mono > low_val], na.rm = TRUE)
)
en_pos <- suppressWarnings(
max(timing_x[timing_pred_mono < 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)
}
}
}
rate <- dmgf_growth_rate(timing_pred)
rate_win <- dmgf_rate_window(
rate = rate,
x_pred = timing_x,
threshold_fraction = rate_threshold_fraction
)
params$peak_rate <- rate_win$peak_rate
if (is.finite(rate_win$start_day)) {
params$rate_start_season_day <- rate_win$start_day
if (!is.na(season_start_date)) {
params$rate_start_date <- season_start_date + (rate_win$start_day - 1)
params$rate_start_day <- lubridate::yday(params$rate_start_date)
}
}
if (is.finite(rate_win$end_day)) {
params$rate_end_season_day <- rate_win$end_day
if (!is.na(season_start_date)) {
params$rate_end_date <- season_start_date + (rate_win$end_day - 1)
params$rate_end_day <- lubridate::yday(params$rate_end_date)
}
}
}
list(
pred = pred_full,
params = params
)
}
# S3 methods -------------------------------------------------------------------
#' Print a dm_growth_fit object
#'
#' @description
#' Prints a compact overview of an object returned by [dm.growth.fit()] or
#' [dm.growth.fit.double()].
#'
#' @param x An object of class \code{"dm_growth_fit"}.
#' @param ... Further arguments passed to or from other methods.
#'
#' @return The input object, invisibly.
#'
#' @method print dm_growth_fit
#' @export
print.dm_growth_fit <- function(x, ...) {
fit_parameters <- x$fit_parameters
fit_statistics <- x$fit_statistics
if (is.null(fit_parameters) || nrow(fit_parameters) == 0) {
cat("<dm_growth_fit>\n")
cat("No fit results available.\n")
return(invisible(x))
}
methods_used <- if ("method" %in% names(fit_parameters)) {
unique(stats::na.omit(fit_parameters$method))
} else {
character(0)
}
n_series <- if ("series" %in% names(fit_parameters)) {
length(unique(stats::na.omit(fit_parameters$series)))
} else {
NA_integer_
}
n_fits <- nrow(fit_parameters)
converged_count <- if ("converged" %in% names(fit_parameters)) {
sum(fit_parameters$converged %in% TRUE, na.rm = TRUE)
} else {
NA_integer_
}
fixed_a_count <- if ("fixed_a_used" %in% names(fit_parameters)) {
sum(fit_parameters$fixed_a_used %in% TRUE, na.rm = TRUE)
} else {
0L
}
fallback_count <- if ("fallback_used" %in% names(fit_parameters)) {
sum(fit_parameters$fallback_used %in% TRUE, na.rm = TRUE)
} else {
0L
}
two_pulse_count <- if (!is.null(fit_statistics) &&
"two_pulse_detected" %in% names(fit_statistics)) {
sum(fit_statistics$two_pulse_detected %in% TRUE, na.rm = TRUE)
} else {
0L
}
cat("Call:\n")
print(x$call)
cat("\n")
cat("<dm_growth_fit>\n")
cat("Series:", n_series, "\n")
cat("Fits:", n_fits, "\n")
cat("Methods used:", paste(methods_used, collapse = ", "), "\n")
cat("Converged fits:", converged_count, "\n")
cat("Fixed-a fits:", fixed_a_count, "\n")
if ("fallback_used" %in% names(fit_parameters)) {
cat("Fallback to single curve:", fallback_count, "\n")
}
if (!is.null(fit_statistics) &&
"two_pulse_detected" %in% names(fit_statistics)) {
cat("True two-pulse fits:", two_pulse_count, "\n")
}
invisible(x)
}
#' Summarize a dm_growth_fit object
#'
#' @description
#' Summarizes an object returned by [dm.growth.fit()] or
#' [dm.growth.fit.double()].
#'
#' @param object An object of class \code{"dm_growth_fit"}.
#' @param ... Further arguments passed to or from other methods.
#'
#' @return An object of class \code{"summary.dm_growth_fit"}.
#'
#' @method summary dm_growth_fit
#' @export
summary.dm_growth_fit <- function(object, ...) {
fit_parameters <- object$fit_parameters
fit_statistics <- object$fit_statistics
season_table <- object$season_table
if (is.null(fit_parameters)) {
fit_parameters <- tibble::tibble()
}
if (is.null(fit_statistics)) {
fit_statistics <- tibble::tibble()
}
if (is.null(season_table)) {
season_table <- tibble::tibble()
}
methods_used <- if ("method" %in% names(fit_parameters)) {
unique(stats::na.omit(fit_parameters$method))
} else {
character(0)
}
fixed_a_count <- if ("fixed_a_used" %in% names(fit_parameters)) {
sum(fit_parameters$fixed_a_used %in% TRUE, na.rm = TRUE)
} else {
0L
}
fallback_count <- if ("fallback_used" %in% names(fit_parameters)) {
sum(fit_parameters$fallback_used %in% TRUE, na.rm = TRUE)
} else {
0L
}
nonfallback_count <- if ("fallback_used" %in% names(fit_parameters)) {
sum(fit_parameters$fallback_used %in% FALSE, na.rm = TRUE)
} else {
NA_integer_
}
two_pulse_count <- if ("two_pulse_detected" %in% names(fit_statistics)) {
sum(fit_statistics$two_pulse_detected %in% TRUE, na.rm = TRUE)
} else {
0L
}
not_two_pulse_count <- if ("two_pulse_detected" %in% names(fit_statistics)) {
sum(fit_statistics$two_pulse_detected %in% FALSE, na.rm = TRUE)
} else {
NA_integer_
}
out <- list(
call = object$call,
n_series = if ("series" %in% names(fit_parameters)) {
length(unique(stats::na.omit(fit_parameters$series)))
} else {
0L
},
n_fits = nrow(fit_parameters),
methods_used = methods_used,
site_mode = if ("site_mode" %in% names(fit_parameters)) {
unique(stats::na.omit(fit_parameters$site_mode))
} else {
character(0)
},
converged = if ("converged" %in% names(fit_parameters)) {
sum(fit_parameters$converged %in% TRUE, na.rm = TRUE)
} else {
0L
},
failed = if ("converged" %in% names(fit_parameters)) {
sum(fit_parameters$converged %in% FALSE, na.rm = TRUE)
} else {
0L
},
fixed_a_used = fixed_a_count,
fallback_used = fallback_count,
no_fallback = nonfallback_count,
true_two_pulse_fits = two_pulse_count,
not_two_pulse_fits = not_two_pulse_count,
season_table = season_table,
statistics_head = utils::head(fit_statistics, 10),
parameter_head = utils::head(fit_parameters, 10),
fit_statistics = fit_statistics,
fit_parameters = fit_parameters
)
class(out) <- "summary.dm_growth_fit"
out
}
#' Print a summary.dm_growth_fit object
#'
#' @description
#' Prints a formatted summary of a \code{"summary.dm_growth_fit"} object.
#'
#' @param x An object of class \code{"summary.dm_growth_fit"}.
#' @param ... Further arguments passed to or from other methods.
#'
#' @return The input object, invisibly.
#'
#' @method print summary.dm_growth_fit
#' @export
print.summary.dm_growth_fit <- function(x, ...) {
cat("Call:\n")
print(x$call)
cat("\n")
cat("dm.growth.fit summary\n")
cat("---------------------\n")
cat("Number of fitted series:", x$n_series, "\n")
cat("Number of fits:", x$n_fits, "\n")
cat("Methods used:", paste(x$methods_used, collapse = ", "), "\n")
cat("Site mode:", paste(x$site_mode, collapse = ", "), "\n")
cat("Converged fits:", x$converged, "\n")
cat("Failed fits:", x$failed, "\n")
cat("Fixed-a fits:", x$fixed_a_used, "\n")
if (!is.null(x$fallback_used)) {
cat("Fits using fallback to single curve:", x$fallback_used, "\n")
}
if (!is.null(x$true_two_pulse_fits)) {
cat("True two-pulse fits:", x$true_two_pulse_fits, "\n")
}
cat("\nSeason table:\n")
print(x$season_table)
cat("\nFit statistics preview:\n")
print(x$statistics_head)
cat("\nParameter overview:\n")
print(x$parameter_head)
invisible(x)
}
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.