### R code from vignette source 'CIMT.rnw'

###################################################
### code chunk number 1: pr-met.rnw:10-14
###################################################
library( lme4 )
library( Epi )
load( file="./data/AD.Rda" )
AD$subjid <- factor( AD$subjid )


###################################################
### code chunk number 2: mod sep base
###################################################
m0 <- lmer( fimtavg ~ -1  + visit + visit:grp + over.65 + pre.ins + sdc +(1|subjid),
            data = AD )
summary( m0 )
round( ee <- ci.lin( m0, subset="visit" ), 4 )
C0 <- rbind(rbind(cbind(diag(2),0,0),
                  cbind(diag(2),diag(2))),
                    c(-1,1,-1,1),
                    c(-1,1, 0,0),
                    c( 0,0,-1,1) )
row.names(C0) <- c( paste( gN[c(1,1,2,2,2,1)],
                           c(c("v1","v7")[c(1,2,1,2)],
                             rep("18m - baseline",2) ) ),
                    paste( gN[2], "vs.", gN[1] ) )
colnames(C0) <- row.names(ee)
C0
e0 <- ci.lin( m0, subset="visit", ctr.mat=C0 )
round( e0, 4 )


###################################################
### code chunk number 3: fu-code
###################################################
mm <- model.matrix( ~ visit:grp-1, data = AD )
m7 <- mm[,xx <- grep("v7",colnames(mm))]
head( m7 )


###################################################
### code chunk number 4: mod ident base
###################################################
m1 <- lmer( fimtavg ~ m7 + over.65 + pre.ins + sdc + (1|subjid),
            data = AD )
summary( m1 )
round( ee <- ci.lin( m1, subset=1:3 ), 4 )
C1 <- rbind( c(1,0,0),
             c(1,1,0),
             c(1,0,0),
             c(1,0,1),
             c(0,0,1),
             c(0,1,0),
             c(0,-1,1) )
colnames(C1) <- row.names(ee)
C1
round( e1 <- ci.lin( m1, subset=1:3, ctr.mat=C1 ), 4 )
rownames( e1 ) <- rownames( e0 )
round( e1, 4 )


###################################################
### code chunk number 5: cmp base
###################################################
round( cbind( e0[,1:2], e1[,1:2] ), 4 )


###################################################
### code chunk number 6: pr-met.rnw:144-147
###################################################
source( "cnr.R" )
cnr
ppclr <- rgb( t(col2rgb(clr)+255)/2, max=255 )


###################################################
### code chunk number 7: effects
###################################################
par( mar=c(3,1,1,1), mgp=c(3,1,0)/1.6 )
plotEst( e0[5:7,c(1,5,6)], lwd=4, cex=1.5,
         xlab="", col=ppclr, vref=0, y=3:1+0.1,
         restore.par=FALSE )
linesEst( e1[5:7,c(1,5,6)], lwd=4, lty=3, cex=1.5,
          xlab="", col=clr, vref=0, y=3:1-0.1,
          restore.par=FALSE )
axis( side=1 )
mtext( "Mean carotid IMT (mm)", side=1, line=3/1.6, col=clr[4] )


###################################################
### code chunk number 8: forestwin
###################################################
win.metafile( "./results/forest1.emf", width=10, height=5, pointsize=24 )
par( mar=c(3,1,1,1), mgp=c(3,1,0)/1.6,
     bg=#"black",
        "transparent",
     col.axis=clr[4], col.lab=clr[4] )
plotEst( e0[5:7,c(1,5,6)], lwd=7,cex=1.5,
         xlab="", col.txt="transparent", col=clr,
         restore.par=FALSE )
for( i in 1:3 ) axis( side=2, at=4-i, labels=rownames(e0)[4+i], col="transparent",
                      col.axis= clr[i], las=1 )
abline( v=0, col=clr[4] )
axis( side=1, col=clr[4] )
mtext( "Carotid intima-media thickness (mm)", side=1, line=3/1.6, col=clr[4] )
dev.off()


###################################################
### code chunk number 9: reshape
###################################################
names( AD )
wimt <- reshape( AD[AD$visit %in% c("v1","v7"),
                    c("fimtavg","visit","subjid","grp","over.65","pre.ins","sdc")],
                 direction = "wide",
                   v.names = "fimtavg",
                   timevar = "visit",
                     idvar = "subjid" )
subset( AD[,c("fimtavg","visit","subjid","grp","over.65","pre.ins","sdc")],
        subjid=="10002" )
subset( wimt, subjid=="10002" )


###################################################
### code chunk number 10: lm7
###################################################
mf <- lm( fimtavg.v7 ~ fimtavg.v1 + grp + over.65 + pre.ins + sdc,
          data = wimt )
summary( mf )
mm <- lm( fimtavg.v7 ~ fimtavg.v1 + grp,
          data = wimt )
summary( mm )
round( cf <- ci.lin( mf, subset=c("fimt","Met") ), 4 )
round( cm <- ci.lin( mm, subset=c("fimt","Met") ), 4 )


###################################################
### code chunk number 11: pr-met.rnw:284-288
###################################################
round( rbind( "Base-simple" = cm[2,],
              "Base-contrl" = cf[2,],
               "RanEff-dif" = e0[7,],
               "RanEff-eql" = e1[7,] ), 4 )


###################################################
### code chunk number 12: effects-x
###################################################
par( mar=c(3,1,1,1), mgp=c(3,1,0)/1.6 )
plotEst( e0[5:7,c(1,5,6)], lwd=4, cex=1.5,
         xlab="", col=ppclr, vref=0, y=3:1-0.1, ylim=c(0,3.5),
         restore.par=FALSE )
linesEst( e1[5:7,c(1,5,6)], lwd=4, lty=3, cex=1.5,
          xlab="", col=clr, vref=0, y=3:1+0.1 )
linesEst( cf[2,c(1,5,6),drop=FALSE], lwd=4, lty=3, cex=1.5,
          xlab="", col=clr[4], vref=0, y=1+0.3 )
axis( side=1 )
mtext( "Mean carotid IMT (mm)", side=1, line=3/1.6, col=clr[4] )


###################################################
### code chunk number 13: imp
###################################################
library( mice )
load( file="./data/imp" )
class( imp )
mf <- with( imp, lm( fimtavg.v7 ~ fimtavg.v1 + grp + over.65 + pre.ins + sdc ) )
round( summary( pool( mf) ), 4 )
round( smf <- ci.lin( pool( mf), subset="grp" ), 4 )
rownames( smf ) <- paste( gN[2],"vs.",gN[1] )


###################################################
### code chunk number 14: forest-met-imp
###################################################
round( rbind( e1[5:6,c(1,5,6)],
             smf[   ,c(1,5,6)]), 3 )
tmpl <-
function()
{
par( mar=c(3,1,1,1), mgp=c(3,1,0)/1.6 )
plotEst( rbind(e1[5:6,c(1,5,6)],smf[,c(1,5,6)]),
         ylim=c(0,3.5),
         txt=rownames(e0)[5:7], lwd=3, cex=1.5, col=clr,
         restore.par=FALSE )
abline( v=0 )
axis( side=1 )
text( -0.03, 0.2, "Improvement"  , adj=0, col=clr[4], cex=0.8 )
text(  0.03, 0.2, "Deterioration", adj=1, col=clr[4], cex=0.8 )
mtext( "Mean carotid IMT (mm)", side=1, line=3/1.6 )
}
tmpl()


###################################################
### code chunk number 15: pr-met.rnw:357-360
###################################################
win.metafile( "./results/forest-met-imp.emf", width=10, height=6, pointsize=20 )
tmpl()
dev.off()


###################################################
### code chunk number 16: pr-ins.rnw:10-13
###################################################
library( Epi )
library( lme4 )
load( file="./data/AD.Rda" )


###################################################
### code chunk number 17: mod sep base ins
###################################################
m0 <- lmer( fimtavg ~ -1  + visit + visit:igr + over.65 + pre.ins + sdc +(1|subjid),
            data = AD )
summary( m0 )
round( ee <- ci.lin( m0, subset="visit" ), 4 )
C0 <- rbind( cbind(diag(2),    0,0,   0,0 ),
             cbind(diag(2),diag(2),   0,0 ),
             cbind(diag(2),    0,0,diag(2)) )
C0 <- rbind( C0, C0[2,]-C0[1,],
                 C0[4,]-C0[3,],
                 C0[6,]-C0[5,] )
C0 <- rbind( C0, C0[8,]-C0[7,],
                 C0[9,]-C0[7,],
                 C0[9,]-C0[8,] )
row.names(C0) <- c(
as.vector( outer( c("v1","v7"), levels(AD$igr), function(x,y) paste(y,x) ) ),
paste( levels(AD$igr), "v7-v1" ),
paste(   "d", levels(AD$igr)[c(2,3,3)],
       "- d", levels(AD$igr)[c(1,1,2)] ) )
colnames(C0) <- row.names(ee)
C0
e0 <- ci.lin( m0, subset="visit", ctr.mat=C0 )
inames <- c("Biphasic","Asp+Det","Detemir")
cbind( levels( AD$igr ), inames )
elab <- c( paste(inames,"18m - baseline"),
           paste(inames[c(2,3,3)], "vs.",
                 inames[c(1,1,2)], "change") )
rownames( e0 )[7:12] <- elab
round( e0[7:12,c(1,5,6,4)], 4 )


