California State University Northridge, genus level, 1992 to present
Olinger LK, Edmunds PJ, Levitan D, Smith TB, Lasker H, Feeley M, Mahoney L, Dahl A. Reproductive mode, but not rarity, influences population trajectories in corals. In Review 2026.
This page reruns the published CSUN analysis on the full monitoring record rather than stopping at the published window. The method, the code, and the site and taxon scoping are identical to the published page; only the year range differs. As new monitoring years are added to the bundled cover file, this page picks them up automatically and shows how the trajectories continue.
The coral-cover data for this analysis comes from the Reef Code benthic cover section: coral genera cover download.
The California State University Northridge program records coral cover at the genus level across 8 sites off St. John, and it holds the longest record in the analysis, starting in 1992. The published analysis covered 1992 through 2023. This page runs the same analysis over the full bundled record, from 1992 through the most recent year available. The full analysis runs below, from the coral-cover file to the figure (Figure 1) and the statistics table (Table 1). Every step is folded. Open any block to read the code and the plain-language notes.
cover_files <- sort(list.files("../../data/cover", pattern = "s2pt4.*\\.csv$", full.names = TRUE))
max_year <- max(readr::read_csv(tail(cover_files, 1), show_col_types = FALSE)$year)
pub_range <- c(1992, max_year); program <- "CSUN"; taxon_level <- "genera"; version <- "updated"suppressWarnings(suppressPackageStartupMessages({
library(tidyverse)
library(sandwich)
library(lmtest)
library(broom)
library(patchwork)
library(kableExtra)
}))
palette_grp <- c(
"Common.Broadcaster" = "#0072B2", "Rare.Broadcaster" = "#56B4E9",
"Common.Brooder" = "#D55E00", "Rare.Brooder" = "#E69F00"
)
label_grp <- c(
"Common.Broadcaster" = "Common broadcaster", "Rare.Broadcaster" = "Rare broadcaster",
"Common.Brooder" = "Common brooder", "Rare.Brooder" = "Rare brooder"
)The analysis begins from the genus-level coral-cover file and the site master table. The cover file holds one row per observation: a site, a year, a survey period, a coral genus, and its percent cover. The site master holds each site’s program, depth, and the year monitoring began. The cover file is matched by pattern in data/cover/, so this step reads whichever genus-level file is bundled, and the analysis grows as that file grows.
cover_files <- sort(list.files("../../data/cover", pattern = "s2pt4.*\\.csv$", full.names = TRUE))
cover_file <- tail(cover_files, 1)
benthiccover <- read.csv(cover_file) |>
select(year, program, site, period, coralGenera, perccover)
sitedat <- read.csv("../../data/site/00_RRS_siteMaster_allSites_data.csv")
c(cover_file = basename(cover_file), cover_rows = nrow(benthiccover), sites_listed = nrow(sitedat)) cover_file
"s2pt4_benthicCoverCoralGenera_49sites_1987_2024.csv"
cover_rows
"234551"
sites_listed
"50" The first move is the explicit publication-range clip: every downstream step runs only on years inside pub_range. At the published range this clip is harmless, since the bundled file does not extend beyond it; as the cover file grows in later monitoring years, the same clip is what keeps this page reproducing the published analysis exactly. After the clip, the analysis holds the rest of the sampling design constant so that trends reflect coral change rather than changes in which reefs or genera were surveyed: the CSUN sites and the coral genera. CSUN identifies corals to genus, so no species-level exclusions apply. One genus name is corrected. The shallow-site rule is applied in the next step, after the year grouping, to match the original analysis order.
# Publication-range clip: the explicit first step of scope.
benthiccover <- benthiccover |> filter(year >= pub_range[1], year <= pub_range[2])
# Sites: CSUN reefs, with a depth category.
sitedat <- sitedat |>
filter(program == "CSUN", yearadded < 2025) |>
mutate(depth_cat = ifelse(depth < 21, "Shallow", "Deep"))
# Cover: CSUN genera, with one genus-name correction.
benthiccover <- benthiccover |>
filter(program == "CSUN") |>
mutate(coralGenera = ifelse(coralGenera == "Isopyhyllastrea", "Isophyllastrea", coralGenera))
c(cover_rows = nrow(benthiccover), genera = n_distinct(benthiccover$coralGenera))cover_rows genera
5502 21 The analysis compares a baseline period to the full record. It bins the survey years into five-year groups and treats the first group as the baseline. CSUN uses its own binning rule for the long record, which drops a trailing group only when it would be shorter than half an interval.
bin_years_csun <- function(df, interval = 5) {
min_yr <- min(df$year, na.rm = TRUE); max_yr <- max(df$year, na.rm = TRUE)
breaks <- seq(min_yr, max_yr, by = interval)
if (length(breaks) > 1 && tail(breaks, 1) > max_yr - (interval / 2)) breaks <- head(breaks, -1)
labels <- paste(breaks, pmin(breaks + (interval - 1), max_yr), sep = "-")
map <- data.frame(year = seq(min_yr, max_yr)) |>
mutate(year_group = cut(year, breaks = c(breaks, max_yr + 1), labels = labels,
right = FALSE, include.lowest = TRUE))
df |> left_join(map, by = "year") |> mutate(year_group = as.factor(year_group))
}
benthiccover <- bin_years_csun(benthiccover, 5)
earliest_year_group <- levels(benthiccover$year_group)[1]
# Keep the shallow CSUN sites.
benthiccover <- benthiccover |>
filter(site %in% sitedat$site) |>
left_join(select(sitedat, site, depth_cat), by = "site") |>
filter(depth_cat == "Shallow")
# earliest_year_groupValidation note: the scoped dataset spans 1992 through 2023, holds 4,809 records across 7 CSUN sites and 21 coral genera. The baseline period is the first five-year group, 1992-1996.
Two presence rules keep the analysis on genera that were part of the community at baseline and at sites within their range. The first keeps genera recorded during the baseline period. The second drops genus-site pairs where the genus was never seen, since a site outside a genus range would enter the model as a structural zero rather than a decline.
genus_yeargroup_observed <- benthiccover |>
group_by(coralGenera, year_group) |>
summarize(seen = max(perccover, na.rm = TRUE) > 0, .groups = "drop") |>
filter(seen)
gen_in_baseline <- genus_yeargroup_observed |>
filter(year_group == earliest_year_group) |>
pull(coralGenera) |> unique()
benthiccover <- benthiccover |> filter(coralGenera %in% gen_in_baseline)
in_range <- benthiccover |>
group_by(coralGenera, site) |>
summarize(seen = max(perccover, na.rm = TRUE) > 0, .groups = "drop") |>
filter(seen) |> mutate(key = paste(coralGenera, site))
benthiccover <- benthiccover |>
filter(paste(coralGenera, site) %in% in_range$key)
benthiccover <- benthiccover |>
group_by(year_group, coralGenera, year, site) |>
summarise(perccover = mean(perccover, na.rm = TRUE), .groups = "drop")
c(rows = nrow(benthiccover), genera = n_distinct(benthiccover$coralGenera)) rows genera
3520 17 Each genus receives its mean percent cover during the baseline period. Ranking the genera by that mean sets the order used in the figure and sets up the commonness split. The cumulative proportion of cover, summed down the ranked list, measures how much of the reef’s baseline cover the top genera account for.
baseline <- benthiccover |>
filter(year_group == earliest_year_group) |>
group_by(coralGenera) |>
summarise(meancov = mean(perccover, na.rm = TRUE), .groups = "drop") |>
filter(meancov > 0) |>
arrange(desc(meancov)) |>
mutate(rank = row_number(),
normcov = meancov / sum(meancov, na.rm = TRUE),
cumsumcov = cumsum(normcov))
# head(baseline, 6)The commonness split follows one rule: the genera that together hold the top 90 percent of baseline cumulative cover are common, and the rest are rare. The manuscript identified 0.90 through a sliding-window analysis as the strictest rarity definition that leaves the temporal slopes stable.
baseline <- baseline |>
mutate(commonness = ifelse(cumsumcov <= 0.90, "Common", "Rare"))
table(baseline$commonness)
Common Rare
3 14 Each genus carries a reproductive mode, brooder or broadcaster. CSUN reports corals at the genus level, so the assignment reads a genus-level trait table from Edmunds, corrects one label, drops unknown entries, and joins the mode onto the baseline genera. Genera that lack a mode leave the analysis, because the model compares brooders to broadcasters and cannot place an unassigned genus.
repro <- read.csv("../../data/collab/reproMode_coralGenera_commonRare_Edmunds.csv") |>
rename(coralGenera = Genus, Reproductive_mode = LH.Strategy) |>
mutate(coralGenera = ifelse(coralGenera == "Montastrea cav", "Montastraea", coralGenera)) |>
filter(Reproductive_mode %in% c("Brooder", "Broadcaster")) |>
distinct(coralGenera, .keep_all = TRUE)
baseline <- baseline |>
left_join(repro |> select(coralGenera, Reproductive_mode), by = "coralGenera") |>
filter(Reproductive_mode %in% c("Brooder", "Broadcaster"))
# with(baseline, table(Commonness = commonness, `Reproductive mode` = Reproductive_mode))The model works on group means rather than individual genera, so that a single abundant genus cannot dominate the trend. This step joins the commonness and reproductive-mode labels onto the full cover series, averages cover within each group at each site and year, adds a year index that starts at zero, and takes the log10 of cover. Log10 cover is the response, because coral loss acts proportionally, and a proportional change is a straight line on the log scale.
analysis_data <- benthiccover |>
inner_join(baseline |> select(coralGenera, commonness, Reproductive_mode), by = "coralGenera") |>
group_by(category = commonness, Reproductive_mode, year, site) |>
summarise(perccover = mean(perccover, na.rm = TRUE), .groups = "drop") |>
mutate(yearind = year - min(year), perccoverLog = log10(perccover)) |>
filter(perccover > 0 & is.finite(perccoverLog)) |>
mutate(category = factor(category, levels = c("Common", "Rare")),
Reproductive_mode = factor(Reproductive_mode, levels = c("Broadcaster", "Brooder")))
c(model_rows = nrow(analysis_data), groups = nlevels(interaction(analysis_data$category, analysis_data$Reproductive_mode)))model_rows groups
854 4 The model predicts log10 percent cover from year, commonness, reproductive mode, and every interaction among them. Common and broadcaster are the reference levels, so the year term is the trend for common broadcasters, and the interactions measure how the other groups differ from it.
model_full <- lm(perccoverLog ~ yearind * category * Reproductive_mode, data = analysis_data)The residual variance is larger for the abundant groups than for the sparse groups, so ordinary standard errors would misstate the uncertainty. The analysis replaces them with HC3 heteroscedasticity-consistent standard errors, computed from a sandwich covariance matrix. Every standard error, p-value, and confidence interval below uses this robust covariance.
vcov_hc3 <- sandwich::vcovHC(model_full, type = "HC3")
coef_robust <- broom::tidy(lmtest::coeftest(model_full, vcov = vcov_hc3))The four group trends are linear combinations of the model coefficients. The common-broadcaster slope is the year term. Each other group adds the interactions that apply to it. The standard error of each slope comes from the robust covariance, so the uncertainty carries through. Annual percent change back-transforms the log10 slope into a yearly rate.
group_slope <- function(terms) {
L <- setNames(rep(0, length(coef(model_full))), names(coef(model_full)))
L[terms] <- 1
est <- sum(L * coef(model_full))
se <- sqrt(as.numeric(t(L) %*% vcov_hc3 %*% L))
tibble(Slope = est, SE = se, `Annual % change` = 100 * (10^est - 1),
`P-value` = 2 * pnorm(-abs(est / se)))
}
group_slopes <- bind_rows(
"Common broadcaster" = group_slope("yearind"),
"Rare broadcaster" = group_slope(c("yearind", "yearind:categoryRare")),
"Common brooder" = group_slope(c("yearind", "yearind:Reproductive_modeBrooder")),
"Rare brooder" = group_slope(c("yearind", "yearind:categoryRare",
"yearind:Reproductive_modeBrooder",
"yearind:categoryRare:Reproductive_modeBrooder")),
.id = "Group"
)The figure pairs the baseline abundance distribution with the temporal trends. The left panel ranks genera by baseline cover and colors each by group. The right panel shows site-level group mean cover over time on a log scale, with a linear fit per group.
sad_df <- baseline |>
filter(meancov > 0) |>
mutate(coralGenera = factor(coralGenera, levels = rev(baseline$coralGenera[order(baseline$rank)])),
grp = interaction(commonness, Reproductive_mode))
p_sad <- ggplot(sad_df, aes(coralGenera, meancov, fill = grp)) +
geom_col(width = 0.8) +
scale_fill_manual(values = palette_grp, labels = label_grp, name = NULL, drop = FALSE) +
coord_flip() +
labs(x = NULL, y = "Mean baseline percent cover", title = "Baseline abundance") +
theme_classic(base_size = 11) +
theme(legend.position = "top", axis.text.y = element_text(face = "italic", size = 8))
trend_df <- analysis_data |> mutate(grp = interaction(category, Reproductive_mode))
p_trend <- ggplot(trend_df, aes(year, perccover, color = grp, fill = grp)) +
geom_jitter(shape = 21, color = "black", alpha = 0.2, width = 0.25, height = 0, show.legend = FALSE) +
geom_smooth(method = "lm", formula = y ~ x, se = TRUE, alpha = 0.15, linewidth = 0.8) +
scale_color_manual(values = palette_grp, labels = label_grp, name = NULL, drop = FALSE) +
scale_fill_manual(values = palette_grp, labels = label_grp, name = NULL, drop = FALSE) +
scale_y_log10() + annotation_logticks(sides = "l", color = "grey70") +
labs(x = "Year", y = "Percent cover (log scale)", title = "Temporal trends") +
theme_classic(base_size = 11) +
theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))
p_sad + p_trend
The table reports the model coefficients with HC3 robust standard errors and p-values. Bold marks terms significant at p less than 0.05.
term_labels <- c(
"(Intercept)" = "Intercept", "yearind" = "Year",
"categoryRare" = "Rarity", "Reproductive_modeBrooder" = "Repro. mode (brooder)",
"yearind:categoryRare" = "Year x Rarity", "yearind:Reproductive_modeBrooder" = "Year x Repro. mode",
"categoryRare:Reproductive_modeBrooder" = "Rarity x Repro. mode",
"yearind:categoryRare:Reproductive_modeBrooder" = "Year x Rarity x Repro. mode"
)
coef_tbl <- coef_robust |>
transmute(Predictor = term_labels[term],
`Estimate (SE)` = sprintf("%.3f (%.3f)", estimate, std.error),
`P-value` = ifelse(p.value < 0.001, "< 0.001", sprintf("%.3f", p.value)),
.sig = p.value < 0.05)
kbl(coef_tbl |> select(-.sig), align = c("l", "r", "r"),
caption = "Table 1 (CSUN). Linear-model coefficients with HC3 robust standard errors.") |>
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
row_spec(which(coef_tbl$.sig), bold = TRUE)| Predictor | Estimate (SE) | P-value |
|---|---|---|
| Intercept | 0.350 (0.081) | < 0.001 |
| Year | -0.011 (0.004) | 0.006 |
| Rarity | -1.431 (0.118) | < 0.001 |
| Repro. mode (brooder) | -0.611 (0.097) | < 0.001 |
| Year x Rarity | -0.001 (0.006) | 0.919 |
| Year x Repro. mode | 0.019 (0.005) | < 0.001 |
| Rarity x Repro. mode | 0.463 (0.154) | 0.003 |
| Year x Rarity x Repro. mode | -0.016 (0.008) | 0.049 |
group_slopes |>
mutate(Slope = sprintf("%.4f", Slope), SE = sprintf("%.4f", SE),
`Annual % change` = sprintf("%.2f", `Annual % change`),
`P-value` = ifelse(`P-value` < 0.001, "< 0.001", sprintf("%.3f", `P-value`))) |>
kbl(align = c("l", "r", "r", "r", "r"),
caption = "Group-specific temporal slopes estimated by Wald tests on linear combinations of model coefficients. The annual percent change is calculated as (10^slope - 1) x 100. P-values test the null hypothesis that the total temporal slope equals zero (H0: slope = 0). Common Broadcasters are the reference group; their slope and p-value come directly from the Year coefficient.") |>
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| Group | Slope | SE | Annual % change | P-value |
|---|---|---|---|---|
| Common broadcaster | -0.0115 | 0.0041 | -2.61 | 0.006 |
| Rare broadcaster | -0.0121 | 0.0045 | -2.75 | 0.008 |
| Common brooder | 0.0080 | 0.0028 | 1.86 | 0.004 |
| Rare brooder | -0.0083 | 0.0042 | -1.89 | 0.050 |
The coefficient table reports the fitted trends. This table steps back to the genera themselves and reports, for every genus in the baseline community, how abundant and how widespread it was at the start of the record and how that changed by the end. Three complementary measures separate the different ways a genus can decline. Cover is the mean percent cover across the genus range. Occupancy is the proportion of surveyed sites where the genus was present, so it tracks range contraction independent of local abundance. Extent is the mean cover at the sites where the genus was actually present, so it tracks thinning within the occupied range. A genus can hold its occupancy while its extent falls, or contract its range while staying dense where it persists; reporting both separates the two. Start values average over the baseline five-year group, and end values average over the last three years of the record. Occupancy counts a site only from the year it entered monitoring, so the denominator follows the site-by-site survey history rather than assuming every reef was watched from the first year.
source("../../_includes/_analysis_helpers.R")
species_table <- compute_change_metrics(benthiccover, baseline, sitedat, earliest_year_group, "coralGenera")
.end_yrs <- (max(benthiccover$year) - 2):max(benthiccover$year)
.desc_tbl <- species_table |>
arrange(rank_start) |>
transmute(
Genus = taxon,
`Repro. mode` = mode,
Commonness = commonness,
`Start rank` = rank_start,
`End rank` = ifelse(is.na(rank_end), "-", as.character(rank_end)),
`Start cover (%)` = sprintf("%.3f", cover_start),
`End cover (%)` = ifelse(is.na(cover_end), "-", sprintf("%.3f", cover_end)),
`Start occ. (%)` = sprintf("%.1f", occ_start),
`End occ. (%)` = sprintf("%.1f", occ_end),
`Start extent (%)` = ifelse(is.na(ext_start), "-", sprintf("%.2f", ext_start)),
`End extent (%)` = ifelse(is.na(ext_end), "-", sprintf("%.2f", ext_end)),
.common = commonness == "Common"
)
kbl(.desc_tbl |> select(-.common), align = c("l", "l", "l", rep("r", 8)),
caption = sprintf(
"Descriptive statistics for coral genera in the CSUN monitoring program. Genera are ordered by their rank in the baseline period (%s). Reproductive Mode ('Repro. Mode') indicates reproductive strategy. Rows in bold signify common genera, determined by the top 90%% cumulative cover threshold during the baseline period. Start Cover (%%), Start Occupancy (Start Occup. (%%)), Start Extent (%%), and Start Rank are metrics from the baseline period. End Cover (%%), End Occupancy (End Occup. (%%)), End Extent (%%), and End Rank are metrics for the recent period (%d-%d).",
earliest_year_group, min(.end_yrs), max(.end_yrs))) |>
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
column_spec(1, italic = TRUE) |>
row_spec(which(.desc_tbl$.common), bold = TRUE)| Genus | Repro. mode | Commonness | Start rank | End rank | Start cover (%) | End cover (%) | Start occ. (%) | End occ. (%) | Start extent (%) | End extent (%) | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Orbicella | Orbicella | Broadcaster | Common | 1 | 1 | 10.203 | 2.347 | 82.9 | 81.0 | 12.31 | 2.90 |
| Siderastrea | Siderastrea | Broadcaster | Common | 2 | 3 | 1.243 | 0.511 | 97.1 | 100.0 | 1.28 | 0.51 |
| Porites | Porites | Brooder | Common | 3 | 2 | 0.693 | 1.298 | 97.1 | 100.0 | 0.71 | 1.30 |
| Agaricia | Agaricia | Brooder | Rare | 4 | 5 | 0.610 | 0.182 | 82.9 | 90.5 | 0.74 | 0.20 |
| Montastraea | Montastraea | Broadcaster | Rare | 5 | 4 | 0.324 | 0.183 | 68.6 | 76.2 | 0.47 | 0.24 |
| Diploria | Diploria | Broadcaster | Rare | 6 | 6 | 0.297 | 0.042 | 65.7 | 52.4 | 0.45 | 0.08 |
| Colpophyllia | Colpophyllia | Broadcaster | Rare | 7 | 8 | 0.287 | 0.037 | 25.7 | 38.1 | 1.12 | 0.10 |
| Meandrina | Meandrina | Broadcaster | Rare | 8 | - | 0.068 | 0.000 | 31.4 | 0.0 | 0.21 | - |
| Dichocoenia | Dichocoenia | Broadcaster | Rare | 9 | 13 | 0.047 | 0.002 | 20.0 | 14.3 | 0.24 | 0.02 |
| Stephanocoenia | Stephanocoenia | Broadcaster | Rare | 10 | 7 | 0.037 | 0.041 | 34.3 | 61.9 | 0.11 | 0.07 |
| Acropora | Acropora | Broadcaster | Rare | 11 | 10 | 0.029 | 0.010 | 11.4 | 4.8 | 0.18 | 0.21 |
| Dendrogyra | Dendrogyra | Broadcaster | Rare | 12 | - | 0.025 | 0.000 | 14.3 | 0.0 | 0.15 | - |
| Eusmilia | Eusmilia | Broadcaster | Rare | 13 | 12 | 0.009 | 0.005 | 14.3 | 23.8 | 0.06 | 0.02 |
| Madracis | Madracis | Brooder | Rare | 14 | 9 | 0.003 | 0.019 | 8.6 | 28.6 | 0.03 | 0.07 |
| Mussa | Mussa | Brooder | Rare | 15 | - | 0.002 | 0.000 | 2.9 | 0.0 | 0.06 | - |
| Favia | Favia | Brooder | Rare | 16 | - | 0.002 | 0.000 | 8.6 | 0.0 | 0.02 | - |
| Manicina | Manicina | Brooder | Rare | 17 | 11 | 0.002 | 0.008 | 5.7 | 19.0 | 0.03 | 0.04 |
The analysis objects are written to a timestamped RData file under data/rdata/, tagged with the program and the version (as published or updated) so the two ranges never overwrite each other. The manuscript-values page loads the most recent file for each program and version and reads every derived number from these saved objects, rather than repeating the computation. The derived CSVs written elsewhere are left unchanged.
The trend model runs only on genera that were part of the baseline community, so that a modeled slope reflects a genus present at the outset rather than one that arrived mid-record. Some genera were observed at least once during the record but not during the baseline period, or could not be assigned a reproductive mode; those are held out of the trend analysis and listed here with the year each was first observed.
.excl <- excluded_taxa_table(genus_yeargroup_observed, baseline$coralGenera, "coralGenera")
kbl(.excl, align = c("l", "r"),
caption = "Taxa excluded from trend analyses because they were not recorded during the baseline period (2001–2005 for TCRMP and VINPS; 1992–1996 for CSUN). For each taxon, the program and year of first observation at shallow monitoring sites are shown. CSUN taxa are at the genus level; TCRMP and VINPS taxa are at the species level.") |>
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) |>
column_spec(1, italic = TRUE)| Taxon | Year first observed |
|---|---|
| Helioseris | 2002 |
| Isophyllia | 1997 |
| Scolymia | 1997 |
| Solenastrea | 2007 |
The coefficient table uses HC3 robust standard errors because the residual variance is not constant across groups. These diagnostic panels show why. Residuals versus fitted values reveal the unequal spread, the normal Q-Q plot shows the tails, and the residual distribution summarizes the whole. The robust standard errors used throughout the analysis are the response to exactly this structure.
model_diagnostics_plot(model_full)