16/11/18 12:03:04.69 Bq/wJd57.net
>>610
# debugged
RR2P=function(a,N1,b,N0,...){
SE=sqrt(1/a-1/N1+1/b-1/N0)
RR=(a/N1)/(b/N0)
f=function(x){
pl=(1-pnorm(log(x/RR)/SE))*2
pu=(1-pnorm(-log(x/RR)/SE))*2
p.value=ifelse(pl<=1,pl,pu)
return(p.value)
}
xx=seq(0,10,by=0.01)
xx=xx[-1]
plot(xx,sapply(xx,f),type="l",lwd=2,xlab="RR",ylab="p.value",...)
abline(h=0.05,lty=3,col=2)
abline(v=RR,lty=2,col=4)
legend("topleft",bty="n",legend=c("p=0.05","estimated RR"),col=c(2,4),lty=3)
names(RR)="Estimated RR"
return(RR)
}
RR2P(4,386,4,1250,xlim=c(0.01,10))
RR2P(4,386,4,1250,xlim=c(0.01,10),log="x")
RR2P(4,386,4,1250,xlim=c(0.01,10),ylim=c(0.001,1),log="y")
RR2P(4,386,4,1250,xlim=c(0.01,10),ylim=c(0.001,1),log="xy")