2017-02-02 14 views
0
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)) 

figure 1をどのように私は、水平赤線に対応する縦線を追加する

次にggplotでYに与えられた値に対応するx値における垂直線を追加し、別の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") 

ただし、エラーが報告されました。コードを変更するにはどうすればよいですか?

答えて

1

solveは、複素数を返します。実際の部分がほしいのであれば、その周りにReを入れてください。

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

enter image description here

関連する問題