### R code from vignette source 'C:/Bendix/teach/NSCE/2022/pracs/cohort-sol.rnw'

###################################################
### code chunk number 1: cohort-sol.rnw:33-43
###################################################
t <- 1
r <- 150/100000
rr <- 2
G <- factor( c("ctr","X") )
N <- c(8000,2000)
Y <- N * t
E <- Y * c(1,rr) * r
D <- rpois( 2, E )
# and print the results nicely
data.frame( G, N, Y, E, D )


###################################################
### code chunk number 2: cohort-sol.rnw:47-51
###################################################
rates <- D/Y
RR <- rates[2]/rates[1]
erf <- exp(1.96 * sqrt(sum(1/D)) )
round( c( RR, RR/erf, RR*erf, erf ), 3 )


###################################################
### code chunk number 3: cohort-sol.rnw:58-60
###################################################
library(Epi)
round( ci.exp( glm( cbind(D,Y) ~ G, family=poisreg ) ), 3 )


###################################################
### code chunk number 4: cohort-sol.rnw:65-80
###################################################
t <- 3
r <- 150/100000
rr <- 2
G <- factor( c("ctr","X") )
N <- c(8000,2000)
Y <- N * t
E <- Y * c(1,rr) * r
D <- rpois( 2, E )
# and print the results nicely
data.frame( G, N, Y, E, D )
rates <- D/Y
RR <- rates[2]/rates[1]
erf <- exp(1.96 * sqrt(sum(1/D)) )
round( c( RR, RR/erf, RR*erf, erf ), 3 )
round( ci.exp( glm( cbind(D,Y) ~ G, family=poisreg ) ), 3 )


###################################################
### code chunk number 5: cohort-sol.rnw:89-104
###################################################
t <- 1
r <- 150/100000
rr <- 2
G <- factor( c("ctr","X") )
N <- c(6000,4000)
Y <- N * t
E <- Y * c(1,rr) * r
D <- rpois( 2, E )
# and print the results nicely
data.frame( G, N, Y, E, D )
rates <- D/Y
RR <- rates[2]/rates[1]
erf <- exp(1.96 * sqrt(sum(1/D)) )
round( c( RR, RR/erf, RR*erf, erf ), 3 )
round( ci.exp( glm( cbind(D,Y) ~ G, family=poisreg ), subset="G" ), 3 )


###################################################
### code chunk number 6: cohort-sol.rnw:118-129
###################################################
sim <- function( t, N )
{
r <- 150/100000
rr <- 2
G <- factor( c("ctr","X") )
Y <- N * t
E <- Y * c(1,rr) * r
D <- rpois( 2, E )
erf <- exp(1.96 * sqrt(sum(1/D)) )
c( ci.exp( glm( cbind(D,Y) ~ G, family=poisreg ), subset="G" ), erf )
}


###################################################
### code chunk number 7: cohort-sol.rnw:133-134
###################################################
sim( t=1, N=c(8000,2000) )


###################################################
### code chunk number 8: cohort-sol.rnw:136-141
###################################################
rbind(
t1=sim( t=1, N=c(8000,2000) ),
t2=sim( t=2, N=c(8000,2000) ),
t3=sim( t=3, N=c(8000,2000) ),
t4=sim( t=4, N=c(8000,2000) ) )


###################################################
### code chunk number 9: cohort-sol.rnw:147-152
###################################################
rbind(
N8.2=sim( t=2, N=c(8000,2000) ),
N7.3=sim( t=2, N=c(7000,3000) ),
N6.4=sim( t=2, N=c(6000,4000) ),
N5.5=sim( t=2, N=c(5000,5000) ) )