###################################################
### code chunk number 18: forest-ins
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
plotEst( e0[7:12,c(1,5,6)], lwd=3, vref=0, y=rev(c(1:3,5:7))+1,
         xlab="", ylim=c(0,8), xlim=c(-0.03,0.03),
                  col.lines =iclr[c(1:3,1,1,2)],
                  col.points=iclr[c(1:3,2,3,3)], cex=1.5 )
axis( side=1 )
text( -0.03, 0.2, "Improvement"  , adj=0, col=clr[4], cex=0.8 )
text(  0.03, 0.2, "Deterioration", adj=1, col=clr[4], cex=0.8 )
mtext( "Mean carotid IMT (mm)", side=1, line=3/1.6 )


###################################################
### code chunk number 19: fu-code
###################################################
mm <- model.matrix( ~ visit:igr-1, data = AD )
m7 <- mm[,xx <- grep("v7",colnames(mm))]
head( m7 )


###################################################
### code chunk number 20: mod ident base
###################################################
m1 <- lmer( fimtavg ~ m7 + over.65 + pre.ins + sdc + (1|subjid),
            data = AD )
summary( m1 )
round( ee <- ci.lin( m1, subset=1:4 ), 4 )
C1 <- rbind( c(1,0,0,0),
             c(1,1,0,0),
             c(1,0,0,0),
             c(1,0,1,0),
             c(1,0,0,0),
             c(1,0,0,1) )
C1 <- rbind( C1, C1[2,]-C1[1,],
                 C1[4,]-C1[3,],
                 C1[6,]-C1[5,] )
C1 <- rbind( C1, C1[8,]-C1[7,],
                 C1[9,]-C1[7,],
                 C1[9,]-C1[8,] )
rownames(C1) <- row.names(C0)
colnames(C1) <- row.names(ee)
C1
round( e1 <- ci.lin( m1, subset=1:4, ctr.mat=C1 ), 4 )
rownames( e1 ) <- rownames( e0 )


###################################################
### code chunk number 21: cmp base
###################################################
round( cbind( e0[,1:2], e1[,1:2] ), 4 )


###################################################
### code chunk number 22: imp
###################################################
library( mice )
load( file="./data/imp" )
class( imp )
mif <- with( imp, lm( fimtavg.v7 ~ fimtavg.v1 + igr + over.65 + pre.ins + sdc ) )
round( smif <- summary( pool( mif) ), 4 )
mim <- with( imp, lm( fimtavg.v7 ~ fimtavg.v1 + igr ) )
round( smim <- summary( pool( mim ) ), 4 )
CI <- rbind( c(1,0), c(0,1), c(-1,1) )
rownames( CI ) <- elab[4:6]
CI
round( i0 <- ci.lin( pool( mif ), subset="igr", ctr.mat=CI ), 4 )


###################################################
### code chunk number 23: forest-impx
###################################################
round( rbind( e1[7:9,c(1,5,6)],
              i0[   ,c(1,5,6)] ), 3 )
tmpl <-
function(xtra=FALSE)
{
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
piclr <- rgb(t(col2rgb( iclr ))/3 + 255*2/3,max=255)
plotEst(  rbind( e1[7:9,c(1,5,6)],
                 i0[   ,c(1,5,6)] ),
          xlab="", ylim=c(0.3,8), xlim=c(-0.03,0.03),
          lwd=3, vref=0, y=c(8:6,4:2),
          col.lines =c(iclr[c(1:3,1,1,2)]),
          col.points=c(iclr[c(1:3,2,3,3)]),
          cex=1.5 )
if( xtra )
  {
linesEst(  e0[10:12,c(1,5,6)],
          lwd=3, vref=0, y=c(4:2-0.1),
          col.lines =piclr[c(1,1,2)],
          col.points=piclr[c(2,3,3)],
          cex=1.5 )
linesEst(  i0[,c(1,5,6)],
          lwd=3, vref=0, y=c(4:2),
          col.lines =iclr[c(1,1,2)],
          col.points=iclr[c(2,3,3)],
          cex=1.5 )
   }
axis( side=1 )
text( -0.03, 0.4, "Improvement"  , adj=0, col=clr[4], cex=0.8 )
text(  0.03, 0.4, "Deterioration", adj=1, col=clr[4], cex=0.8 )
mtext( "Mean carotid IMT (mm)", side=1, line=3/1.6 )
}
tmpl(TRUE)


###################################################
### code chunk number 24: forest-imp
###################################################
tmpl(FALSE)


###################################################
### code chunk number 25: pr-ins.rnw:179-182
###################################################
win.metafile( "./results/forest-ins-imp.emf", width=7, height=4, pointsize=20 )
tmpl(FALSE)
dev.off()


###################################################
### code chunk number 26: pr-int.rnw:6-10
###################################################
options( width=90 )
library( Epi )
library( lme4 )
load( file="./data/AD.Rda" )


###################################################
### code chunk number 27: pr-int.rnw:18-18
###################################################



###################################################
### code chunk number 28: mod sep base ins
###################################################
AD$subjid <- factor( AD$subjid )
m0 <- lmer( fimtavg ~ -1 + visit +
                      + over.65 + pre.ins + sdc +(1|subjid), data = AD )
mM <- lmer( fimtavg ~ -1 + visit + visit:grp
                      + over.65 + pre.ins + sdc +(1|subjid), data = AD )
mI <- lmer( fimtavg ~ -1 + visit + visit:igr
                      + over.65 + pre.ins + sdc +(1|subjid), data = AD )
mA <- lmer( fimtavg ~ -1 + visit + visit:grp + visit:igr
                      + over.65 + pre.ins + sdc +(1|subjid), data = AD )
mi <- lmer( fimtavg ~ -1 + visit + visit:interaction(igr,grp)
                      + over.65 + pre.ins + sdc +(1|subjid), data = AD )
round( ci.lin( m0 ), 3 )
round( ci.lin( mi ), 3 )
mod <- rbind( anova( m0, mM, mA     )[-1,],
              anova( m0, mI, mA, mi )[-1,] )
rownames( mod ) <- c("Met|0","Ins|Met","Ins|0","Met|Ins","Interact")
attr( mod, "heading" ) <- NULL
round( mod, 4 )[c(1,3,4,2,5),]


###################################################
### code chunk number 29: pr-int.rnw:46-67
###################################################
round( ee <- ci.lin( mi, subset="visit" ), 4 )
C0 <- diag(12)
C0[,1] <- rep(1:0,6)
C0[,2] <- rep(0:1,6)
igrps <- as.vector( outer( iN, gN, FUN=function(x,y) paste(x,y,sep=", ") ) )
for( i in 1:6 ) C0 <- rbind( C0, C0[2*i,]-C0[2*i-1,] )
rownames( C0 ) <- c( as.vector( outer( levels( AD$visit )[c(1,7)],
                                       igrps,
                                       paste ) ),
                     paste( igrps, "18m - baseline") )
C0
dim( C0 )
C0 <- rbind( C0, C0[15+1:3     ,]-C0[12+1:3     ,],
                 C0[12+c(2,3,3),]-C0[12+c(1,1,2),],
                 C0[15+c(2,3,3),]-C0[15+c(1,1,2),] )
rownames( C0 )[18+1:9] <- c( paste( iN, ":", gN[2], "vs", gN[1] ),
                             paste( rep(gN,each=3), ":",
                                    iN[c(2,3,3)], "vs", iN[c(1,1,2)] ) )
C0
e0 <- ci.lin( mi, subset="visit", ctr.mat=C0 )
round( e0, 4 )


###################################################
### code chunk number 30: forest-int
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
plotEst( e0[13:27,c(1,5,6)], lwd=3, vref=0,
         y = 18-c(1:3,
                  1:3+3.5,
                  1:3+7.5,
                  1:3+11.1,
                  1:3+14.5),
         xlab="", ylim=c(0,17), xlim=c(-0.06,0.06) )
abline( h=18-7.5, col=gray(0.7), lty=2 )
axis( side=1 )
text( -0.06, -0.2, "Improvement"  , adj=0, cex=0.8 )
text(  0.06, -0.2, "Deterioration", adj=1, cex=0.8 )
mtext( "Mean carotid IMT (mm)", side=1, line=3/1.6 )


###################################################
### code chunk number 31: imp-int
###################################################
library( mice )
load( file="./data/imp" )
class( imp )
mii <- with( imp, lm( fimtavg.v7 ~ fimtavg.v1 + interaction(igr,grp) + over.65 + pre.ins + sdc ) )
round( smii <- summary( pool( mii) ), 4 )
ci.lin( pool( mii ), subset="igr" )
( CI <- -C0[19:27,1:5*2+1] )
round( i0 <- ci.lin( pool( mii ), subset="igr", ctr.mat=CI ), 4 )


###################################################
### code chunk number 32: forest-int-impx
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
plotEst( e0[13:27,c(1,5,6)], lwd=3, vref=0,
         y = 18-c(1:3,
                  1:3+ 3.5,
                  1:3+ 7.5-0.1,
                  1:3+11.1-0.1,
                  1:3+14.5-0.1),
         xlab="", ylim=c(0,17), xlim=c(-0.06,0.06) )
linesEst( i0[,c(1,5,6)], lwd=3, vref=0,
         y = 18-c(1:3+ 7.5+0.1,
                  1:3+11.1+0.1,
                  1:3+14.5+0.1),
         col=gray(0.6) )
