Plotting results of logistic regression with binom

2019-06-14 16:55发布

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:

enter image description here

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:

enter image description here

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?

1条回答
做个烂人
2楼-- · 2019-06-14 17:04

You can use the ggeffects-package for this, either with ggpredict() or ggeffect() (see ?ggpredict for the difference for these two functions, the first calls predict(), the latter effects::Effect()).

library(ggeffects)
library(sjmisc)
library(lme4)
data(mtcars)

mtcars <- std(mtcars, wt)
mtcars$am <- as.factor(mtcars$am)

m <- glmer(vs ~ wt_z + am + (1|carb), mtcars, family = binomial, na.action = "na.fail")

# Note the use of the "all"-tag here, see help for details
ggpredict(m, "wt_z [all]") %>% plot()

enter image description here

ggpredict(m, "am") %>% plot()

enter image description here

查看更多
登录 后发表回答