### R code from vignette source 'DM-Scot-s.rnw'

###################################################
### code chunk number 1: DM-Scot-data.rnw:2-8
###################################################
options( width=100,
         prompt=" ", continue=" ",
         SweaveHooks=list( fig=function()
         par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,bty="n",las=1) ) )
library( Epi )
library( splines )


###################################################
### code chunk number 2: DM-Scot-data.rnw:18-20
###################################################
library( Epi )
sessionInfo()


###################################################
### code chunk number 3: DM-Scot-data.rnw:31-40
###################################################
pop <- read.csv( "../data/PopulationSIMD2009.csv" )
names( pop ) <- tolower( names(pop) )
names( pop )
names( pop )[c(1,4:6)] <- c("per","simd","D","N")
pop$sex <- factor( pop$sex, labels=c("M","F") )
str( pop )
head( pop )
tail( pop, 25 )
summary( pop )


###################################################
### code chunk number 4: DM-Scot-data.rnw:47-57
###################################################
with( pop, ftable( per, sex, simd ) )
round( ftable( xtabs( N/1000 ~
                      sex + per + simd,
                      data=pop ) ), 1 )
odd <- function( x ) x[length(x)]/sum(x) * 1000
ftable( xtabs( D ~ sex + simd + per, data=pop ) )
round( ftable( addmargins( xtabs( D ~ sex + simd + per, data=pop ),
                           margin = 3:2,
                           FUN = list( list("2005-12"=sum),
                                       list(sum,"'11' (per 1000)"=odd) ) ) ) )


###################################################
### code chunk number 5: DM-Scot-data.rnw:68-70
###################################################
pop <- subset( pop, simd < 11 )
pop$D <- pmax( 0, pop$D, na.rm=TRUE )


###################################################
### code chunk number 6: readDM
###################################################
system.time(
DM <- read.csv(
     # "http://bendixcarstensen.com/AdvCoh/Scot-2014/data/dm-data.csv" )
       "../data/dm-data.csv" )
            )
str( DM )
for( i in 4:6 ) DM[,i] <- cal.yr( as.Date( DM[,i], format="%d/%m/%Y" ) )
names( DM )[1] <- "simd"
levels( DM$sex ) <- c("F","M")
head( DM )
summary( DM )


###################################################
### code chunk number 7: DM-Scot-data.rnw:97-98
###################################################
subset( DM, doDM < 1935 )


###################################################
### code chunk number 8: DM-Scot-data.rnw:106-111
###################################################
with( DM, ftable( addmargins( table( sex,
                              doDM = floor( pmax(doDM,1999.999) ),
                                     simd,
                                     useNA="ifany"),
                              margin = 2:3 ) ) )


###################################################
### code chunk number 9: DM-Scot-data.rnw:122-129
###################################################
with( DM, ftable( BB= !is.na(dob),
                  DM= !is.na(doDM),
                  DD= !is.na(dod),
              "b<dm"= dob<doDM,
              "dm<d"= doDM<dod,
                 useNA="ifany",
                 row.vars=1:3 ) )


###################################################
### code chunk number 10: DM-Scot-data.rnw:141-143
###################################################
subset( DM, dob>doDM )
subset( DM, doDM>dod )


###################################################
### code chunk number 11: DM-Scot-data.rnw:150-153
###################################################
nrow( DM )
DM <- subset( DM, dob<doDM & pmin(0,dod-doDM,na.rm=T)>=0 )
nrow(DM)