abline( h=18-7.5, col=gray(0.7), lty=2 )
axis( side=1 )
text( -0.06, -0.2, "Improvement"  , adj=0, cex=0.8 )
text(  0.06, -0.2, "Deterioration", adj=1, cex=0.8 )
mtext( "Mean carotid IMT (mm)", side=1, line=3/1.6 )


###################################################
### code chunk number 33: forest-int-imp
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 )
plotEst( e0[13:27,c(1,5,6)], lwd=3, vref=0,
         y = 18-c(1:3,
                  1:3+ 3.5,
                  1:3+ 7.5,
                  1:3+11.1,
                  1:3+14.5), col=rep(c("black","transparent"),
                                     c(6,9) ),
         xlab="", ylim=c(0,17), xlim=c(-0.06,0.06) )
linesEst( i0[,c(1,5,6)], lwd=3, vref=0,
         y = 18-c(1:3+ 7.5,
                  1:3+11.1,
                  1:3+14.5),
         col=gray(0.0) )
abline( h=18-7.5, col=gray(0.0), lty=2 )
axis( side=1 )
text( -0.06, -0.2, "Improvement"  , adj=0, cex=0.8 )
text(  0.06, -0.2, "Deterioration", adj=1, cex=0.8 )
mtext( "Mean carotid IMT (mm)", side=1, line=3/1.6 )


###################################################
### code chunk number 34: sec-met.rnw:8-18
###################################################
library( Epi )
library( lme4 )
library( abind )
load( file="./data/AD.Rda" )
with( AD, ftable( addmargins( table(grp,igr,PP=pp),
                              2:3 ) ) )
with( AD, ftable( addmargins( tapply( subjid,
                                      list(grp,igr,PP=pp),
                                      function(x) length(unique(x) ) ),
                              2:3 ) ) )


###################################################
### code chunk number 35: sec-met.rnw:29-33
###################################################
load( file="./data/base.Rda" )
tt <- addmargins( with( base, table( grp, hba1c.b7>7.0 )[2:1,] ) )
cbind( tt, round( 100*tt/tt[,3], 1 ) )[,-6]
twoby2( tt )


###################################################
### code chunk number 36: sec-met.rnw:40-43
###################################################
m1 <- glm( hba1c.b7<=7.0 ~ grp -1, family=binomial(link="log"), data=base )
round( 100*ci.exp( m1 ), 1 )
anova( m1, update( m1, . ~ 1 ), test="Chisq" )


###################################################
### code chunk number 37: sec-met.rnw:49-53
###################################################
mi <- update( m1, . ~ . + 1 )
mx <- update( mi, . ~ . + igr )
round( rbind( ci.exp( mi ),
              ci.exp( mx ) ), 2 )


###################################################
### code chunk number 38: lmer 1
###################################################
ITa <- lmer( hba1c ~ grp + grp:factor(visit) - 1 +
                     sdc + over.65 + pre.ins +
                     (1|subjid),
                     data = AD )
round( ci.lin( ITa ), 3 )
# Per protocol analysis
PPa <- update( ITa, data = subset( AD, pp ) )


###################################################
### code chunk number 39: sec-met.rnw:78-80
###################################################
ITr <- update( ITa, . ~ grp + grp:factor(visit) - 1 + (1|subjid) )
PPr <- update( ITr, data = subset( AD, pp ) )


###################################################
### code chunk number 40: sec-met.rnw:84-86
###################################################
ITc <- update( ITa, . ~ . + sex + statin + gad.pos + cvd )
PPc <- update( ITc, data = subset( AD, pp ) )


###################################################
### code chunk number 41: sec-met.rnw:99-102
###################################################
eM <- rbind( cbind(1,rbind(0,diag(6))), rep(0:1,c(6,1)) )
rownames( eM ) <- c( paste( "Vis", 1:7, sep="" ), "v7-v1" )
eM


###################################################
### code chunk number 42: sec-met.rnw:106-110
###################################################
round( ci.exp( ITa, Exp=F ), 3 )
round( ci.exp( ITa, subset="Met", Exp=F ), 3 )
round( ci.exp( ITa, subset="Plc", Exp=F ), 3 )
round( ci.lin( ITa, subint=c("Met","v7"), Exp=F ), 3 )


###################################################
### code chunk number 43: sec-met.rnw:112-117
###################################################
c.plc <- ci.exp( ITa, ctr.mat=eM, subset="Plc", Exp=F )
c.met <- ci.exp( ITa, ctr.mat=eM, subset="Met", Exp=F )
c.dif <- ci.exp( ITa, ctr.mat=cbind(eM,-eM),
                 subset=c("Met","Plc"), Exp=F )
round( cbind( c.met, c.plc, c.dif), 3 )


###################################################
### code chunk number 44: sec-met.rnw:124-135
###################################################
test <- rbind( Wald( ITa, subset=c("Met","Plc"), ctr.mat=cbind(-eM,eM)[-8,] ),
               Wald( ITa, subset=c("Met","Plc"), ctr.mat=cbind(-eM,eM)[7,,drop=F] ),
               Wald( ITa, subset=c("Met","Plc"), ctr.mat=cbind(-eM,eM)[8,,drop=F] ),
               Wald( ITa, subint=c("Met","v7" ) ),
               Wald( ITa, subint=c("Plc","v7" ) ) )
rownames( test ) <- c("    All equal",
                      "Visit 7 equal",
                      " Change equal",
                      "  Met chg = 0",
                      "  Plc chg = 0")
round( test, 4 )


###################################################
### code chunk number 45: resfun
###################################################
resfun <-
function( ITa )
{
c.plc <- ci.exp( ITa, ctr.mat=eM, subset="Plc", Exp=F )
c.met <- ci.exp( ITa, ctr.mat=eM, subset="Met", Exp=F )
c.dif <- ci.exp( ITa, subset=c("Met","Plc"), Exp=F, ctr.mat=cbind(eM,-eM) )

test <- rbind( Wald( ITa, subset=c("Met","Plc"), ctr.mat=cbind(-eM,eM)[-8,] ),
               Wald( ITa, subset=c("Met","Plc"), ctr.mat=cbind(-eM,eM)[7,,drop=F] ),
               Wald( ITa, subset=c("Met","Plc"), ctr.mat=cbind(-eM,eM)[8,,drop=F] ),
               Wald( ITa, subint=c("Met","v7" ) ),
               Wald( ITa, subint=c("Plc","v7" ) ) )
rownames( test ) <- c("    All equal",
                      "Visit 7 equal",
                      " Change equal",
                      "  Met chg = 0",
                      "  Plc chg = 0")
eff <- cbind( c.met, c.plc, c.dif )
colnames( eff )[1+0:2*3] <- c("Met","Plc","M-P")
print( round( eff, 3 ) )
print( round( test, 4 ) )
invisible( list( eff=eff, test=test ) )
}


###################################################
### code chunk number 46: sec-met.rnw:165-177
###################################################
cat( "\nIntention to treat, stratum variables:\n----------------------\n" )
hbr <- resfun( ITa )
cat( "\nIntention to treat, no stratum variables:\n----------------------\n" )
resfun( ITr )
cat( "\nIntention to treat, stratum variables + confounders:\n----------------------\n" )
resfun( ITc )
cat( "\nPer protocol, stratum variables:\n----------------------\n" )
resfun( PPa )
cat( "\nPer protocol, no stratum variables:\n----------------------\n" )
resfun( PPr )
cat( "\nPer protocol, stratum variables + confounders:\n----------------------\n" )
resfun( PPc )


###################################################
### code chunk number 47: hba1c
###################################################
plt <-
function()
{
    matplot( 0:6*3, hbr$eff[1:7,],
         xlab="Months since trial entry", xaxt="n", xlim=c(0,19.5),
         ylim=c(7,9.5),  yaxs="i", yaxt="n",
         ylab=expression( "Hb"*A[1][c]*" (%)  [Mean (95% CI)]" ),
         type="n", lwd=c(4,1,1), lty=1, col=rep(clr[1:2],each=3) )
abline( h=seq(7,9.5,0.2), col=gray(0.8) )
matlines( 0:6*3, hbr$eff[1:7,],
         type="l", lwd=c(5,1,1), lty=1, col=rep(clr[1:2],each=3) )
axis( side=1, at=0:6*3, col=clr[4] )
axis( side=2, at=7:9, col=clr[4] )
text( c(20,20), c(9.3,9.1), c("Metformin","Placebo"), col=clr[1:2],
     font=2, adj=1, cex=1 )
text( 20, 8.9, substitute( "Equal change in Hb"*A[1][c]*" 0-18 mth: P ="*pval,
                            list(pval = formatC( hbr$test[3,3], format="f", digits=4 )) ),
      col=clr[4], adj=1 )
text( 19.9, 8.7, "HbA1c change\n from baseline", col=clr[4], adj=c(1,0.5) )
text( c(20,20), hbr$eff[7,c(1,4)],
      paste( formatC( hbr$eff[8,c(1,4)],
                      format="f", digits=1 ), "%", sep="" ),
      col=clr[1:2], adj=1 )
}
oldpar <- par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1, bty="n" )
plt()


###################################################
### code chunk number 48: ana-fun-def
###################################################
var <- "hba1c"


###################################################
### code chunk number 49: sec-met.rnw:238-240
###################################################
# Names of tests performed
( tnames <- rownames( test ) )


