# Predict effects (direct, total)
m <- shipley.sem
e <- shipley.sem.eff
dir <- getDirEff(e)
tot <- getTotEff(e)
f.dir <- predEff(m, effects = dir, type = "response")
f.tot <- predEff(m, effects = tot, type = "response")
f.dir$Live[1:10]
f.tot$Live[1:10]
# Using new data for predictors
d <- na.omit(shipley)
xn <- c("lat", "DD", "Date", "Growth")
seq100 <- function(x) seq(min(x), max(x), length = 100)
nd <- data.frame(sapply(d[xn], seq100))
f.dir <- predEff(m, nd, dir, type = "response")
f.tot <- predEff(m, nd, tot, type = "response")
f.dir$Live[1:10]
f.tot$Live[1:10]
# Add CIs
# dir.b <- getDirEff(e, "boot")
# tot.b <- getTotEff(e, "boot")
# f.dir <- predEff(m, nd, dir, dir.b, type = "response")
# f.tot <- predEff(m, nd, tot, tot.b, type = "response")
# Predict an interactive effect (e.g. Live ~ Growth * DD)
xn <- c("Growth", "DD")
d[xn] <- scale(d[xn]) # scale predictors (improves fit)
m <- lme4::glmer(Live ~ Growth * DD + (1 | site) + (1 | tree),
family = binomial, data = d)
nd <- with(d, expand.grid(
Growth = seq100(Growth),
DD = mean(DD) + c(-sd(DD), sd(DD)) # two levels for DD
))
f <- predEff(m, nd, type = "response", interaction = "Growth:DD")
f$fit[1:10]
f$interaction
# Add CIs (need to bootstrap model...)
# system.time(B <- bootEff(m, R = 1000, ran.eff = "site"))
# f <- predEff(m, nd, B, type = "response", interaction = "Growth:DD")
# Model-averaged predictions (several approaches)
m <- shipley.growth # candidate models (list)
w <- runif(length(m), 0, 1) # weights
e <- stdEff(m, w) # averaged effects
f1 <- predEff(m[[1]], effects = e) # pass avg. effects
f2 <- predEff(m, weights = w) # pass weights argument
f3 <- avgEst(predEff(m), w) # use avgEst function
stopifnot(all.equal(f1, f2))
stopifnot(all.equal(f2, f3))
# Compare model fitted values: predEff() vs. fitted()
m <- shipley.sem$Live
f1 <- predEff(m, unique.eff = FALSE, re.form = NULL, type = "response")
f2 <- fitted(m)
stopifnot(all.equal(f1, f2))
# Compare predictions using standardised vs. raw effects (same)
f1 <- predEff(m)
f2 <- predEff(m, use.raw = TRUE)
stopifnot(all.equal(f1, f2))
Run the code above in your browser using DataLab