Common errors & warnings encountered in DImulti()"

In DImodelsMulti: Fit Multivariate Diversity-Interactions Models with Repeated Measures

```{css styling, echo=FALSE} span.R { font-family: Courier New; }

span.bad { color: red; }

```r
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

options(crayon.enabled = TRUE)

ansi_aware_handler <- function(x, options)
{
  paste0(
    "<pre class=\"r-output\"><code>",
    fansi::sgr_to_html(x = x, warn = FALSE, term.cap = "256"),
    "</code></pre>"
  )
}

old_hooks <- fansi::set_knit_hooks(knitr::knit_hooks,
                                 which = c("output", "message", "error", "warning"))
knitr::knit_hooks$set(
  output = ansi_aware_handler,
  message = ansi_aware_handler,
  warning = ansi_aware_handler,
  error = ansi_aware_handler
)

library(DImodelsMulti)
data("dataBEL"); data("dataSWE"); data("simMV"); data("simRM"); data("simMVRM")

Incorrect argument specifications

The DImulti() function has a long list of possible parameters, which can make it difficult to ensure that all arguments supplied are in the requested format. See below for an example of each parameter and the included datasets, dataBEL, dataSWE, simMV, simRM and simMVRM, for examples of dataset formatting.

y

This is a required argument. It can be passed in three forms, using either column names or indices:

• Univariate

In this case, we pass a single name or index of a numerical column containing the ecosystem function response value. See ?dataSWE for an example of a univariate dataset suitable for use with DImulti().

DImulti(y = "YIELD", ..., data = dataSWE)
DImulti(y = 9, ..., data = dataSWE)

dataSWE[1:6, "YIELD", drop = FALSE]

• Multivariate stacked

In this case, we pass a single name or index of a numerical column containing the ecosystem function response value. See ?dataBEL for an example of a multivariate stacked dataset suitable for use with DImulti().

DImulti(y = "Y", ..., data = dataBEL)
DImulti(y = 9, ..., data = dataBEL)

dataBEL[1:6, "Y", drop = FALSE]

• Multivariate wide

In this case, we pass multiple names or indices of numerical columns containing the ecosystem function response values. See ?simMV for an example of a multivariate wide dataset suitable for use with DImulti().

DImulti(y = c("Y1", "Y2", "Y3", "Y4"), ..., data = simMV)
DImulti(y = 9:12, ..., data = simMV)

simMV[1:6, c("Y1", "Y2", "Y3", "Y4"), drop = FALSE]

eco_func

This is a required argument for multivariate data and excluded for univariate data. It can be passed in two forms:

• Multivariate stacked

In this case, we pass a vector of size two. The first index contains a single column name for the ecosystem function response type, a factor, the second contains the autocorrelation structure desired, either "UN" or "CS". See ?dataBEL for an example of a multivariate stacked dataset suitable for use with DImulti().

DImulti(eco_func = c("Var", "UN"), ..., data = dataBEL)

dataBEL[1:6, "Var", drop = FALSE]

• Multivariate wide

In this case, we pass a vector of size two. The first index contains the string "NA", the second contains the autocorrelation structure desired, either "UN" or "CS". See ?simMV for an example of a multivariate stacked dataset suitable for use with DImulti().

DImulti(eco_func = c("NA", "UN"), ..., data = simMV)

time

This is a required argument for repeated measures data and excluded for data taken at a single time point. It can be passed in the following form:

• Repeated Measures

In this case, we pass a vector of size two. The first index contains a single column name for the time point reference, a factor, the second contains the autocorrelation structure desired, either "UN", "CS", or "AR1". See ?dataSWE for an example of a repeated measures dataset suitable for use with DImulti().

DImulti(time = c("YEARN", "AR1"), ..., data = dataSWE)

dataSWE[c(1, 49, 97, 2, 50, 98), "YEARN", drop = FALSE]

unit_IDs

This is a required argument for any data.

It can be passed a column name or index, referencing a unique reference for each experimental unit, used for grouping responses from different ecosystem functions or time points.

DImulti(unit_IDs = "plot", ..., data = simMVRM)
DImulti(unit_IDs = 1, ..., data = simMVRM)

simMVRM[1:6, "plot", drop = FALSE]

prop

This is a required argument for any data.

It can be passed as a vector of column names or indices, referencing the (numerical) initial proportions of any species included in the experimental design. In the event that a species is included as a factor as opposed to a numerical value, e.g. if it only appears at values 0 or 1, only the numerical species should be included here, while the factor species is passed through the extra_fixed parameter. A warning will be printed to the console warning the user of not meeting the simplex requirement.

DImulti(unit_IDs = paste0("p", 1:4), ..., data = simMVRM)
DImulti(unit_IDs = 2:5, ..., data = simMVRM)

simMVRM[1:6, 2:5, drop = FALSE]

data

This is a required argument for any data.

This argument should be a dataframe or tibble which contains all columns referenced in the DImulti() call. There are some restrictions on columns names allowed in this dataset. An error will be printed to the console to inform a user of name changes required.

DImulti(..., data = simMVRM)
DImulti(..., data = simMVRM)

simMVRM[1:6, , drop = FALSE]

DImodel

This is a required argument for any data.

This argument should be a single string selected from the list of included interaction structures: "STR", "ID", "FULL", "E", "AV", "ADD", "FG".

DImulti(DImodel = "FULL", ...)
DImulti(DImodel = "AV", ...)

FG

This is a required argument if the argument "FG" is passed through DImodel and can be excluded otherwise.

This argument should be a vector of strings, of the same length as prop, which maps each species in the experiment to a functional grouping.

DImulti(prop = c("G1", "G2", "L1", "L2"), DImodel = "FG", 
        FG = c("Grass", "Grass", "Legume", "Legume"), ..., data = dataBEL)

