Habitat Sensitivity, Rates and Seasonality
Heidi K. Hirsh
2025-10-15
1 Set up
1.1 Settings
doc_dir <- knitr::opts_knit$get("root.dir")
if (is.null(doc_dir) || !nzchar(doc_dir)) doc_dir <- dirname(rstudioapi::getActiveDocumentContext()$path)
project_dir <- normalizePath(doc_dir)
# project_dir <- "/Users/heidi.k.hirsh/Documents/GitHub/sensitivity_repo"
slopes_dir <- "Tables/Slopes"
cc_path_rel <- "InputData/CCc_plusArea.csv"
# --- params with safe fallbacks when knitting interactively ---
ndays_focus <- if (exists("params") && !is.null(params$ndays_focus)) params$ndays_focus else 5L
zone_focus <- if (exists("params") && !is.null(params$zone_focus)) params$zone_focus else "Inshore"
hab_focus <- if (exists("params") && !is.null(params$hab_focus)) params$hab_focus else "SG"
y_label_area <- switch(
toupper(hab_focus),
"SG" = "Seagrass Area (ha)",
"ALG" = "Algae Area (ha)",
"CALC" = "Calcifier Area (ha)",
"Area (ha)"
)
# define readable names
hab_labels <- c(
"SG" = "Seagrass",
"CALC" = "Calcifiers",
"ALG" = "Noncalcifying Algae"
)
metric_labels <- c(
"DIC" = "Dissolved Inorganic Carbon",
"TA" = "Total Alkalinity"
)
# pick the right names
hab_name <- hab_labels[params$hab_focus]
metric_name <- metric_labels[params$metric]
zone_safe <- gsub("[^A-Za-z0-9]+", "", zone_focus)
loess_span <- 0.5
# loess_span <- 0.3
window_days <- 1
save_figs <- TRUE
set.seed(42)
# --- timestamp (define BEFORE using in tag helpers) ---
TS <- function() format(Sys.time(), "%Y%m%d_%H%M")
stamp <- TS()
# --- habitat + metric controls ---
HAB_LIST <- c("SG","ALG","CALC")
METRIC <- if (exists("params") && !is.null(params$metric)) {
match.arg(toupper(params$metric), c("DIC","TA"))
} else {"DIC"}
TARGET_YR <- 5.35 # µmol kg^-1 yr^-1 to offset
# --- filename helpers add habitat + metric + ndz/zone ---
tag_hab <- function(h) paste0("_", toupper(h), "_", METRIC)
tag_ndz <- paste0("_nd", ndays_focus, "_", zone_safe, "_", stamp)
# In Settings (knobs) — set this ABOVE the auto-picker block or remove the auto-picker
# --- fixed slopes file path (explicit, no guessing) ---
# slopes_file <- "/Users/heidi.k.hirsh/Documents/GitHub/sensitivity_repo/Tables/Slopes/slopes_per_visit_all_ndays_ALLHABS_20251007_2018.csv"
slopes_file <- params$slopes_file
if (!file.exists(slopes_file)) {
stop("slopes_file not found at:\n", slopes_file,
"\nPlease confirm the path and filename.")
} else {
message("Using fixed slopes_file: ", slopes_file)
}
# --- output dirs ---
out_dir <- file.path(project_dir, "Tables", "AnnualRates")
fig_dir <- file.path(project_dir, "Figures")
dir.create("Reports", recursive = TRUE, showWarnings = FALSE) # <-- root-level Reports
dir.create(file.path(fig_dir, "JdayPlots"), recursive = TRUE, showWarnings = FALSE)
dir.create(file.path(fig_dir, "AnnualRates"), recursive = TRUE, showWarnings = FALSE)
dir.create(file.path(fig_dir, "SeasonalRates"), recursive = TRUE, showWarnings = FALSE)
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
# --- guardrail message ---
message("Knit settings → ",
"Habitat = ", hab_focus,
" | Metric = ", METRIC,
" | ndays = ", ndays_focus,
" | Zone = ", zone_focus,
" | Save figs = ", save_figs,
" | Timestamp = ", stamp)
message("Working dir: ", getwd())
message("HTML/PDF will be written under: ", normalizePath("Reports", mustWork = FALSE))1.2 Settings summary
Parameters:
Habitat: CALC
Metric: DIC
ndays: 5
Zone: Inshore
Save figures: TRUE
Timestamp: 20251015_1125
Slopes file: slopes_per_visit_all_ndays_ALLHABS_20251007_2018.csv
Habitat: CALC
Metric: DIC
ndays: 5
Zone: Inshore
Save figures: TRUE
Timestamp: 20251015_1125
Slopes file: slopes_per_visit_all_ndays_ALLHABS_20251007_2018.csv
1.3 Load libraries
suppressPackageStartupMessages({
library(tidyverse); library(readr); library(lubridate); library(scales); library(ggplot2);
library(ggbreak); library(purrr); library(tidyr); library(forcats); library(patchwork); library(dplyr)})1.4 Order factors
# ---- factor orders ----
site_levels <- c("1","2","3","EK_IN","EK_MID","EK_OFF","4","5","5.5","6","6.5","MR","UK_IN","UK_MID","UK_OFF","7","8","9","9.5","10","11","12","13","14","15","15.5","16","17","18","19","20","21","LK","21.5","WS","24","23","22","22.5")
season_levels <- c("Winter","Spring","Summer","Fall")
subreg_levels <- c("BB","UK","MK","LK")
zone_levels <- c("Inshore","Mid channel","Offshore","Oceanic")1.5 Plotting settings
# Title suffix used across all plots/files
title_suffix <- paste0("(ndays = ", ndays_focus, ", zone = ", zone_focus, ")")
# Convenience: wrap a base title with the suffix
title_ndz <- function(txt) paste(txt, title_suffix)
season_pal <- c(
"Winter" = "#2C7BB6", # blue
"Spring" = "#1FA187", # green
"Summer" = "#FDB863", # warm yellow/orange
"Fall" = "#D7191C" # red
)
sub_region_pal <- c(
"Biscayne Bay" = "#1b9e77",
"Upper Keys" = "#d95f02",
"Middle Keys" = "#7570b3",
"Lower Keys" = "#e7298a"
)
# SiteID often has many levels → generate on demand, but stable for an explicit set
# make_site_pal <- function(levels) {grDevices::setNames(grDevices::hcl.colors(length(levels), palette = "Dark3"), levels)}
SITE_PAL_MASTER <- setNames(hcl.colors(length(site_levels), palette = "Dark3"), site_levels)1.6 Helper functions
## ---- helpers ---------------------------------------------------------------
## These utilities keep plotting consistent across data frames, provide
## common scales/palettes, and support seasonality + LOESS smoothing.
# -----------------------------------------------------------------------------
# build_common_x_levels()
# -----------------------------------------------------------------------------
#' Build a canonical order of x-axis categories across many data frames.
#'
#' @param dfs List of data frames that each contain an x column (e.g., SiteID).
#' @param x_col Name of the x column as a string (default "SiteID").
#' @param canonical Optional vector with a preferred global order
#' (e.g., site_levels). If given, we keep this order and drop
#' any entries not present in the data; otherwise we sort unique().
#' @return Character vector of x-levels to use everywhere.
build_common_x_levels <- function(dfs, x_col = "SiteID", canonical = NULL) {
# All unique values present across the provided data frames
vals <- Reduce(union, lapply(dfs, function(d) as.character(d[[x_col]])))
# Keep canonical order if provided; otherwise alphabetical
if (!is.null(canonical)) {
vals <- intersect(canonical, vals)
} else {
vals <- sort(vals)
}
vals
}
# -----------------------------------------------------------------------------
# apply_x_levels()
# -----------------------------------------------------------------------------
#' Apply a common set of factor levels to one data frame.
#'
#' @param df Data frame.
#' @param x_col Name of the x column as a string.
#' @param levels_vec Vector of levels in the desired order.
#' @return Data frame with x_col coerced to an ordered factor.
apply_x_levels <- function(df, x_col = "SiteID", levels_vec) {
df[[x_col]] <- factor(df[[x_col]], levels = levels_vec)
df
}
# -----------------------------------------------------------------------------
# scale_x_locked()
# -----------------------------------------------------------------------------
#' Discrete x scale with fixed limits across plots (prevents reordering/dropping).
#'
#' @param levels_vec Levels returned by build_common_x_levels().
#' @param expand_add Additive expansion (left/right padding in discrete units).
#' @param drop Whether to drop unused levels (default FALSE to preserve gaps).
#' @return A ggplot2 discrete x scale.
scale_x_locked <- function(levels_vec, expand_add = 0.6, drop = FALSE) {
ggplot2::scale_x_discrete(
limits = levels_vec,
drop = drop,
expand = ggplot2::expansion(add = expand_add)
)
}
# -----------------------------------------------------------------------------
# Reusable color/fill scales
# -----------------------------------------------------------------------------
# Keep plotting code terse and consistent with central palettes.
## Season scales
scale_color_season <- function(...) {
ggplot2::scale_color_manual(values = season_pal, limits = names(season_pal), ...)
}
scale_fill_season <- function(...) {
ggplot2::scale_fill_manual(values = season_pal, limits = names(season_pal), ...)
}
## Subregion scales
scale_color_subregion <- function(...) {
ggplot2::scale_color_manual(values = sub_region_pal, limits = names(sub_region_pal), ...)
}
scale_fill_subregion <- function(...) {
ggplot2::scale_fill_manual(values = sub_region_pal, limits = names(sub_region_pal), ...)
}
## Site scales (auto-subset to sites present in the data, while preserving global order)
# Helper to grab just the sites present in a given data frame, but keeping the global order from `site_levels`
present_sites <- function(data) {
x <- as.character(droplevels(data$SiteID))
present <- unique(stats::na.omit(x))
intersect(site_levels, present) # keep global order, drop unused
}
# Use the master site palette but restrict to present sites (preserves color stability across plots)
scale_color_site <- function(data, ...) {
lvls <- present_sites(data)
vals <- SITE_PAL_MASTER[lvls]
scale_color_manual(values = vals, limits = lvls, ...)
}
scale_fill_site <- function(data, ...) {
lvls <- present_sites(data)
vals <- SITE_PAL_MASTER[lvls]
scale_fill_manual(values = vals, limits = lvls, ...)
}
# -----------------------------------------------------------------------------
# season_from_jday()
# -----------------------------------------------------------------------------
#' Map Julian day (1–365) to meteorological seasons using a non-leap anchor year.
#'
#' @details Anchor date is "2023-01-01" to avoid 366 creeping in.
#' Winter: Dec–Feb, Spring: Mar–May, Summer: Jun–Aug, Fall: Sep–Nov.
#' @param j Integer Julian day (1–365). Values outside are clamped elsewhere.
#' @return Character season name ("Winter", "Spring", "Summer", "Fall").
season_from_jday <- function(j) {
d <- ymd("2023-01-01") + days(j - 1)
m <- month(d)
case_when(
m %in% c(12, 1, 2) ~ "Winter",
m %in% 3:5 ~ "Spring",
m %in% 6:8 ~ "Summer",
TRUE ~ "Fall"
)
}
# -----------------------------------------------------------------------------
# fill_edges()
# -----------------------------------------------------------------------------
#' Forward/backward-fill the leading/trailing NAs of a numeric vector.
#'
#' @details Useful after predict(loess) where edges can be NA; prevents
#' blank strips at the start/end of the year when plotting.
#' Internal NAs are left as-is.
#' @param x Numeric vector (possibly with NAs).
#' @return Numeric vector with only the leading/trailing NAs filled.
fill_edges <- function(x) {
if (length(x) == 0 || all(is.na(x))) return(x)
i1 <- which(!is.na(x))[1]; i2 <- tail(which(!is.na(x)), 1)
if (!is.na(i1) && i1 > 1) x[1:(i1 - 1)] <- x[i1]
if (!is.na(i2) && i2 < length(x)) x[(i2 + 1):length(x)] <- x[i2]
x
}
# -----------------------------------------------------------------------------
# fit_loess_and_predict_daily()
# -----------------------------------------------------------------------------
#' Fit a site-level LOESS of slope vs Julian day, predict daily values.
#'
#' @details
#' • Requires ≥ 5 observations per site to fit (else returns empty tibble).
#' • Clamps jday to [1, 366] but forces a 365-day prediction grid to avoid
#' leap-year issues.
#' • If your raw slopes are computed over rolling windows >1 day, pass the
#' window size via `window_days` to convert to per-day rates.
#' • Assumes the slope column is named `slope`. If you are using a
#' metric-specific column (e.g., `slope_use`), rename before calling,
#' or adapt this function to use that column.
#'
#' @param df_site Data frame for a single site (already grouped/split).
#' @param span LOESS span (default 0.3).
#' @param window_days Size of the slope window (1 = already per day).
#' @return Tibble with columns: jday (1–365), pred_per_day (numeric).
fit_loess_and_predict_daily <- function(df_site, span = 0.3, window_days = 1) {
d <- df_site %>%
tidyr::drop_na(jday, slope) %>% # drop rows missing jday or slope
dplyr::mutate(
jday = as.integer(jday) |> pmin(366L) |> pmax(1L), # clamp to [1, 366]
slope_per_day = slope / window_days # convert to per-day units if needed
)
# Require at least 5 points to fit a smoother
if (nrow(d) < 5) return(tibble(jday = integer(), pred_per_day = numeric()))
grid_days <- 1:365 # force 365-day grid
newd <- tibble(jday = grid_days)
# Fit smoothed trend (LOESS, degree=1 is locally linear)
mod <- loess(slope_per_day ~ jday, data = d, span = span, degree = 1,
control = loess.control(surface = "direct"))
# Predict and fill leading/trailing NAs so edges aren’t blank
preds <- as.numeric(predict(mod, newdata = newd)) %>% fill_edges()
newd %>% mutate(pred_per_day = preds)
}
# -----------------------------------------------------------------------------
# norm_hab_col()
# -----------------------------------------------------------------------------
#' Ensure a standardized Habitat factor exists with levels SG / ALG / CALC.
#'
#' @details
#' Looks for one of the following columns: "Habitat", "habitat", or "habitat_col".
#' If none exist or if values are missing (NA or ""), all rows are filled with
#' a default habitat (by default "SG") so existing single-habitat pipelines
#' continue to work seamlessly.
#'
#' @param df Data frame that may or may not include a habitat column.
#' @param default Character. Habitat to assign when values are missing (default = "SG").
#' @return Data frame with a standardized `Habitat` factor with levels c("SG", "ALG", "CALC").
# Ensure a standardized Habitat factor exists with levels SG/ALG/CALC,
# and fill missing values with a default (SG by default).
# norm_hab_col <- function(df, default = "SG") {
# cand <- intersect(names(df), c("Habitat","habitat","habitat_col"))
# if (length(cand)) {
# df <- df %>% dplyr::rename(Habitat = !!cand[1])
# } else {
# df <- df %>% dplyr::mutate(Habitat = NA_character_)
# }
# df %>%
# dplyr::mutate(
# Habitat = toupper(trimws(as.character(Habitat))),
# Habitat = dplyr::if_else(is.na(Habitat) | Habitat == "", default, Habitat),
# Habitat = factor(Habitat, levels = c("SG","ALG","CALC"))
# )
# }
# Standardize Habitat from habitat_col ("SGi_m2", "ALGi_m2", "CALC_m2")
norm_hab_col <- function(df) {
if (!"habitat_col" %in% names(df)) {
warning("No 'habitat_col' column found — defaulting Habitat = 'SG'")
df$Habitat <- factor("SG", levels = c("SG","ALG","CALC"))
return(df)
}
df <- df %>%
dplyr::mutate(
Habitat = dplyr::case_when(
habitat_col == "SGi_m2" ~ "SG",
habitat_col == "ALGi_m2" ~ "ALG",
habitat_col == "CALC_m2" ~ "CALC",
TRUE ~ NA_character_
),
Habitat = factor(Habitat, levels = c("SG","ALG","CALC"))
)
df
}
# -----------------------------------------------------------------------------
# build_focus()
# -----------------------------------------------------------------------------
#' Slice CCc_with_slopes by Habitat, ndays, and Zone.
#'
#' @details Creates jdate for plotting and renames slope_use → slope
#' so older code (LOESS, plots, etc.) continues to work.
#'
#' @param hab Character, one of "SG", "ALG", "CALC".
#' @return Filtered data frame for that habitat/zone.
build_focus <- function(hab = "SG") {
hab <- toupper(hab)
out <- CCc_with_slopes %>%
dplyr::filter(Habitat == hab,
ndays == ndays_focus,
Zone == zone_focus) %>%
dplyr::mutate(
jday = pmin(pmax(as.integer(jday), 1L), 365L),
jdate = as.Date("2024-01-01") + (jday - 1),
slope = slope_use # expose for legacy plotting code
)
message("build_focus: ", nrow(out), " rows for Habitat=", hab,
", ndays=", ndays_focus, ", Zone=", zone_focus)
out
}2 Load Data
2.1 Sensitivity slopes
## ---- Step 2: Load & normalize slopes (use slopes_file from Settings) ----
# DO NOT force a stop here anymore — the Settings chunk already validated it.
message("Reading slopes from: ", slopes_file)
slopes_all <- readr::read_csv(slopes_file, show_col_types = FALSE) %>%
dplyr::mutate(
visitID = trimws(as.character(visitID)),
ndays = as.integer(ndays),
Zone = factor(trimws(as.character(Zone)), levels = zone_levels),
Season = factor(Season, levels = season_levels),
Sub_region = factor(Sub_region), # levels applied later
SiteID = factor(as.character(SiteID), levels = site_levels)
) %>%
norm_hab_col() %>%
{ if ("chemistry" %in% names(.)) dplyr::filter(., toupper(chemistry) == toupper(METRIC)) else . } %>%
dplyr::mutate(
slope_use = dplyr::coalesce(
!!!dplyr::select(., tidyselect::any_of(c("slope","slope_use","slope_DIC","slope_TA")))
)
)
if (!("slope_use" %in% names(slopes_all)) || !any(is.finite(slopes_all$slope_use))) {
stop("No usable slope values found (checked slope, slope_use, slope_DIC, slope_TA).")
}
message("Habitats present in slopes: ",
paste(unique(as.character(slopes_all$Habitat)), collapse = ", "))
message("Param slopes_file seen by Rmd: ", slopes_file)2.2 Carbonate chemistry
## ---- Step 3: Load CCc, join to slopes by (visitID, ndays) ----
# 1) Load CCc
CCc <- readr::read_csv(file.path(project_dir, cc_path_rel), show_col_types = FALSE) %>%
dplyr::mutate(
visitID = as.character(visitID),
ndays = as.integer(ndays),
Sub_region = factor(Sub_region,
levels = subreg_levels,
labels = c("Biscayne Bay","Upper Keys","Middle Keys","Lower Keys")),
Season = factor(Season, levels = season_levels),
SiteID = factor(as.character(SiteID), levels = site_levels),
Zone = factor(as.character(Zone), levels = zone_levels)
)
# 2) Enforce uniqueness of slopes per (visitID, ndays, Habitat)
# (If duplicates exist, keep the first or choose a rule you prefer.)
slopes_keyed <- slopes_all %>%
dplyr::arrange(visitID, ndays, Habitat) %>% # or arrange by quality (e.g., desc(n_points), asc(slope_se))
dplyr::distinct(visitID, ndays, Habitat, .keep_all = TRUE) %>%
dplyr::select(visitID, ndays, Habitat, slope_use)
# # Count duplicate combinations of (visitID, ndays, Habitat)
# dupes <- slopes_all %>%
# dplyr::count(visitID, ndays, Habitat) %>%
# dplyr::filter(n > 1)
#
# # View how many there are and which they are
# nrow(dupes) # total number of duplicated combos
# head(dupes, 10) # peek at the first few
# 3) Join by the minimal keys
CCc_with_slopes <- CCc %>%
dplyr::left_join(slopes_keyed, by = c("visitID","ndays"))
# 4) Sanity checks
if (!any(is.finite(CCc_with_slopes$slope_use))) {
n_overlap <- dplyr::inner_join(
dplyr::distinct(CCc, visitID, ndays),
dplyr::distinct(slopes_keyed, visitID, ndays),
by = c("visitID","ndays")
) %>% nrow()
stop("After join: no usable slope values. Check that visitID/ndays match between files.\n",
"Overlapping (visitID, ndays) pairs: ", n_overlap)
}
message("Joined rows: ", nrow(CCc_with_slopes))
# nrow(CCc) # total rows in CCc (13776)
# length(unique(CCc$visitID)) # how many distinct visitIDs (984)
# length(unique(CCc$ndays)) # how many ndays values (14)
# length(unique(slopes_all$Habitat)) # usually 3 (3)
message("Habitats present post-join: ",
paste(unique(na.omit(as.character(CCc_with_slopes$Habitat))), collapse = ", "))2.3 Subset: Calcifiers, Dissolved Inorganic Carbon
# IN <- CCc_with_slopes %>%
# filter(ndays == ndays_focus, Zone == zone_focus) %>%
# mutate(jdate = ymd("2024-01-01") + days(jday - 1))
# IN <- build_focus("SG") # generalize to ALG/CALC later
IN <- build_focus(hab_focus)
# # 1. What habitats are actually in the joined table?
# unique(CCc_with_slopes$Habitat)
#
# # 2. What ndays values are present?
# unique(CCc_with_slopes$ndays)
#
# # 3. What zone names exist?
# unique(CCc_with_slopes$Zone)## ---- jday-preview (same views, clearer names) -------------------------
current_hab <- as.character(unique(IN$Habitat))[1]
month_letters <- function(x) substr(as.character(lubridate::month(x, label = TRUE, abbr = TRUE)), 1, 1)
# 1) Site facets • points+ribbon by Sub_region (your old "B")
p_site_facets__subregion_loess <- ggplot(IN, aes(x = jdate, y = slope, color = Sub_region, fill = Sub_region)) +
geom_point(alpha = 0.5) +
geom_hline(yintercept = 0, color = "black", linewidth = 0.3) +
geom_smooth(method = "loess", se = TRUE, span = loess_span, alpha = 0.2) +
scale_color_subregion() +
scale_fill_subregion() +
scale_x_date(date_breaks = "1 month", labels = month_letters) +
labs(
title = title_ndz(paste0("Figure 1. Slope of Δ", METRIC, " vs jday — ", current_hab)),
x = "Month",
y = bquote("Slope (" * Delta * .(METRIC) ~ mu * "mol kg"^{-1} ~ " km"^{-2} * ")")
) +
facet_wrap(~ SiteID) +
theme_bw()
p_site_facets__subregion_loess
# 2) Site facets • points colored by Season; smooth lines by Sub_region (your old "C")
p_site_facets__season_points_subregion_lines <- ggplot(IN, aes(x = jdate, y = slope)) +
geom_point(aes(color = Season), alpha = 0.5) +
geom_hline(yintercept = 0, color = "black", linewidth = 0.3) +
geom_smooth(aes(color = Sub_region), method = "loess", se = FALSE, span = loess_span) +
scale_color_season() +
scale_x_date(date_breaks = "1 month", labels = month_letters) +
labs(
title = title_ndz(paste0("Figure 2. Slope of Δ", METRIC, " vs jday — ", current_hab)),
x = "Month",
y = bquote("Slope (" * Delta * .(METRIC) ~ mu * "mol kg"^{-1} ~ " km"^{-2} * ")")
) +
facet_wrap(~ SiteID, scales = "free_y") + # same behavior as your original "C"
theme_bw()
p_site_facets__season_points_subregion_lines
# 3) Subregion facets • x=jday with month labels; color/fill by SiteID (your old "p_jday")
p_subregion_facets__site_colors <- ggplot(IN, aes(x = jday, y = slope, color = SiteID, fill = SiteID)) +
geom_hline(yintercept = 0, linewidth = 0.3, color = "black") +
geom_point(alpha = 0.7) +
geom_smooth(method = "loess", se = TRUE, span = loess_span, alpha = 0.2) +
scale_color_site(IN) +
scale_fill_site(IN) +
scale_x_continuous(
breaks = c(15,46,75,105,135,166,196,227,258,288,319,349),
labels = month.abb
) +
labs(
title = title_ndz(paste0("Figure 3. Slope vs Julian day — ", current_hab)),
x = "Month",
y = bquote(Delta * .(METRIC) ~ " slope per km"^2)
) +
facet_wrap(~ Sub_region, ncol = 1) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_subregion_facets__site_colors
#---- save previews ------------------------------------
if (save_figs) {
ggsave(file.path(fig_dir, "JdayPlots",
paste0("Fig1.jday_site-facets_subregion-loess", tag_hab(current_hab), tag_ndz, ".png")),
p_site_facets__subregion_loess, width = 11, height = 8, dpi = 300)
ggsave(file.path(fig_dir, "JdayPlots",
paste0("Fig2.jday_site-facets_season-points_subregion-lines", tag_hab(current_hab), tag_ndz, ".png")),
p_site_facets__season_points_subregion_lines, width = 11, height = 8, dpi = 300)
ggsave(file.path(fig_dir, "JdayPlots",
paste0("Fig3.jday_subregion-facets_site-colors", tag_hab(current_hab), tag_ndz, ".png")),
p_subregion_facets__site_colors, width = 9, height = 8, dpi = 300)
}3 Seasonal and Annual Rates
3.1 Quarterly Approach
This section uses seasonal medians (Winter, Spring, Summer, Fall)
based on Season in CCc.
## ---- get-rates-quarterly-method, message=FALSE ------------------------
## ---- quarterly-rates (habitat-aware, uses IN$slope) -------------------
# NOTE: IN already has:
# - Habitat filtered (e.g., "SG")
# - ndays == ndays_focus, Zone == zone_focus
# - slope column set (from slope_use)
# - Season/Sub_region/SiteID factors ready
current_hab <- as.character(unique(IN$Habitat))[1] # e.g., "SG"
season_days <- 365.25 / 4 # = 91.3125
# 1) Seasonal median daily slope, then multiply by seasonal days
SSm <- IN %>%
dplyr::group_by(Habitat, Sub_region, SiteID, Season) %>%
dplyr::summarise(mdd = median(slope, na.rm = TRUE), .groups = "drop") %>%
dplyr::mutate(
season_total_km2 = mdd * season_days, # per km^2
season_total_ha = season_total_km2 / 100
)
# 2) Annual totals per site (sum seasons), then required area to offset TARGET_YR
nAnn <- SSm %>%
dplyr::group_by(Habitat, Sub_region, SiteID) %>%
dplyr::summarise(annual_ha = sum(season_total_ha, na.rm = TRUE), .groups = "drop") %>%
dplyr::mutate(
annual_m2 = annual_ha / 1e4, # per m^2
HA_offOA = -TARGET_YR / annual_m2 / 1e4, # hectares to offset TARGET_YR
HA_offOA = dplyr::if_else(is.finite(HA_offOA), HA_offOA, NA_real_)
)
## 3) Diagnostic scatter (per m2 vs hectares needed)
# diag_quarterly <- ggplot(nAnn, aes(x = annual_m2, y = HA_offOA, color = Sub_region)) +
# geom_hline(yintercept = 0, linewidth = 0.3) +
# geom_vline(xintercept = 0, linewidth = 0.3, linetype = "dashed") +
# geom_point(size = 3, alpha = 0.8, na.rm = TRUE) +
# scale_color_subregion() +
# scale_x_continuous(labels = scales::scientific) +
# scale_y_continuous(labels = scales::comma) +
# labs(
# title = title_ndz(paste0("Diagnostic (Quarterly, ", current_hab, "): required ha vs annual per m\u00B2 [", METRIC, "]")),
# x = paste0("Annual \u0394", METRIC, " per m\u00B2"),
# y = "Required area (ha)"
# ) +
# theme_bw()
# 4) Annual bar: \u0394 per ha (quarterly)
p_quarterly_annual <- ggplot(nAnn, aes(x = SiteID, y = annual_ha, fill = Sub_region)) +
geom_hline(yintercept = 0, linewidth = 0.3) +
geom_col(width = 0.7) +
scale_fill_subregion() +
labs(
title = title_ndz(paste0("Figure 4. Annual \u0394 per ha [", METRIC, "] — Quarterly (", current_hab, ")")),
x = "SiteID",
y = paste0("Annual \u0394", METRIC, " (\u03BCmol kg^-1) per ha")
) +
scale_y_continuous(labels = scales::label_number(big.mark = ",")) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
geom_text(aes(label = round(annual_ha, 2),
vjust = ifelse(annual_ha >= 0, -0.3, 1.2)),
size = 3)
# # 5) Area (ha) required to offset TARGET_YR — Quarterly
# p_quarterly_area <- ggplot(nAnn, aes(x = SiteID, y = HA_offOA, fill = Sub_region)) +
# geom_col(width = 0.7) +
# scale_fill_subregion() +
# ggbreak::scale_y_break(c(40, 600)) + # tweak break as needed
# geom_text(aes(label = scales::comma(HA_offOA, accuracy = 0.1)),
# vjust = 1.2, color = "white", size = 3, na.rm = TRUE) +
# labs(
# title = bquote(Area~(ha)~required~to~offset~.(TARGET_YR)~mu~mol~kg^{-1}~yr^{-1}~
# "[" * .(METRIC) * "] — Quarterly (" * .(current_hab) * ")"),
# x = "SiteID",
# y = "Area (ha)"
# ) +
# theme_bw() +
# theme(axis.text.x = element_text(angle = 45, hjust = 1))
# 6) Stacked seasonal per ha — Quarterly + net labels at the top
quarterly_season <- SSm %>%
mutate(
Season = factor(Season, levels = c("Winter","Spring","Summer","Fall")),
METRIC_season_ha = season_total_ha
)
# ---- labels: place above positive stack, show NET value ----
q_labels <- quarterly_season %>%
dplyr::group_by(SiteID) %>%
dplyr::summarise(
pos_top = sum(pmax(METRIC_season_ha, 0), na.rm = TRUE), # height of positive stack
net = sum(METRIC_season_ha, na.rm = TRUE), # NET (pos + neg)
.groups = "drop"
) %>%
dplyr::mutate(
SiteID = factor(SiteID, levels = levels(quarterly_season$SiteID)),
y_off = 0.05 * max(abs(c(pos_top, net)), na.rm = TRUE),
y_lab = ifelse(pos_top > 0, pos_top + y_off, 0 + y_off), # placement
lab_val = net # text = NET
)
p_quarterly_stacked <- ggplot(quarterly_season,
aes(SiteID, METRIC_season_ha, fill = Season)) +
geom_hline(yintercept = 0, linewidth = 0.3) +
geom_col(position = "stack") +
scale_fill_season() +
labs(
title = title_ndz(paste0("Figure 5. Seasonal \u0394", METRIC, " per ha — Quarterly (", current_hab, ")")),
x = "SiteID",
y = paste0("\u0394", METRIC, " per ha")
) +
geom_text(
data = q_labels,
mapping = aes(x = SiteID, y = y_lab, label = round(lab_val, 2)),
inherit.aes = FALSE,
size = 3,
fontface = "bold"
) +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.12))) +
coord_cartesian(clip = "off") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# # print
# p_quarterly_annual
# p_quarterly_area
# p_quarterly_stacked
# diag_quarterly
# 7) Save (habitat-aware filenames)
if (save_figs) {
ggsave(file.path(fig_dir, "AnnualRates",
paste0("Fig4.annual_quarterly", tag_hab(current_hab), tag_ndz, ".png")),
p_quarterly_annual, width = 9, height = 5, dpi = 300)
# ggsave(file.path(fig_dir, "AnnualRates",
# paste0("area_quarterly", tag_hab(current_hab), tag_ndz, ".png")),
# p_quarterly_area, width = 9, height = 5, dpi = 300)
# ggsave(file.path(fig_dir, "AnnualRates",
# paste0("diagnostic_quarterly", tag_hab(current_hab), tag_ndz, ".png")),
# diag_quarterly, width = 8, height = 5, dpi = 300)
ggsave(file.path(fig_dir, "SeasonalRates",
paste0("Fig5.quarterly_seasonal_stacked", tag_hab(current_hab), tag_ndz, ".png")),
p_quarterly_stacked, width = 11, height = 6, dpi = 300)
}
# 8) Save tables (habitat-aware filenames)
readr::write_csv(SSm,
file.path(out_dir, paste0("seasonal_medians", tag_hab(current_hab), tag_ndz, ".csv")))
readr::write_csv(nAnn,
file.path(out_dir, paste0("annual_quarterly", tag_hab(current_hab), tag_ndz, ".csv")))
# print
p_quarterly_annual
# p_quarterly_area
p_quarterly_stacked
# diag_quarterly3.2 Julian Day Approach
stopifnot(exists("IN"))
current_hab <- as.character(unique(IN$Habitat))[1]
# 1) Deterministic season calendar (1–365 → Winter/Spring/Summer/Fall)
season_calendar <- tibble(jday = 1:365) %>%
mutate(
Season_cal = cut(
jday,
breaks = c(0, 90, 182, 273, 365), # 1–90, 91–182, 183–273, 274–365
labels = c("Winter","Spring","Summer","Fall"),
include.lowest = TRUE, right = TRUE
),
Season_cal = factor(Season_cal, levels = season_levels)
)
# 2) LOESS predictions per site → daily per-ha contributions
daily_preds <- IN %>%
dplyr::transmute(Habitat, Sub_region, SiteID, jday = as.integer(jday), slope) %>%
dplyr::group_by(Habitat, Sub_region, SiteID) %>%
dplyr::group_modify(~ fit_loess_and_predict_daily(.x, span = loess_span, window_days = window_days)) %>%
dplyr::ungroup() %>%
dplyr::left_join(season_calendar, by = "jday") %>%
dplyr::mutate(contrib_ha = pred_per_day / 100)
# 3) Seasonal sums (per ha) using Season_cal; then annual sums
seasonal_jday <- daily_preds %>%
dplyr::group_by(Habitat, Sub_region, SiteID, Season_cal) %>%
dplyr::summarise(METRIC_season_ha = sum(contrib_ha, na.rm = TRUE), .groups = "drop")
annual_jday <- seasonal_jday %>%
dplyr::group_by(Habitat, Sub_region, SiteID) %>%
dplyr::summarise(annual_ha = sum(METRIC_season_ha, na.rm = TRUE), .groups = "drop")
# 4) Area required to offset TARGET_YR (ha), with safe handling of zeros
annual_jday_area <- annual_jday %>%
dplyr::mutate(
annual_m2 = annual_ha / 1e4,
HA_offOA = -TARGET_YR / annual_m2 / 1e4,
HA_offOA = dplyr::if_else(is.finite(HA_offOA), HA_offOA, NA_real_)
)
# 5) Save tables (habitat-aware filenames)
readr::write_csv(seasonal_jday,
file.path(out_dir, paste0("seasonal_jday", tag_hab(current_hab), tag_ndz, ".csv")))
readr::write_csv(annual_jday,
file.path(out_dir, paste0("annual_jday", tag_hab(current_hab), tag_ndz, ".csv")))
readr::write_csv(annual_jday_area,
file.path(out_dir, paste0("annual_jday_area", tag_hab(current_hab), tag_ndz, ".csv")))
# 6) (Optional) Save LOESS models for this habitat as RDS
loess_models <- IN %>%
dplyr::group_by(Sub_region, SiteID) %>%
dplyr::group_split() %>%
purrr::map(~ stats::loess(slope ~ jday, data = ., span = loess_span, degree = 1,
control = loess.control(surface = "direct")))
names(loess_models) <- IN %>% dplyr::distinct(SiteID) %>% dplyr::pull() %>% as.character()
saveRDS(loess_models,
file.path(out_dir, paste0("loess_models", tag_hab(current_hab), tag_ndz, ".rds")))
# ---- JDAY PLOTS (goes at end of rates_jday_tables chunk) ---------------
# Annual (jday) Δ per ha
p_jday_annual <- ggplot(annual_jday,
aes(SiteID, annual_ha, fill = Sub_region)) +
geom_hline(yintercept = 0, linewidth = 0.3) +
geom_col(width = 0.7) +
scale_fill_subregion() +
labs(title = title_ndz(paste0("Figure 6. Annual (jday) Δ", METRIC, " per ha — ", current_hab)),
x = "SiteID",
y = paste0("Δ", METRIC, " (\u03BCmol kg^-1) per ha")) +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
geom_text(aes(label = round(annual_ha, 2),
vjust = ifelse(annual_ha >= 0, -0.3, 1.2)),
size = 3)
p_jday_annual
# Area (ha) required to offset TARGET_YR, Jday method
nAnn_jday <- annual_jday %>%
mutate(
annual_m2 = annual_ha / 1e4,
HA_offOA = -TARGET_YR / annual_m2 / 1e4,
HA_offOA = if_else(is.finite(HA_offOA), HA_offOA, NA_real_)
)
# p_jday_area <- ggplot(nAnn_jday,
# aes(SiteID, HA_offOA, fill = Sub_region)) +
# geom_col(width = 0.7) +
# scale_fill_subregion() +
# ggbreak::scale_y_break(c(40, 200)) + # adjust break to taste
# geom_text(aes(label = scales::comma(HA_offOA, accuracy = 0.1)),
# vjust = 1.2, color = "white", size = 3, na.rm = TRUE) +
# labs(title = title_ndz(paste0("Area (ha) required to offset ", TARGET_YR,
# " \u03BCmol kg^-1 yr^-1 — Jday")),
# x = "SiteID", y = y_label_area) +
# theme_bw() +
# theme(axis.text.x = element_text(angle = 45, hjust = 1))
# p_jday_area
# Seasonal (jday) stacked bars + net annual labels
jday_season <- seasonal_jday %>%
rename(Season = Season_cal) %>%
mutate(Season = factor(Season, levels = c("Winter","Spring","Summer","Fall")))
# Compute per-site positive stack height and net value
jday_labels <- jday_season %>%
dplyr::group_by(SiteID) %>%
dplyr::summarise(
pos_top = sum(pmax(METRIC_season_ha, 0), na.rm = TRUE), # positive stack height
net = sum(METRIC_season_ha, na.rm = TRUE), # NET value
.groups = "drop"
) %>%
dplyr::mutate(
SiteID = factor(SiteID, levels = levels(jday_season$SiteID)),
y_off = 0.05 * max(abs(c(pos_top, net)), na.rm = TRUE),
y_lab = ifelse(pos_top > 0, pos_top + y_off, 0 + y_off), # always ABOVE bar
lab_val = net
)
p_jday_stacked <- ggplot(jday_season,
aes(SiteID, METRIC_season_ha, fill = Season)) +
geom_hline(yintercept = 0, linewidth = 0.3) +
geom_col(position = "stack") +
scale_fill_season() +
labs(
title = title_ndz(paste0("Figure 7. Jday: seasonal \u0394", METRIC, " per ha — ", current_hab)),
x = "SiteID", y = paste0("\u0394", METRIC, " per ha")
) +
geom_text(
data = jday_labels,
mapping = aes(x = SiteID, y = y_lab, label = round(lab_val, 2)),
inherit.aes = FALSE,
size = 3,
fontface = "bold"
) +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.12))) +
coord_cartesian(clip = "off") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p_jday_stacked
# Optional saves
if (save_figs) {
ggsave(file.path(fig_dir, "AnnualRates",
paste0("Fig6.jday_annual", tag_hab(current_hab), tag_ndz, ".png")),
p_jday_annual, width = 9, height = 5, dpi = 300)
# ggsave(file.path(fig_dir, "AnnualRates",
# paste0("jday_area", tag_hab(current_hab), tag_ndz, ".png")),
# p_jday_area, width = 9, height = 5, dpi = 300)
ggsave(file.path(fig_dir, "SeasonalRates",
paste0("Fig7.jday_seasonal_stacked", tag_hab(current_hab), tag_ndz, ".png")),
p_jday_stacked, width = 11, height = 6, dpi = 300)
}
message("Jday tables saved for ", current_hab, " (", METRIC, ").")4 Comparison Plots
## ---- legacy-plots (overlay + side-by-side), message=FALSE -------------
stopifnot(exists("IN"), exists("SSm"), exists("nAnn"),
exists("seasonal_jday"), exists("annual_jday"))
current_hab <- as.character(unique(IN$Habitat))[1]
month_letters <- function(x) substr(as.character(lubridate::month(x, label = TRUE, abbr = TRUE)), 1, 1)
# ---------- Prep seasonal-per-ha tables for both methods ----------
quarterly_season <- SSm %>%
dplyr::mutate(
Season = factor(Season, levels = c("Winter","Spring","Summer","Fall")),
METRIC_season_ha = season_total_ha
)
jday_season <- seasonal_jday %>%
dplyr::rename(Season = Season_cal) %>%
dplyr::mutate(Season = factor(Season, levels = c("Winter","Spring","Summer","Fall")))
# Align SiteID levels across both tables so bars line up perfectly
COMMON_X_QJ <- build_common_x_levels(
list(quarterly_season, jday_season),
x_col = "SiteID",
canonical = site_levels
)
quarterly_season <- apply_x_levels(quarterly_season, "SiteID", COMMON_X_QJ)
jday_season <- apply_x_levels(jday_season, "SiteID", COMMON_X_QJ)
# ---- Determine true stacked y-limits (no clipping of negatives) ----
stack_limits <- function(df, y = "METRIC_season_ha", x = "SiteID") {
df %>%
dplyr::group_by(.data[[x]]) %>%
dplyr::summarise(
ymin = sum(pmin(.data[[y]], 0), na.rm = TRUE),
ymax = sum(pmax(.data[[y]], 0), na.rm = TRUE),
.groups = "drop"
) %>%
summarise(
ymin = min(ymin, na.rm = TRUE),
ymax = max(ymax, na.rm = TRUE)
) %>%
as.numeric()
}
# compute limits across BOTH methods so the two panels share the same frame
lim_q <- stack_limits(quarterly_season, y = "METRIC_season_ha", x = "SiteID")
lim_j <- stack_limits(jday_season, y = "METRIC_season_ha", x = "SiteID")
lim_both <- c(min(lim_q[1], lim_j[1]), max(lim_q[2], lim_j[2]))
# pad a bit so labels/lines don’t clip
pad_stack <- function(lim, lower = 0.08, upper = 0.08) {
span <- diff(lim)
c(lim[1] - lower*span, lim[2] + upper*span)
}
y_range_stack <- pad_stack(lim_both, lower = 0.08, upper = 0.08)
# Optional: dynamic main titles
main_title_stack <- sprintf("Figure 9. Seasonal totals per ha — %s — %s (STACKED)", METRIC, current_hab)
main_title_group <- sprintf("Figure 10. Seasonal totals per ha — %s — %s (GROUPED)", METRIC, current_hab)
# ---------- SIDE-BY-SIDE: STACKED ----------
p_quarterly_stacked <- ggplot(quarterly_season,
aes(SiteID, METRIC_season_ha, fill = Season)) +
geom_hline(yintercept = 0, linewidth = 0.3) +
geom_col(position = "stack") +
scale_fill_season() +
scale_x_locked(COMMON_X_QJ) +
labs(title = "Quarterly method (Season × 91.3 days)",
x = "SiteID", y = bquote(Delta*.(METRIC)~" per ha")) +
theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
coord_cartesian(ylim = y_range_stack)
p_jday_stacked <- ggplot(jday_season,
aes(SiteID, METRIC_season_ha, fill = Season)) +
geom_hline(yintercept = 0, linewidth = 0.3) +
geom_col(position = "stack") +
scale_fill_season() +
scale_x_locked(COMMON_X_QJ) +
labs(title = "Jday method (daily to seasonal sum)",
x = "SiteID", y = bquote(Delta*.(METRIC)~" per ha")) +
theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
coord_cartesian(ylim = y_range_stack)
p_side_by_side_stacked <- (
(p_quarterly_stacked | p_jday_stacked) +
patchwork::plot_layout(guides = "collect") +
patchwork::plot_annotation(
title = main_title_stack,
theme = ggplot2::theme(
plot.title = ggplot2::element_text(size = 14, face = "bold", hjust = 0.5)
)
)
) & ggplot2::theme(legend.position = "right")
# ---------- SIDE-BY-SIDE: GROUPED (wider bars, tighter spacing) ----------
p_quarterly_grouped <- ggplot(quarterly_season,
aes(SiteID, METRIC_season_ha, fill = Season)) +
geom_hline(yintercept = 0, linewidth = 0.3) +
geom_col(position = position_dodge(width = 0.85), width = 0.8) +
scale_fill_season() +
scale_x_locked(COMMON_X_QJ, expand_add = 0.15) +
labs(title = "Quarterly method (Season × 91.3 days)",
x = "SiteID", y = bquote(Delta*.(METRIC)~" per ha")) +
theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_y_continuous(expand = expansion(mult = c(0.04, 0.06)))
p_jday_grouped <- ggplot(jday_season,
aes(SiteID, METRIC_season_ha, fill = Season)) +
geom_hline(yintercept = 0, linewidth = 0.3) +
geom_col(position = position_dodge(width = 0.85), width = 0.8) +
scale_fill_season() +
scale_x_locked(COMMON_X_QJ, expand_add = 0.15) +
labs(title = "Jday method (daily to season sum)",
x = "SiteID", y = bquote(Delta*.(METRIC)~" per ha")) +
theme_bw() + theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
scale_y_continuous(expand = expansion(mult = c(0.04, 0.06)))
p_side_by_side_grouped <- (
(p_quarterly_grouped | p_jday_grouped) +
patchwork::plot_layout(guides = "collect") +
patchwork::plot_annotation(
title = main_title_group,
theme = ggplot2::theme(
plot.title = ggplot2::element_text(size = 14, face = "bold", hjust = 0.5)
)
)
) & ggplot2::theme(legend.position = "right")
# ---------- OVERLAY: quarterly “flat” vs Jday LOESS + raw points ----------
SSm_perday <- SSm %>% dplyr::transmute(Sub_region, SiteID, Season, rate_per_day = mdd)
season_calendar <- tibble::tibble(jday = 1:365) %>%
dplyr::mutate(
Season = cut(jday, breaks = c(0, 90, 182, 273, 365),
labels = c("Winter","Spring","Summer","Fall"),
include_lowest = TRUE, right = TRUE)
)
quarterly_daily <- season_calendar %>%
dplyr::mutate(jdate = as.Date("2024-01-01") + (jday - 1)) %>%
dplyr::inner_join(SSm_perday, by = "Season", relationship = "many-to-many")
daily_preds_plot <- (IN %>%
dplyr::transmute(Sub_region, SiteID, jday = as.integer(jday), slope) %>%
dplyr::group_by(Sub_region, SiteID) %>%
dplyr::group_modify(~ fit_loess_and_predict_daily(.x, span = loess_span, window_days = window_days)) %>%
dplyr::ungroup()) %>%
dplyr::mutate(jdate = as.Date("2024-01-01") + (jday - 1),
rate_per_day = pred_per_day)
p_overlay <- ggplot() +
geom_hline(yintercept = 0, color = "black", linewidth = 0.3) +
geom_col(data = quarterly_daily,
aes(jdate, rate_per_day, fill = Season),
width = 1, alpha = 0.25, show.legend = TRUE) +
scale_fill_season() +
geom_line(data = daily_preds_plot,
aes(jdate, rate_per_day, group = SiteID),
linewidth = 0.7, color = "black") +
geom_point(data = IN,
aes(jdate, slope, color = Season), alpha = 0.5, size = 1.8) +
scale_color_season() +
scale_x_date(date_breaks = "1 month", labels = month_letters) +
labs(title = title_ndz(paste0("Figure 8. Quarterly (flat) vs Jday (LOESS) + raw points — ", current_hab)),
x = "Month", y = paste0("Slope (Δ", METRIC, " per day per km\u00B2)")) +
facet_wrap(~ SiteID, scales = "free_y") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.title = element_blank())
# ---------- print ----------
p_overlay
p_side_by_side_stacked
p_side_by_side_grouped
# ---------- optional saves ----------
if (save_figs) {
ggsave(file.path(fig_dir, "JdayPlots",
paste0("Fig8.overlay_quarterly-vs-jday", tag_hab(current_hab), tag_ndz, ".png")),
p_overlay, width = 11, height = 8, dpi = 300)
ggsave(file.path(fig_dir, "SeasonalRates",
paste0("Fig9.side-by-side_stacked_quarterly_vs_jday", tag_hab(current_hab), tag_ndz, ".png")),
p_side_by_side_stacked, width = 11, height = 6, dpi = 300)
ggsave(file.path(fig_dir, "SeasonalRates",
paste0("Fig10.side-by-side_grouped_quarterly_vs_jday", tag_hab(current_hab), tag_ndz, ".png")),
p_side_by_side_grouped, width = 11, height = 6, dpi = 300)
}# ---- session info ----
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Ventura 13.5.2
##
## Matrix products: default
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] fs_1.6.3 glue_1.8.0 rmarkdown_2.30 patchwork_1.2.0 ggbreak_0.1.6 scales_1.3.0 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2
## [12] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] here_1.0.2 lattice_0.22-6 rprojroot_2.1.1 digest_0.6.37 utf8_1.2.4 R6_2.6.1 evaluate_1.0.5 pillar_1.9.0
## [9] ggfun_0.2.0 yulab.utils_0.2.1 rlang_1.1.2 rstudioapi_0.17.1 jquerylib_0.1.4 Matrix_1.6-4 textshaping_0.4.0 labeling_0.4.3
## [17] splines_4.2.1 bit_4.0.5 munsell_0.5.1 compiler_4.2.1 xfun_0.53 systemfonts_1.1.0 pkgconfig_2.0.3 gridGraphics_0.5-1
## [25] mgcv_1.8-40 htmltools_0.5.8.1 tidyselect_1.2.1 fansi_1.0.6 crayon_1.5.3 tzdb_0.4.0 withr_3.0.2 rappdirs_0.3.3
## [33] grid_4.2.1 nlme_3.1-157 jsonlite_1.8.8 gtable_0.3.5 lifecycle_1.0.4 magrittr_2.0.3 cli_3.6.2 stringi_1.8.3
## [41] vroom_1.6.5 cachem_1.1.0 farver_2.1.2 bslib_0.8.0 ragg_1.3.2 generics_0.1.3 vctrs_0.6.5 tools_4.2.1
## [49] bit64_4.0.5 ggplotify_0.1.2 hms_1.1.3 rsconnect_1.5.1.9000 parallel_4.2.1 fastmap_1.2.0 yaml_2.3.10 timechange_0.3.0
## [57] colorspace_2.1-2 aplot_0.2.8 knitr_1.50 sass_0.4.10#rsconnect::deployDoc("path/to/your.Rmd")
# or for RPubs:
# rpubs::rpubsUpload("My title", "path/to/your.html")