###################################################################################################
###################################################################################################
# \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(simRM)
# We call DImulti() to fit a series of models, with increasing complexity, and test whether the
# additional terms are worth keeping.
# We begin with an ID DI model, ensuring to use method = "ML" as we will be comparing fixed effects
RMmodel <- DImulti(y = "Y", time = c("time", "AR1"), unit_IDs = "plot",
prop = paste0("p", 1:4), data = simRM, DImodel = "ID",
method = "ML")
print(RMmodel)
# We can simplify the above model by grouping the ID effect terms into two categories, G1 and G2
RMmodel_ID <- DImulti(y = "Y", time = c("time", "AR1"), unit_IDs = "plot",
prop = paste0("p", 1:4), data = simRM, DImodel = "ID",
ID = c("G1", "G1", "G2", "G2"),
method = "ML")
# We can compare the two models by calling the anova function to perform a likelihood ratio test
anova(RMmodel, RMmodel_ID)
# As the p-value < an alpha value of 0.05, we reject the null hypthesis that the extra terms are
# equal to zero, therefore we continue with the larger model.
# We can confirm this result by selecting the model with the lower AIC value.
AIC(RMmodel); AIC(RMmodel_ID)
# Next, we include the simplest interaction structure available in this package, "AV", which adds
# a single extra term per ecosystem function
RMmodel_AV <- DImulti(y = "Y", time = c("time", "AR1"), unit_IDs = "plot",
prop = paste0("p", 1:4), data = simRM, DImodel = "AV",
method = "ML")
anova(RMmodel, RMmodel_AV)
# Once again, we select the more complex model.
#
# We can continue increasing the complexity of the interaction structure in the same fashion, this
# time we elect to use the additive interaction structure
RMmodel_ADD <- DImulti(y = "Y", time = c("time", "AR1"), unit_IDs = "plot",
prop = paste0("p", 1:4), data = simRM, DImodel = "ADD",
method = "ML")
anova(RMmodel_AV, RMmodel_ADD)
# We continue with the additive interactions structure, the next interaction structure to test
# functional group interactions, is not nested with our additive model, so we compare th two using
# information criterion. In this case we choose to use AICc, to better penalise extra terms, as AIC
# becomes biased towards the more complex model as parameters numbers increase.
# The functional group strcuture requires an additional parameter, FG, to specify groups.
RMmodel_FG <- DImulti(y = "Y", time = c("time", "AR1"), unit_IDs = "plot",
prop = paste0("p", 1:4), data = simRM, DImodel = "FG",
FG = c("G1", "G1", "G2", "G2"),
method = "ML")
AICc(RMmodel_ADD); AICc(RMmodel_FG)
# In this case, we choose the lower valued model, which is the additive structure.
#
# The last interaction structure to test is the full pairwise model. Which we can see is not worth
# the extra terms. So our final interaction structure chosen is additive.
RMmodel_FULL <- DImulti(y = "Y", time = c("time", "AR1"), unit_IDs = "plot",
prop = paste0("p", 1:4), data = simRM, DImodel = "FULL",
method = "ML")
anova(RMmodel_ADD, RMmodel_FULL)
# Finally, we can also increase the model complexity via the inclusion of the non-linear parameter
# theta, which we can estimate, or select a value for. We also choose to estimate using the "REML"
# method as we will do no further fixed effect model comparisons
RMmodel_theta <- DImulti(y = "Y", time = c("time", "AR1"), unit_IDs = "plot",
prop = paste0("p", 1:4), data = simRM, DImodel = "ADD",
estimate_theta = TRUE, method = "REML")
print(RMmodel_theta)
# We can however test changes to the correlation structure, testing the AR(1) structure we have been
# using against the simpler CS strcuture.
RMmodel_CS <- DImulti(y = "Y", time = c("time", "CS"), unit_IDs = "plot",
prop = paste0("p", 1:4), data = simRM, DImodel = "ADD",
theta = 0.9991, method = "REML")
AICc(RMmodel_theta); AICc(RMmodel_CS)
# We see that they are practically identical, but the AR(1) strcuture is preferred.
##################################################################################################
#
##################################################################################################
## Code to simulate data
# \donttest{
set.seed(523)
props <- data.frame(plot = integer(),
p1 = integer(),
p2 = integer(),
p3 = integer(),
p4 = integer())
index <- 1 #row number
#Monocultures
for(i in 1:4) #4 species
{
for(j in 1:2) #2 technical reps
{
props[index, i+1] <- 1
index <- index + 1
}
}
#Equal Mixtures
for(rich in sort(rep(2:4, 3))) #3 reps at each richness level
{
sp <- sample(1:4, rich) #randomly pick species from pool
for(j in 1:2) #2 technical reps
{
for(i in sp)
{
props[index, i+1] <- 1/rich #equal proportions
}
index <- index + 1
}
}
#Unequal Mixtures
for(rich in sort(rep(c(2, 3, 4), 15))) #15 reps at each richness level
{
sp <- sample(1:4, rich, replace = TRUE) #randomly pick species from pool
for(j in 1:2) #2 technical reps
{
for(i in 1:4)
{
props[index, i+1] <- sum(sp==i)/rich #equal proportions
}
index <- index + 1
}
}
props[is.na(props)] <- 0
mySimData <- props
mySimData$Y <- NA
ADDs <- DImodels::DI_data(prop=2:5, what=c("ADD"), data=mySimData)
mySimData <- cbind(mySimData, ADDs)
E_AV <- DImodels::DI_data(prop=2:5, what=c("E", "AV"), data=mySimData)
mySimData <- cbind(mySimData, E_AV)
mySimData$plot <- 1:nrow(mySimData)
mySimData$time <- 1
mySimDataT1 <- mySimData
mySimDataT2 <- mySimData
mySimDataT3 <- mySimData
mySimDataT2$time <- 2
mySimDataT3$time <- 3
n <- 3 #Number of Ys
p <- qr.Q(qr(matrix(stats::rnorm(n^2), n))) #Principal Components (make sure it's positive definite)
S <- crossprod(p, p*(n:1)) #Sigma
m <- stats::runif(n, -0.25, 1.5)
#runif(11, -1, 7) #decide on betas randomly
for(i in 1:nrow(mySimData))
{
#Within subject error
error <- MASS::mvrnorm(n=1, mu=m, Sigma=S)
mySimDataT1$Y[i] <- 4.1*mySimDataT1$p1[i] + 2.1*mySimDataT1$p2[i] + 3.6*mySimDataT1$p3[i] +
4.8*mySimDataT1$p4[i] + 3.3*mySimDataT1$p1_add[i] + 4.0*mySimDataT1$p2_add[i] +
1.3*mySimDataT1$p3_add[i] + 5.2*mySimDataT1$p4_add[i] + error[1]
mySimDataT2$Y[i] <- 2.3*mySimDataT2$p1[i] + 2.4*mySimDataT2$p2[i] + 5.1*mySimDataT2$p3[i] +
5.0*mySimDataT2$p4[i] + 3.6*mySimDataT2$p1_add[i] + 4.5*mySimDataT2$p2_add[i] +
0.5*mySimDataT2$p3_add[i] + 6.5*mySimDataT2$p4_add[i] + error[2]
mySimDataT3$Y[i] <- 0.7*mySimDataT3$p1[i] + 2.3*mySimDataT3$p2[i] + 6.5*mySimDataT3$p3[i] +
5.9*mySimDataT3$p4[i] + 4.5*mySimDataT3$p1_add[i] + 5.2*mySimDataT3$p2_add[i] +
0.6*mySimDataT3$p3_add[i] + 7.8*mySimDataT3$p4_add[i] + error[3]
}
mySimData <- rbind(mySimDataT1, mySimDataT2)
mySimData <- rbind(mySimData, mySimDataT3)
mySimData$time <- as.factor(mySimData$time)
mySimData$plot <- as.factor(mySimData$plot)
# }
Run the code above in your browser using DataLab