Search code examples
rplotglmsjplot

How to plot two glms in one plot with double axis?

Imagine that I have the sample data below:

RespVar1 <- runif(n=18, min=55, max=120)
RespVar2 <- runif(n=18, min=0.3, max=0.5)
PredVar <- c(-2, -1, 0, 1, 2, 3)
df <- data.frame(RespVar1, RespVar2, PredVar)

Where I run these models:

M1 <- glm(RespVar1 ~ PredVar, data=df, family=gaussian())
M2 <- glm(RespVar2 ~ PredVar, data=df, family=gaussian())

My question is: How can I plot the two model plots below (including the data and the model line) in a double-axis plot in R? (for instance, RespVar1 on the left axis, RespVar2 on the right axis, and PredVar on the x axis)

library(sjPlot)    
plot_model(M1, type="pred", terms="PredVar", show.data=T)
plot_model(M2, type="pred", terms="PredVar", show.data=T)

Thanks!

Solution

  • This technically does want you want, you will have adjust how you plot your lines so that the axes make more sense. But in essence this is it. You could replace the lines() calls for the CI's with polygon() calls and a col=c() setting that matches your needs. You will also probably need to adjust the label values for the second y axis to your specific needs.

    FWIW, a much wiser person than me said: "if you need 2 Y axes, you're doing it wrong"

    Your data

    RespVar1 <- runif(n=18, min=55, max=120)
RespVar2 <- runif(n=18, min=0.3, max=0.5)
PredVar <- c(-2, -1, 0, 1, 2, 3)
df <- data.frame(RespVar1, RespVar2, PredVar)

    The models

    M1 <- glm(RespVar1 ~ PredVar, data=df, family=gaussian())
M2 <- glm(RespVar2 ~ PredVar, data=df, family=gaussian())

    'New data' for predictions, a lazy way to do this.

    new.data<-data.frame(PredVar=PredVar)

    Predictions

    pred1<-predict(M1, newdata = new.data, se=TRUE)
pred2<-predict(M2, newdata = new.data, se=TRUE)

    Confidence intervals for plotting

    ci_lwr1 <- with(pred1, fit + qnorm(0.025)*se.fit)
ci_upr1 <- with(pred1, fit + qnorm(0.975)*se.fit)

ci_lwr2 <- with(pred2, fit + qnorm(0.025)*se.fit)
ci_upr2 <- with(pred2, fit + qnorm(0.975)*se.fit)

    The plot:

    plot(pred1$fit~new.data$PredVar, type="l", ylim=c(0,120), col='red')
lines(ci_lwr1~new.data$PredVar, col="red")
lines(ci_upr1~new.data$PredVar, col="red")

lines(pred2$fit~new.data$PredVar, col="blue")
lines(ci_lwr2~new.data$PredVar, col="blue") # CIs are hard to see
lines(ci_upr2~new.data$PredVar, col="blue") # CIs are hard to see

points(RespVar1~PredVar, data=df)
points(RespVar2~PredVar, data=df)
axis(4, at=c(0,20,40,60,80, 100, 120), labels=round(seq(0,1,length=7),2))

    enter image description here