####################################
## Code to simulate the sim3 dataset
# \donttest{
## Simulate dataset sim3 with 9 species, three functional groups and two levels of a treatment.
## The species 1-5 are FG1, species 6-7 are FG2 and species 8-9 are FG3.
## Assume ID effects and the FG interaction model, with a treatment (factor with two levels).
## Set up proportions
data("design_a")
sim3a <- design_a
# Replicate the design over two treatments
sim3b <- sim3a[rep(seq_len(nrow(sim3a)), each = 2), ]
sim3c <- data.frame(treatment = factor(rep(c("A","B"), times = 206)))
sim3 <- data.frame(sim3b[,1:2], sim3c, sim3b[,3:11])
row.names(sim3) <- NULL
## Create the functional group interaction variables
FG_matrix <- DI_data(prop = 4:12, FG = c("FG1","FG1","FG1","FG1","FG1","FG2","FG2","FG3","FG3"),
data = sim3, what = "FG")
sim3 <- data.frame(sim3, FG_matrix)
## To simulate the response, first create a matrix of predictors that includes p1-p9, the treatment
## and the interaction variables.
X <- model.matrix(~ p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 + treatment
+ bfg_FG1_FG2 + bfg_FG1_FG3 + bfg_FG2_FG3
+ wfg_FG1 + wfg_FG2 + wfg_FG3 -1, data=sim3)
## Create a vector of 'known' parameter values for simulating the response.
## The first nine are the p1-p9 parameters, and the second set of two are the treatment effects
## and the third set of six are the interaction parameters.
sim3_coeff <- c(10,9,8,7,11, 6,5, 8,9, 3,0, 4,9,3, 2,3,1)
## Create response and add normally distributed error
sim3$response <- as.numeric(X %*% sim3_coeff)
set.seed(1657914)
r <- rnorm(n = 412, mean = 0, sd = 1.2)
sim3$response <- round(sim3$response + r, digits = 3)
sim3[,13:18] <- NULL
# }
###########################
## Analyse the sim3 dataset
## Load the sim3 data
data(sim3)
## View the first few entries
head(sim3)
## Explore the variables in sim3
str(sim3)
## Check characteristics of sim3
hist(sim3$response)
summary(sim3$response)
plot(sim3$richness, sim3$response)
plot(sim3$p1, sim3$response)
plot(sim3$p2, sim3$response)
plot(sim3$p3, sim3$response)
plot(sim3$p4, sim3$response)
plot(sim3$p5, sim3$response)
plot(sim3$p6, sim3$response)
plot(sim3$p7, sim3$response)
plot(sim3$p8, sim3$response)
plot(sim3$p9, sim3$response)
# \donttest{
## What model fits best? Selection using F-test in autoDI
auto1 <- autoDI(y = "response", prop = 4:12, treat = "treatment",
FG = c("FG1","FG1","FG1","FG1","FG1","FG2","FG2","FG3","FG3"), data = sim3,
selection = "Ftest")
summary(auto1)
# }
## Fit the functional group model, with treatment, using DI and the FG tag
m1 <- DI(y = "response", prop = 4:12,
FG = c("FG1","FG1","FG1","FG1","FG1","FG2","FG2","FG3","FG3"), treat = "treatment",
DImodel = "FG", data = sim3)
summary(m1)
plot(m1)
# \donttest{
## Check goodness-of-fit using a half-normal plot with a simulated envelope
library(hnp)
hnp(m1)
# }
## Create the functional group interactions and store them in a new dataset
FG_matrix <- DI_data(prop = 4:12, FG = c("FG1","FG1","FG1","FG1","FG1","FG2","FG2","FG3","FG3"),
data = sim3, what = "FG")
sim3a <- data.frame(sim3, FG_matrix)
## Test if the FG interaction variables interact with treatment using 'extra_formula'
m2 <- DI(y = "response", prop = 4:12,
FG = c("FG1","FG1","FG1","FG1","FG1","FG2","FG2","FG3","FG3"),
treat = "treatment", DImodel = "FG", extra_formula = ~ bfg_FG1_FG2:treatment
+ bfg_FG1_FG3:treatment + bfg_FG2_FG3:treatment + wfg_FG1:treatment + wfg_FG2:treatment
+ wfg_FG3:treatment, data = sim3a)
summary(m2)
## Fit the functional group model using DI and custom_formula
## Set up a dummy variable for treatment first (required).
sim3a$treatmentA <- as.numeric(sim3a$treatment=="A")
m3 <- DI(y = "response", custom_formula = response ~ 0 + p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8
+ p9 + treatmentA + bfg_FG1_FG2 + bfg_FG1_FG3 + bfg_FG2_FG3 + wfg_FG1 + wfg_FG2
+ wfg_FG3, data = sim3a)
summary(m3)
Run the code above in your browser using DataLab