####################################################################################################
####################################################################################################
# \dontshow{
## Set up for R markdown for crayon and cli output if user has packages installed
if(requireNamespace("fansi", quietly = TRUE) &
requireNamespace("crayon", quietly = TRUE) &
requireNamespace("cli", quietly = TRUE))
{
options(crayon.enabled = TRUE)
ansi_aware_handler <- function(x, options)
{
paste0(
"",
fansi::sgr_to_html(x = x, warn = FALSE, term.cap = "256"),
""
)
}
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
)
}
# }
## Modelling Example
# For a more thorough example of the workflow of this package, please see vignette
# DImulti_workflow using the following code:
# vignette("DImulti_workflow")
head(dataSWE)
# We begin by transforming the time identifier column to a factor, to better group coefficients
dataSWE$YEARN <- as.factor(dataSWE$YEARN)
# We fit a repeated measures DI model using the main function DImulti().
# We begin by passing the initial species proportions column names through 'prop' and the
# response value column name through 'y'.
# As the data contains multiple time points, we include the 'time' parameter, through which we
# specify the column name containing the repeated measure indicator and the autocorrelation
# structure we want to use, in this case, we use compound symmetry (CS).
# The experimental unit ID is stored in the column "PLOT" and so we pass this to 'unit_IDs'.
# We request that the method fits an average interaction structure and we include the treatments
# DENS and TREAT as additional fixed effects.
# We opt to use the ML estimation method to allow for model comparisons with different fixed
# effects.
# Finally, we specify the data object, dataSWE, through 'data'.
SWEmodel <- DImulti(y = "YIELD", time = c("YEARN", "CS"), unit_IDs = "PLOT",
prop = 5:8, data = dataSWE, DImodel = "AV",
extra_fixed = ~ 1:(DENS:TREAT),
method = "ML")
print(SWEmodel)
# We can now make any changes to the model's fixed effects and use a series of anovas (if the
# models are nested) or information criteria such as AIC, AICc, BIC, or BICc to select a final model
# We choose to change the interactions structure to functional groups, grouping the grasses and
# legumes, and simplify the extra_fixed terms, as an example, by only crossing the ID effect of G1
# with the treatments, and no longer crossing them with each other.
SWEmodel2 <- DImulti(y = "YIELD", time = c("YEARN", "CS"), unit_IDs = "PLOT",
prop = 5:8, data = dataSWE, DImodel = "FG", FG = c("G", "G", "L", "L"),
extra_fixed = ~ G1_ID:(DENS+TREAT),
method = "ML")
BICc(SWEmodel); BICc(SWEmodel2)
# We select the model with the lower information criteria, which is our functional group model, and
# refit it using the "REML" estimation method, for unbiased estimates.
SWEmodel2 <- DImulti(y = "YIELD", time = c("YEARN", "CS"), unit_IDs = "PLOT",
prop = 5:8, data = dataSWE, DImodel = "FG", FG = c("G", "G", "L", "L"),
extra_fixed = ~ G1_ID:(DENS+TREAT),
method = "REML")
# With this model, we can examine the coefficients
coef(SWEmodel2)
# or the variance and correlation structures
SWEmodel2$corr
nlme::getVarCov(SWEmodel2)
Run the code above in your browser using DataLab