Fit linear models to plankton time series data

Source: R/utils_analysis.R

Fit linear models to time series data for trend analysis, accounting for both long-term trends and seasonal cycles. This function is typically used before extracting coefficients or plotting trends with model predictions.

pr_model_data(dat)

Arguments

dat

A dataframe from pr_get_EOVs(), pr_get_Indices(), or similar functions containing time series data. Must contain only one parameter at a time.

Value

The input dataframe with model objects stored as a "Model" attribute

Details

The function fits linear models of the form:

Values ~ Year + sin(Month) + cos(Month)

Where:

  • Year term captures long-term linear trends (increasing or decreasing over time)

  • Harmonic terms capture seasonal cycles using Fourier basis functions (k=1, allowing for one complete seasonal cycle per year)

The model is fitted separately for each station (NRS/LTM) or bioregion (CPR). For LTM data, values are first averaged across depths to create a single depth- integrated estimate.

Model Interpretation

  • Positive Year coefficient: Variable increasing over time

  • Negative Year coefficient: Variable decreasing over time

  • Significant harmonic terms: Strong seasonal cycle present

The fitted model objects are stored as an attribute of the dataframe and can be extracted using pr_get_model() and analysed using pr_get_coeffs().

Limitations

  • Assumes linear trends (non-linear trends may be better captured with GAMs)

  • Assumes single seasonal cycle (some variables may have sub-annual variation)

  • Only handles one parameter at a time

See also

pr_get_model() to extract model objects, pr_get_coeffs() to extract tidy coefficients, pr_remove_outliers() to remove outliers before modelling

Examples

# Fit models to zooplankton biomass
dat <- pr_get_EOVs(Survey = "NRS") %>%
  dplyr::filter(Parameters == "Biomass_mgm3") %>%
  pr_remove_outliers(2) %>%
  pr_model_data()

# Extract and view model coefficients
models <- pr_get_model(dat)
coeffs <- pr_get_coeffs(models)
print(coeffs)
#> # A tibble: 40 × 7
#>    Station         term              estimate std.error statistic p.value signif
#>    <chr>           <chr>                <dbl>     <dbl>     <dbl>   <dbl> <chr> 
#>  1 Darwin          (Intercept)       -4.51e+3  2388.       -1.89  6.11e-2 ""    
#>  2 Darwin          Year_Local         2.27e+0     1.18      1.91  5.79e-2 ""    
#>  3 Darwin          pr_harmonic(Mont… -5.88e+0     6.53     -0.901 3.69e-1 ""    
#>  4 Darwin          pr_harmonic(Mont… -2.22e+0     7.20     -0.308 7.58e-1 ""    
#>  5 Esperance       (Intercept)        1.65e+2  1177.        0.140 8.91e-1 ""    
#>  6 Esperance       Year_Local        -7.79e-2     0.585    -0.133 8.97e-1 ""    
#>  7 Esperance       pr_harmonic(Mont… -1.31e+0     1.01     -1.29  2.24e-1 ""    
#>  8 Esperance       pr_harmonic(Mont…  2.87e+0     1.18      2.44  3.29e-2 "*"   
#>  9 Kangaroo Island (Intercept)       -1.24e+3   254.       -4.87  8.76e-6 "***" 
#> 10 Kangaroo Island Year_Local         6.17e-1     0.126     4.90  7.95e-6 "***" 
#> # ℹ 30 more rows

# Fit models to CPR phytoplankton diversity
dat_cpr <- pr_get_EOVs(Survey = "CPR") %>%
  dplyr::filter(Parameters == "ShannonPhytoDiversity") %>%
  pr_remove_outliers(2) %>%
  pr_model_data()