### 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/PopulationSIMD2009V2.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
###################################################
DM  <- read.csv( "../data/dm_data.csv" )
names( DM )
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:94-95
###################################################
subset( DM, doDM < 1935 )


###################################################
### code chunk number 8: DM-Scot-data.rnw:103-108
###################################################
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:119-126
###################################################
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:138-140
###################################################
subset( DM, dob>doDM )
subset( DM, doDM>dod )


###################################################
### code chunk number 11: DM-Scot-data.rnw:147-150
###################################################
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:213-214
###################################################
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-FU.rnw:83-92
###################################################
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 29: DM-Scot-FU.rnw:102-114
###################################################
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 30: DM-Scot-FU.rnw:119-122
###################################################
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 31: DM-Scot-FU.rnw:131-136
###################################################
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 32: 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 33: DM-Scot-inc.rnw:12-21
###################################################
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 )
cbind( round(tt[,,1]), round(tt[,,2],1) )


###################################################
### code chunk number 34: DM-Scot-inc.rnw:27-30
###################################################
Iana <- transform( subset( Atab, A<2011 ),
                   A = A + 0.5,
                   p = P + 0.5 )


###################################################
### code chunk number 35: DM-Scot-inc.rnw:39-46
###################################################
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 36: DM-Scot-inc.rnw:48-49
###################################################
round( cbind( ci.exp( im1 ), ci.exp( if1 ) ), 3 )


###################################################
### code chunk number 37: DM-Scot-inc.rnw:64-70
###################################################
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 38: 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 39: DM-Scot-inc.rnw:90-94
###################################################
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 40: anova
###################################################
anova( imi, im1, test="Chisq" )
anova( ifi, if1, test="Chisq" )


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


###################################################
### code chunk number 42: DM-Scot-inc.rnw:121-126
###################################################
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 43: 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 44: 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 45: DM-Scot-mort.rnw:14-23
###################################################
library( Epi )
library( splines )
load( file="../data/Atab" )
summary( 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 46: DM-Scot-mort.rnw:31-34
###################################################
Mana <- transform( subset( Atab, A<2012 & Y.dm>0 ),
                   A = A + 0.5,
                   p = P + 0.5 )


###################################################
### code chunk number 47: DM-Scot-mort.rnw:38-45
###################################################
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 48: DM-Scot-mort.rnw:47-48
###################################################
round( cbind( ci.exp( mm1 ), ci.exp( mf1 ) ), 3 )


###################################################
### code chunk number 49: DM-Scot-mort.rnw:58-64
###################################################
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 50: inc
###################################################
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 51: DM-Scot-mort.rnw:84-88
###################################################
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 52: anova
###################################################
anova( mmi, mm1, test="Chisq" )
anova( mfi, mf1, test="Chisq" )


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


###################################################
### code chunk number 54: DM-Scot-mort.rnw:115-120
###################################################
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 55: mort-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, mi[,"M",],
         lwd=2:3, col=gray(1:10/14), lty=1, type="l",
         log="y", xlab="Age (years)", ylab=" ", ylim=c(1/100,250) )
matplot( nd$A, mi[,"F",],
         lwd=2:3, col=gray(1:10/14), lty=1, type="l",
         log="y", xlab="Age (years)", ylab=" ", ylim=c(1/100,250) )
mtext( "Mortality rate among DM pateints (per 1000 PY)",
       side=2, line=2.5, las=0, outer=TRUE )


###################################################
### code chunk number 56: mort-int-lin
###################################################
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, mi[,"M",],
         lwd=2:3, col=gray(1:10/14), lty=1, type="l",
         xlab="Age (years)", ylab=" ", ylim=c(0,250) )
matplot( nd$A, mi[,"F",],
         lwd=2:3, col=gray(1:10/14), lty=1, type="l",
         xlab="Age (years)", ylab=" ", ylim=c(0,250) )
mtext( "Mortality rate among DM pateints (per 1000 PY)",
       side=2, line=2.5, las=0, outer=TRUE )