### R code from vignette source 'EASD.rnw' ################################################### ### code chunk number 1: EASD.rnw:31-39 ################################################### imt <- read.csv2("./data/imt.csv") str( imt ) names( imt )[c(5,9,12)] for( i in c(5,9,12) ) imt[,i] <- as.numeric(as.character(imt[,i])) imt$datescanned <- as.Date( substr(imt$datescanned,1,10) ) str(imt) addmargins( with( imt, table(visit,cca_side) ) ) with( imt, table( table( subjid ) ) ) ################################################### ### code chunk number 2: EASD.rnw:42-44 ################################################### rnd <- read.csv2("./data/PTNR_MetforminCode_InsulinTypeCode.csv" ) str( rnd ) ################################################### ### code chunk number 3: EASD.rnw:48-53 ################################################### dim( imt ) rimt <- merge( imt, rnd, by.x="subjid", by.y="PTNR", all.y=TRUE ) dim( rimt ) str( rimt ) with( rimt, table(table(subjid)) ) ################################################### ### code chunk number 4: EASD.rnw:57-61 ################################################### tt <- with( imt, table( subjid ) ) has4i <- names( tt[tt==4] ) isrnd <- rnd$PTNR setdiff( has4i, isrnd ) ################################################### ### code chunk number 5: pairs ################################################### with( rimt, pairs( cbind( avg= as.numeric(as.character(mean_fimtavg)), min=mean_fimtmin, max=mean_fimtmax, gns=(mean_fimtmin+mean_fimtmax)/2 ), gap=0, pch=16, cex=0.5 ) ) ################################################### ### code chunk number 6: EASD.rnw:75-94 ################################################### what <- 5:22 names( rimt )[what] mimt <- aggregate( rimt[,what], list( subjid=rimt$subjid, visit=rimt$visit), FUN=mean, na.rm=TRUE ) wh <- c(1,2,5:12,18:20) names( imt )[wh] Limt <- subset( rimt, cca_side=="L" )[,wh] Rimt <- subset( rimt, cca_side=="R" )[,wh] names( Limt ) <- gsub( "mean_", "", names( Limt ) ) names( Rimt ) <- gsub( "mean_", "", names( Rimt ) ) names( Limt )[-(1:2)] <- paste( "left_" , names( Limt )[-(1:2)], sep="" ) names( Rimt )[-(1:2)] <- paste( "right_", names( Rimt )[-(1:2)], sep="" ) mimt <- merge( merge( mimt, Limt, all.x=TRUE ), Rimt, all.x=TRUE ) head( mimt ) subset( rimt, subjid==10001 ) subset( mimt, subjid==10001 ) ################################################### ### code chunk number 7: EASD.rnw:97-100 ################################################### st <- read.csv2( "./data/strata.csv" )[,-2] names( st )[2:4] <- c("o.65","pre.ins","SDC") mimt <- merge( mimt, st, by.x="subjid", by.y="ptnr" ) ################################################### ### code chunk number 8: EASD.rnw:104-105 ################################################### with( mimt, ftable(addmargins(table(Met=MetforminCode, Ins=InsulinTypeCode, visit ))) ) ################################################### ### code chunk number 9: EASD.rnw:109-114 ################################################### save( mimt, file="./data/mimt.Rda" ) write.foreign( mimt, datafile = "./sas/mimt.dat", codefile = "./sas/mimt.sas", package = "SAS" ) ################################################### ### code chunk number 10: EASD.rnw:122-130 ################################################### library( lme4 ) library( Epi ) load( file="./data/mimt.Rda" ) mimt <- transform( mimt, subjid = factor(subjid), met = factor(MetforminCode,labels=c("Plc","Met")), ins = factor(InsulinTypeCode,labels=c("NR+Lv","Lv","NM30")) ) with( mimt, table(MetforminCode, met )) with( mimt, table(InsulinTypeCode,ins)) ################################################### ### code chunk number 11: EASD.rnw:169-189 ################################################### m0 <- lmer( mean_fimtavg ~ -1 + visit:met + o.65 + pre.ins + SDC +(1|subjid), data = mimt ) summary( m0 ) round( ee <- ci.lin( m0, subset="visit" ), 4 ) CM <- rbind(diag(4),c(0,0,-1,1), c(-1,1,0,0), c(1,-1,-1,1)) row.names(CM) <- c("Plc 1a","Plc 7a", "Met 1a","Met 7a", "Met 7a-1a", "Plc 7a-1a", "dMet - dPlc" ) colnames(CM) <- row.names(ee) CM round( e0 <- ci.lin( m0, subset="visit", ctr.mat=CM ), 4 ) elab <- c("Metformin 18m - baseline", "Placebo 18m - baseline", "Metformin vs. Placebo change") rownames( e0 )[5:7] <- elab round( e0, 4 ) ################################################### ### code chunk number 12: EASD.rnw:191-205 ################################################### eclr <- c("green","yellow","orange") win.metafile( "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="white", col.lab="white" ) plotEst( e0[5:7,c(1,5,6)], lwd=7,cex=1.5, xlab="", col.txt="transparent", col=eclr, restore.par=FALSE ) for( i in 1:3 ) axis( side=2, at=4-i, labels=elab[i], col="transparent", col.axis= eclr[i], las=1 ) abline(v=0,col="white") axis( side=1, col="white" ) mtext( "Carotid intima-media thickness (mm)", side=1, line=3/1.6, col="white" ) dev.off() ################################################### ### code chunk number 13: forest1 ################################################### par( mar=c(3,1,1,1), mgp=c(3,1,0)/1.6 ) plotEst( e0[5:7,c(1,5,6)], lwd=3, cex=1.5, vref=0, xlab="Carotid intima-media thickness (mm)", txt=elab ) ################################################### ### code chunk number 14: EASD.rnw:245-256 ################################################### names( mimt ) wimt <- reshape( mimt[,c("mean_fimtavg","visit","subjid", "MetforminCode","o.65","pre.ins","SDC")], direction = "wide", v.names = "mean_fimtavg", timevar = "visit", idvar = "subjid" ) subset( mimt[,c("mean_fimtavg","visit","subjid", "MetforminCode","o.65","pre.ins","SDC")], subjid=="10002" ) subset( wimt, subjid=="10002" ) ################################################### ### code chunk number 15: EASD.rnw:260-268 ################################################### mf <- lm( mean_fimtavg.7a ~ mean_fimtavg.1a + MetforminCode + o.65 + pre.ins + SDC, data = wimt ) summary( mf ) mm <- lm( mean_fimtavg.7a ~ mean_fimtavg.1a + MetforminCode, 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 16: EASD.rnw:271-272 ################################################### round( e0, 4 ) ################################################### ### code chunk number 17: EASD.rnw:276-280 ################################################### round( rbind( "Simple"=cf[2,], "Contrl"=cm[2,], "RanEff"=e0[7,] ), 4 ) ################################################### ### code chunk number 18: summ ################################################### options( width=90 ) summary( wimt ) ################################################### ### code chunk number 19: EASD.rnw:306-332 ################################################### library( mice ) library( foreign ) llch <- read.ssd( "./data", "mi" ) names( llch ) <- tolower( gsub("_",".",names(llch)) ) mimt <- merge( wimt, llch ) names( mimt ) <- gsub( "mean_", "", names(mimt) ) mimt <- transform( mimt, sex = factor( sex , labels=c("F","M") ), # code 0/1 ! ami = factor( ami , labels=c("No","Yes") ), heartsur = factor( heartsur, labels=c("No","Yes") ), apopl = factor( apopl , labels=c("No","Yes") ), tci = factor( tci , labels=c("No","Yes") ), sdc = factor( sdc , labels=c("no","yes") ), o.65 = factor( o.65 , labels=c("no","yes") ), pre.ins = factor( pre.ins , labels=c("no","yes") ), metformi = factor( metformi, labels=c("Plc","Met") ) ) mimt <- mimt[,c(1,42,43,41,3,4,6:40,44:47)] cbind( 1:ncol(mimt), t( apply( mimt, 2, function(x) c(sum(is.na(x)), length(unique(x))) ) ) ) ftable( addmargins( xtabs( cbind( "Baseline"=!is.na(fimtavg.1a), "FU18mth"=!is.na(fimtavg.7a), "Diff."=!is.na(fimtavg.1a)- !is.na(fimtavg.7a), "MissingFU"= is.na(fimtavg.7a) ) ~ sex + metformi, data=mimt ), margin=1:2 ) ) ################################################### ### code chunk number 20: EASD.rnw:336-338 ################################################### zz <- md.pattern( mimt ) t( zz ) ################################################### ### code chunk number 21: pairs1 ################################################### pairs( mimt[,c(7,8,15,16:24)], gap=0, pch=16, cex=0.2 ) ################################################### ### code chunk number 22: pairs2 ################################################### pairs( mimt[,25:36], gap=0, pch=16, cex=0.2 ) ################################################### ### code chunk number 23: pairs3 ################################################### pairs( mimt[,37:45], gap=0, pch=16, cex=0.2 ) ################################################### ### code chunk number 24: mimimt ################################################### mimimt <- transform( mimt, l.alb = log( avgnatua ), l.nsev = log( nonsever ) )[,-c(14,16,26:30,42:45)][c(1:15,36,16:35)] names( mimimt ) save( mimimt, file="./data/mimimt.Rda" ) ################################################### ### code chunk number 25: ppairs1 ################################################### pairs( mimimt[,7:21], gap=0, pch="." ) ################################################### ### code chunk number 26: ppairs2 ################################################### pairs( mimimt[,22:36], gap=0, pch="." ) ################################################### ### code chunk number 27: mice ################################################### set.seed( 876528358 ) system.time( imp <- mice( mimimt[,2:36], m=50 ) ) save( imp, file="./data/imp" ) ################################################### ### code chunk number 28: imp ################################################### class( imp ) i.mf <- with( imp, lm( fimtavg.7a ~ fimtavg.1a + metformi + o.65 + pre.ins + sdc ) ) round( smf <- summary( pool(i.mf) ), 4 ) i.mm <- with( imp, lm( fimtavg.7a ~ fimtavg.1a + metformi ) ) round( smm <- summary( pool( i.mm ) ), 4 ) ################################################### ### code chunk number 29: forest2 ################################################### par( mar=c(3,1,1,1), mgp=c(3,1,0)/1.6 ) plotEst( rbind(e0[5:6,c(1,5,6)],smm[3,c("est","lo 95","hi 95")]), txt=rownames(e0)[5:7], lwd=7,cex=1.5, restore.par=FALSE ) abline( v=0 ) axis( side=1 ) mtext( "Carotid intima-media thickness (mm)", side=1, line=3/1.6 ) ################################################### ### code chunk number 30: EASD.rnw:407-429 ################################################### eclr <- c("green","yellow","orange") win.metafile( "forest2.emf", width=10, height=6, pointsize=20 ) par( mar=c(3,1,1,1), mgp=c(3,1,0)/1.6, bg="transparent", col.axis="white", col.lab="white" ) plotEst( rbind(e0[5:6,c(1,5,6)],smm[3,c("est","lo 95","hi 95")]), ylim=c(0,3), txt=rownames(e0)[5:7], lwd=7,cex=1.5, xlab="", col.txt="transparent", col=eclr, restore.par=FALSE ) for( i in 1:3 ) axis( side=2, at=4-i, labels=elab[i], col="transparent", col.axis= eclr[i], las=1 ) abline(v=0,col="white") axis( side=1, col="white" ) text( c(e0[5:6,1],smm[3,"est"]), 0.7+2:0, paste( "P =", formatC( c(e0[5:6,"P"],smm[3,"Pr(>|t|)"]), format="f", digits=3 ) ), col="white" ) text( -0.03, 0.2, "Improvement" , adj=0, col="white", cex=0.8 ) text( 0.03, 0.2, "Deterioration", adj=1, col="white", cex=0.8 ) mtext( "Carotid intima-media thickness (mm)", side=1, line=3/1.6, col="white" ) dev.off() ################################################### ### code chunk number 31: EASD.rnw:445-456 ################################################### options( width=110 ) library( Epi ) library( foreign ) library( lme4 ) hba <- read.xport( "./data/hba.xpt" ) names( hba ) <- tolower( names( hba ) ) names( hba ) <- gsub( "_", ".", names( hba ) ) hba$grp <- factor( hba$metformi, labels=c("plc","met") ) names( hba )[2] <- "inspw" str( hba ) summary( hba ) ################################################### ### code chunk number 32: EASD.rnw:459-461 ################################################### ( tt <- with( hba, table(table(subjid)) ) ) round( 100*tt/sum(tt), 1 ) ################################################### ### code chunk number 33: EASD.rnw:466-472 ################################################### vv <- with( hba, table(subjid) ) cplv <- names( vv[vv==7] ) hbc <- subset( hba, subjid %in% cplv ) with( hba, table(table(subjid)) ) with( hbc, table(table(subjid)) ) names( hbc ) ################################################### ### code chunk number 34: lmer 1 ################################################### mc <- lmer( hba1c ~ grp + grp:factor(visit) - 1 + factor(sdc) + factor(above.si) + factor(prior.in) + (1|subjid), data=hbc ) ma <- update( mc, data=hba ) ( eM <- cbind(1,rbind(0,diag(6))) ) ci.lin( mc ) c.plc <- ci.lin( mc, ctr.mat=eM, subset="plc" )[,-(2:4)] c.met <- ci.lin( mc, ctr.mat=eM, subset="met" )[,-(2:4)] c.dif <- ci.lin( mc, subset=c("met","plc"), ctr.mat=cbind(-eM,eM) ) round( c.dif, 3 ) ci.lin( ma ) a.plc <- ci.lin( ma, ctr.mat=cbind(1,rbind(0,diag(6))), subset="plc" )[,-(2:4)] a.met <- ci.lin( ma, ctr.mat=cbind(1,rbind(0,diag(6))), subset="met" )[,-(2:4)] a.dif <- ci.lin( ma, subset=c("met","plc"), ctr.mat=cbind(-eM,eM) ) round( a.dif, 3 ) P.all <- Wald( ma, H0=rep(0,7), subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) P.last<- Wald( ma, H0=0, subset=c("met","plc"), ctr.mat=cbind(eM,-eM)[7,,drop=F] ) win.metafile( "hba1c.emf", width=10, height=6, pointsize=20 ) eclr <- c("yellow","green","orange") par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1, bg=#"black", "transparent", col.axis="white", col.lab="white", bty="n" ) matplot( 0:6*3, cbind( a.plc, a.met ), 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(eclr[1:2],each=3) ) # abline( h=seq(7,9.5,0.5), col=gray(0.5) ) matlines( 0:6*3, cbind( a.plc, a.met ), type="l", lwd=c(5,1,1), lty=1, col=rep(eclr[1:2],each=3) ) axis( side=1, at=0:6*3, col="white" ) axis( side=2, at=7:9, col="white" ) text( c(20,20), c(9.1,9.3), c("Placebo","Metformin"), col=eclr[1:2], font=2, adj=1, cex=1 ) text( 20, 8.8, substitute( "Equal Hb"*A[1][c]*" 0-18 mth: P ="*pval, list(pval = formatC( P.all[3], format="f", digits=3 )) ), col="white", adj=1 ) text( 20, 8.6,substitute( "Equal Hb"*A[1][c]*" at end: P ="*pval, list(pval = formatC( P.last[3], format="f", digits=3 )) ), col="white", adj=1 ) # text( 19.2, 8.55, "HbA1c change\n from baseline", col="white", adj=c(1,1) ) text( c(20,20), c(a.plc[7,1],a.met[7,1]), paste( formatC( c(a.plc[7,1],a.met[7,1])-c(a.plc[1,1],a.met[1,1]), format="f", digits=1 ), "%", sep="" ), col=eclr[1:2], adj=1 ) dev.off() ################################################### ### code chunk number 35: EASD.rnw:533-541 ################################################### vv <- with( hba, table(subjid) ) cplv <- names( vv[vv==7] ) hbc <- subset( hba, subjid %in% cplv ) with( hba, table(table(subjid)) ) with( hbc, table(table(subjid)) ) names( hbc ) with( hba, tapply( inspw, visit, function(x) sum( !is.na(x) ) ) ) with( hba, tapply( inspw, visit, function(x) sum( is.na(x) ) ) ) ################################################### ### code chunk number 36: lmer 1 ################################################### mc <- lmer( inspw ~ grp + grp:factor(visit) - 1 + factor(sdc) + factor(above.si) + factor(prior.in) + (1|subjid), data = subset( hba, visit!=7) ) round( ci.lin( mc ), 3 ) ( eM <- cbind(1,rbind(0,diag(5))) ) c.plc <- ci.lin( mc, ctr.mat=eM, subset="plc" )[,-(2:4)] c.met <- ci.lin( mc, ctr.mat=eM, subset="met" )[,-(2:4)] c.dif <- ci.lin( mc, subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) round( c.dif, 3 ) P.all <- Wald( mc, H0=rep(0,6), subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) P.last<- Wald( mc, H0=0, subset=c("met","plc"), ctr.mat=cbind(eM,-eM)[6,,drop=F] ) ################################################### ### code chunk number 37: EASD.rnw:562-563 ################################################### ( IL <- with( hba, tapply( insdose2, list( metformi, visit), mean ) )[,6] ) ################################################### ### code chunk number 38: EASD.rnw:565-592 ################################################### eclr <- c("yellow","green","orange") win.metafile( "inspw.emf", width=10, height=6, pointsize=20 ) par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1, bg=#"black", "transparent", col.axis="white", col.lab="white", bty="n" ) matplot( 0:5*3, cbind( c.plc, c.met ), xlab="Months since trial entry", xaxt="n", xlim=c(0,18), ylim=c(0,1.3), yaxs="i", yaxt="n", ylab=expression( "Insulin (IU/kg/day) [Mean (95% CI)]" ), type="n", lwd=c(4,1,1), lty=1, col=rep(eclr[1:2],each=3) ) abline( h=seq(7,9.5,0.5), col=gray(0.5) ) matlines( 0:5*3, cbind( c.plc, c.met ), type="l", lwd=c(5,1,1), lty=1, col=rep(eclr[1:2],each=3) ) axis( side=1, at=0:6*3, col="white" ) axis( side=2, at=seq(0,1.2,0.2), col="white" ) text( c(18,18), c(0.35,0.45), c("Placebo","Metformin"), col=eclr[1:2], font=2, adj=1, cex=1 ) text( 18, 0.2, paste( "Equal insulin dose 0-15 mth: P =", formatC( P.all[3], format="f", digits=3 ) ), col="white", adj=1 ) text( 18, 0.1, paste( "Equal insulin dose at end: P =", formatC( P.last[3], format="f", digits=3 ) ), col="white", adj=1 ) text( c(18,18), c(c.plc[6,1],c.met[6,1]), paste( round(IL), "IU/day"), col=eclr[1:2], adj=1 ) dev.off() ################################################### ### code chunk number 39: EASD.rnw:599-602 ################################################### names( hba ) with( hba, tapply( weight, visit, function(x) sum( !is.na(x) ) ) ) with( hba, tapply( weight, visit, function(x) sum( is.na(x) ) ) ) ################################################### ### code chunk number 40: lmer 1 ################################################### mc <- lmer( weight ~ grp + grp:factor(visit) - 1 + factor(sdc) + factor(above.si) + factor(prior.in) + (1|subjid), data = hba ) round( ci.lin( mc ), 3 ) ( eM <- cbind(1,rbind(0,diag(6))) ) c.plc <- ci.lin( mc, ctr.mat=eM, subset="plc" )[,-(2:4)] c.met <- ci.lin( mc, ctr.mat=eM, subset="met" )[,-(2:4)] c.dif <- ci.lin( mc, subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) round( c.dif, 3 ) P.all <- Wald( mc, H0=rep(0,7), subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) P.last<- Wald( mc, H0=0, subset=c("met","plc"), ctr.mat=cbind(eM,-eM)[7,,drop=F] ) ################################################### ### code chunk number 41: EASD.rnw:624-650 ################################################### win.metafile( "weight.emf", width=10, height=6, pointsize=24 ) eclr <- c("yellow","green","orange") par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1, bg=#"black", "transparent", col.axis="white", col.lab="white", bty="n" ) matplot( 0:6*3, cbind( c.plc, c.met ), xlab="Months since trial entry", xaxt="n", xlim=c(0,18), ylim=c(90,110), yaxs="i", yaxt="n", ylab=expression( "Weight (kg) [Mean (95% CI)]" ), type="n", lwd=c(4,1,1), lty=1, col=rep(eclr[1:2],each=3) ) # abline( h=seq(7,9.5,0.5), col=gray(0.5) ) matlines( 0:6*3, cbind( c.plc, c.met ), type="l", lwd=c(5,1,1), lty=1, col=rep(eclr[1:2],each=3) ) axis( side=1, at=0:6*3, col="white" ) axis( side=2, at=seq(90,110,5), col="white" ) text( c(17,17), c(0.25,0.35), c("Placebo","Metformin"), col=eclr[1:2], font=2, adj=1, cex=1 ) text( 17, 0.15, paste( "Equal insulin dose 0-15 mth: P =", formatC( P.all[3], format="f", digits=3 ) ), col="white", adj=1 ) text( 17, 0.1, paste( "Equal insulin dose at end: P =", formatC( P.last[3], format="f", digits=3 ) ), col="white", adj=1 ) dev.off() ################################################### ### code chunk number 42: EASD.rnw:658-667 ################################################### round( ci.lin( mc ), 3 ) ( eM <- cbind(0,rbind(0,diag(6))) ) c.plc <- ci.lin( mc, ctr.mat=eM, subset="plc" )[,-(2:4)] c.met <- ci.lin( mc, ctr.mat=eM, subset="met" )[,-(2:4)] c.dif <- ci.lin( mc, subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) round( c.dif, 3 ) P.all <- Wald( mc, H0=rep(0,6), subset=c("met","plc"), ctr.mat=cbind(eM,-eM)[-1,] ) P.last<- Wald( mc, H0=0, subset=c("met","plc"), ctr.mat=cbind(eM,-eM)[7,,drop=F] ) ################################################### ### code chunk number 43: EASD.rnw:669-694 ################################################### win.metafile( "wchg.emf", width=10, height=6, pointsize=20 ) eclr <- c("yellow","green","orange") par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1, bg=#"black", "transparent", col.axis="white", col.lab="white", bty="n" ) matplot( 0:6*3, cbind( c.plc, c.met ), xlab="Months since trial entry", xaxt="n", xlim=c(0,18), ylim=c(0,6.2), yaxt="n", ylab=expression( "Weight change (kg) [Mean (95% CI)]" ), type="n", lwd=c(4,1,1), lty=1, col=rep(eclr[1:2],each=3) ) # abline( h=seq(7,9.5,0.5), col=gray(0.5) ) matlines( 0:6*3, cbind( c.plc, c.met ), type="l", lwd=c(5,1,1), lty=1, col=rep(eclr[1:2],each=3) ) axis( side=1, at=0:6*3, col="white" ) axis( side=2, at=0:6, col="white" ) text( c(0.5,0.5), c(5.7,6.2), c("Placebo","Metformin"), col=eclr[1:2], font=2, adj=0, cex=1 ) text( 0.5, 5.2, paste( "Equal weight change 0-18 mth: P =", formatC( P.all[3], format="f", digits=3 ) ), col="white", adj=0 ) text( 0.5, 4.7, paste( "Equal weight change end: P =", formatC( P.last[3], format="f", digits=3 ) ), col="white", adj=0 ) dev.off()