代码之家 › 专栏 › 技术社区 › Leprechault

用poisson-glm方法绘制二次曲线

quadratic-curve poisson glm plot r

0

Leprechault · 技术社区 · 7 年前

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

#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的曲线。有什么解决办法吗?

1 回复 | 直到 7 年前

1

eipi10 7 年前

在下面的代码中,我们使用 poly 函数生成二次模型,而不需要在数据帧中创建额外的列。此外,我们还创建了一个预测数据帧,以生成跨范围的模型预测 condition 每个级别的 mating . 这个 predict type="response" 在结果的范围内生成预测,而不是在线性预测范围内生成预测,线性预测范围是默认值。而且,我们改变了 200 到 100 在创建 交配 为了避免在每一个水平上都有完全相同的结果数据 .

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