Created
July 24, 2018 22:30
-
-
Save conjugateprior/844ca2899277bfbaaf150f4a8c3562a0 to your computer and use it in GitHub Desktop.
Predicted probabilities for logistic regression models using R and ggplot2
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# convenience function for logit models | |
invlogit <- function(x){ 1 / (1 + exp(-x)) } | |
# some data cribbed from the R help pages | |
dd <- data.frame(ldose = rep(0:5, 2), | |
dead = c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16), | |
sex = rep(c("M", "F"), each = 6)) | |
# 30 trials in each row of which 'dead' beasties died | |
# fit a model | |
mod <- glm(cbind(dead, 20 - dead) ~ sex * ldose, family = binomial, data = dd) | |
# cases / range that we want predicted probabilities for | |
ndata <- expand.grid(ldose = seq(0, 5, 0.1), sex = c("M", "F")) | |
# get predicted probabilities | |
preds <- predict(mod, newdata = ndata, se.fit = TRUE) | |
# construct intervals by transforming predictions and 2*SEs | |
# and bundle with ndata for ease of plotting together | |
plottable <- with(preds, | |
data.frame(ndata, | |
fit = invlogit(fit), | |
lwr = invlogit(fit + 2*se.fit), | |
upr = invlogit(fit - 2*se.fit)) ) | |
# plot the lot | |
library(ggplot2) | |
theme_set(theme_minimal()) # replace default 'mouldy waffle' background | |
ggplot(plottable, aes(x = ldose, y = fit, col = sex, fill = sex)) + | |
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.4, colour = NA) + | |
geom_line(size = 1, alpha = 0.4) + | |
scale_color_brewer(palette="Set1") + # replace default colors and fill | |
scale_fill_brewer(palette="Set1") + | |
labs(x = "Dose", y = "Probability of death", title = "Dead Beetles") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Many thanks for sharing the code. There are some issues for me about the code. I just copy-pasted the code to RStudio and run it.
Here are my issues:
ggplot shows Male in Pink and Female in Blue. Blue is the traditional color to represent Male, and Pink is the traditional color to represent Female in world. So, is there an error in the code while labelling the gender in legend of the plot? Or labelling was done without caring their traditional coloring?
You say, " 30 trials in each row of which 'dead' beasties died"
Is it 30 or 12? Could you please explain the experiment design and problem you deal with this code a bit further?
I couldn't grasp the problem that this code solved.
Best and warmest regards.