Построение квадратичных кривых с пуассоновским glm с взаимодействиями в категориальных / числовых переменных

Question

Я хочу знать, возможно ли построить квадратичные кривые с помощью Пуассона glm с взаимодействиями в категориальных / числовых переменных.В моем случае:

##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)

В ситуации, когда в модели значимы спаривание, условие и условие2, я делаю:

#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")
#

Не работает и у меня нетмои желательные изгибы.Я хотел бы кривую для T1 и другие для T2 в переменной спаривания.Есть какое-то решение для этого?

Answer 1

В приведенном ниже коде мы используем функцию poly для генерации квадратичной модели без необходимости создания дополнительного столбца во фрейме данных.Кроме того, мы создаем фрейм данных прогноза для генерации прогнозов модели в диапазоне значений condition и для каждого уровня mating.Функция predict с type="response" создает прогнозы по шкале результата, а не по шкале линейного предиктора, которая используется по умолчанию.Кроме того, мы изменили 200 на 100 при создании данных для mating, чтобы избежать получения одинаковых итоговых данных для каждого уровня 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")

Теперь строим график с помощью ggplot2и с основанием 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"))