###################################################
### code chunk number 50: sec-met.rnw:246-271
###################################################
CM <-
function( mod, eff )
{
# Contrast matrix to fish out the estimates, resp tests
cM <- cbind(1,rbind(0,diag(6)))
rownames( cM ) <- paste( "Vis", 1:7 )
# Find missing visits for the variable, make two versions of the matrix
tp <- sapply( strsplit(rownames(ci.exp(mod,subset="Met")), "" ),
              function(x) x[length(x)] )
tp[1] <- 1
tp <- as.numeric( tp )
# The contrast row should be between the last existing measurement and
# the first:
cM <- rbind( cM, cM[tp[length(tp)],]-cM[1,] )
# Are there any missing timepoints in this analysis?
mv <- setdiff(1:7,tp)
if( length(mv) > 0 )
  {
  # If used to extract effects keep rows for all timepoints, but
  # return NAs for timepoints with no measurements:
  cM[mv,] <- NA
  cM <- if(eff) cM[,-mv] else cM[-mv,-mv]
  }
cM
}


###################################################
### code chunk number 51: sec-met.rnw:274-284
###################################################
eFUN <-
function( mod )
{
eM <- CM( mod, eff=TRUE )
# Insert some NA-rows to get results in right dimension
met <- ci.exp( mod, ctr.mat=      eM,      subset=  "Met",        Exp=F )
plc <- ci.exp( mod, ctr.mat=      eM,      subset=        "Plc",  Exp=F )
dif <- ci.exp( mod, ctr.mat=cbind(eM,-eM), subset=c("Met","Plc"), Exp=F )
abind( list(Met=met, Plc=plc, Dif=dif), rev.along=0 )
}


###################################################
### code chunk number 52: sec-met.rnw:287-295
###################################################
sFUN <-
function( mod )
{
res <- c( attr( VarCorr( mod )$subjid, "stddev" ),
          attr( VarCorr( mod ), "sc" ) )
names( res ) <- c("Btw","Res")
res
}


###################################################
### code chunk number 53: sec-met.rnw:298-313
###################################################
tFUN <-
function( mod )
{
eX <- CM( mod, eff=FALSE )
np <- nrow( eX )
lastM <- grep( "Met", rownames(ci.exp(mod)) ) ; lastM <- lastM[length(lastM)]
lastP <- grep( "Plc", rownames(ci.exp(mod)) ) ; lastP <- lastP[length(lastP)]
test <- rbind( Wald( mod, subset=c("Met","Plc"), ctr.mat=cbind(-eX,eX)[-np,] ),
               Wald( mod, subset=c("Met","Plc"), ctr.mat=cbind(-eX,eX)[ np-1,,drop=F] ),
               Wald( mod, subset=c("Met","Plc"), ctr.mat=cbind(-eX,eX)[ np  ,,drop=F] ),
               Wald( mod, subset=lastM ),
               Wald( mod, subset=lastP ) )
rownames( test ) <- tnames
test
}


###################################################
### code chunk number 54: sec-met.rnw:317-323
###################################################
# Names of variable to be used in repeated measures analyses
vnam <- c("hba1c","weight","bmi","whr","gluc","ins","cpep",
          "idos","ipkg","sys","dia","pulse",
          "chol","ldl","hdl","vldl","trig",
          "fimtavg","fimtmax","iem","csc2","imtareal","n.pl")
match( vnam, names(AD) )


###################################################
### code chunk number 55: pairs
###################################################
pairs( AD[,c("visit",vnam)], gap=0, pch=16, cex=0.3 )


###################################################
### code chunk number 56: hist
###################################################
par( mfrow=c(6,4), mar=c(2,1,1,1), mgp=c(3,1,0)/1.6 )
for( vv in vnam ) hist( AD[,vv], breaks=50, border="black",
                        col="black", main=vv, xlab="",ylab="",yaxt="n" )
hist( log10(AD[,"ins"]), breaks=50, border="black",
                         col="black", main="log10(ins)", xlab="",ylab="",yaxt="n" )


###################################################
### code chunk number 57: sec-met.rnw:351-359
###################################################
ltrf <- rep( FALSE, length(vnam) )
ltrf[c(6:9,16,17,23)] <- TRUE
cbind( vnam, ltrf )
for( vv in vnam )
   {
cat( "-----------------------------------\n",vv )
print( table( !is.na(AD[,vv]), AD$visit ) )
   }


###################################################
### code chunk number 58: sec-met.rnw:367-387
###################################################
# Arrays to hold results over variables
# Variance components
Std <- NArray( list( var = paste( vnam, ifelse(ltrf,"(l)",""), sep="" ),
                     ana = c("IT","PP"),
                     mod = c("Prim","Rest","Conf"),
                     std = c("Btw","Res") ) )
# Estimates
Eff <- NArray( c( dimnames(Std)[-4],
                  list( par = c(paste("Vis",1:7),"V7-V1"),
                              c("Est","lo","hi"),
                              c("Met","Plc","Dif") ) ) )
# Tests
Tst <- NArray( c( dimnames(Std)[-4],
                  list( test = tnames,
                        what = c("Chisq","df","Pval") ) ) )
# Smaller versions of the arrys to be used inside the function that
# does the calculations for each response variable
std <- NArray( dimnames(Std)[-1] )
eff <- NArray( dimnames(Eff)[-1] )
tst <- NArray( dimnames(Tst)[-1] )


###################################################
### code chunk number 59: sec-met.rnw:395-467
###################################################
# Function to do PPplots for residuals
PPplot <-
function( x, xl )
{
n <- length(x)
plot( pnorm( sort(x/sd(x)) ), (1:n-0.5)/n,
      pch=16, cex=0.4,
      xlim=0:1, ylim=0:1,
      xlab="", ylab="",
      xaxt="n", yaxt="n" )
text( 0, 1, xl, adj=c(0,1) )
abline( 0, 1 )
}

# A function that returns NA instead of -Inf for 0 argument
logI <- function(x) ifelse(x>0,log(x),NA)

# Function that makes all analyses for one response variable
ana.fun <-
function( var, log.tr=FALSE, respl=FALSE )
{
# We need a single name of the response variable
AD$Y <- if(log.tr) logI(AD[,var]) else AD[,var]
# Fit models to total and PP part of dataset:
# r: Reduced model
# a: Primary model
# c: Confounder-expanded model
ITr <- lmer( Y ~ grp + grp:factor(visit) - 1 + (1|subjid), data = AD )
ITa <- lmer( Y ~ grp + grp:factor(visit) - 1 + sdc + over.65 + pre.ins
                                             + (1|subjid), data = AD )
ITc <- lmer( Y ~ grp + grp:factor(visit) - 1 + sdc + over.65 + pre.ins
                                             + sex + statin + gad.pos + cvd
                                             + (1|subjid), data = AD )
# ITa <- update( ITr, . ~ . + sdc + over.65 + pre.ins )
# ITc <- update( ITa, . ~ . + sex + statin + gad.pos + cvd )
PPa <- update( ITa, data = subset( AD, pp ) )
PPr <- update( ITr, data = subset( AD, pp ) )
PPc <- update( ITc, data = subset( AD, pp ) )

# Plots if specified - but only for ITa analysis - note the S4 way of
# extracting slots from the objects, does not work for medMod objects
# if( respl )
# {
# PPplot( ITa@ranef, paste(var,if(log.tr)"(log-tr)","\nptt") )
# PPplot( ITa@resid, paste(var,if(log.tr)"(log-tr)","\nres") )
# }

# Fixed effects
eff["IT","Prim",,,] <- eFUN( ITa )
eff["IT","Rest",,,] <- eFUN( ITr )
eff["IT","Conf",,,] <- eFUN( ITc )
eff["PP","Prim",,,] <- eFUN( PPa )
eff["PP","Rest",,,] <- eFUN( PPr )
eff["PP","Conf",,,] <- eFUN( PPc )
if( log.tr ) eff <- exp( eff )
# Variance components
std["IT","Prim",] <- sFUN( ITa )
std["IT","Rest",] <- sFUN( ITr )
std["IT","Conf",] <- sFUN( ITc )
std["PP","Prim",] <- sFUN( PPa )
std["PP","Rest",] <- sFUN( PPr )
std["PP","Conf",] <- sFUN( PPc )
# Tests
tst["IT","Prim",,] <- tFUN( ITa )
tst["IT","Rest",,] <- tFUN( ITr )
tst["IT","Conf",,] <- tFUN( ITc )
tst["PP","Prim",,] <- tFUN( PPa )
tst["PP","Rest",,] <- tFUN( PPr )
tst["PP","Conf",,] <- tFUN( PPc )
# Return all in a list
list( eff=eff, std=std, tst=tst )
}


###################################################
### code chunk number 60: PP
###################################################
par( mfrow=c(8,6), mar=c(0,0,0,0), omi=c(3,3,1,1)/4 )
for( i in 1:length(vnam) )
   {
cat( vnam[i], "\n" )
   res <- ana.fun( vnam[i], ltrf[i], respl=TRUE )
   Eff[i,,,,,] <- res$eff
   Std[i,,,]   <- res$std
   Tst[i,,,,]  <- res$tst
#   mtext( expression(Phi^-1*"[sort(std. res.)]"),
#          side=1, line=2, outer=TRUE )
#   mtext( "Uniform[0,1]", side=2, line=2, outer=TRUE )
   }
