Modeling of indicator responses to single pressures with GAMMs

Source: R/model_gamm.R

model_gamm accounts for temporal autocorrelation (TAC) in the time series by fitting Generalized Additive Mixed Models (GAMMs) that include AR or ARMA correlation structures using the gamm function. The GAMMs are applied to all IND~pressure combinations provided as input or only those with significant TAC in the GAM residuals (using filter argument).

model_gamm(
  init_tbl,
  k = 5,
  family = stats::gaussian(),
  excl_outlier = NULL,
  filter = NULL,
  lme_control = NULL
)

Arguments

init_tbl

The output tibble of the ind_init function.
k

Choice of knots (for the smoothing function s); the default is 5.
family

A description of the error distribution and link to be used in the GAM. This needs to be defined as a family function (see also family). All standard family functions can be used as well some of the distribution families in the mgcv package (see family.mgcv; e.g.negbin). Note that nb, which estimates theta parameter, cannot be used for gamm.
excl_outlier

A list of values identified as outliers in specific IND~pressure GAMMs, which should be excluded in this modeling step (the output tibble of this function includes the variable `pres_outlier`, which is a column-list containing all indices of values with cook`s distance > 1 (see below). The function can be re-run again, then excluding all these outliers provided in $pres_outlier from the the first run (see example)).
filter

logical; a filter used to select specific rows in init_tbl (row gets selected if value TRUE). That could be the tac column in the model_gam output tibble which indicates whether the model residuals show TAC.
lme_control

A specification list of the control values for the lme fit. The default is NULL so the control settings in gamm are used (list(niterEM=0,optimMethod="L-BFGS-B",returnObject=TRUE)). To overwrite these values use best the lmeControl function (see example).

Value

Returns a model output tibble that contains for each filtered IND~pressure pair 6 rows with the individual GAMM outputs. The structure remains the same as in model_gam except for the explained deviance, which is not computed by the gamm function. The selection of the final correlation structure for each IND~pressure model can be done manually on this tibble or with an automatized routine using select_model.

Details

Modeling first-differenced indicator time series can be an alternative solution to avoid temporal dependence between observations. However, this approach does often not help reducing the significant auto-correlation while GAMMs do as found in Otto et al. (2018). Such an extension implies that the single elements of the response variable are not independent anymore and that the correlation between the residuals at time t1 and t2 only depends on their time difference t1 – t2 (Wood, 2006).

In model_gamm six GAMMs are computed for each filtered IND~pressure pair, i.e.

  1. no correlation structure (for AIC comparison)

  2. auto-regressive error structure of order p=1 (AR1)

  3. auto-regressive error structure of order p=2 (AR2)

  4. auto-regressive moving average of order p=1 and q=1 (ARMA11)

  5. auto-regressive moving average of order p=1 and q=2 (ARMA12)

  6. auto-regressive moving average of order p=2 and q=1 (ARMA21)

References

Otto, S.A., Kadin, M., Casini, M., Torres, M.A., Blenckner, T. (2018) A quantitative framework for selecting and validating food web indicators. Ecological Indicators, 84: 619-631, doi: https://doi.org/10.1016/j.ecolind.2017.05.045

Wood, S.N. (2006) Generalized Additive Models: An Introduction with R. Chapman and Hall/CRC Press

See also

gamm for more information on GAMMs and plot_diagnostics for assessing the model diagnostics

Other IND~pressure modeling functions: find_id(), ind_init(), model_gam(), plot_diagnostics(), plot_model(), scoring(), select_model(), test_interaction()

Examples

# Using the Baltic Sea demo data in this package
dat_init <- ind_init(
  ind_tbl = data.frame(Cod = ind_ex$Sprat),
  press_tbl = press_ex[, c("Fsprat", "Fher")],
  time = ind_ex[ ,1])
gam_tbl <- model_gam(dat_init)
# Any temporal autocorrelation
gam_tbl$tac

#> [1] TRUE TRUE
# Applying model_gamm function and passing the $tac variable as filter
gamm_tbl <- model_gamm(dat_init, filter = gam_tbl$tac)
# Use different control values for lme fit:
 gamm_tbl <- model_gamm(
 dat_init, filter = gam_tbl$tac,
 lme_control = nlme::lmeControl(niterEM = 100, msMaxIter = 200)
 )
# \donttest{
 # run analysis will full demo dataset
 gam_tbl <- model_gam(ind_init_ex)
 gamm_tbl <- model_gamm(ind_init_ex, filter = gam_tbl$tac)
 # now run the same analysis but exclude all outliers in the GAMMs
 gamm_tbl2 <- model_gamm(ind_init_ex, filter = gam_tbl$tac,
   excl_outlier = gamm_tbl$pres_outlier)
# }