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

###################################################
### code chunk number 1: case-con-MI-sol.rnw:7-11
###################################################
D1 <- c(141, 49)
D0 <- c(144, 32)
C1 <- c(208, 58)
C0 <- c(112, 45)


###################################################
### code chunk number 2: case-con-MI-sol.rnw:14-18
###################################################
names(D1) <-
names(D0) <-
names(C1) <-
names(C0) <- c("M","F")


###################################################
### code chunk number 3: case-con-MI-sol.rnw:23-27
###################################################
EOR <- (D1/D0)/(C1/C0)
SE.leor <- sqrt(1/D1 + 1/D0 + 1/C1 + 1/C0)
EOR.95low <- EOR / exp(1.96*SE.leor)
EOR.95up  <- EOR * exp(1.96*SE.leor)


###################################################
### code chunk number 4: case-con-MI-sol.rnw:32-34
###################################################
strata <- cbind(EOR, SE.leor, EOR.95low, EOR.95up)
round( strata, 3 )


###################################################
### code chunk number 5: case-con-MI-sol.rnw:52-57
###################################################
EOR.ratio <- EOR[1] / EOR[2]
V.logEOR.ratio <- SE.leor[1]^2 + SE.leor[2]^2 
Wald <- log(EOR.ratio)^2 / V.logEOR.ratio
P.Wald <- 1 - pchisq(Wald, df=1)
round( cbind(Wald,P.Wald), 4 )


###################################################
### code chunk number 6: case-con-MI-sol.rnw:72-77
###################################################
EOR.crude <- (sum(D1)/sum(D0)) / (sum(C1)/sum(C0))
SE.lc <- sqrt( 1/sum(D1) + 1/sum(D0) + 1/sum(C1) + 1/sum(C0) )
EOR.c95low <- EOR.crude / exp(1.96*SE.lc )
EOR.c95up <- EOR.crude * exp(1.96*SE.lc )
cbind(EOR.crude, EOR.c95low, EOR.c95up)


###################################################
### code chunk number 7: case-con-MI-sol.rnw:132-138
###################################################
library(Epi)
y <- cbind( D=c(D0,D1), C=c(C0,C1) )
sex <- factor( rep(c("M","F"),2) )
phys <- factor( rep(c("N","Y"),each=2) )
data.frame( y, sex, phys ) 
cbind( y, sex, phys ) 


###################################################
### code chunk number 8: case-con-MI-sol.rnw:146-149
###################################################
y
mamod <- glm( y ~ sex + phys, family=binomial )
summary( mamod )


###################################################
### code chunk number 9: case-con-MI-sol.rnw:154-157
###################################################
mcmod <- glm( y ~ phys, family=binomial )
round( rbind( ci.exp(mamod,subset="phys"),
              ci.exp(mcmod,subset="phys") ), 3 )


###################################################
### code chunk number 10: case-con-MI-sol.rnw:164-166
###################################################
mimod <- glm( y ~ sex + sex:phys, family=binomial )
summary( mimod )


###################################################
### code chunk number 11: case-con-MI-sol.rnw:174-176
###################################################
ci.exp( mimod )
round( ci.exp( mimod, subset="phys" ), 3 )


###################################################
### code chunk number 12: case-con-MI-sol.rnw:180-183
###################################################
mcmod <- glm( y ~ phys, family=binomial )
summary( mcmod )
round( ci.exp( mcmod, subset="phys" ), 3 )


###################################################
### code chunk number 13: case-con-MI-sol.rnw:191-192
###################################################
anova( mimod, mamod, mcmod, test="Chisq" )