0

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!

jpsmith
  • 11,023
  • 5
  • 15
  • 36
Teresa
  • 15
  • 3
  • look here to start? https://stackoverflow.com/questions/55978360/how-to-plot-multiple-glmer-models-into-one-single-plot – CAWA Aug 31 '23 at 16:55
  • Thanks for your comment, @CAWA. However, this solution doesn´t include both models in the same plotting region. – Teresa Sep 01 '23 at 08:10

1 Answers1

1

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

CAWA
  • 123
  • 10