Statistical tests for evaluating trends are available in this
package. These methods are considered “secondary” analyses that use
results from a fitted GAM to evaluate trends or changes over time. In
particular, significance of changes over time are evaluated using
mixed-effect meta-analysis (Sera et al.
2019) applied to the GAM results to allow for full propagation of
uncertainty between methods. Each test includes a plotting method to
view the results.
Evaluating seasonal changes over time
The anlz_metseason(), anlz_mixmeta(), and
show_metseason() functions evaluate seasonal metrics (e.g.,
mean, max, etc.) between years, including an assessment of the trend for
selected years using mixed-effects meta-analysis modelling. These
functions require inputs for the seasonal ranges to evaluate
(doyend, doystr) and years for assessing the
trend in the seasonal averages/metrics (yrstr,
yrend).
The anlz_metseason() function estimates the seasonal
metrics (including uncertainty as standard error) for results from the
GAM fit. The seasonal metric can be any summary function available in R,
such as seasonal maxima (max), minima (min),
variance (var), or others. The function uses repeated
resampling of the GAM model coefficients to simulate multiple time
series as an estimate of uncertainty for the summary parameter.
The inputs for anlz_metseason() include the seasonal
range as day of year using start (doystr) and end
(doyend) days and the metfun and
nsim arguments to specify the summary function and number
of simulations, respectively. Here we show the estimate for the maximum
chlorophyll in each season, using a relatively low number of
simulations. Repeating this function will produce similar but slightly
different results because the estimates are stochastic. In practice, a
large value for nsim should be used to produce accurate
results (e.g., nsim = 1e5).
metseason <- anlz_metseason(mod, metfun = max, doystr = 90, doyend = 180, nsim = 100)
metseason
#> # A tibble: 29 × 7
#> yr met se bt_lwr bt_upr bt_met dispersion
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1991 0.264 0.486 0.230 18.5 2.06 0.0434
#> 2 1992 0.825 0.0856 5.09 11.0 7.49 0.0434
#> 3 1993 1.41 0.0929 18.9 43.6 28.7 0.0434
#> 4 1994 1.12 0.107 9.02 23.8 14.6 0.0434
#> 5 1995 1.27 0.0878 14.2 31.3 21.0 0.0434
#> 6 1996 1.31 0.108 13.9 37.1 22.7 0.0434
#> 7 1997 1.41 0.100 18.2 45.0 28.6 0.0434
#> 8 1998 1.96 0.112 61.1 168. 101. 0.0434
#> 9 1999 1.58 0.136 23.0 78.4 42.5 0.0434
#> 10 2000 1.32 0.125 13.3 40.8 23.3 0.0434
#> # ℹ 19 more rows
The anlz_mixmeta() function uses results from the
anlz_metseason() to estimate the trend in the seasonal
metric over a selected year range. Here, we evaluate the seasonal trend
from 2006 to 2017 for the seasonal estimate of the model results
above.
anlz_mixmeta(metseason, yrstr = 2006, yrend = 2017)
#> Call: mixmeta::mixmeta(formula = met ~ yr, S = S, data = totrnd, random = ~1 |
#> yr, method = "reml")
#>
#> Fixed-effects coefficients:
#> (Intercept) yr
#> 70.1911 -0.0343
#>
#> 12 units, 1 outcome, 12 observations, 2 fixed and 1 random-effects parameters
#> logLik AIC BIC
#> 8.4389 -10.8778 -9.9701
The show_metseason() function plots the seasonal metrics
and trends over time. The anlz_metseason() and
anlz_mixmeta() functions are used internally to get the
predictions. The same arguments for these functions are used for
show_metseason, with the mean as the default metric.
show_metseason(mod, doystr = 90, doyend = 180, yrstr = 2006, yrend = 2017, ylab = "Chlorophyll-a (ug/L)")
To plot only the seasonal metrics, the regression line showing trends
over time can be suppressed by setting one or both of yrstr
and yrend as NULL.
show_metseason(mod, doystr = 90, doyend = 180, yrstr = NULL, yrend = NULL, ylab = "Chlorophyll-a (ug/L)")
Adding an argument for metfun to
show_metseason() will plot results and trends for a metric
other than the average. Note the use of nsim in this
example. In practice, a much higher value should be used (e.g.,
nsim = 1e5)
show_metseason(mod, metfun = max, nsim = 100, doystr = 90, doyend = 180, yrstr = 2006, yrend = 2017, ylab = "Chlorophyll-a (ug/L)")
For convenience, the anlz_sumstats() function returns a
list of summary statistics for the GAM and associated mixed-effect
meta-analysis model. This function can be useful for creating tabular
results of the models. The list output includes mixmet as a
mixmeta object of the fitted mixed-effects meta-analysis
trend model, metseason as a tibble object of the fitted
seasonal metrics as returned by anlz_metseason() or
anlz_avgseason(), summary of the
mixmet object, and coeffs as a tibble object
of the slope estimate coefficients from mixmet. An
approximately linear slope estimate will be included as
slope.approx in coeffs if
trans = 'log10' for the GAM used in mod.
anlz_sumstats(mod, metfun = mean, doystr = 90, doyend = 180, yrstr = 2006, yrend = 2017)
#> $mixmet
#> Call: mixmeta::mixmeta(formula = met ~ yr, S = S, data = totrnd, random = ~1 |
#> yr, method = "reml")
#>
#> Fixed-effects coefficients:
#> (Intercept) yr
#> 46.4436 -0.0226
#>
#> 12 units, 1 outcome, 12 observations, 2 fixed and 1 random-effects parameters
#> logLik AIC BIC
#> 8.6873 -11.3745 -10.4668
#>
#>
#> $metseason
#> # A tibble: 29 × 7
#> yr met se bt_lwr bt_upr bt_met dispersion
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1991 0.115 0.527 0.135 15.8 1.46 0.0434
#> 2 1992 0.684 0.0655 4.03 7.28 5.41 0.0434
#> 3 1993 0.689 0.0976 3.53 8.51 5.48 0.0434
#> 4 1994 0.790 0.0864 4.68 10.2 6.91 0.0434
#> 5 1995 0.728 0.0798 4.18 8.59 5.99 0.0434
#> 6 1996 1.20 0.0914 11.7 26.8 17.7 0.0434
#> 7 1997 0.835 0.0906 5.10 11.6 7.68 0.0434
#> 8 1998 1.02 0.0795 8.28 17.0 11.9 0.0434
#> 9 1999 1.05 0.154 6.25 25.1 12.5 0.0434
#> 10 2000 0.811 0.0981 4.66 11.3 7.25 0.0434
#> # ℹ 19 more rows
#>
#> $summary
#> Call: mixmeta::mixmeta(formula = met ~ yr, S = S, data = totrnd, random = ~1 |
#> yr, method = "reml")
#>
#> Univariate extended random-effects meta-regression
#> Dimension: 1
#> Estimation method: REML
#>
#> Fixed-effects coefficients
#> Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
#> (Intercept) 46.4436 18.6863 2.4854 0.0129 9.8191 83.0681 *
#> yr -0.0226 0.0093 -2.4345 0.0149 -0.0408 -0.0044 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Random-effects (co)variance components
#> Formula: ~1 | yr
#> Structure: General positive-definite
#> Std. Dev
#> 0.0531
#>
#> Univariate Cochran Q-test for residual heterogeneity:
#> Q = 13.3462 (df = 10), p-value = 0.2050
#> I-square statistic = 25.1%
#>
#> 12 units, 1 outcome, 12 observations, 2 fixed and 1 random-effects parameters
#> logLik AIC BIC
#> 8.6873 -11.3745 -10.4668
#>
#>
#> $coeffs
#> # A tibble: 1 × 8
#> slope.approx slope slope.se z p likelihood ci.lb ci.ub
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 -0.538 -0.0226 0.00929 -2.43 0.0149 0.993 -0.0408 -0.00441
The seasonal estimates and mixed-effects meta-analysis regression can
be used to estimate the rate of seasonal change across the time series.
For any given year and seasonal metric, a trend can be estimated within
a specific window (i.e., yrstr and yrend
arguments in show_metseason()). This trend can be estimated
for every year in the period of record to estimate the rate of change
over time for the seasonal estimates.
The anlz_trndseason() function estimates the rate of
change and the show_trndseason() function plots the
results. For both, all inputs required for the
anlz_metseason() function are required, in addition to the
desired window width to evaluate for each year (win) and
the justification for the window as "left",
"right", or "center" from each year
(justify).
It’s important to note the behavior of the centering for window
widths (win argument) if choosing even or odd values. For
left and right windows, the exact number of years in win is
used. For example, a left-centered window for 1990 of ten years will
include exactly ten years from 1990, 1991, … , 1999. The same applies to
a right-centered window, e.g., for 1990 it would include 1981, 1982, …,
1990 (if those years have data). However, for a centered window, picking
an even number of years for the window width will create a slightly
off-centered window because it is impossible to center on an even number
of years. For example, if win = 8 and
justify = 'center', the estimate for 2000 will be centered
on 1997 to 2004 (three years left, four years right, eight years total).
Centering for window widths with an odd number of years will always
create a symmetrical window, i.e., if win = 7 and
justify = 'center', the estimate for 2000 will be centered
on 1997 and 2003 (three years left, three years right, seven years
total).
trndseason <- anlz_trndseason(mod, doystr = 90, doyend = 180, justify = 'left', win = 5)
head(trndseason)
#> # A tibble: 6 × 12
#> yr met se bt_lwr bt_upr bt_met dispersion yrcoef pval appr_yrcoef
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1991 0.115 0.529 0.134 15.9 1.46 0.0434 0.0275 0.399 0.362
#> 2 1992 0.684 0.0651 4.04 7.26 5.41 0.0434 0.103 0.0384 1.80
#> 3 1993 0.689 0.0971 3.53 8.49 5.48 0.0434 0.0701 0.263 1.29
#> 4 1994 0.790 0.0856 4.70 10.2 6.91 0.0434 0.0583 0.337 1.25
#> 5 1995 0.728 0.0805 4.17 8.62 5.99 0.0434 0.0490 0.463 1.18
#> 6 1996 1.20 0.0918 11.7 26.8 17.7 0.0434 -0.0614 0.210 -1.51
#> # ℹ 2 more variables: yrcoef_lwr <dbl>, yrcoef_upr <dbl>
The show_trndseason() function can be used to plot the
results directly, one model at a time.
show_trndseason(mod, doystr = 90, doyend = 180, justify = 'left', win = 5, ylab = 'Chl. change/yr, average')
As before, adding an argument for metfun to
show_trndseason() will plot results and trends for a metric
other than the average. Note the use of nsim in this
example. In practice, a much higher value should be used (e.g.,
nsim = 1e5)
show_trndseason(mod, metfun = max, nsim = 100, doystr = 90, doyend = 180, justify = 'left', win = 5, ylab = 'Chl. change/yr, maximum')
The results supplied by show_trndseason() can be
extended to multiple window widths by stacking the results into a single
plot. Below, results for window widths from 5 to 15 years are shown
using the show_sumtrndseason() function for a selected
seasonal range using a left-justified window. This function only works
with average seasonal metrics due to long processing times with other
metrics. To retrieve the results in tabular form, use
anlz_sumtrndseason().
show_sumtrndseason(mod, doystr = 90, doyend = 180, justify = 'left', win = 5:15)
Lastly, the plots returned by show_metseason() and
show_trndseason() can be combined using the
show_mettrndseason() function. This plot will show the
seasonal metrics from the GAM as in show_metseason() with
the colors of the points for the seasonal metrics colored by the
significance of the moving window trends shown in
show_trndseason(). The four colors indicate increasing,
decreasing, no trend, or no estimate (i.e., too few points for the
window). Most of the arguments for show_metseason() and
show_trndseason() apply to
show_mettrndseason().
show_mettrndseason(mod, metfun = mean, doystr = 90, doyend = 180, ylab = "Chlorophyll-a (ug/L)", win = 5, justify = 'left')
Four colors are used to define increasing, decreasing, no trend, or
no estimate. The cmbn argument can be used to combine the
no trend and no estimate colors into one color and label. Although this
may be desired for aesthetic reasons, the colors and labels may be
misleading with the default names since no trend is shown for points
where no estimates were made.
show_mettrndseason(mod, metfun = mean, doystr = 90, doyend = 180, ylab = "Chlorophyll-a (ug/L)", win = 5, justify = 'left', cmbn = T)