apply( Eff, 1:2, FUN=function(x) mean(is.na(x)) )


###################################################
### code chunk number 61: sec-met.rnw:513-526
###################################################
round( ftable( Eff[,"IT","Prim",,,], col.vars=4:3 ), 3 )
round( ftable( Eff[,"IT","Prim",,,], col.vars=4:3 ), 3 )
round( ftable( Eff[,"IT","Rest",,,], col.vars=4:3 ), 3 )
round( ftable( Eff[,"IT","Rest",,,], col.vars=4:3 ), 3 )
round( ftable( Eff[,"IT","Conf",,,], col.vars=4:3 ), 3 )
round( ftable( Eff[,"IT","Conf",,,], col.vars=4:3 ), 3 )
#for( vv in 1:(dim(Eff)[1]) )
#   {
#   cat( "--------------------------------\n",  dimnames(Eff)[[1]][vv], "\n" )
#   dec <- 3 - ( max(Eff[vv,,,,,1:3])>100  ) -
#              ( max(Eff[vv,,,,,1:3])>1000 )
#   print( round( ftable( Eff[vv,,,,,], col.vars=5:4 ), dec ) )
#   }


###################################################
### code chunk number 62: sec-met.rnw:529-534
###################################################
for( vv in dimnames(Tst)[[1]] )
   {
   cat( "--------------------------------\n", vv, "\n" )
   print( round( ftable( Tst[vv,,,,"Pval"], col.vars=3 ), 4 ) )
   }


###################################################
### code chunk number 63: sec-met.rnw:542-620
###################################################
source( "cnr.R" )
cnr
pleff <-
function( var,
         ytxt = var,
           yl = NULL, # yaxis limits
           yt = NULL, # yaxis major ticks
           ym = NULL, # yaxis minor ticks
           hl = ym )
#         ltyp = if ltrf[match(var,vnam)] "pch" else "ach",
#                # Type of labelling of curves:
#                # "pch" - percent change
#               # "ach" - absolute change
#                # Any other text will be printed
#         lval = NULL ) # Values plotted at end of curves )
{
# Get the estimates, but only for the times where they are, otherwise
# the the points will not be connected (a feature)
var <- dimnames(Eff)$var[grep(var,vnam <- dimnames(Eff)$var)]
eff <- Eff[var,"IT","Prim",,,]
pef <- cbind( eff[-8,,"Met"],
              eff[-8,,"Plc"] )
pef <- pef[complete.cases(pef),]
tim <- (as.numeric( gsub("Vis ","",rownames(pef) ) )-1)*3

if( is.null(yl) ) yl <- range(pef)*c(0.9,1.2)
matplot( tim, pef, type="n",
         xlab="Months since trial entry", xaxt="n", xlim=c(0,19.5),
         ylab=ytxt, ylim=yl, yaxs="i", yaxt="n", bty="n" )
if( !is.null(hl) ) abline( h=hl, col=gray(0.8) )
axis( side=1, at=0:6*3, col=clr[4] )
axis( side=2, if(!is.null(yt)) at=yt, col=clr[4] )
axis( side=2, if(!is.null(ym)) at=ym, labels=rep("",length(ym)), col=clr[4] )
text( cnr(c(98,98),c(92,98)), gN, col=clr[2:1],
      font=2, adj=c(1,1), cex=1.2 )
if( FALSE )
{
text( cnr(98,74)[[1]], cnr(98,74)[[2]],
      substitute( "Equal over 0-18 mth: P = "*pval,
                  list( pval = formatC( Tst[var,"IT","Prim","    All equal","Pval"],
                        format="f", digits=4 )) ),
      col=clr[4], adj=1 )
text( cnr(98,78)[[1]], cnr(98,78)[[2]],
      substitute( "Equal at end: P = "*pval,
                  list( pval = formatC( Tst[var,"IT","Prim","Visit 7 equal","Pval"],
                        format="f", digits=4 )) ),
      col=clr[4], adj=1 )
}
text( cnr(98,82)[[1]], cnr(98,82)[[2]],
      substitute( "Equal change: P = "*pval,
                  list( pval = formatC( Tst[var,"IT","Prim"," Change equal","Pval"],
                        format="f", digits=4 )) ),
      col=clr[4], adj=1 )
if( ltrf[match(var,vnam)] ) { text( cnr(98,74)[[1]], cnr(98,74)[[2]],
                                  "% change\nfrom baseline", col=clr[4], adj=c(1,0.5) )
} else {                       text( cnr(98,74)[[1]], cnr(98,74)[[2]],
                                    "Change\nfrom baseline", col=clr[4], adj=c(1,0.5) ) }
####################################################Experimental
#if( ltyp=="pch" ) text( cnr(98,74)[[1]], cnr(98,74)[[2]],
#                        "% change\nfrom baseline", col=clr[4], adj=c(1,0.5) )
#else
#if( ltyp=="ach" ) text( cnr(98,74)[[1]], cnr(98,74)[[2]],
#                        "Change\nfrom baseline", col=clr[4], adj=c(1,0.5) )
#else text( cnr(98,74)[[1]], cnr(98,74)[[2]],
#           ltyp, col=clr[4], adj=c(1,0.5) )
EF <- eff[8,"Est",c("Met","Plc")]
# It might not be visit 7 that is the last non-missing, so find that
ends <- eff[1:7,"Est",c("Met","Plc")]
ends <- ends[max(which(!is.na(ends[,1]))),]
dv <- ends-mean(ends) ; dv <- dv/abs(dv)
if( abs(diff(ends)) < diff(yl)/20 ) ends <- mean(ends)+dv*diff(yl)/30
if( ltrf[match(var,vnam)] ) EF <- (EF-1)*100
print( c(EF, ends ) )
text( c(20,20), ends,
      paste( formatC( EF, format="f", digits=1 ), sep="" ),
      col=clr[1:2], adj=1, font=2 )
matlines( tim, pef, type="l", lwd=c(5,1,1), lty=1, col=rep(clr[1:2],each=3) )
}


###################################################
### code chunk number 64: tr-hba
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 )
pleff( "hba1c", expression("Hb"*A[1][c]*" (%)  [Mean (95% CI)]"),
       yl = c(7,10),
       yt = 7:10,
       ym = seq(7,10,0.5) )


###################################################
### code chunk number 65: tr-trig
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 )
pleff("trig","Triglycerides (mmol/l)  [Mean (95% CI)]",
              yl = c(1,2.5),       # yaxis limits
              yt = 0:4,            # yaxis major ticks
              ym = seq(1,3,0.25) ) # yaxis minor ticks


###################################################
### code chunk number 66: first4
###################################################
tmpl <-
function()
{
par( mfrow=c(2,2),  mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 )
# HbA1c:
pleff( "hba1c",
       expression("Hb"*A[1][c]*" (%)  [Mean (95% CI)]"),
       yl = c(7,10),
       yt = 7:10,
       ym = seq(7,10,0.5) )
text( cnr(2,98), "a", adj=c(0,1), font=2, cex=1.5 )
# Gluc
pleff("gluc","Fasting p-glucose (mmol/L) [Mean (95% CI)]",
              yl = c(7,12),      # yaxis limits
              yt = seq(7,13,2),  # yaxis major ticks
              ym = seq(7,13,1) ) # yaxis minor ticks
text( cnr(2,98), "b", adj=c(0,1), font=2, cex=1.5 )
# Insulin dose
pleff("ipkg","Insulin dose (IU/day/kg)  [Mean (95% CI)]",
              yl = c(0.0,2.4),        # yaxis limits
              yt = seq(0.0,2.0,1.0),  # yaxis major ticks
              ym = seq(0.0,2.4,0.2) ) # yaxis minor ticks
text( cnr(2,98), "c", adj=c(0,1), font=2, cex=1.5 )
# Weight:
pleff("weight","Weight (kg)  [Mean (95% CI)]",
              yl = c(90,120),      # yaxis limits
              yt = seq(90,120,10), # yaxis major ticks
              ym = seq(90,120,5) ) # yaxis minor ticks
text( cnr(2,98), "d", adj=c(0,1), font=2, cex=1.5 )
}
tmpl()


###################################################
### code chunk number 67: sec-met.rnw:684-687
###################################################
win.metafile( "./results/first4-met.emf", height=9, width=8 )
tmpl()
dev.off()


###################################################
### code chunk number 68: sec-ins.rnw:8-18
###################################################
library( Epi )
library( lme4 )
library( abind )
load( file="./data/AD.Rda" )
with( AD, ftable( addmargins( table(grp,igr,PP=pp),
                              2:3 ) ) )
with( AD, ftable( addmargins( tapply( subjid,
                                      list(grp,igr,PP=pp),
                                      function(x) length(unique(x) ) ),
                              2:3 ) ) )


###################################################
### code chunk number 69: sec-ins.rnw:29-32
###################################################
load( file="./data/base.Rda" )
tt <- addmargins( with( base, table( igr, hba1c.b7>7.0 ) ) )
cbind( tt, round( 100*tt/tt[,3], 1 ) )[,-6]


###################################################
### code chunk number 70: sec-ins.rnw:37-40
###################################################
m1 <- glm( hba1c.b7<=7.0 ~ igr -1, family=binomial(link="log"), data=base )
round( 100*ci.exp( m1 ), 1 )
anova( m1, update( m1, . ~ 1 ), test="Chisq" )


