### R code from vignette source 'pwrprc.rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: pwrprc.rnw:62-63
###################################################
c(20,30) * 7.94


###################################################
### code chunk number 2: pwrprc.rnw:72-73
###################################################
exp( 1.96/sqrt(c(159,238)) )


###################################################
### code chunk number 3: pwrprc.rnw:94-102
###################################################
library( Epi )
nsim <- 10000
prpw <- NArray( list( RR = seq(1.05,1.5,0.05),
                      pr = seq(5,50,5),
                      PY = 8:13,
                     inc = c(20,30),
                    what = c("RR","Erf","Pwr") ) )
str( prpw )


###################################################
### code chunk number 4: pwrprc.rnw:106-132
###################################################
for( ir in dimnames(prpw)[["RR"]] )
for( ip in dimnames(prpw)[["pr"]] )
for( iy in dimnames(prpw)[["PY"]] )
for( ii in dimnames(prpw)[["inc"]] )
{
nr <- as.numeric(ir)
np <- as.numeric(ip)/100
ny <- as.numeric(iy)
ni <- as.numeric(ii)
# Cumulative incidence in exposed and non-exposed:
cuminc <- ny * c(np,1-np) * ni * c(nr,1)
# Events in exposed:
nx <- rpois( nsim, cuminc[1] )
# Events in unexposed:
nr <- rpois( nsim, cuminc[2] )
# Rate Ratio
rr <-  (nx/np) / (nr/(1-np))
# Error factor
erf <- exp( 1.96*sqrt(1/nx+1/nr) )
# Significance
sgn <- abs(log(rr)) > log(erf)
# Store results
prpw[ir,ip,iy,ii,"RR" ] <- median( rr )
prpw[ir,ip,iy,ii,"Erf"] <- median( erf )
prpw[ir,ip,iy,ii,"Pwr"] <- mean( sgn )
}


###################################################
### code chunk number 5: pwrprc.rnw:136-138
###################################################
round( ftable( prpw[,,"12",,"RR" ],
               col.vars=1, row.vars=c(3,2) ), 2 )


###################################################
### code chunk number 6: pwrprc.rnw:143-145
###################################################
round( ftable( prpw[,,c("8","12"),,"Erf"],
               col.vars=1, row.vars=c(3,4,2) ), 2 )


###################################################
### code chunk number 7: pwrprc.rnw:149-151
###################################################
round( ftable( prpw[,,c("8","12"),,"Pwr"]*100,
               col.vars=1, row.vars=c(3,4,2) ), 1 )


###################################################
### code chunk number 8: pwrprc.rnw:154-157
###################################################
print( ftable( (prpw[,,c("8","12"),,"Pwr"]>0.8)*1,
               col.vars=1, row.vars=c(3,4,2) ),
       zero.print="." )