ID

This is an optional argument.

This argument should be a vector of strings, of the same length as prop, which maps each species in the experiment to a identity grouping, assuming functional redundancy between within-group species for their ID effects, this grouping does not affect interactions.

DImulti(prop = c("G1", "G2", "L1", "L2"), ID = c("Grass", "Grass", "Legume", "Legume"), 
        ..., data = dataBEL)

extra_fixed

This is an optional argument.

This parameter can be passed a formula of any additional fixed effects to be included, e.g. treatments. These effects can be easily crossed with the automatic DI model formula by prefacing the formula passed with 1: or 1*. See ?dataSWE for a dataset containing treatments suitable for use with DImulti().

DImulti(extra_fixed = ~DENS, ..., data = dataSWE)
DImulti(extra_fixed = ~DENS + TREAT, ..., data = dataSWE)
DImulti(extra_fixed = ~1*DENS, ..., data = dataSWE)
DImulti(extra_fixed = ~1:(DENS+TREAT), ..., data = dataSWE)
DImulti(extra_fixed = ~1:DENS:TREAT, ..., data = dataSWE)

dataSWE[1:6, c("DENS", "TREAT"), drop = FALSE]

estimate_theta

This is an optional argument.

This parameter can be passed a boolean value (FALSE), indicating whether the user wants the function to use profile likelihood to estimate values for the non-linear parameter $\theta$. Defaults to "STR" or "ID" was passed through the parameter DImodel, then setting estimate_theta = TRUE will cause a warning to be printed to the console, as $\theta$ only applies to interaction terms, which do not exist for those options.

DImulti(estimate_theta = TRUE, ...)
DImulti(estimate_theta = FALSE, ...)

theta

This is an optional argument.

This argument should contain numerical values representing the non-linear parameter $\theta$, which will be applied as a power to the products of pairs of species in the interaction terms of the model. A single value may be passed, which will be applied to each ecosystem function, or a different value may be passed for each ecosystem function. If "STR" or "ID" was passed through the parameter DImodel, then setting a non-1 value for theta will cause a warning to be printed to the console, as $\theta$ only applies to interaction terms, which do not exist for those options.

DImulti(theta = 1, ...)
DImulti(theta = c(0.5, 1, 1.2), ...)

method

This is an optional argument.

This parameter is used to change the estimation method used by DImulti(). The options are "REML", referring to restricted maximum likelihood, and "ML", referring to maximum likelihood. Defaults to "REML".

DImulti(method = "REML", ...)
DImulti(method = "ML", ...)

Singular fit

A singular fit error, thrown by nlme::gls() when fitting the model, occurs when an element with value of exactly zero exists in the fixed effect variance covariance matrix of the model. This usually caused by rank deficiency in the model, i.e. two terms explaining the same variance. The fix is usually to simplify/change the fixed effect terms, be it the interactions structure, treatment terms, or theta values. To aid in this change, the fixed effect formula and any estimated theta values are printed to the console with the error.

As an example, we use the dataset dataSWE</span, which is included in this package.

XXDO NOT RUNXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

DImulti(prop = 5:8, y = "YIELD", time = c("YEARN", "CS"), unit_IDs = "PLOT", data = dataSWE, 
        DImodel = "ID", extra_fixed = ~DENS+TREAT)

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX As we include no theta values or interaction terms, and use the simple structure "CS", compound symmetry, for our repeated measure, the only complication in this model is the inclusion of the two treatment terms, DENS and TREAT. When multiple treatments are included, it is usually good to check that one cannot be inferred from the other, e.g. if DENS = High then we know that TREAT = 1. Luckily this is not the case here.

dataSWE[c(1, 16, 31, 40), c("DENS", "TREAT"), drop = FALSE]

Instead, we change the formatting of our formula, from creating intercepts to crossing the treatment with our ID effects, but not each other.

DImulti(prop = 5:8, y = "YIELD", time = c("YEARN", "CS"), unit_IDs = "PLOT", data = dataSWE, 
        DImodel = "ID", extra_fixed = ~1:(DENS+TREAT))

Which has resolved our issue.

We show another example using the multivariate dataset dataBEL and estimating $\theta$. We would first set theta = 1 and use model selection to choose an appropriate interaction structure, which gives us the following model, with a functional group structure:

model1 <- DImulti(prop = 2:5, y = "Y", eco_func = c("Var", "UN"), unit_IDs = 1, data = dataBEL,
                  FG = c("Grass", "Grass", "Leg", "Leg"), DImodel = "FG", extra_fixed = ~Density, 
                  method = "ML", theta = 1)

Comparing models: ML vs REML

The optional parameter method in the function DImulti() allows a user to change the estimation method used, between the default "REML", restricted maximum likelihood, and "ML", maximum likelihood.

"REML" was chosen as the default as it provides unbiased estimates, as opposed to "ML", however, "REML" is not appropriate for use with model comparison where fixed effects vary between models, as we can see from the following warning message:

XXDO NOT RUNXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

model1 <- DImulti(prop = 2:5, y = "Y", eco_func = c("Var", "UN"), unit_IDs = 1, data = dataBEL,
                  FG = c("Grass", "Grass", "Leg", "Leg"), DImodel = "FG", 
                  method = "REML")

model2 <- DImulti(prop = 2:5, y = "Y", eco_func = c("Var", "UN"), unit_IDs = 1, data = dataBEL,
                  FG = c("Grass", "Grass", "Leg", "Leg"), DImodel = "FG", extra_fixed = ~Density,
                  method = "REML")

anova(model1, model2)

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

In this case, we should use "ML" to fit the models for comparison, then refit the chosen model using "REML", for the unbiased estimates.

DImodelsMulti documentation built on May 29, 2024, 2:15 a.m.