###################################################
### code chunk number 71: sec-ins.rnw:47-56
###################################################
CM = rbind( c(-1,1,0), c(-1,0,1), c(0,-1,1) )
rownames( CM ) <- paste( levels(base$igr)[c(2,3,3)], "vs.",
                         levels(base$igr)[c(1,1,2)] )
CM
round( ci.exp( m1 ), 2 )
cbind(
round( ci.exp( m1, ctr.mat= CM ), 2 ),
round( ci.exp( m1, ctr.mat=-CM ), 2 ),
round( ci.lin( m1, ctr.mat= CM )[,4,drop=F], 3 ) )


###################################################
### code chunk number 72: lmer 1
###################################################
MM <- model.matrix( ~ igr:factor(visit), data=AD )[,-(1:4)]
colnames( MM )
IT0 <- lmer( hba1c ~ MM +
                     sdc + over.65 + pre.ins +
                     (1|subjid),
                     data = AD )
round( ci.lin( IT0 ), 3 )
ITa <- lmer( hba1c ~ igr + igr:factor(visit) - 1 +
                     sdc + over.65 + pre.ins +
                     (1|subjid),
                     data = AD )
round( ci.lin( ITa ), 3 )
# Per protocol analysis
PPa <- update( ITa, data = subset( AD, pp ) )


###################################################
### code chunk number 73: sec-ins.rnw:88-90
###################################################
ITr <- update( ITa, . ~ igr + igr:factor(visit) - 1 + (1|subjid) )
PPr <- update( ITr, data = subset( AD, pp ) )


###################################################
### code chunk number 74: sec-ins.rnw:94-96
###################################################
ITc <- update( ITa, . ~ . + sex + statin + gad.pos + cvd )
PPc <- update( ITc, data = subset( AD, pp ) )


###################################################
### code chunk number 75: sec-ins.rnw:109-112
###################################################
eM <- rbind( cbind(1,rbind(0,diag(6))), rep(0:1,c(6,1)) )
rownames( eM ) <- c( paste( "Vis", 1:7, sep="" ), "v7-v1" )
eM


###################################################
### code chunk number 76: sec-ins.rnw:116-121
###################################################
round( ci.exp( ITa, Exp=F ), 3 )
round( ci.exp( ITa, subset="Biph", Exp=F ), 3 )
round( ci.exp( ITa, subset="AspD", Exp=F ), 3 )
round( ci.exp( ITa, subset="Detm", Exp=F ), 3 )
round( ci.exp( ITa, subint=c("Detm","v7"), Exp=F ), 3 )


###################################################
### code chunk number 77: sec-ins.rnw:123-137
###################################################
( ilev <- levels( AD$igr ) )
c.1 <- ci.exp( ITa, ctr.mat=eM, subset=ilev[1], Exp=F )
c.2 <- ci.exp( ITa, ctr.mat=eM, subset=ilev[2], Exp=F )
c.3 <- ci.exp( ITa, ctr.mat=eM, subset=ilev[3], Exp=F )
c.21 <- ci.exp( ITa, ctr.mat=cbind(eM,-eM), subset=ilev[c(2,1)], Exp=F )
c.31 <- ci.exp( ITa, ctr.mat=cbind(eM,-eM), subset=ilev[c(3,1)], Exp=F )
c.32 <- ci.exp( ITa, ctr.mat=cbind(eM,-eM), subset=ilev[c(3,2)], Exp=F )
eff <- cbind( c.1, c.2, c.3 )
colnames(eff)[1+0:2*3] <- ilev
dff <- cbind( c.21, c.31, c.32)
colnames(dff)[1+0:2*3] <- paste( ilev[c(2,3,3)], "-",
                                 ilev[c(1,1,2)], sep="" )
round( eff, 3 )
round( dff, 3 )


###################################################
### code chunk number 78: sec-ins.rnw:144-165
###################################################
test <- rbind( Wald( ITa, subset=ilev[c(2,1)], ctr.mat=cbind(-eM,eM)[-8,] ),
               Wald( ITa, subset=ilev[c(3,1)], ctr.mat=cbind(-eM,eM)[-8,] ),
               Wald( ITa, subset=ilev[c(3,2)], ctr.mat=cbind(-eM,eM)[-8,] ),
               Wald( ITa, subset=ilev[c(2,1)], ctr.mat=cbind(-eM,eM)[7,,drop=F] ),
               Wald( ITa, subset=ilev[c(3,1)], ctr.mat=cbind(-eM,eM)[7,,drop=F] ),
               Wald( ITa, subset=ilev[c(3,2)], ctr.mat=cbind(-eM,eM)[7,,drop=F] ),
               Wald( ITa, subset=ilev[c(2,1)], ctr.mat=cbind(-eM,eM)[8,,drop=F] ),
               Wald( ITa, subset=ilev[c(3,1)], ctr.mat=cbind(-eM,eM)[8,,drop=F] ),
               Wald( ITa, subset=ilev[c(3,2)], ctr.mat=cbind(-eM,eM)[8,,drop=F] ),
               Wald( ITa, subint=c(ilev[1],"v7" ) ),
               Wald( ITa, subint=c(ilev[2],"v7" ) ),
               Wald( ITa, subint=c(ilev[3],"v7" ) ) )
tnames <-
rownames( test ) <- c( paste( rep(c("    All equal",
                                    "Visit 7 equal",
                                    " Change equal"),each=3),
                              rep( paste( ilev[c(2,3,3)], "vs.",
                                          ilev[c(1,1,2)] ), 3 ) ),
                       paste( "       ", ilev, "change=0" ) )
round( test, 4 )
round( test[c(1,4,7,2,5,8,3,6,9),], 4 )


###################################################
### code chunk number 79: resfun
###################################################
ilev <- levels( AD$igr )
resfun <-
function( ITa )
{
c.1 <- ci.exp( ITa, ctr.mat=eM, subset=ilev[1], Exp=F )
c.2 <- ci.exp( ITa, ctr.mat=eM, subset=ilev[2], Exp=F )
c.3 <- ci.exp( ITa, ctr.mat=eM, subset=ilev[3], Exp=F )
c.21 <- ci.exp( ITa, ctr.mat=cbind(eM,-eM), subset=ilev[c(2,1)], Exp=F )
c.31 <- ci.exp( ITa, ctr.mat=cbind(eM,-eM), subset=ilev[c(3,1)], Exp=F )
c.32 <- ci.exp( ITa, ctr.mat=cbind(eM,-eM), subset=ilev[c(3,2)], Exp=F )
eff <- cbind( c.1, c.2, c.3 )
colnames(eff)[1+0:2*3] <- ilev
dff <- cbind( c.21, c.31, c.32)
colnames(dff)[1+0:2*3] <- paste( ilev[c(2,3,3)], "-",
                                 ilev[c(1,1,2)], sep="" )
test <- rbind( Wald( ITa, subset=ilev[c(2,1)], ctr.mat=cbind(-eM,eM)[-8,] ),
               Wald( ITa, subset=ilev[c(3,1)], ctr.mat=cbind(-eM,eM)[-8,] ),
               Wald( ITa, subset=ilev[c(3,2)], ctr.mat=cbind(-eM,eM)[-8,] ),
               Wald( ITa, subset=ilev[c(2,1)], ctr.mat=cbind(-eM,eM)[7,,drop=F] ),
               Wald( ITa, subset=ilev[c(3,1)], ctr.mat=cbind(-eM,eM)[7,,drop=F] ),
               Wald( ITa, subset=ilev[c(3,2)], ctr.mat=cbind(-eM,eM)[7,,drop=F] ),
               Wald( ITa, subset=ilev[c(2,1)], ctr.mat=cbind(-eM,eM)[8,,drop=F] ),
               Wald( ITa, subset=ilev[c(3,1)], ctr.mat=cbind(-eM,eM)[8,,drop=F] ),
               Wald( ITa, subset=ilev[c(3,2)], ctr.mat=cbind(-eM,eM)[8,,drop=F] ),
               Wald( ITa, subint=c(ilev[1],"v7" ) ),
               Wald( ITa, subint=c(ilev[2],"v7" ) ),
               Wald( ITa, subint=c(ilev[3],"v7" ) ) )
rownames( test ) <- tnames
print( round( eff, 3 ) )
print( round( dff, 3 ) )
print( round( test, 4 ) )
invisible( list( eff=eff, dff=dff, test=test ) )
}


###################################################
### code chunk number 80: sec-ins.rnw:206-212
###################################################
hbr <- resfun( ITa )
resfun( ITr )
resfun( ITc )
resfun( PPa )
resfun( PPr )
resfun( PPc )


