7

I would like to plot the results of a multivariate logistic regression analysis (GLM) for a specific independent variables adjusted (i.e. independent of the confounders included in the model) relationship with the outcome (binary).

I have seen posts that recommend the following method using the predict command followed by curve, here's an example;

x     <- data.frame(binary.outcome, cont.exposure)
model <- glm(binary.outcome ~ cont.exposure, family=binomial, data=x)
plot(cont.exposure, binary.outcome, xlab="Temperature",ylab="Probability of Response") 
curve(predict(model, data.frame(cont.exposure=x), type="resp"), add=TRUE, col="red")

However this does not seem to work for multivariate regression models. I get the following error when I add 'age' (arbitrary - could be any variable of same length) as a confounding variable;

> x     <- data.frame(binary.outcome, cont.exposure, age)
> model <- glm(binary.outcome ~ cont.exposure + age, family=binomial, data=x)
> plot(cont.exposure, binary.outcome, xlab="Temperature",ylab="Probability of Response") 
> curve(predict(model, data.frame(cont.exposure=x), type="resp"), add=TRUE, col="red")
Error in model.frame.default(Terms, newdata, na.action = na.action, xlev = object$xlevels) : 
  variable lengths differ (found for 'age')
In addition: Warning message:
  'newdata' had 101 rows but variable(s) found have 698 rows

The above model is a simplified version of the models I'd like to run, but the principle is the same; I would like to plot the relationship between a binary outcome variable and a continuous exposure, independent of confounding factors..

It would be great to get either a workaround for the above, or an alternative way to view the relationship I am interested in. Many thanks.

Luke
  • 643
  • 1
  • 10
  • 17
  • 1
    You could have a look at the `crPlots` function in the `car` package. – BenBarnes Jul 02 '12 at 10:36
  • @BenBarnes thanks for that. I've had a look and a quick play with the data, and the function doesn't recognize that I am doing logistic regressions. However, the if I use a linear regression (ie. .my exposure is now my outcome, my binary variable an independent variable) then I do get exactly what I want. Will you post this as an answer for me to accept, or shall I? – Luke Jul 02 '12 at 11:01
  • I'm going to upvote Thierry's answer! – BenBarnes Jul 02 '12 at 11:49
  • @Luke did you ever find a good solution to this question? The accepted answer sort of missed the intent of the question, as you point out in your comment. – colin Mar 17 '16 at 18:34
  • @colin sadly no, and I have not seen anything in the last few years to do what I'm after really. Do you have any suggestions? – Luke Oct 11 '16 at 11:05
  • @Luke I've resorted to plotting the mean probability and its 95% credible interval over the range of the dependent variable, holding all other predictors constant at their mean. Looks nice, but would be nicer if there were data points on the figure that seemed in any way related. – colin Oct 20 '16 at 14:00

1 Answers1

8
set.seed(12345)
dataset <- expand.grid(Temp = rnorm(30), Age = runif(10))
dataset$Truth <- with(dataset, plogis(2 * Temp - 3 * Age))
dataset$Sample <- rbinom(nrow(dataset), size = 1, prob = dataset$Truth)
model <- glm(Sample ~ Temp + Age, data = dataset, family = binomial)
newdata <- expand.grid(
  Temp = pretty(dataset$Temp, 20), 
  Age = pretty(dataset$Age, 5))
newdata$Sample <- predict(model, newdata = newdata, type = "response")
library(ggplot2)
ggplot(newdata, aes(x = Temp, y = Sample)) + geom_line() + facet_wrap(~Age)

enter image description here

ggplot(newdata, aes(x = Temp, y = Sample, colour = Age, group = Age)) + 
  geom_line()

enter image description here

Thierry
  • 18,049
  • 5
  • 48
  • 66
  • 1
    Thanks for the response - I should have specified in my question, but I provided a simplified model in my example. There are actually numerous confounding variables (some of which are continuous, others factors), and I would like to vizualize the relationship between my outcome (binary) and an exposure (continuous), independent of the other variables I include in the model. Do you know whether ggplot can do this? Thanks again – Luke Jul 02 '12 at 13:23
  • @BenBarnes does provide a good method for doing this with continuous outcomes; by running a linear regression with my binary variable as a exposure I can get a nice plot, but this is not the original logistic framework I was using, so an alternative would also be great. – Luke Jul 02 '12 at 13:27
  • @LukeTheDuke: in that case set a range of values for the covariate of interest and fix all other covariates at a relevant level. – Thierry Jul 02 '12 at 17:53