### R code from vignette source 'C:/Bendix/teach/NSCE/2022/pracs/SMR-sol.rnw'

###################################################
### code chunk number 1: SMR-sol.rnw:6-17
###################################################
library( Epi )
D1 <- c(11, 15, 10)
Y1 <- c(10, 6, 2)
D0 <- c(15, 60, 150)
Y0 <- c(30, 50, 70)
names(D1) <-
names(Y1) <-
names(D0) <-
names(Y0) <-
c("30-39", "40-49", "50-59")
cbind( D1, Y1, D0, Y0 )


###################################################
### code chunk number 2: SMR-sol.rnw:23-27
###################################################
I1 <- D1/Y1
I0 <- D0/Y0
IR <- I1/I0
round(cbind(I1, I0, IR),2 )


###################################################
### code chunk number 3: SMR-sol.rnw:32-34
###################################################
EF <- exp( 1.96 * sqrt(1/D1+1/D0) )
round( cbind(I1, I0, IR, EF), 2 )


###################################################
### code chunk number 4: SMR-sol.rnw:40-44
###################################################
I1.c <- sum(D1) / sum(Y1)
I0.c <- sum(D0) / sum(Y0)
IR.c <- I1.c / I0.c
round( cbind(I1.c, I0.c, IR.c), 2)


###################################################
### code chunk number 5: SMR-sol.rnw:54-59
###################################################
D <- c(D1,D0)
Y <- c(Y1,Y0)
A <- factor( rep(c('30-39','40-49','50-59'),2) )
G <- factor( rep(c("Wrk","Pop"),each=3), levels=c("Pop","Wrk") )
data.frame( D, Y, A, G )


###################################################
### code chunk number 6: SMR-sol.rnw:62-66
###################################################
mc <- glm( cbind(D,Y) ~     G, family=poisreg )
ma <- glm( cbind(D,Y) ~ A + G, family=poisreg )
round( rbind( ci.exp( mc, subset="G" ),
              ci.exp( ma, subset="G" ) ), 3 )


###################################################
### code chunk number 7: SMR-sol.rnw:81-85
###################################################
Obs <- sum( D1 )
Exp <- sum( I0 * Y1)
SMR <- Obs / Exp
round( cbind(Obs, Exp, SMR), 2)


###################################################
### code chunk number 8: SMR-sol.rnw:92-96
###################################################
I1.s <- sum( Y0*I1 ) / sum( Y0 )
I0.s <- sum( Y0*I0 ) / sum( Y0 )
IR.s <- I1.s / I0.s
round( cbind(I1.s, I0.s, IR.s), 2 )


###################################################
### code chunk number 9: SMR-sol.rnw:102-106
###################################################
I1.x <- sum( Y1*I1 ) / sum( Y1 )
I0.x <- sum( Y1*I0 ) / sum( Y1 )
IR.x <- I1.x / I0.x
round( cbind(I1.x, I0.x, IR.x), 2)


###################################################
### code chunk number 10: SMR-sol.rnw:122-127
###################################################
D <- c(D1,D0)
Y <- c(Y1,Y0)
A <- factor( rep(c('30-39','40-49','50-59'),2) )
G <- factor( rep(c("Wrk","Pop"),each=3), levels=c("Pop","Wrk") )
data.frame( D, Y, A, G )


###################################################
### code chunk number 11: SMR-sol.rnw:132-134
###################################################
mc <- glm( cbind(D,Y) ~ G - 1, family=poisreg )
round( ci.exp( mc ), 2)


###################################################
### code chunk number 12: SMR-sol.rnw:139-141
###################################################
mc <- glm( cbind(D,Y) ~ G, family=poisreg )
round( ci.exp( mc ), 2)


###################################################
### code chunk number 13: SMR-sol.rnw:147-151
###################################################
mi <- glm( cbind(D,Y) ~ A:G -1, family=poisreg )
round( ci.exp( mi ), 2)
mi <- glm( cbind(D,Y) ~ A-1 + A:G, family=poisreg )
round( ci.exp( mi ), 2)


###################################################
### code chunk number 14: SMR-sol.rnw:155-158
###################################################
ms <- glm( cbind(D,Y) ~ A + G, family=poisreg )
ms
round( ci.exp( ms ), 2 )


###################################################
### code chunk number 15: SMR-sol.rnw:165-166
###################################################
anova( ms, mi, test="Chisq" )