###################################################
### code chunk number 81: hba1c
###################################################
plt <-
function()
{
matplot( 0:6*3, hbr$eff[1:7,],
         xlab="Months since trial entry", xaxt="n", xlim=c(0,19.5),
         ylim=c(7,9.5),  yaxs="i", yaxt="n",
         ylab=expression( "Hb"*A[1][c]*" (%)  [Mean (95% CI)]" ),
         type="n", lwd=c(4,1,1), lty=1, col=rep(iclr[1:3],each=3) )
abline( h=seq(7,9.5,0.2), col=gray(0.8) )
matlines( 0:6*3, hbr$eff[1:7,],
          type="l", lwd=c(5,1,1), lty=1, col=rep(iclr[1:3],each=3) )
axis( side=1, at=0:6*3, col=clr[4] )
axis( side=2, at=7:9, col=clr[4] )
text( c(20,20,20), c(9.3,9.1,8.9), iN, col=iclr[1:3],
      font=2, adj=1, cex=1 )
#text( 20, 8.9, substitute( "Equal change in Hb"*A[1][c]*" 0-18 mth: P ="*pval,
#                            list(pval = formatC( hbr$test[3,3], format="f", digits=4 )) ),
#      col=clr[4], adj=1 )
#text( 19.9, 8.7, "HbA1c change\n from baseline", col=clr[4], adj=c(1,0.5) )
#text( c(20,20), hbr$eff[7,c(1,4)],
#      paste( formatC( hbr$eff[8,c(1,4)],
#                      format="f", digits=1 ), "%", sep="" ),
#      col=clr[1:2], adj=1 )
}
oldpar <- par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1, bty="n" )
plt()


###################################################
### code chunk number 82: ana-fun-def
###################################################
var <- "hba1c"
mod <- ITa


###################################################
### code chunk number 83: sec-ins.rnw:275-300
###################################################
CM <-
function( mod, eff )
{
# Contrast matrix to fish out the estimates, resp tests
cM <- cbind(1,rbind(0,diag(6)))
rownames( cM ) <- paste( "Vis", 1:7 )
# Find missing visits for the variable, make two versions of the matrix
tp <- sapply( strsplit(rownames(ci.exp(mod,subset=ilev[1])), "" ),
              function(x) x[length(x)] )
tp[1] <- 1
tp <- as.numeric( tp )
# The contrast row should be between the last existing measurement and
# the first:
cM <- rbind( cM, cM[tp[length(tp)],]-cM[1,] )
# Are there any missing timepoints in this analysis?
mv <- setdiff(1:7,tp)
if( length(mv) > 0 )
  {
  # If used to extract effects keep rows for all timepoints, but
  # return NAs for timepoints with no measurements:
  cM[mv,] <- NA
  cM <- if(eff) cM[,-mv] else cM[-mv,-mv]
  }
cM
}


###################################################
### code chunk number 84: sec-ins.rnw:303-319
###################################################
eFUN <-
function( mod )
{
eM <- CM( mod, eff=TRUE )
# Insert some NA-rows to get results in right dimension
c.1 <- ci.exp( mod, ctr.mat=eM, subset=ilev[1], Exp=F )
c.2 <- ci.exp( mod, ctr.mat=eM, subset=ilev[2], Exp=F )
c.3 <- ci.exp( mod, ctr.mat=eM, subset=ilev[3], Exp=F )
c.21 <- ci.exp( mod, ctr.mat=cbind(eM,-eM), subset=ilev[c(2,1)], Exp=F )
c.31 <- ci.exp( mod, ctr.mat=cbind(eM,-eM), subset=ilev[c(3,1)], Exp=F )
c.32 <- ci.exp( mod, ctr.mat=cbind(eM,-eM), subset=ilev[c(3,2)], Exp=F )
aa <- abind( list( c.1, c.2, c.3, c.21, c.31, c.32 ), rev.along=0 )
dimnames( aa )[[3]] <- c( ilev, paste( ilev[c(2,3,3)], "-",
                                       ilev[c(1,1,2)], sep="" ) )
aa
}


###################################################
### code chunk number 85: sec-ins.rnw:322-330
###################################################
sFUN <-
function( mod )
{
res <- c( attr( VarCorr( mod )$subjid, "stddev" ),
          attr( VarCorr( mod ), "sc" ) )
names( res ) <- c("Btw","Res")
res
}


###################################################
### code chunk number 86: sec-ins.rnw:333-357
###################################################
tFUN <-
function( mod )
{
eX <- CM( mod, eff=FALSE )
np <- nrow( eX )
# Find the position of the last parameter
last1 <- grep( ilev[1], rownames(ci.exp(mod)) ) ; last1 <- last1[length(last1)]
last2 <- grep( ilev[2], rownames(ci.exp(mod)) ) ; last2 <- last2[length(last2)]
last3 <- grep( ilev[3], rownames(ci.exp(mod)) ) ; last3 <- last3[length(last3)]
test <- rbind( Wald( mod, subset=ilev[c(2,1)], ctr.mat=cbind(-eX,eX)[-np,] ),
               Wald( mod, subset=ilev[c(3,1)], ctr.mat=cbind(-eX,eX)[-np,] ),
               Wald( mod, subset=ilev[c(3,2)], ctr.mat=cbind(-eX,eX)[-np,] ),
               Wald( mod, subset=ilev[c(2,1)], ctr.mat=cbind(-eX,eX)[np-1,,drop=F] ),
               Wald( mod, subset=ilev[c(3,1)], ctr.mat=cbind(-eX,eX)[np-1,,drop=F] ),
               Wald( mod, subset=ilev[c(3,2)], ctr.mat=cbind(-eX,eX)[np-1,,drop=F] ),
               Wald( mod, subset=ilev[c(2,1)], ctr.mat=cbind(-eX,eX)[np,,drop=F] ),
               Wald( mod, subset=ilev[c(3,1)], ctr.mat=cbind(-eX,eX)[np,,drop=F] ),
               Wald( mod, subset=ilev[c(3,2)], ctr.mat=cbind(-eX,eX)[np,,drop=F] ),
               Wald( mod, subset=last1 ),
               Wald( mod, subset=last2 ),
               Wald( mod, subset=last3 ) )
rownames( test ) <- tnames
test
}


###################################################
### code chunk number 87: sec-ins.rnw:361-367
###################################################
# Names of variable to be used in repeated measures analyses
vnam <- c("hba1c","weight","bmi","whr","gluc","ins","cpep",
          "idos","ipkg","sys","dia","pulse",
          "chol","ldl","hdl","vldl","trig",
          "fimtavg","fimtmax","iem","csc2","imtareal","n.pl")
match( vnam, names(AD) )


###################################################
### code chunk number 88: sec-ins.rnw:372-375
###################################################
ltrf <- rep( FALSE, length(vnam) )
ltrf[c(6:9,16,17,23)] <- TRUE
vnam[ltrf]


###################################################
### code chunk number 89: sec-ins.rnw:383-405
###################################################
# Arrays to hold results over variables
# Variance components
Std <- NArray( list( var = paste( vnam, ifelse(ltrf,"(l)",""), sep="" ),
                     ana = c("IT","PP"),
                     mod = c("Prim","Rest","Conf"),
                     std = c("Btw","Res") ) )
# Estimates
Eff <- NArray( c( dimnames(Std)[-4],
                  list( par = c(paste("Vis",1:7),"V7-V1"),
                              c("Est","lo","hi"),
                              c( levels(AD$igr),
                          paste( levels(AD$igr)[c(2,3,3)], "vs.",
                                 levels(AD$igr)[c(1,1,2)] ) ) ) ) )
# Tests
Tst <- NArray( c( dimnames(Std)[-4],
                  list( test = tnames,
                        what = c("Chisq","df","Pval") ) ) )
# Smaller versions of the arrys to be used inside the function that
# does the calculations for each response variable
std <- NArray( dimnames(Std)[-1] )
eff <- NArray( dimnames(Eff)[-1] )
tst <- NArray( dimnames(Tst)[-1] )


###################################################
### code chunk number 90: sec-ins.rnw:413-483
###################################################
# Function to do PPplots for residuals
PPplot <-
function( x, xl )
{
n <- length(x)
plot( pnorm( sort(x/sd(x)) ), (1:n-0.5)/n,
      pch=16, cex=0.4,
      xlim=0:1, ylim=0:1,
      xlab="", ylab="",
      xaxt="n", yaxt="n" )
text( 0, 1, xl, adj=c(0,1) )
abline( 0, 1 )
}

# A function that returns NA instead of -Inf for 0 argument
logI <- function(x) ifelse(x>0,log(x),NA)

# Function that makes all analyses for one response variable
ana.fun <-
function( var, log.tr=FALSE, respl=FALSE )
{
# We need a single name of the response variable
AD$Y <- if(log.tr) logI(AD[,var]) else AD[,var]
# Fit models to total and PP part of dataset:
# r: Reduced model
# a: Primary model
# c: Confounder-expanded model
ITr <- lmer( Y ~ igr + igr:factor(visit) - 1 + (1|subjid), data = AD )
ITa <- lmer( Y ~ igr + igr:factor(visit) - 1 + sdc + over.65 + pre.ins
                                             + (1|subjid), data = AD )
ITc <- lmer( Y ~ igr + igr:factor(visit) - 1 + sdc + over.65 + pre.ins
                                             + sex + statin + gad.pos + cvd
                                             + (1|subjid), data = AD )
PPa <- update( ITa, data = subset( AD, pp ) )
PPr <- update( ITr, data = subset( AD, pp ) )
PPc <- update( ITc, data = subset( AD, pp ) )

# Plots if specified - but only for ITa analysis - note the S4 way of
# extracting slots from the objects, does not work for medMod objects
if( respl )
{
PPplot( ranef(ITa), paste(var,if(log.tr)"(log-tr)","\nptt") )
PPplot( resid(ITa), paste(var,if(log.tr)"(log-tr)","\nres") )
}

# Fixed effects
eff["IT","Prim",,,] <- eFUN( ITa )
eff["IT","Rest",,,] <- eFUN( ITr )
eff["IT","Conf",,,] <- eFUN( ITc )
eff["PP","Prim",,,] <- eFUN( PPa )
eff["PP","Rest",,,] <- eFUN( PPr )
eff["PP","Conf",,,] <- eFUN( PPc )
if( log.tr ) eff <- exp( eff )
# Variance components
std["IT","Prim",] <- sFUN( ITa )
std["IT","Rest",] <- sFUN( ITr )
std["IT","Conf",] <- sFUN( ITc )
std["PP","Prim",] <- sFUN( PPa )
std["PP","Rest",] <- sFUN( PPr )
std["PP","Conf",] <- sFUN( PPc )
# Tests
tst["IT","Prim",,] <- tFUN( ITa )
tst["IT","Rest",,] <- tFUN( ITr )
tst["IT","Conf",,] <- tFUN( ITc )
tst["PP","Prim",,] <- tFUN( PPa )
tst["PP","Rest",,] <- tFUN( PPr )
tst["PP","Conf",,] <- tFUN( PPc )
# Return all in a list
invisible( list( eff=eff, std=std, tst=tst ) )
}


