###################################################################################################
###################################################################################################
# \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")
# We use head() to view the dataset
head(simMV)
# We can use the function DImulti() to fit a multivariate DI model, with an intercept for "treat"
# for each ecosystem function. The dataset is wide, so the Y columns are all entered through 'y'
# and the first index of eco_func is "NA". We fit the average interaction structure and use "ML"
# so that we can perform model comparisons with varying fixed effects
MVmodel <- DImulti(y = paste0("Y", 1:4), eco_func = c("NA", "UN"), unit_IDs = 1,
prop = paste0("p", 1:6), data = simMV, DImodel = "AV", extra_fixed = ~ treat,
method = "ML")
print(MVmodel)
# We can adjust the previous model to now cross "treat" with each other ID effect. We also specify
# different values of theta for each ecosystem function and use the "REML" method to get unbiased
# estimates.
# \donttest{
MVmodel_theta <- DImulti(y = paste0("Y", 1:4), eco_func = c("NA", "UN"), unit_IDs = 1,
prop = paste0("p", 1:6), data = simMV, DImodel = "AV", extra_fixed = ~ 1:treat,
theta = c(1, 0.5, 0.8, 0.6), method = "REML")
summary(MVmodel_theta)
# We can now use any S3 method compatible with a gls object, for example, predict()
predict(MVmodel_theta, newdata = simMV[which(simMV$plot == 1), ])
# }
##################################################################################################
#
##################################################################################################
## Code to simulate data
#
#
# \donttest{
set.seed(412)
props <- data.frame(plot = integer(),
p1 = integer(),
p2 = integer(),
p3 = integer(),
p4 = integer(),
p5 = integer(),
p6 = integer())
index <- 1 #row number
#Monocultures
for(i in 1:6) #6 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:6, 3))) #3 reps at each richness level
{
sp <- sample(1:6, 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, 5, 6), 15))) #15 reps at each richness level
{
sp <- sample(1:6, rich, replace = TRUE) #randomly pick species from pool
for(j in 1:2) #2 technical reps
{
for(i in 1:6)
{
props[index, i+1] <- sum(sp==i)/rich #equal proportions
}
index <- index + 1
}
}
props[is.na(props)] <- 0
mySimData <- props
mySimData$treat <- 0
mySimDataDupe <- mySimData
mySimDataDupe$treat <- 1
mySimData <- rbind(mySimData, mySimDataDupe)
mySimData$plot <- 1:nrow(mySimData)
mySimData$Y1 <- NA
mySimData$Y2 <- NA
mySimData$Y3 <- NA
mySimData$Y4 <- NA
ADDs <- DImodels::DI_data(prop=2:7, what=c("ADD"), data=mySimData)
mySimData <- cbind(mySimData, ADDs)
E_AV <- DImodels::DI_data(prop=2:7, what=c("E", "AV"), data=mySimData)
mySimData <- cbind(mySimData, E_AV)
n <- 4 #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(8, -1, 7) #decide on betas randomly
for(i in 1:nrow(mySimData))
{
#Within subject error
error <- MASS::mvrnorm(n=1, mu=m, Sigma=S)
mySimData$Y1[i] <- 6.9*mySimData$p1[i] + -0.3*mySimData$p2[i] + 6.6*mySimData$p3[i] +
1.7*mySimData$p4[i] + -0.8*mySimData$p5[i] + 4.3*mySimData$p6[i] + 3.5*mySimData$treat[i] +
1.8*mySimData$AV[i] + error[1]
mySimData$Y2[i] <- 4.9*mySimData$p1[i] + 3.6*mySimData$p2[i] + 4.4*mySimData$p3[i] +
2.3*mySimData$p4[i] + 4.3*mySimData$p5[i] + 6.6*mySimData$p6[i] + -0.3*mySimData$treat[i] +
6.8*mySimData$AV[i] + error[2]
mySimData$Y3[i] <- -0.3*mySimData$p1[i] + 4.6*mySimData$p2[i] + 1.2*mySimData$p3[i] +
6.8*mySimData$p4[i] + 1.4*mySimData$p5[i] + 6.9*mySimData$p6[i] + 5.5*mySimData$treat[i] +
1.4*mySimData$AV[i] + error[3]
mySimData$Y4[i] <- 4.1*mySimData$p1[i] + 4.2*mySimData$p2[i] + -0.5*mySimData$p3[i] +
4.9*mySimData$p4[i] + 6.7*mySimData$p5[i] + -0.9*mySimData$p6[i] + -1.0*mySimData$treat[i] +
0.3*mySimData$AV[i] + error[4]
}
mySimData$treat <- as.factor(mySimData$treat)
# }
Run the code above in your browser using DataLab