Plotting quadratic curves with poisson glm with interactions in categorical/numeric variables

Question

I want to know if it's possible to plot quadratic curves with Poisson glm with interactions in categorical/numeric variables. In my case:

##Data set artificial
set.seed(20)
d <- data.frame(
  behv = c(rpois(100,10),rpois(100,100)),
  mating=sort(rep(c("T1","T2"), 200)),
  condition = scale(rnorm(200,5))
) 

#Condition quadratic
d$condition2<-(d$condition)^2

#Binomial GLM ajusted
md<-glm(behv ~ mating + condition + condition2, data=d, family=poisson)
summary(md)

In a situation where mating, condition and condition2 are significant in the model, I make:

#Create x's vaiues
x<-d$condition## 
x2<-(d$condition)^2 

# T1 estimation
y1<-exp(md$coefficients[1]+md$coefficients[3]*x+md$coefficients[4]*x2)
#
# T2 estimation
y2<-exp(md$coefficients[1]+md$coefficients[2]+md$coefficients[3]*x+md$coefficients[4]*x2)
#
#
#Separete data set
d_T1<-d[d[,2]!="T2",] 
d_T2<-d[d[,2]!="T1",] 

#Plot
plot(d_T1$condition,d_T1$behv,main="", xlab="condition", ylab="behv", 
xlim=c(-4,3), ylim=c(0,200), col= "black")
points(d_T2$condition,d_T2$behv, col="gray")
lines(x,y1,col="black")
lines(x,y2,col="grey")
#

Doesn't work and I don't have my desirable curves. I'd like a curve for T1 and other for T2 in mating variable. There are any solution for this?

Answer 1

In the code below, we use the poly function to generate a quadratic model without needing to create an extra column in the data frame. In addition, we create a prediction data frame to generate model predictions across the range of condition values and for each level of mating. The predict function with type="response" generates predictions on the scale of the outcome, rather than on the linear predictor scale, which is the default. Also, we change 200 to 100 in creating the data for mating in order to avoid having the exact same outcome data for each level of mating.

library(ggplot2)

# Fake data
set.seed(20)
d <- data.frame(
  behv = c(rpois(100,10),rpois(100,100)),
  mating=sort(rep(c("T1","T2"), 100)),   # Changed from 200 to 100
  condition = scale(rnorm(200,5))
)

# Model with quadratic condition
md <- glm(behv ~ mating + poly(condition, 2, raw=TRUE), data=d, family=poisson)
#summary(md)

# Get predictions at range of condition values
pred.data = data.frame(condition = rep(seq(min(d$condition), max(d$condition), length=50), 2),
                       mating = rep(c("T1","T2"), each=50))
pred.data$behv = predict(md, newdata=pred.data, type="response")

Now plot with ggplot2 and with base R:

ggplot(d, aes(condition, behv, colour=mating)) +
  geom_point() +
  geom_line(data=pred.data)

plot(NULL, xlim=range(d$condition), ylim=range(d$behv),
     xlab="Condition", ylab="behv")
with(subset(d, mating=="T1"), points(condition, behv, col="red"))
with(subset(d, mating=="T2"), points(condition, behv, col="blue"))
with(subset(pred.data, mating=="T1"), lines(condition, behv, col="red"))
with(subset(pred.data, mating=="T2"), lines(condition, behv, col="blue"))
legend(-3, 70, title="Mating", legend=c("T1","T2"), pch=1, col=c("blue", "red"))