###################################################
### code chunk number 12: hists
###################################################
par( mfrow=c(2,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
hist( DM$dob, breaks=seq(1900,2013,1),
      col="black", main="", xlab="Date of birth" )
abline(v=seq(1900,2010,10),col="red")
hist( DM$doDM, breaks=seq(1915,2013,1),
      col="black", main="", xlab="Date of DM diagnosis" )
abline(v=seq(1920,2010,10),col="red")
hist( DM$doDM[DM$doDM>2000], breaks=seq(2000,2013,1/12),
      col="black", main="", xlab="Date of DM diagnosis" )
abline(v=seq(2000,2013,1),col="red")
abline(v=2005,col="limegreen")
hist( DM$dod, breaks=seq(2005,2013,1/12),
      col="black", main="", xlab="Date of death" )
abline(v=2004:2013,col="red")


###################################################
### code chunk number 13: DM-Scot-data.rnw:216-217
###################################################
save( pop, DM, file="../data/DM.Rda" )


###################################################
### code chunk number 14: DM-Scot-prev.rnw:2-9
###################################################
options( width=100,
#         prompt=" ", continue=" ",
         SweaveHooks=list( fig=function()
         par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,bty="n",las=1) ) )
library( Epi )
library( splines )
load( file="../data/DM.Rda" )


###################################################
### code chunk number 15: DM-Scot-prev.rnw:20-21
###################################################
cal.yr( as.Date("2005-07-01") )


###################################################
### code chunk number 16: DM-Scot-prev.rnw:24-29
###################################################
p2005 <- subset( DM, doDM<2005.5 & (dod>2005.5 | is.na(dod) ) )
p2005 <- with( p2005, as.data.frame( table( sex, simd,
                                            age=floor(2005.5-dob) ) ) )
str( p2005 )
head( p2005 )


###################################################
### code chunk number 17: DM-Scot-prev.rnw:37-54
###################################################
prv <- NULL
for( y in 2005:2011 )
   {
   my <- y + 0.5 # Mid-year date
   sb <- subset( DM, doDM < my & ( dod > my | is.na(dod) ) )
   prv <- rbind( prv,
                 cbind( per = y,
                        with( sb, as.data.frame( table( sex,
                                                        simd,
                                                    age=floor(my-dob) ) ) ) ) )
   }
str( prv )
names( prv )[5] <- "X"
prv <- transform( prv, age = as.numeric( as.character( age  ) ),
                      simd = as.numeric( as.character( simd ) ) )
str( prv )
str( pop )


###################################################
### code chunk number 18: DM-Scot-prev.rnw:60-65
###################################################
prv <- merge( subset( prv, age<90 ),
              subset( pop, age<90 & per<2012 ) )
summary( prv )
save( prv, file="../data/prv.Rda" )
load( file="../data/prv.Rda" )


###################################################
### code chunk number 19: DM-Scot-prev.rnw:85-95
###################################################
library( splines )
a.kn <- seq( 5, 85, 10 )
m0 <- glm( cbind(X,N-X) ~ Ns( age+0.5, knots=a.kn ) +
                            factor( per ) +
                            factor( simd ),
           family = binomial(link="log"),
             data = subset(prv,sex=="M") )
f0 <- update( m0, data = subset(prv,sex=="F") )
round( cbind(  ci.exp( m0 ),
               ci.exp( f0 ) ), 3 )


###################################################
### code chunk number 20: grpests
###################################################
mests <- ci.exp( m0, subset="fact" )
fests <- ci.exp( f0, subset="fact" )
rownames( mests ) <- gsub( "factor\\(per\\)", "", rownames( mests ) )
rownames( mests ) <- gsub( "factor\\(simd\\)", "", rownames( mests ) )
mests <- rbind( 1, mests[1:6,], 1, mests[7:15,] )
fests <- rbind( 1, fests[1:6,], 1, fests[7:15,] )
rownames( mests )[c(1,8)] <- c("Year 2005", "Social class 1")
mests
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
 plotEst( mests, y=17:1+0.1, lwd=3, cex=0.5, xlog=TRUE, vref=1, col="blue",
                 xtic=c(c(5,7)/10,1,1.1,1.3), grid=5:13/10,
                 xlab="Relative prevalence" )
linesEst( fests, y=17:1-0.1, lwd=3, cex=0.5, col="red" )


###################################################
### code chunk number 21: DM-Scot-prev.rnw:130-138
###################################################
prv <- subset( prv, per<2011 )
m0 <- update( m0, data = subset(prv,sex=="M") )
f0 <- update( f0, data = subset(prv,sex=="F") )
ml <- update( m0, . ~ . - factor(simd) - factor(per)
                        + simd + I(per-2008) )
fl <- update( ml, data = subset(prv,sex=="F") )
round( (cbind( ci.exp( ml, subset=c("simd","per") ),
               ci.exp( fl, subset=c("simd","per") ) ) - 1 )*100, 3 )


###################################################
### code chunk number 22: DM-Scot-prev.rnw:153-163
###################################################
mpl <- update( m0, . ~ . +  Ns( age+0.5, knots=a.kn ):I(per-2005) )
mpi <- update( m0, . ~ . +  Ns( age+0.5, knots=a.kn ):factor(per) )
msl <- update( m0, . ~ . +  Ns( age+0.5, knots=a.kn ):I(simd-2005) )
msi <- update( m0, . ~ . +  Ns( age+0.5, knots=a.kn ):factor(simd) )
fpl <- update( f0, . ~ . +  Ns( age+0.5, knots=a.kn ):I(per-2005) )
fpi <- update( f0, . ~ . +  Ns( age+0.5, knots=a.kn ):factor(per) )
fsl <- update( f0, . ~ . +  Ns( age+0.5, knots=a.kn ):I(simd-2005) )
fsi <- update( f0, . ~ . +  Ns( age+0.5, knots=a.kn ):factor(simd) )
anova( mpi, mpl, m0, msl, msi, test="Chisq" )
anova( fpi, fpl, f0, fsl, fsi, test="Chisq" )


###################################################
### code chunk number 23: DM-Scot-prev.rnw:167-171
###################################################
mspl <- update( msl, . ~ . + Ns( age+0.5, knots=a.kn ):I(per-2005) )
fspl <- update( fsl, . ~ . + Ns( age+0.5, knots=a.kn ):I(per-2005) )
anova( mpl, mspl, msl, test="Chisq" )
anova( fpl, fspl, fsl, test="Chisq" )


###################################################
### code chunk number 24: prev2x2
###################################################
nd <- data.frame( age = 0:90,
                  per = 2008,
                 simd = 5 )
mpr2008 <- fpr2008 <- NULL
for( sc in 1:10 )
   {
   mpr2008 <- cbind( mpr2008, ci.pred( mspl, newdata = transform(nd,simd=sc) ) )
   fpr2008 <- cbind( fpr2008, ci.pred( fspl, newdata = transform(nd,simd=sc) ) )
   }
dim( mpr2008 )
mprcl5 <- fprcl5 <- NULL
for( yy in 2005+0:5 )
   {
   mprcl5 <- cbind( mprcl5, ci.pred( mspl, newdata = transform(nd,per=yy) ) )
   fprcl5 <- cbind( fprcl5, ci.pred( fspl, newdata = transform(nd,per=yy) ) )
   }
dim( mprcl5 )
par( mfrow=c(2,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, bty="n" )
matplot( nd$age, mpr2008[,1+0:9*3]*100,
         ylim=c(0,19), xlab="Age",
         ylab="Prevalence at 1.1.2008",
         type="l", lty=1, lwd=3, col=gray(3:12/15) )
matplot( nd$age, fpr2008[,1+0:9*3]*100,
         ylim=c(0,19), xlab="Age",
         ylab="Prevalence at 1.1.2008",
         type="l", lty=1, lwd=3, col=gray(3:12/15) )
matplot( nd$age, mprcl5[,1+0:5*3]*100,
         ylim=c(0,19), xlab="Age",
         ylab="Prevalence in class 5",
         type="l", lty=1, lwd=3, col=gray(8:3/11) )
matplot( nd$age, fprcl5[,1+0:5*3]*100,
         ylim=c(0,19), xlab="Age",
         ylab="Prevalence in class 5",
        type="l", lty=1, lwd=3, col=gray(8:3/11) )


###################################################
### code chunk number 25: DM-Scot-FU.rnw:2-9
###################################################
options( width=100,
#         prompt=" ", continue=" ",
         SweaveHooks=list( fig=function()
         par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,bty="n",las=1) ) )
library( Epi )
library( splines )
load( file="../data/DM.Rda" )


###################################################
### code chunk number 26: DM-Scot-FU.rnw:54-64
###################################################
options( width=120 )
summary( DM )
( dox <- cal.yr( as.Date("2012-05-18") ) )
LD <- Lexis( entry = list( per = pmax( 2005, doDM ),
                           age = pmax( 2005, doDM ) - dob ),
              exit = list( per = pmin( dox, dod, na.rm=TRUE ) ),
       exit.status = pmin( dod, dox, na.rm=TRUE ) < dox,
            states = c("ALive","Dead"),
              data = DM )
summary( LD )


###################################################
### code chunk number 27: DM-Scot-FU.rnw:72-73
###################################################
save( LD, file="../data/LD.Rda" )


###################################################
### code chunk number 28: DM-Scot-timesplit.rnw:3-5
###################################################
library( Epi )
load( file="../data/LD.Rda " )


###################################################
### code chunk number 29: DM-Scot-timesplit.rnw:22-27
###################################################
table( entry(LD,"per")<2011, useNA="ifany" )
LD <- subset( LD, entry(LD,"per")<2011 )
nch <- 20
( ll <- round( seq(0,nrow(LD),,nch+1) ) )
diff( ll )


###################################################
### code chunk number 30: DM-Scot-timesplit.rnw:35-37
###################################################
i <- 1
range( whr <- (ll[i]+1):ll[i+1] )


###################################################
### code chunk number 31: DM-Scot-timesplit.rnw:44-46
###################################################
sl <- splitLexis( LD[whr,], 1990:2015, "per" )
sl <- splitLexis( sl      ,    0:150 , "age" )


###################################################
### code chunk number 32: DM-Scot-timesplit.rnw:51-59
###################################################
agg <- with( sl, aggregate( cbind( y = lex.dur,
                                   d = (lex.Xst=="Dead") ),
                      by = list( sex = sex,
                                   A = floor(age),
                                   P = floor(per),
                                  sC = simd ),
                                 FUN = sum ) )
str( agg )


###################################################
### code chunk number 33: DM-Scot-timesplit.rnw:70-72
###################################################
DMtab <- cbind( agg, Y.dm=0, D.dm=0 )
str( DMtab )


###################################################
### code chunk number 34: DM-Scot-timesplit.rnw:86-92
###################################################
DMtab <- transform( DMtab, Y.dm = pmax( Y.dm, 0, na.rm=TRUE ) +
                                  pmax( y   , 0, na.rm=TRUE ),
                           D.dm = pmax( D.dm, 0, na.rm=TRUE ) +
                                  pmax( d   , 0, na.rm=TRUE ) )[
         c("sex","A","P","sC","Y.dm","D.dm") ]
str( DMtab )


###################################################
### code chunk number 35: DM-Scot-timesplit.rnw:95-119
###################################################
for( i in 2:nch )
{
whr <- (ll[i]+1):ll[i+1]
sl <- splitLexis( LD[whr,], 1990:2015, "per" )
sl <- splitLexis( sl      ,    0:150 , "age" )
agg <- with( sl, aggregate( cbind( y = lex.dur,
                                   d = (lex.Xst=="Dead") ),
                           list( sex = sex,
                                   A = floor(age),
                                   P = floor(per),
                                  sC = simd ),
                                 FUN = sum ) )
DMtab <- merge( DMtab, agg, all=TRUE )
DMtab <- transform( DMtab, Y.dm = pmax( Y.dm, 0, na.rm=TRUE ) +
                                  pmax( y   , 0, na.rm=TRUE ),
                           D.dm = pmax( D.dm, 0, na.rm=TRUE ) +
                                  pmax( d   , 0, na.rm=TRUE ) )[
                           c("sex","A","P","sC","Y.dm","D.dm") ]
cat( "Merged in chunk", i, "now", nrow(DMtab), "rows, at",
     format(Sys.time(),format="%Y-%m-%d %H:%M:%S"), "\n" )
flush.console()
}
str( DMtab )
save( DMtab, file="../data/DMtab.Rda" )


###################################################
### code chunk number 36: DM-Scot-FU.rnw:84-93
###################################################
load( file="../data/DMtab.Rda" )
str( DMtab )
str( pop )
Atab <- merge( subset( DMtab, A<90 ),
               subset( pop  , age<90 & simd<11 ),
               by.x = c("sex","sC"  ,"A"  ,"P"  ),
               by.y = c("sex","simd","age","per"),
               all = TRUE )
summary( Atab )


###################################################
### code chunk number 37: DM-Scot-FU.rnw:103-115
###################################################
head( LD )
DMinc <- with( subset(LD,entry(LD,"per")<2011),
               aggregate( !is.na(doDM),
                          list( A = floor(age),
                                P = floor(per),
                              sex = sex,
                               sC = simd ),
                              FUN = sum ) )
names( DMinc )[5] <- "I.dm"
str( DMinc )
Atab <- merge( Atab, subset( DMinc, A < 90 & P < 2012 ), all=TRUE )
head( Atab )


###################################################
### code chunk number 38: DM-Scot-FU.rnw:120-123
###################################################
Atab <- transform( Atab, Y.dm = pmax( 0, Y.dm, na.rm=TRUE ),
                         D.dm = pmax( 0, D.dm, na.rm=TRUE ),
                         I.dm = pmax( 0, I.dm, na.rm=TRUE ) )


###################################################
### code chunk number 39: DM-Scot-FU.rnw:132-137
###################################################
Atab <- transform( Atab, Y.nd = pmax( 0, N-Y.dm, na.rm=TRUE ),
                         D.nd = pmax( 0, D-D.dm, na.rm=TRUE ) )
str( Atab )
summary( Atab )
save( Atab, file="../data/Atab" )


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


###################################################
### code chunk number 41: DM-Scot-inc.rnw:12-22
###################################################
library( Epi )
library( splines )
load( file="../data/Atab" )
summary( Atab )
tt <- addmargins( xtabs( cbind(I.dm,Y.nd/1000) ~ A + P,
                         data=Atab ),
                  margin = 1 )
str( tt )
      tt[,,1]
round(tt[,,2],1)


###################################################
### code chunk number 42: DM-Scot-inc.rnw:28-31
###################################################
Iana <- transform( subset( Atab, A<2011 ),
                   A = A + 0.5,
                   P = P + 0.5 )


###################################################
### code chunk number 43: DM-Scot-inc.rnw:40-47
###################################################
a.kn <- seq(5,85,,10)
p.kn <- c(2006.5,2008,2009.5)
im1 <- glm( I.dm ~ Ns(A,kn=a.kn) + Ns(P,kn=p.kn) + factor(sC),
                   offset = log(Y.nd),
                   family = poisson,
                     data = subset(Iana,sex=="M") )
if1 <- update( im1,  data = subset(Iana,sex=="F") )


###################################################
### code chunk number 44: DM-Scot-inc.rnw:49-50
###################################################
round( cbind( ci.exp( im1 ), ci.exp( if1 ) ), 3 )


###################################################
### code chunk number 45: DM-Scot-inc.rnw:65-71
###################################################
nd <- data.frame( A = 0:90,
                  P = 2008,
                 sC = 5,
               Y.nd = 1000 )
minc2008 <- ci.pred( im1, newdata = nd )
finc2008 <- ci.pred( if1, newdata = nd )


###################################################
### code chunk number 46: inc
###################################################
par( mar=c(3,4,1,1) )
matplot( nd$A, cbind( minc2008,
                      finc2008 ),
         lwd=c(3,1,1), col=rep(c("blue","red"),each=3), lty=1, type="l",
         log="y", xlab="Age (years)", ylab=" ")
mtext( "DM incidence rate (per 1000 PY)", side=2, line=2.5, las=0 )


###################################################
### code chunk number 47: DM-Scot-inc.rnw:93-97
###################################################
r.kn <- seq(2,88,,4)
imi <- update( im1, . ~ . + factor(sC):Ns(A,knots=r.kn) )
ifi <- update( if1, . ~ . + factor(sC):Ns(A,knots=r.kn) )
summary( ifi )


###################################################
### code chunk number 48: anova
###################################################
anova( imi, im1, test="Chisq" )
anova( ifi, if1, test="Chisq" )


###################################################
### code chunk number 49: DM-Scot-inc.rnw:117-121
###################################################
ii <- NArray( list( A = nd$A,
                  sex = c("M","F"),
                   sC = 1:10 ) )
str( ii )


###################################################
### code chunk number 50: DM-Scot-inc.rnw:124-129
###################################################
for( sc in 1:10 )
   {
 ii[,"M",sc] <- ci.pred( imi, newdata = transform( nd, sC=sc ) )[,1]
 ii[,"F",sc] <- ci.pred( ifi, newdata = transform( nd, sC=sc ) )[,1]
   }


###################################################
### code chunk number 51: inc-int
###################################################
par( mfrow=c(1,2), mar=c(3,1,1,1), oma=c(0,4,0,0), mgp=c(3,1,0)/1.6,
     las=1, bty="n" )
matplot( nd$A, ii[,"M",],
         lwd=2:3, col=gray(1:10/14), lty=1, type="l",
         log="y", xlab="Age (years)", ylab=" ", ylim=c(5/199,10) )
matplot( nd$A, ii[,"F",],
         lwd=2:3, col=gray(1:10/14), lty=1, type="l",
         log="y", xlab="Age (years)", ylab=" ", ylim=c(5/199,10) )
mtext( "DM incidence rate (per 1000 PY)", side=2, line=2.5, las=0, outer=TRUE )


###################################################
### code chunk number 52: DM-Scot-mort.rnw:2-8
###################################################
options( width=100,
#         prompt=" ", continue=" ",
         SweaveHooks=list( fig=function()
         par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n") ) )
library( Epi )
library( splines )


###################################################
### code chunk number 53: DM-Scot-mort.rnw:14-21
###################################################
load( file="../data/Atab" )
str( Atab )
tt <- addmargins( xtabs( cbind(D.dm,Y.dm/1000) ~ A + P,
                         data=Atab ),
                  margin = 1 )
str( tt )
cbind( round(tt[,,1]), round(tt[,,2],1) )


###################################################
### code chunk number 54: DM-Scot-mort.rnw:29-32
###################################################
Mana <- transform( subset( Atab, A<2012 & Y.dm>0 ),
                   A = A + 0.5,
                   p = P + 0.5 )


###################################################
### code chunk number 55: DM-Scot-mort.rnw:36-43
###################################################
a.kn <- c(10,20,40,seq(50,85,,5))
p.kn <- c(2006.5,2008.5,2010.5)
mm1 <- glm( D.dm ~ Ns(A,kn=a.kn) + Ns(P,kn=p.kn) + factor(sC),
                   offset = log(Y.dm),
                   family = poisson,
                     data = subset(Mana,sex=="M") )
mf1 <- update( mm1,  data = subset(Mana,sex=="F") )


###################################################
### code chunk number 56: DM-Scot-mort.rnw:45-46
###################################################
round( cbind( ci.exp( mm1 ), ci.exp( mf1 ) ), 3 )


###################################################
### code chunk number 57: DM-Scot-mort.rnw:56-62
###################################################
nd <- data.frame( A = 0:90,
                  P = 2008,
                 sC = 5,
               Y.dm = 1000 )
mmort2008 <- ci.pred( mm1, newdata = nd )
fmort2008 <- ci.pred( mf1, newdata = nd )


###################################################
### code chunk number 58: mort-MF
###################################################
par( mar=c(3,4,1,1) )
matplot( nd$A, cbind( mmort2008,
                      fmort2008 ),
         lwd=c(3,1,1), col=rep(c("blue","red"),each=3), lty=1, type="l",
         log="y", xlab="Age (years)", ylab=" ",ylim=c(1,300))
mtext( "Mortalty rate in DM patients (per 1000 PY)", side=2, line=2.5, las=0 )


###################################################
### code chunk number 59: DM-Scot-mort.rnw:92-96
###################################################
( r.kn <- seq(30,85,,4) )
mmi <- update( mm1, . ~ . + factor(sC):Ns(A,knots=r.kn) )
mfi <- update( mf1, . ~ . + factor(sC):Ns(A,knots=r.kn) )
summary( mfi )


###################################################
### code chunk number 60: anova
###################################################
anova( mm1, mmi, test="Chisq" )
anova( mf1, mfi, test="Chisq" )


###################################################
### code chunk number 61: DM-Scot-mort.rnw:116-120
###################################################
mi <- NArray( list( A = nd$A,
                  sex = c("M","F"),
                   sC = 1:10 ) )
str( mi )


###################################################
### code chunk number 62: DM-Scot-mort.rnw:124-129
###################################################
for( sc in 1:10 )
   {
 mi[,"M",sc] <- ci.pred( mmi, newdata = transform( nd, sC=sc ) )[,1]
 mi[,"F",sc] <- ci.pred( mfi, newdata = transform( nd, sC=sc ) )[,1]
   }


###################################################
### code chunk number 63: mort-int
###################################################
par( mfcol=c(2,2), mar=c(0,0,0,0), oma=c(4,5,0,0), mgp=c(3,1,0)/1.6,
     las=1, bty="n" )
for( sx in c("M","F") )
for( ax in c("y","") )
   {
plot( NA, NA, log=ax,
      xaxt=if( ax=="y" ) "n" else "s",
      yaxt=if( sx=="F" ) "n" else "s",
      xlab="", ylab="", xlim=c(0,90), ylim=c(1/500,250) )
abline( h= if( ax=="y" ) outer(c(1:9),-3:3,function(x,y) x*10^y)
                    else seq(0,260,10),
        v=seq(0,90,10), col=gray(0.9) )
matlines( nd$A, mi[,sx,],
          lwd=2:3, col=gray(1:10/14),
          lty=1, type="l" )
if( ax=="y" ) text( 15, sqrt(20*50), sx, font=2, cex=1.2, adj=c(0.5,0.5) )
   }
mtext( "Age at follow-up", side=1, line=2.5, las=0, outer=TRUE )
mtext( "Mortality rate among DM patients (per 1000 PY)",
       side=2, line=3.0, las=0, outer=TRUE )


###################################################
### code chunk number 64: mosaic
###################################################
dd <- xtabs( D.dm ~ A + sC + sex, data=Atab )
str(dd)
par(mfrow=c(1,2), mar=c(0,0,1,0))
mosaicplot(dd[,,"M"],col=heat.colors(10),off=0,border="transparent",main="")
text( 0.5, 0.95, "M", col="white", adj=c(0.5,1), font=2, cex=2 )
mosaicplot(dd[,,"F"],col=heat.colors(10),off=0,border="transparent",main="")
text( 0.5, 0.95, "F", col="white", adj=c(0.5,1), font=2, cex=2 )


###################################################
### code chunk number 65: DM-Scot-smr.rnw:2-8
###################################################
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") ) )
library( Epi )
library( splines )


###################################################
### code chunk number 66: DM-Scot-smr.rnw:15-17
###################################################
load( file="../data/Atab" )
str( Atab )


###################################################
### code chunk number 67: DM-Scot-smr.rnw:24-29
###################################################
Sana <- transform( subset( Atab, A<2012 & Y.dm>0 ),
                   A = A + 0.5,
                   P = P + 0.5,
                   E = Y.dm * (D.nd/Y.nd) )
Sana <- subset( Sana, E>0 )


###################################################
### code chunk number 68: DM-Scot-smr.rnw:38-45
###################################################
( a.kn <- c(10,30,seq(50,85,,5)) )
p.kn <- c(2006.5,2008.5,2010.5)
sm1 <- glm( D.dm ~ Ns(A,kn=a.kn) + Ns(P,kn=p.kn) + factor(sC),
                   offset = log(E),
                   family = poisson,
                     data = subset(Sana,sex=="M") )
sf1 <- update( sm1,  data = subset(Sana,sex=="F") )


###################################################
### code chunk number 69: DM-Scot-smr.rnw:47-48
###################################################
round( cbind( ci.exp( sm1 ), ci.exp( sf1 ) ), 3 )


###################################################
### code chunk number 70: DM-Scot-smr.rnw:60-66
###################################################
nd <- data.frame( A = 10:90,
                  P = 2008,
                 sC = 5,
                  E = 1 )
msmr2008 <- ci.pred( sm1, newdata = nd )
fsmr2008 <- ci.pred( sf1, newdata = nd )


###################################################
### code chunk number 71: SMR-MF
###################################################
par( mar=c(3,4,1,1) )
matplot( nd$A, cbind( msmr2008,
                      fsmr2008 ),
         lwd=c(3,1,1), col=rep(c("blue","red"),each=3), lty=1, type="l",
         log="y", xlab="Age (years)", ylab=" ",ylim=c(0.8,20) )
mtext( "SMR of DM patients relative to the non-DM population", side=2, line=2.5, las=0 )


###################################################
### code chunk number 72: DM-Scot-smr.rnw:91-95
###################################################
( r.kn <- seq(30,85,,4) )
smi <- update( sm1, . ~ . + factor(sC):Ns(A,knots=r.kn) )
sfi <- update( sf1, . ~ . + factor(sC):Ns(A,knots=r.kn) )
ci.exp( sfi )


###################################################
### code chunk number 73: anova
###################################################
anova( sm1, smi, test="Chisq" )
anova( sf1, sfi, test="Chisq" )


###################################################
### code chunk number 74: DM-Scot-smr.rnw:115-120
###################################################
si <- NArray( list( A = nd$A,
                  sex = c("M","F"),
                   sC = 1:10,
                   CI = c("SMR","lo","hi") ) )
str( si )


###################################################
### code chunk number 75: DM-Scot-smr.rnw:124-129
###################################################
for( sc in 1:10 )
   {
 si[,"M",sc,] <- ci.pred( smi, newdata = transform( nd, sC=sc ) )
 si[,"F",sc,] <- ci.pred( sfi, newdata = transform( nd, sC=sc ) )
   }


###################################################
### code chunk number 76: SMR-int
###################################################
par( mfcol=c(1,2), mar=c(0,0,0,0), oma=c(4,5,0,0) )
#   , mgp=c(3,1,0)/1.6, las=1, bty="n" )
for( sx in c("M","F") )
   {
plot( NA, NA, log="y",
      yaxt=if( sx=="F" ) "n" else "s",
      xlab="", ylab="", xlim=c(30,90), ylim=c(0.8,20) )
abline( h=outer(c(1:9),-3:3,function(x,y) x*10^y), v=seq(0,90,10), col=gray(0.9) )
matlines( nd$A, si[,sx,,1],
          lwd=2:1*2+1, col=gray(10:1/14),
          lty=1, type="l" )
abline( h=1 )
text( 65, 4.5, sx, font=2, cex=1.2, adj=c(0.5,0.5) )
   }
mtext( "Age at follow-up", side=1, line=2.5, las=0, outer=TRUE )
mtext( "SMR among DM patients vs. non-DM population",
       side=2, line=3.0, las=0, outer=TRUE )


###################################################
### code chunk number 77: DM-Scot-smr.rnw:200-208
###################################################
tt <- xtabs( D.nd ~ floor(A/5) + sex + sC + P, data=Atab )
str( tt )
data.frame( gtN = 1:5*100,
            pct = round( rbind( mean( tt>100 ),
                                mean( tt>200 ),
                                mean( tt>300 ),
                                mean( tt>400 ),
                                mean( tt>500 ) )*100, 1 ) )


###################################################
### code chunk number 78: DM-Scot-smr.rnw:221-236
###################################################
Rdm <- transform( subset( Atab, A<2012 & Y.dm>0 ),
                  A = A + 0.5,
                  P = P + 0.5,
                 DM = "DM",
                  D = D.dm,
                  Y = Y.dm )
Rnd <- transform( subset( Atab, A<2012 & Y.nd>0 ),
                  A = A + 0.5,
                  P = P + 0.5,
                 DM = "nD",
                  D = D.nd,
                  Y = Y.nd )
Rana <- rbind( Rnd, Rdm )[,c("sex","A","P","sC","DM","D","Y")]
str( Rana )
head( Rana )


###################################################
### code chunk number 79: DM-Scot-smr.rnw:273-280 (eval = FALSE)
###################################################
## rmx <- glm( D ~ Ns(A,kn=a.kn) +
##                               factor(sC) +
##                 Ns(A,kn=r.kn):factor(sC) +
##                 Ns(P,kn=p.kn):factor(sC),
##                 offset = log(Y),
##                 family = poisson,
##                   data = subset(Rana,sex=="M") )


###################################################
### code chunk number 80: RR0
###################################################
rm0 <- glm( D ~ DM + Ns(A,kn=a.kn) +
                                   factor(sC) +
                     Ns(A,kn=r.kn):factor(sC) +
                     Ns(P,kn=p.kn):factor(sC),
                  offset = log(Y/1000),
                  family = poisson,
                    data = subset(Rana,sex=="M") )
rf0 <- update( rm0, data = subset(Rana,sex=="F") )
round( ci.exp(rm0), 3 )


###################################################
### code chunk number 81: DM-Scot-smr.rnw:300-302
###################################################
round( rbind( ci.exp(rm0,subset="DM"),
              ci.exp(rf0,subset="DM") ), 3 )


###################################################
### code chunk number 82: RR1
###################################################
rm1 <- update( rm0, . ~ . + DM:Ns(A,kn=a.kn) + DM:Ns(P,kn=p.kn) )
rf1 <- update( rm1, data = subset(Rana,sex=="F") )


###################################################
### code chunk number 83: Cmat1
###################################################
round( cbind( ci.exp( rm1, subset="DMDM" ),
              ci.exp( rf1, subset="DMDM" ) ), 3 )
CA <- Ns( nd$A, knots=a.kn )
Cp <- Ns( nd$P, knots=p.kn )
m1.rr <- ci.exp( rm1, subset="DMDM", ctr.mat=cbind(1,CA,Cp) )
f1.rr <- ci.exp( rf1, subset="DMDM", ctr.mat=cbind(1,CA,Cp) )


###################################################
### code chunk number 84: RR1-MF
###################################################
par( mar=c(3,4,1,1) )
matplot( nd$A, cbind( m1.rr,
                      f1.rr ),
         lwd=c(3,1,1), col=rep(c("blue","red"),each=3), lty=1, type="l",
         log="y", xlab="Age (years)", ylab=" ",ylim=c(0.8,20) )
matlines( nd$A, cbind( msmr2008, fsmr2008 ),
         lwd=c(3,1,1), col=rep(c("blue","red"),each=3), lty=3, type="l" )
mtext( "RR of DM patients relative to the non-DM population",
       side=2, line=2.5, las=0 )


###################################################
### code chunk number 85: RRi
###################################################
rmi <- update( rm1, . ~ . + factor(sC):DM
                          + DM:Ns(A,kn=r.kn):factor(sC) )
rfi <- update( rmi, data = subset(Rana,sex=="F") )
round( ci.exp( rmi, subset="DM" ), 3 )


###################################################
### code chunk number 86: DM-Scot-smr.rnw:386-391
###################################################
ri <- NArray( list( A = nd$A,
                  sex = c("M","F"),
                   sC = 1:10,
                   CI = c("RR","lo","hi") ) )
str( ri )


###################################################
### code chunk number 87: DM-Scot-smr.rnw:401-416
###################################################
a.pt <- nd$A
N <- length( a.pt )
p.rf <- 2008
CA <- Ns(     a.pt   , knots=a.kn )
RA <- Ns(     a.pt   , knots=r.kn )
Cp <- Ns( rep(p.rf,N), knots=p.kn )
Rx <- RA[,rep(c(1,1:ncol(RA)),each=10)]*0
for( sc in 1:10 )
   {
   Rx[,sc+(1:(ncol(RA)+1)-1)*10] <- cbind(1,RA)
   CM <- cbind(1,CA,Cp,Rx[,-1])
   ri[,"M",sc,] <- ci.exp( rmi, subset="DMDM", ctr.mat=CM )
   ri[,"F",sc,] <- ci.exp( rfi, subset="DMDM", ctr.mat=CM )
   Rx <- Rx*0
   }


###################################################
### code chunk number 88: RR-int
###################################################
par( mfcol=c(1,2), mar=c(0,0,0,0), oma=c(4,5,0,0) )
# , mgp=c(3,1,0)/1.6, las=1, bty="n" )
for( sx in c("M","F") )
   {
plot( NA, NA, log="y",
      yaxt=if( sx=="F" ) "n" else "s",
      xlab="", ylab="", xlim=c(20,90), ylim=c(0.8,20) )
abline( h=outer(c(1:19/2),-3:3,function(x,y) x*10^y),
        v=seq(0,90,10), col=gray(0.9) )
matlines( nd$A, ri[,sx,,1],
          lwd=2:1*2+1, col=gray(10:1/14),
          lty=1, type="l" )
abline( h=1 )
text( 65, 4.5, sx, font=2, cex=1.2, adj=c(0.5,0.5) )
   }
mtext( "Age at follow-up", side=1, line=2.5, las=0, outer=TRUE )
mtext( "Relative mortality among DM patients vs. non-DM population",
       side=2, line=3.0, las=0, outer=TRUE )


###################################################
### code chunk number 89: RR-int-ci
###################################################
par( mfcol=c(1,2), mar=c(0,0,0,0), oma=c(4,5,0,0) )
# , mgp=c(3,1,0)/1.6, las=1, bty="n" )
for( sx in c("M","F") )
   {
plot( NA, NA, log="y",
      yaxt=if( sx=="F" ) "n" else "s",
      xlab="", ylab="", xlim=c(20,90), ylim=c(0.8,20) )
abline( h=outer(c(1:19/2),-3:3,function(x,y) x*10^y),
        v=seq(0,90,10), col=gray(0.9) )
matlines( nd$A, cbind( ri[,sx,,1], ri[,sx,,2], ri[,sx,,3] ),
          lwd=c(rep(2:1*2+1,5),rep(1,20)), col=gray(10:1/14),
          lty=1, type="l" )
abline( h=1 )
text( 75, 6.25, sx, font=2, cex=1.2, adj=c(0.5,0.5) )
   }
mtext( "Age at follow-up", side=1, line=2.5, las=0, outer=TRUE )
mtext( "Relative mortality among DM patients vs. non-DM population",
       side=2, line=3.0, las=0, outer=TRUE )