I'm trying to display the results of a logistic regression. My model was fit using glmer() from the lme4 package, I then used MuMIn for model averaging.
Simplified version of my model using the mtcars
dataset:
glmer(vs ~ wt + am + (1|carb), database, family = binomial, na.action = "na.fail")
My desired output is two plots that show the predicted probability that vs
=1, one for wt
, which is continuous, one for am
, which is binomial.
UPDATED:
I got this much working after comments from @KamilBartoń:
database <- mtcars
# Scale data
database$wt <- scale(mtcars$wt)
database$am <- scale(mtcars$am)
# Make global model
model.1 <- glmer(vs ~ wt + am + (1|carb), database, family = binomial, na.action = "na.fail")
# Model selection
model.1.set <- dredge(model.1, rank = "AICc")
# Get models with <10 delta AICc
top.models.1 <- get.models(model.1.set,subset = delta<10)
# Model averaging
model.1.avg <- model.avg(top.models.1)
# make dataframe with all values set to their mean
xweight <- as.data.frame(lapply(lapply(database[, -1], mean), rep, 100))
# add new sequence of wt to xweight along range of data
xweight$wt <- (wt = seq(min(database$wt), max(database$wt), length = 100))
# predict new values
yweight <- predict(model.1.avg, newdata = xweight, type="response", re.form=NA)
# Make plot
plot(database$wt, database$vs, pch = 20, xlab = "WEIGHT (g)", ylab = "VS")
# Add predicted line
lines(xweight$wt, yweight)
Produces:
The remaining issue is that the data are scaled and centred around 0, meaning interpretation of the graph is impossible. I'm able to unscale the data using an answer from @BenBolker to this question but this does not display correctly:
## Ben Bolker's unscale function:
## scale variable x using center/scale attributes of variable y
scfun <- function(x,y) {
scale(x,
center=attr(y,"scaled:center"),
scale=attr(y,"scaled:scale"))
}
## scale prediction frame with scale values of original data -- for all variables
xweight_sc <- transform(xweight,
wt = scfun(wt, database$wt),
am = scfun(am, database$am))
# predict new values
yweight <- predict(model.1.avg, newdata = xweight_sc, type="response", re.form=NA)
# Make plot
plot(mtcars$wt, mtcars$vs, pch = 20, xlab = "WEIGHT (g)", ylab = "VS")
# Add predicted line
lines(xweight$wt, yweight)
Produces:
I can see the plot line is there but it's in the wrong place. I've tried this a few different ways but can't work out what the problem is. What have I done wrong?
Also, another remaining issue: How do I make a binomial plot for am
?
You can use the ggeffects-package for this, either with
ggpredict()
orggeffect()
(see?ggpredict
for the difference for these two functions, the first callspredict()
, the lattereffects::Effect()
).