###################################################
### code chunk number 91: PP
###################################################
# par( mfrow=c(8,6), mar=c(0,0,0,0), omi=c(3,3,1,1)/4 )
for( i in 1:length(vnam) )
   {
   res <- ana.fun( vnam[i], ltrf[i], respl=FALSE )
cat( vnam[i], " done\n" )
flush.console()
   Eff[i,,,,,] <- res$eff
   Std[i,,,]   <- res$std
   Tst[i,,,,]  <- res$tst
#   mtext( expression(Phi^-1*"[sort(std. res.)]"),
#          side=1, line=2, outer=TRUE )
#   mtext( "Uniform[0,1]", side=2, line=2, outer=TRUE )
   }


###################################################
### code chunk number 92: sec-ins.rnw:528-549
###################################################
dimnames( Eff )[[6]][4:6] <- c("Asp-Det","Asp-Bip","Det-Bip")
round( ftable( Eff[,"IT","Prim",,,1:3  ], col.vars=4:3 ), 3 )
round( ftable( Eff[,"IT","Prim",,,1:3+3], col.vars=4:3 ), 3 )
round( ftable( Eff[,"IT","Rest",,,1:3  ], col.vars=4:3 ), 3 )
round( ftable( Eff[,"IT","Rest",,,1:3+3], col.vars=4:3 ), 3 )
round( ftable( Eff[,"IT","Conf",,,1:3  ], col.vars=4:3 ), 3 )
round( ftable( Eff[,"IT","Conf",,,1:3+3], col.vars=4:3 ), 3 )
round( ftable( Eff[,"PP","Prim",,,1:3  ], col.vars=4:3 ), 3 )
round( ftable( Eff[,"PP","Prim",,,1:3+3], col.vars=4:3 ), 3 )
round( ftable( Eff[,"PP","Rest",,,1:3  ], col.vars=4:3 ), 3 )
round( ftable( Eff[,"PP","Rest",,,1:3+3], col.vars=4:3 ), 3 )
round( ftable( Eff[,"PP","Conf",,,1:3  ], col.vars=4:3 ), 3 )
round( ftable( Eff[,"PP","Conf",,,1:3+3], col.vars=4:3 ), 3 )
#for( vv in dimnames(Eff)[[1]] )
#   {
#   cat( "-------------------------------------\n", vv, "\n" )
#   dec <- 3 - ( max(Eff[vv,,,,,1:3])>100  ) -
#              ( max(Eff[vv,,,,,1:3])>1000 )
#   print( round( ftable( Eff[vv,,,,,1:3  ], col.vars=5:4 ), dec ) )
#   print( round( ftable( Eff[vv,,,,,1:3+3], col.vars=5:4 ), dec ) )
#   }


###################################################
### code chunk number 93: sec-ins.rnw:552-557
###################################################
for( vv in dimnames(Tst)[[1]] )
   {
   cat( "--------------------------------\n", vv, "\n" )
   print( round( ftable( Tst[vv,,,,"Pval"], row.vars=3 ), 4 ) )
   }


###################################################
### code chunk number 94: sec-ins.rnw:565-618
###################################################
source( "cnr.R" )
cnr
pleff <-
function( var,
         ytxt = var,
           yl = NULL, # yaxis limits
           yt = NULL, # yaxis major ticks
           ym = NULL, # yaxis minor ticks
           hl = ym )
{
# Get the estimates, but only for the times where they are, otherwise
# the the points will not be connected (a feature)
var <- dimnames(Eff)$var[grep(var,vnam <- dimnames(Eff)$var)]
eff <- Eff[var,"IT","Prim",,,]
pef <- cbind( eff[-8,,ilev[1]],
              eff[-8,,ilev[2]],
              eff[-8,,ilev[3]] )
pef <- pef[complete.cases(pef),]
tim <- (as.numeric( gsub("Vis ","",rownames(pef) ) )-1)*3

if( is.null(yl) ) yl <- range(pef)*c(0.9,1.2)
matplot( tim, pef, type="n",
         xlab="Months since trial entry", xaxt="n", xlim=c(0,19.5),
         ylab=ytxt, ylim=yl, yaxs="i", yaxt="n", bty="n" )
if( !is.null(hl) ) abline( h=hl, col=gray(0.8) )
axis( side=1, at=0:6*3, col=clr[4] )
axis( side=2, if(!is.null(yt)) at=yt, col=clr[4] )
axis( side=2, if(!is.null(ym)) at=ym, labels=rep("",length(ym)), col=clr[4] )
text( cnr(c(98,98,98),c(98,92,86)), iN, col=iclr[1:3],
      font=2, adj=c(1,1), cex=1.2 )
if( ltrf[match(var,vnam)] ) { text( cnr(98,74)[[1]], cnr(98,74)[[2]],
                                  "% change\nfrom baseline", col=clr[4], adj=c(1,0.5) )
} else {                       text( cnr(98,74)[[1]], cnr(98,74)[[2]],
                                    "Change\nfrom baseline", col=clr[4], adj=c(1,0.5) ) }
EF <- eff[8,"Est",ilev]
# It might not be vist 7 that is the last non-missing, so find that
ends <- eff[1:7,"Est",ilev]
ends <- ends[max(which(!is.na(ends[,1]))),]
# dv <- ends-mean(ends) ; dv <- dv/abs(dv)
oe <- order( ends )
eo <- c(mean(ends),ends[oe])
de <- diff(eo)
nd <- de[-1]/abs(de[-1])*pmax( abs(de[-1]), diff(yl)/30 )
de <- c(de[1]+nd[1]-de[2],nd)
eo <- eo[1] + cumsum(de)
ends <- eo[order(oe)]
if( ltrf[match(var,vnam)] ) EF <- (EF-1)*100
# print( c(EF, ends ) )
text( c(20,20,20), ends,
      paste( formatC( EF, format="f", digits=1 ), sep="" ),
      col=iclr[1:3], adj=1, font=2 )
matlines( tim, pef, type="l", lwd=c(5,1,1), lty=1, col=rep(iclr[1:3],each=3) )
}


###################################################
### code chunk number 95: tr-hba
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 )
pleff( "hba1c", expression("Hb"*A[1][c]*" (%)  [Mean (95% CI)]"),
       yl = c(7,10),
       yt = 7:10,
       ym = seq(7,10,0.5) )


###################################################
### code chunk number 96: tr-trig
###################################################
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 )
pleff("trig","Triglycerides (mmol/l)  [Mean (95% CI)]",
              yl = c(1,2.5),       # yaxis limits
              yt = 0:4,            # yaxis major ticks
              ym = seq(1,3,0.25) ) # yaxis minor ticks


###################################################
### code chunk number 97: first4
###################################################
tmpl <-
function()
{
par( mfrow=c(2,2),  mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 )
# HbA1c:
pleff( "hba1c",
       expression("Hb"*A[1][c]*" (%)  [Mean (95% CI)]"),
       yl = c(7,10),
       yt = 7:10,
       ym = seq(7,10,0.5) )
text( cnr(2,98), "a", adj=c(0,1), font=2, cex=1.5 )
# Gluc
pleff("gluc","Fasting p-glucose (mmol/L) [Mean (95% CI)]",
              yl = c(7,12),      # yaxis limits
              yt = seq(7,13,2),  # yaxis major ticks
              ym = seq(7,13,1) ) # yaxis minor ticks
text( cnr(2,98), "b", adj=c(0,1), font=2, cex=1.5 )
# Insulin dose
pleff("ipkg","Insulin dose (IU/day/kg)  [Mean (95% CI)]",
              yl = c(0.0,2.4),        # yaxis limits
              yt = seq(0.0,2.0,1.0),  # yaxis major ticks
              ym = seq(0.0,2.4,0.2) ) # yaxis minor ticks
text( cnr(2,98), "c", adj=c(0,1), font=2, cex=1.5 )
# Weight:
pleff("weight","Weight (kg)  [Mean (95% CI)]",
              yl = c(90,120),      # yaxis limits
              yt = seq(90,120,10), # yaxis major ticks
              ym = seq(90,120,5) ) # yaxis minor ticks
text( cnr(2,98), "d", adj=c(0,1), font=2, cex=1.5 )
}
tmpl()


###################################################
### code chunk number 98: sec-ins.rnw:682-685
###################################################
win.metafile( "./results/first4-ins.emf", height=9, width=8 )
tmpl()
dev.off()