####################################
## Code to simulate the sim1 dataset
# \donttest{
## Simulate dataset sim1 with species identity effects and block effects, but no interaction effect
## Use the proportions from the first fifteen plots in Switzerland
data(Switzerland)
## Repeat the 15 plots over four blocks.
## Give each community type a unique (community) number.
sim1 <- data.frame(community = rep(1:15, each = 4),
block = factor(rep(1:4, times = 15)),
p1 = rep(Switzerland$p1[1:15], each = 4),
p2 = rep(Switzerland$p2[1:15], each = 4),
p3 = rep(Switzerland$p3[1:15], each = 4),
p4 = rep(Switzerland$p4[1:15], each = 4))
## To simulate the response, first create a matrix of predictors that includes
## p1-p4 and the four block dummy variables.
X <- model.matrix(~ p1 + p2 + p3 + p4 + block -1, data = sim1)
## Create a vector of 'known' parameter values for simulating the response.
## The first four are the p1-p4 parameters, the second four are the block effects.
sim1_coeff <- c(10,9,8,7, 1,1.5,2,0)
## Create response and add normally distributed error
sim1$response <- as.numeric(X %*% sim1_coeff)
set.seed(2020)
r <- rnorm(n = 60, mean = 0, sd = 1.1)
sim1$response <- round(sim1$response + r, digits = 3)
# }
###########################
## Analyse the sim1 dataset
## Load the sim1 data
data(sim1)
## View the first few entries
head(sim1)
## Explore the variables in sim1
str(sim1)
## Check that the proportions sum to 1 (required for DI models)
## p1 to p4 are in the 3rd to 6th columns in sim1
sim1sums <- rowSums(sim1[3:6])
summary(sim1sums)
## Check characteristics of sim1
hist(sim1$response)
summary(sim1$response)
plot(sim1$p1, sim1$response)
plot(sim1$p2, sim1$response)
plot(sim1$p3, sim1$response)
plot(sim1$p4, sim1$response)
## Find the best DI model using autoDI and F-test selection
auto1 <- autoDI(y = "response", prop = c("p1","p2","p3","p4"), block = "block", data = sim1,
selection = "Ftest")
summary(auto1)
## Fit the identity model using DI and the ID tag
m1 <- DI(y = "response", prop = c("p1","p2","p3","p4"), block = "block", DImodel = "ID",
data = sim1)
summary(m1)
plot(m1)
# \donttest{
## Check goodness-of-fit using a half-normal plot with a simulated envelope
library(hnp)
hnp(m1)
# }
## Fit the identity model using DI and custom_formula
m2 <- DI(y = "response", custom_formula = response ~ 0 + p1 + p2 + p3 + p4 + block, data = sim1)
summary(m2)
Run the code above in your browser using DataLab