# NOT RUN {
## Predict effects (direct, total)
m <- Shipley.SEM
e <- Shipley.SEM.Eff
dir <- dirEff(e)
tot <- totEff(e)
f.dir <- predEff(m, effects = dir, type = "response")
f.tot <- predEff(m, effects = tot, type = "response")
## 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")
## Add CI's
# }
# NOT RUN {
dir.b <- dirEff(e, "boot")
tot.b <- totEff(e, "boot")
f.dir <- predEff(m, nd, dir, dir.b, type = "response")
f.tot <- predEff(m, nd, tot, tot.b, type = "response")
# }
# NOT RUN {
## Predict an interactive effect (e.g. Live ~ Growth * DD)
xn <- c("Growth", "DD")
d[xn] <- scale(d[xn]) # standardise 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")
## Add CI's (need to bootstrap model - will take a while)
# }
# NOT RUN {
system.time(B <- bootEff(m, ran.eff = "site", R = 1000))
est <- B$t0; est.b <- B$t # estimates
f <- predEff(m, nd, est, est.b, type = "response", interaction = "Growth:DD")
# }
# NOT RUN {
## Model-averaged predictions (several approaches)
m <- Shipley.Growth # candidate models (list)
w <- runif(length(m), 0, 1) # weights
e <- stdCoeff(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. predict
m <- Shipley.SEM$Live
f1 <- predEff(m, unique.x = FALSE, re.form = NULL)
f2 <- predict(m)
stopifnot(all.equal(f1, f2))
# }
Run the code above in your browser using DataLab