Reputation: 757
library(polynom)
set.seed(12345)
x<-as.numeric(seq(1,5,0.01))
lp<-rnorm(401,-0.7,1)-4*x+0.9*x^2
link_lp <- exp(lp)/(1 + exp(lp))
y<-ifelse((runif(401) < link_lp),0,1)
df<-data.frame(y=y,x=x)
f<-glm(y~poly(x,degree=2,raw=TRUE),family="binomial")
ymark<-c(0.2,0.3,0.35,0.4,0.45,0.5,0.6)
ggplot(df, aes(x=x, y=y))+
geom_jitter(size=1, alpha=0.2,height=0.05)+
stat_smooth(method="glm",
formula=y~poly(x,degree=2,raw=TRUE),
method.args =list(family = "binomial"),
colour="blue", size=1.5)+
xlab("x")+
ylab("Probability of event")+theme_bw()+
geom_hline(yintercept=ymark,col="red")+
scale_y_continuous(breaks=c(0,ymark,1))
Next I want to add vertical lines corresponding to the horizontal red lines and added another geom_vline()
ggplot(df, aes(x=x, y=y))+
geom_jitter(size=1, alpha=0.2,height=0.05)+
stat_smooth(method="glm",
formula=y~poly(x,degree=2,raw=TRUE),
method.args =list(family = "binomial"),
colour="blue", size=1.5)+
xlab("x")+
ylab("Probability of event")+theme_bw()+
geom_hline(yintercept=ymark,col="red")+
scale_y_continuous(breaks=c(0,ymark,1))+
geom_vline(xintercept=solve(polynomial(coef = coef(f)),
b=log(ymark/(1-ymark))),col="green")
However, it reported error, how can I change my code?
Upvotes: 1
Views: 173
Reputation: 1948
solve
returns a complex number. If you just want the real part put Re
around it.
ggplot(df, aes(x=x, y=y))+
geom_jitter(size=1, alpha=0.2,height=0.05)+
stat_smooth(method="glm",
formula=y~poly(x,degree=2,raw=TRUE),
method.args =list(family = "binomial"),
colour="blue", size=1.5)+
xlab("x")+
ylab("Probability of event")+theme_bw()+
geom_hline(yintercept=ymark,col="red")+
scale_y_continuous(breaks=c(0,ymark,1))+
geom_vline(xintercept=Re(solve(polynomial(coef = coef(f)),
b=log(ymark/(1-ymark)))),col="green")
Upvotes: 2