### R code from vignette source 'cacon'
### Encoding: UTF-8

###################################################
### code chunk number 1: cacon.rnw:2-6
###################################################
options( width=90,
#        prompt=" ", continue=" ",
         SweaveHooks=list( fig=function()
         par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n") ) )


###################################################
### code chunk number 2: cacon.rnw:9-10
###################################################
library(Epi)


###################################################
### code chunk number 3: cacon.rnw:29-35
###################################################
library(Epi)
mm <- cbind( c(20,180), c(10,190) )
colnames(mm) <- c("ca","co")
rownames(mm) <- c("red","oth")
mm
twoby2(mm)


###################################################
### code chunk number 4: cacon.rnw:244-249
###################################################
library( Epi )
data( lep )
( cctab <- with( lep, table( bcg, factor(d,labels=c("Co","Ca")) ) ) )
( do <- cctab[,"Ca"]/cctab[,"Co"] )
( OR <- do["yes"]/do["no"] )


###################################################
### code chunk number 5: cacon.rnw:312-315
###################################################
cctab 
cctab[2:1,2:1]
twoby2( cctab[2:1,2:1] )


###################################################
### code chunk number 6: cacon.rnw:336-340
###################################################
effx( response = d,
          type = 'binary',
      exposure = bcg,
          data = lep )


###################################################
### code chunk number 7: cacon.rnw:388-394
###################################################
x <- 0:1
ca <- cctab[,"Ca"]
co <- cctab[,"Co"]
cbind( ca, co, x )
mg <- glm( cbind(ca,co) ~ x, family=binomial )
ci.exp( mg )


###################################################
### code chunk number 8: cacon.rnw:425-427
###################################################
mi <- glm( d ~ bcg, family=binomial, data=lep )
ci.exp( mi )


###################################################
### code chunk number 9: cacon.rnw:536-546
###################################################
set.seed( 1952 ) # to ensure reproducibility
library(e1071)
str( whcc <- matchControls( d ~ age,  data=lep, case=1, cont=0 ) )
mlep <- subset( lep, !is.na(whcc$factor) )
ca <- subset( lep, d==1 )
co <- subset( lep, d==0 )
co <- co[sample(nrow(co),nrow(ca)),]
rlep <- rbind( ca, co )
cbind( with( mlep, addmargins( table(age,d) ) ),
       with( rlep, addmargins( table(age,d) ) ) )


###################################################
### code chunk number 10: or-est
###################################################
mm.0 <- glm( d ~ bcg      , family=binomial, data=mlep )
mm.a <- glm( d ~ bcg + age, family=binomial, data=mlep )
rm.0 <- glm( d ~ bcg      , family=binomial, data=rlep )
rm.a <- glm( d ~ bcg + age, family=binomial, data=rlep )
or <- rbind( ci.exp(mm.0,subset="bcg"),
             ci.exp(mm.a,subset="bcg"),
             ci.exp(rm.0,subset="bcg"),
             ci.exp(rm.a,subset="bcg") )
rownames( or ) <- c("mm.0","mm.a","rm.0","rm.a")
or
plotEst( or, lwd=2, xlog=TRUE, grid=1:10/10, vref=1 )


###################################################
### code chunk number 11: cacon.rnw:619-621
###################################################
anova( mm.a, mm.0, test="Chisq" )
anova( rm.a, rm.0, test="Chisq" )


###################################################
### code chunk number 12: cacon.rnw:780-788
###################################################
library( Epi )
library( survival )
data( bdendo )
head( bdendo, 4 )
ms <- clogit(    d ~ gall+hyp+est+dur+ strata(set), data=bdendo )
mE <- clogistic( d ~ gall+hyp+est+dur, strata=set , data=bdendo )
round( cbind( ci.exp(ms), ci.exp(mE) ), 4 )
with( bdendo, table(est,dur) )


###################################################
### code chunk number 13: cacon.rnw:814-815
###################################################
with( bdendo, table(est,dur) )


###################################################
### code chunk number 14: clogEst
###################################################
with( bdendo, table( interaction(est,dur) ) )
is <- clogit   ( d ~ gall+hyp+factor(interaction(est,dur))+strata(set), 
                 data=bdendo )
iE <- clogistic( d ~ gall+hyp+interaction(est,dur), strata=set, 
                 data=bdendo )
( RR <- ci.exp(iE,subset="est") )
rownames(RR) <- paste( "Estrogen level", 0:4, "vs. none" )
plotEst( RR, xlog=T, grid=TRUE,
         xlab="RR of endometrial cancer relative to no Estrogen" )


###################################################
### code chunk number 15: cacon.rnw:990-1002
###################################################
data( diet )
set.seed( 1952 ) # random numbers will be reproducible 
ncc <- ccwc( entry = doe, 
              exit = dox, 
              fail = chd, 
            origin = dob, 
          controls = 2, 
             match = job,
           include = list(dob,doe,dox,energy), 
              data = diet, 
            silent = FALSE )
head( ncc )


###################################################
### code chunk number 16: cacon.rnw:1030-1035
###################################################
table( tt <- table(ncc$Map) )
who <- names( tt[tt>1] )
with( subset(ncc,Map %in% who), table(Fail,Map) )
( wh.sets <- subset( ncc, Map %in% c(40,51))$Set )
subset( ncc, Set %in% wh.sets )


###################################################
### code chunk number 17: cacon.rnw:1080-1093
###################################################
set.seed( 1952 ) # random numbers will be reproducible 
diet <- transform( diet, aoe = (doe-dob)/365.25,
                         aox = (dox-dob)/365.25 )
acc <- ccwc( entry = aoe, 
              exit = aox, 
              fail = chd, 
            origin = 0, 
          controls = 2, 
             match = job,
           include = list(dob,doe,dox,energy), 
              data = diet, 
            silent = FALSE )
head( acc )


###################################################
### code chunk number 18: cacon.rnw:1188-1203
###################################################
library( survival )
data( nwtco )
str( nwtco )
nwtco <- transform( nwtco,
                   histol = factor(histol,labels=c("FH","UH")),
                    stage = factor(stage,labels=c("I","II","III","IV")),
                      age = age/12 )
ccoh <- subset( nwtco, in.subcohort==1 | rel==1 ) 
c( nrow(nwtco), nrow(ccoh) )
with( ccoh, table(in.subcohort,rel) )
fit.all <- coxph( Surv(edrel,rel) ~ stage + histol + age, data = nwtco )
fit.ccp <-   cch( Surv(edrel,rel) ~ stage + histol + age, data = ccoh,
                  subcoh = ~ in.subcohort, 
                      id = ~ seqno, 
             cohort.size = nrow(nwtco) )


###################################################
### code chunk number 19: ccoh-forest
###################################################
round( cbind( ci.exp(fit.all), ci.exp(fit.ccp) ), 3 )
 plotEst( ci.exp(fit.ccp), xlog=TRUE, lwd=3, grid=c(5:9/10,1:10), 
          vref=1, xlab="Rate ratio" )
linesEst( ci.exp(fit.all), y=5:1-0.2, col=gray(0.5), lwd=3 )