#############################################################
### Inhibition constant variance ###
#############################################################
# mu1=signal mean, mu2=negative control mean
# sigma=SD,
# inh=percentage of inhibition
# inh=1 means 100% inhibition, i.e. reduces to negative control
# z=Z-factor
getPower.inh.cv <- function(mu1,mu2,sigma,inh,z=NULL){
if(is.null(z))
{ z=1-6*sigma/(mu1-mu2)
pp=pnorm(inh*(mu1-mu2)/sigma-3)}
pp=pnorm(inh*6/(1-z)-3)
x=c(round(z,3),round(pp,3))
names(x)=c("Z",paste(inh*100,"% Inhibition Power",sep=""))
return(x)
}
#############################################################
### Inhibition constant CV ###
#############################################################
# mu1=signal mean, mu2=negative control mean
# tau=CV,
# inh=percentage of inhibition
# inh=1 means 100% inhibition, i.e. reduces to negative control
# z=Z-factor
# sb=S/B
## if mu1 and mu2 are known, then either tau or z is required.
## if mu1 and mu2 are unknown, then at least two items of tau, sb and z are required
getPower.inh.ccv <- function(mu1,mu2,tau=NULL,inh,sb=NULL,z=NULL){
if(is.null(sb) & is.null(z))
{ z=1-3*tau*(mu1+mu2)/(mu1-mu2)
pp=pnorm((inh*(mu1-mu2)-3*tau*mu1)/(tau*(mu2+(1-inh)*(mu1-mu2))))}
if(is.null(sb) & is.null(tau))
{ tau=(1-z)*(mu1-mu2)/(3*(mu1+mu2))
pp=pnorm((inh*(mu1-mu2)-3*tau*mu1)/(tau*(mu2+(1-inh)*(mu1-mu2))))}
if(!is.null(sb) & is.null(z))
{ z=1-3*tau*(sb+1)/(sb-1)
pp=pnorm((inh*(sb-1)-3*tau*sb)/(tau*(inh+(1-inh)*sb)))}
if(!is.null(sb) & !is.null(z))
{ tau=(1-z)*(sb-1)/(3*(sb+1))
pp=pnorm(3*(inh*(sb+1)-sb*(1-z))/((1-z)*(inh+(1-inh)*sb)))}
if(!is.null(tau) & !is.null(z))
{ sb=(1-z+3*tau)/(1-z-3*tau)
pp=pnorm((6*inh-3*(3*tau+1-z))/(1-z+3*tau-6*inh*tau))}
x=c(round(z,3),round(tau,3),round(pp,3))
names(x)=c("Z","CV",paste(inh*100,"% Inhibition Power",sep=""))
return(x)
}
################ Figures
### power vs S/B
makefigure1A <- function(inh=0.5){
getp2 <- function(sb,tau,inh) pnorm((inh*(sb-1)-3*tau*sb)/(tau*(inh+(1-inh)*sb)))
tau=c(0.05,0.1,0.15,0.2,0.3)
sb=seq(0,15,0.1)
plot(sb,sapply(sb,function(x) getp2(x,tau[1],inh)),
ylab="power",xlim=c(1,15),type="l",xlab="Signal-to-background ratio")
lines(sb,sapply(sb,function(x) getp2(x,tau[2],inh)),lty=2)
lines(sb,sapply(sb,function(x) getp2(x,tau[3],inh)),lty=3)
lines(sb,sapply(sb,function(x) getp2(x,tau[4],inh)),lty=4)
lines(sb,sapply(sb,function(x) getp2(x,tau[5],inh)),lty=5)
legend(12,.5,legend=c("CV=.05",".1",".15",".2",".3"),lty=1:5,cex=.6)
title(paste(100*inh,"% inhibition"))
}
### power vs Z-factor
makefigure1B <- function(inh=0.5){
getp <- function(z,tau,inh) pnorm((6*inh-3*(3*tau+1-z))/(1-z+3*tau-6*inh*tau))
tau=.05;z=seq(-1,1-3*tau,len=50)
plot(z,getp(z=z,tau=tau,inh=inh),
ylab="power",ylim=c(0,1),xlim=c(-1,1),type="l",xlab="Z-factor")
tau=.1;z=seq(-1,1-3*tau,len=25)
lines(z,getp(z=z,tau=tau,inh=inh),lty=2)
tau=1/6;z=seq(-1,1-3*tau,len=25)
lines(z,getp(z=z,tau=tau,inh=inh),lty=3)
tau=.2;z=seq(-1,1-3*tau,len=25)
lines(z,getp(z=z,tau=tau,inh=inh),lty=4)
tau=.3;z=seq(-1,1-3*tau,len=25)
lines(z,getp(z=z,tau=tau,inh=inh),lty=5)
legend(-1,1,legend=c("CV=.05",".1",".15",".2",".3"),lty=1:5,cex=.6)
title(paste(100*inh,"% inhibition"))
}