`test_interaction`

tests for each significant GAM(M) whether any other
pressure modifies the IND response to the original pressure using a threshold
formulation of the GAM.

test_interaction( init_tbl, mod_tbl, interactions, sign_level = 0.05, k = 4, a = 0.2, b = 0.8, excl_outlier = FALSE )

init_tbl | The (full) output tibble of the |
---|---|

mod_tbl | A model output tibble from |

interactions | A tibble of all potential pressure combinations to test for interactions for specific INDs. |

sign_level | Significance level for selecting models on which to test for interactions. Only models with a p value <= sign_level will be selected; the default is 0.05. |

k | Choice of knots (for the smoothing function |

a | The lower quantile value of the selected threshold variable, which the estimated threshold is not allowed to exceed; the default is 0.2. |

b | The upper quantile value of the selected threshold variable, which the estimated threshold is not allowed to exceed; the default is 0.8. |

excl_outlier | logical; if TRUE, the outliers excluded in the original models will be also excluded in the interaction models. |

The function returns the input model tibble `mod_tbl` with the following 5 columns added:

`interaction`

logical; if TRUE, at least one thresh_gam performs better than its corresponding gam based on LOOCV value.

`thresh_var`

A list-column with the threshold variables of the better performing thresh_models.

`thresh_models`

A list-column with nested lists containing the better performing thresh_models.

`thresh_error`

A list-column capturing potential error messages that occurred as side effects when fitting each threshold GAMs and performing the LOOCV.

`tac_in_thresh`

logical vector; indicates for every listed thresh_model whether temporal autocorrelation (TAC) was detected in the residuals. TRUE if model residuals show TAC.

**Threshold-GAMs**

To identify potential interactions between pressures (relevant for sub-crit. 10.4),
threshold formulations are applied to every significant IND-pressure GAM.
Threshold-GAMs or TGAMs represent a special type of varying-coefficient models
and were first introduced by Ciannelli *et al.* (2004). In varying-coefficient
models coefficients are allowed to vary as a smooth functions of other variables
(Hastie and Tibshirani, 1993), which allows to detect interactions between two
external pressures. Threshold-GAMs are particularly useful if the response to the
interacting pressure variables is not continuous but rather step-wise. They have
been applied to data from real ecosystems to model population dynamics
(e.g. Otto *et al.*, 2014) up to food web dynamics in response to
climate, nutrient, and fishing interactions (e.g. Llope *et al.*, 2011).
The following threshold formulation is applied to every sign. IND-Pressure GAM:

$$ IND_t = \alpha_1 + s_1(press1_t) + \epsilon_t if press2 <= r$$ $$ IND_t = \alpha_2 + s_2(press1_t) + \epsilon_t if press2 > r$$

where the thin-plate regression spline function *s* can differ depending on
whether pressure 2 is below or above a threshold *r* with possible changes
in the intercept (from *a_1* to *a_2*). The index *t* represents
each observation in the training data. The threshold formulation can be
implemented in the GAM by including two smoothing functions for the same pressure
and using the `by`

argument in the smoothing function `s`

to select specific observations. The threshold itself is estimated from the data and
chosen by minimizing the GCV score (termed "gcvv" in the threshold-GAM object)
over an interval defined by the lower and upper quantiles (see the
`a`

and `b`

arguments respectively) of Pressure 2.

**Selection of threshold-GAMs**

To compare the performance of a GAM without interactions with the threshold-GAM we
implemented a leave-one-out cross-validation (LOOCV) approach as suggested by
Ciannelli *et al.* (2004). LOOCV involves splitting the full set of observations
n (i.e. the training data defined in `ind_init`

) into two parts: a new
training set containing n-1 observations and a test set containing 1 observations.
The threshold-GAM and corresponding GAM are then fit on the new training data
and a prediction is made for the excluded test observation using the corresponding
pressure 1 and pressure 2 values for that time step. A square prediction error
is then calculated for each predicted test observation. Repeating this approach
n times produces n squared errors, MSE_1, . . . , MSE_n. The LOOCV estimate for the
test MSE, also termed (genuine) cross-validatory squared prediction error, is the
root of the average of these n test error estimates. This approach properly accounts
for the estimation of the threshold value as well as the effective degrees of
freedom of all model terms.

**Implementation of threshold modeling**

For each IND~pressure pair, specific pressures to test for interactions can be selected
by creating a tibble containing the IND (termed `ind`), the pressure 1 (termed `press`)
and the pressure 2 (termed `t_var`). The easiest is to use the helper function
`select_interaction`

: it creates all combinations of IND~press pairs and
the threshold variables based on the input model tibble. If specific combinations should
not be modeled simply delete them from this data frame.

`test_interaction`

takes this data frame and the ind_init and model tibble as input
and applies the following procedure:

Filters significant GAMs and the corresponding IND ~ pressure | threshold pressure combination.

Extracts all data needed for the model, excluding outliers if set to TRUE.

Computes the LOOCV for each IND ~ pressure | threshold pressure model (threshold-GAM and GAM).

If the LOOCV estimate of the threshold-GAM is better than its corresponding GAM the threshold-GAM and its model output are saved in the returned model tibble. Note, however, that it is crucial to inspect the model diagnostic plots (i.e. the

`$thresh_plot`

from the`plot_diagnostics`

function! The performance of the final model (in terms of its Generalized Cross-Validation value) might be not much lower than for models with other thresholds indicating a lack of robustness. In this case, the interaction might need to be ignored but that needs to be decided on a case-by-case basis.

There is no threshold-GAMM implemented in this package yet. Until then, threshold-GAMs are
also applied to GAMMs when GAM residuals show temporal autocorrelation (TAC). However, the residuals
of threshold GAMs often show less TAC due to the splitting of observations into a low and high
regime of the threshold variable. In the case of significant TAC (indicated by the output variable
`tac_in_thresh`

) the user can decide whether that interaction should be neglected and
modify the `interaction`

) output variable accordingly.

Ciannelli, L., Chan, K.-S., Bailey, K.M., Stenseth, N.C. (2004) Nonadditive
effects of the environment on the survival of a large marine fish population.
*Ecology* 85, 3418-3427.

Hastie, T., Tibshirani, R. (1993) Varying-Coefficient Models. *Journal
of the Royal Statistical Society. Series B (Methodological)* 55, 757-796.

Llope, M., Daskalov, G.M., Rouyer, T.A., Mihneva, V., Chan, K.-S., Grishin, A.N.,
Stenseth, N.C. (2011) Overfishing of top predators eroded the resilience of the
Black Sea system regardless of the climate and anthropogenic conditions.
*Global Change Biology* 17, 1251-1265.

Otto, S.A., Kornilovs, G., Llope, M., Möllmann, C. (2014) Interactions among
density, climate, and food web effects determine long-term life cycle
dynamics of a key copepod. *Marine Ecology Progress Series* 498, 73-84.

`plot_diagnostics`

for assessing the model diagnostics

Other IND~pressure modeling functions:
`find_id()`

,
`ind_init()`

,
`model_gamm()`

,
`model_gam()`

,
`plot_diagnostics()`

,
`plot_model()`

,
`scoring()`

,
`select_model()`

# \donttest{ # Using some models of the Baltic Sea demo data in this package init_tbl <- ind_init_ex mod_tbl <- merge_models_ex[c(5:7),] interactions <- select_interaction(mod_tbl) test <- test_interaction(init_tbl, mod_tbl, interactions)#>#>#>#># if only one combination should be tested interactions <- select_interaction(mod_tbl)[1:2, ] test <- test_interaction(init_tbl, mod_tbl, interactions)#>#># Determine manually what to test for (e.g. only TZA ~ Fsprat | Pwin) interactions <- tibble::tibble(ind = "TZA", press = "Fsprat", t_var = "Pwin" ) test <- test_interaction(init_tbl, mod_tbl, interactions)#># }