varpred | R Documentation |
Computes predictor effect for generalized linear (mixed) models together with the associated confidence intervals anchored on some values, mostly, model center. It also incorporates proaches for correcting bias in predictions for GL(M)Ms involving nonlinear link functions.
varpred(
mod,
focal_predictors,
x.var = NULL,
type = c("response", "link"),
isolate = TRUE,
isolate.value = NULL,
level = 0.95,
steps = 100,
at = list(),
dfspec = 100,
true.beta = NULL,
vcov. = NULL,
internal = FALSE,
avefun = mean,
offset = NULL,
bias.adjust = c("none", "delta", "quantile", "population"),
sigma = c("mod", "lp", "total"),
include.re = FALSE,
modelname = NULL,
returnall = FALSE,
...
)
mod |
fitted model object. See details for supported class of models. |
focal_predictors |
a character vector of one or more predictors. For models with interaction, the marginal predictions are obtained by specifying the corresponding predictors together. For example |
x.var |
a character specifying the predictor to define the x variable (horizontal axis on the plot). The default is |
type |
a character specifying the desired prediction. |
isolate |
logical. If |
isolate.value |
numeric (default |
level |
desired confidence interval for computing marginal predictions. Default is |
steps |
number of points to evaluate numerical predictors in |
at |
default |
dfspec |
default |
vcov. |
a function or a matrix. If a function, it is used to compute the variance-covariance matrix of the model coefficients. The function should take model as it's first (or maybe only) argument. A matrix of variance-covariance matrix of the estimated coefficient can also be used. Otherwise |
internal |
logical. If |
avefun |
the averaging scheme (function) to be used in conditioning non-focal predictors. Default is |
offset |
a function or a value (FIXME:). |
bias.adjust |
specifies the bias correction method. If "none" (default), no bias correction method is applied; if "delta", delta method is used; if "population", all the values of non-focal predictors are used; otherwise, if "quantile", quantiles of non-focal numerical predictors are use. The options "quantile" and "population" (both EXPERIMENTAL) are used for bias correction in GL(M)M models involving non-linear link functions. |
sigma |
standard deviation used in delta method (only if |
include.re |
logical. If |
returnall |
logical. If |
character |
string naming the predictions. Useful when comparing several predictions. |
Predictor effects computes E(Y|X)
by meaningfully holding the non-focal predictors constant (or averaged in some meaningful way) while varying the focal predictor, with the goal that the response (E(Y|X)
) represents how the model responds to the changes in the focal predictor.
The traditional way to compute variances for predictions is \sigma^2 = \textrm{Diag}(\bX^\star \Sigma \bX^{\star\top})
, so that the confidence intervals are \eta \pm q\sigma
, where q
is an appropriate quantile of Normal or t distribution. This approach incorporates all the uncertainties – including the uncertainty due to non-focal predictors. But what if we are only interested in the uncertainty as a result of the focal predictor, so that the confidence intervals are \eta \pm q \sigma_f
(what we call anchored CIs)? There are two ways to anchor CIs: variance-covariance matrix based which requires properly scaled input predictors prior to model fitting; and centered model matrix which is more general and does not require scaled input predictors prior to model fitting.
Currently, the package supports lm, glm, lme4
and glmmTMB
models.
plot.vareffects
# Set theme for ggplot. Comment out if not needed
varefftheme()
set.seed(101)
N <- 100
x1_min <- 1
x1_max <- 9
b0 <- 0.3
b1 <- 0.1
b2 <- -0.6
b3 <- 0.01
x2_levels <- factor(c("A", "B", "D"))
df <- expand.grid(x1u = runif(n=N, min=x1_min, max=x1_max)
, x2u = x2_levels
)
X <- model.matrix(~x1u + x2u, df)
betas <- c(b0, b1, b2, b3)
df$y <- rnorm(nrow(df), mean= X %*% betas, sd=1)
df2 <- df
df <- transform(df
, x1c = drop(scale(x1u, scale=FALSE))
)
head(df)
# Unscaled model
m1u <- lm(y ~ x1u + x2u, df)
# Predictor rffects of x1u
pred1u <- varpred(m1u, "x1u")
plot(pred1u)
# Scaled model (x1 centered = x1 - mean(x1))
m1c <- lm(y ~ x1c + x2u, df)
# All uncertainities included
pred1c <- varpred(m1c, "x1c")
plot(pred1c)
# Centered predictor effects
# Results similar to m1c
# Using zero_vcov by specifying vcov.
vv <- zero_vcov(m1c, "x1c")
pred2c <- varpred(m1c, "x1c", vcov. = vv)
plot(pred2c)
# Using mean centering (isolate)
pred3c <- varpred(m1u, "x1u", isolate = TRUE)
plot(pred3c)
all.equal(pred2c$pred[,-1], pred3c$pred[,-1], check.attributes = FALSE)
# Compare across groups
pred4c <- varpred(m1c, c("x1c", "x2u"), x.var = "x1c", isolate = TRUE)
plot(pred4c)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.