# Input p under the null hypothesis and the cells of a 2x2 table (in # row order). The output includes pexact$pval (for the one-sided # alternative p_2>p_1) , pexact$rr, the tables with a larger # pooled Z-test statistic than pexact$zobs, the observed test stat. pexact=function(p0,a0,b0,c0,d0){ np10=a0+b0 np110=np10+1 np20=c0+d0 np210=np20+1 n=a0+b0+c0+d0 # Generate a grid of all possible outcomes a=rep(0:np10,rep(np210,np110)) c=rep(0:np20,np110) # Construct pooled Z-test statistics p=(a+c)/n np1=a+c p1=a/np10 p2=c/np20 den=sqrt(p*(1-p)*(1/np10+1/np20)) z=(p2-p1)/den t1=cbind(a=a,c=c,z=z) # Find the observed test statistic in the list of test stats zobs=t1[,3][(t1[,1]==a0) & (t1[,2]==c0)] # Find all outcomes in the rejection region rr=t1[t1[,3]>=zobs,] rr=rr[!is.na(rr[,3]),] # Compute the p-value x=dbinom(rr[,1],np10,p0)*dbinom(rr[,2],np20,p0) pval=sum(x[!is.na(x)]) out=list(pval=pval,rr=rr,zobs=zobs) out }