##Halve WaterRunoff data to reduce time to execute
data(WaterRunoff.dat)
tmp <- subset(WaterRunoff.dat, Date == "05-18")
##Use asreml to get predictions and associated statistics
if (FALSE) {
#Analyse pH
m1.asr <- asreml(fixed = pH ~ Benches + (Sources * (Type + Species)),
random = ~ Benches:MainPlots,
keep.order=TRUE, data= tmp)
current.asrt <- as.asrtests(m1.asr, NULL, NULL)
current.asrt <- as.asrtests(m1.asr)
current.asrt <- rmboundary(current.asrt)
m1.asr <- current.asrt$asreml.obj
#Get predictions and associated statistics
TS.diffs <- predictPlus.asreml(classify = "Sources:Type",
asreml.obj = m1.asr, tables = "none",
wald.tab = current.asrt$wald.tab,
present = c("Type","Species","Sources"))
#Use sort.alldiffs and save order for use with other response variables
TS.diffs.sort <- sort(TS.diffs, sortFactor = "Sources",
sortParallelToCombo = list(Type = "Control"))
sort.order <- attr(TS.diffs.sort, which = "sortOrder")
#Analyse Turbidity
m2.asr <- asreml(fixed = Turbidity ~ Benches + (Sources * (Type + Species)),
random = ~ Benches:MainPlots,
keep.order=TRUE, data= tmp)
current.asrt <- as.asrtests(m2.asr)
#Use pH sort.order to sort Turbidity alldiffs object
diffs2.sort <- predictPlus(m2.asr, classify = "Sources:Type",
pairwise = FALSE, error.intervals = "Stand",
tables = "none", present = c("Type","Species","Sources"),
sortFactor = "Sources",
sortOrder = sort.order)
}
## Use lmeTest and emmmeans to get predictions and associated statistics
if (requireNamespace("lmerTest", quietly = TRUE) &
requireNamespace("emmeans", quietly = TRUE))
{
#Analyse pH
m1.lmer <- lmerTest::lmer(pH ~ Benches + (Sources * (Type + Species)) +
(1|Benches:MainPlots),
data=na.omit(tmp))
TS.emm <- emmeans::emmeans(m1.lmer, specs = ~ Sources:Type)
TS.preds <- summary(TS.emm)
den.df <- min(TS.preds$df, na.rm = TRUE)
## Modify TS.preds to be compatible with a predictions.frame
TS.preds <- as.predictions.frame(TS.preds, predictions = "emmean",
se = "SE", interval.type = "CI",
interval.names = c("lower.CL", "upper.CL"))
## Form an all.diffs object and check its validity
TS.vcov <- vcov(TS.emm)
TS.diffs <- allDifferences(predictions = TS.preds,
classify = "Sources:Type",
vcov = TS.vcov, tdf = den.df)
validAlldiffs(TS.diffs)
#Use sort.alldiffs and save order for use with other response variables
TS.diffs.sort <- sort(TS.diffs, sortFactor = "Sources",
sortParallelToCombo = list(Type = "Control"))
sort.order <- attr(TS.diffs.sort, which = "sortOrder")
#Analyse Turbidity
m2.lmer <- lmerTest::lmer(Turbidity ~ Benches + (Sources * (Type + Species)) +
(1|Benches:MainPlots),
data=na.omit(tmp))
TS.emm <- emmeans::emmeans(m2.lmer, specs = ~ Sources:Type)
TS.preds <- summary(TS.emm)
den.df <- min(TS.preds$df, na.rm = TRUE)
## Modify TS.preds to be compatible with a predictions.frame
TS.preds <- as.predictions.frame(TS.preds, predictions = "emmean",
se = "SE", interval.type = "CI",
interval.names = c("lower.CL", "upper.CL"))
## Form an all.diffs object, sorting it using the pH sort.order and check its validity
TS.vcov <- vcov(TS.emm)
TS.diffs2.sort <- allDifferences(predictions = TS.preds,
classify = "Sources:Type",
vcov = TS.vcov, tdf = den.df,
sortFactor = "Sources",
sortOrder = sort.order)
validAlldiffs(TS.diffs2.sort)
}
Run the code above in your browser using DataLab