1### R code from vignette source 'simLexis' 2### Encoding: UTF-8 3 4################################################### 5### code chunk number 1: simLexis.rnw:23-26 6################################################### 7options( width=90, 8 SweaveHooks=list( fig=function() 9 par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n") ) ) 10 11 12################################################### 13### code chunk number 2: start 14################################################### 15options( width=90 ) 16library( Epi ) 17print( sessionInfo(), l=F ) 18 19 20################################################### 21### code chunk number 3: Lexis 22################################################### 23data(DMlate) 24dml <- Lexis( entry = list(Per=dodm, Age=dodm-dobth, DMdur=0 ), 25 exit = list(Per=dox), 26 exit.status = factor(!is.na(dodth),labels=c("DM","Dead")), 27 data = DMlate ) 28 29 30################################################### 31### code chunk number 4: cut 32################################################### 33dmi <- cutLexis( dml, cut = dml$doins, 34 pre = "DM", 35 new.state = "Ins", 36 new.scale = "t.Ins", 37 split.states = TRUE ) 38summary( dmi, timeScales=T ) 39 40 41################################################### 42### code chunk number 5: boxes 43################################################### 44boxes( dmi, boxpos = list(x=c(20,20,80,80), 45 y=c(80,20,80,20)), 46 scale.R = 1000, show.BE = TRUE ) 47 48 49################################################### 50### code chunk number 6: split 51################################################### 52Si <- splitLexis( dmi, seq(0,20,1/4), "DMdur" ) 53summary( Si ) 54print( subset( Si, lex.id==97 )[,1:10], digits=6 ) 55 56 57################################################### 58### code chunk number 7: knots 59################################################### 60nk <- 5 61( ai.kn <- with( subset(Si,lex.Xst=="Ins" & lex.Cst!=lex.Xst ), 62 quantile( Age+lex.dur , probs=(1:nk-0.5)/nk ) ) ) 63( ad.kn <- with( subset(Si,lex.Xst=="Dead"), 64 quantile( Age+lex.dur , probs=(1:nk-0.5)/nk ) ) ) 65( di.kn <- with( subset(Si,lex.Xst=="Ins" & lex.Cst!=lex.Xst ), 66 c(0,quantile( DMdur+lex.dur, probs=(1:(nk-1))/nk ) )) ) 67( dd.kn <- with( subset(Si,lex.Xst=="Dead"), 68 c(0,quantile( DMdur+lex.dur, probs=(1:(nk-1))/nk ) )) ) 69( ti.kn <- with( subset(Si,lex.Xst=="Dead(Ins)"), 70 c(0,quantile( t.Ins+lex.dur, probs=(1:(nk-1))/nk ) )) ) 71 72 73################################################### 74### code chunk number 8: Poisson 75################################################### 76library( splines ) 77DM.Ins <- glm( (lex.Xst=="Ins") ~ Ns( Age , knots=ai.kn ) + 78 Ns( DMdur, knots=di.kn ) + 79 I(Per-2000) + sex, 80 family=poisson, offset=log(lex.dur), 81 data = subset(Si,lex.Cst=="DM") ) 82ci.exp( DM.Ins ) 83class( DM.Ins ) 84 85 86################################################### 87### code chunk number 9: simLexis.rnw:282-288 88################################################### 89DM.Ins <- glm.Lexis( Si, from = "DM", to = "Ins", 90 formula = ~ Ns( Age , knots=ai.kn ) + 91 Ns( DMdur, knots=di.kn ) + 92 I(Per-2000) + sex ) 93ci.exp( DM.Ins ) 94class( DM.Ins ) 95 96 97################################################### 98### code chunk number 10: simLexis.rnw:293-302 99################################################### 100DM.Dead <- glm.Lexis( Si, from = "DM", to = "Dead", 101 formula = ~ Ns( Age , knots=ad.kn ) + 102 Ns( DMdur, knots=dd.kn ) + 103 I(Per-2000) + sex ) 104Ins.Dead <- glm.Lexis( Si, from = "Ins", 105 formula = ~ Ns( Age , knots=ad.kn ) + 106 Ns( DMdur, knots=dd.kn ) + 107 Ns( t.Ins, knots=ti.kn ) + 108 I(Per-2000) + sex ) 109 110 111################################################### 112### code chunk number 11: prop-haz 113################################################### 114All.Dead <- glm.Lexis( Si, to = c("Dead(Ins)","Dead"), 115 formula = ~ Ns( Age , knots=ad.kn ) + 116 Ns( DMdur, knots=dd.kn ) + 117 lex.Cst + 118 I(Per-2000) + sex ) 119round( ci.exp( All.Dead ), 3 ) 120 121 122################################################### 123### code chunk number 12: get-dev 124################################################### 125what <- c("null.deviance","df.null","deviance","df.residual") 126( rD <- unlist( DM.Dead[what] ) ) 127( rI <- unlist( Ins.Dead[what] ) ) 128( rA <- unlist( All.Dead[what] ) ) 129round( c( dd <- rA-(rI+rD), "pVal"=1-pchisq(dd[3],dd[4]+1) ), 3 ) 130 131 132################################################### 133### code chunk number 13: pr-array 134################################################### 135pr.rates <- NArray( list( DMdur = seq(0,12,0.1), 136 DMage = 4:7*10, 137 r.Ins = c(NA,0,2,5), 138 model = c("DM/Ins","All"), 139 what = c("rate","lo","hi") ) ) 140str( pr.rates ) 141 142 143################################################### 144### code chunk number 14: mknd 145################################################### 146nd <- data.frame( DMdur = as.numeric( dimnames(pr.rates)[[1]] ), 147 lex.Cst = factor( 1, levels=1:4, 148 labels=levels(Si$lex.Cst) ), 149 sex = factor( 1, levels=1:2, labels=c("M","F")) ) 150 151 152################################################### 153### code chunk number 15: make-pred 154################################################### 155for( ia in dimnames(pr.rates)[[2]] ) 156 { 157dnew <- transform( nd, Age = as.numeric(ia)+DMdur, 158 Per = 1998+DMdur ) 159pr.rates[,ia,1,"DM/Ins",] <- ci.pred( DM.Dead, newdata = dnew ) 160pr.rates[,ia,1,"All" ,] <- ci.pred( All.Dead, newdata = dnew ) 161for( ii in dimnames(pr.rates)[[3]][-1] ) 162 { 163dnew = transform( dnew, lex.Cst = factor( 2, levels=1:4, 164 labels=levels(Si$lex.Cst) ), 165 t.Ins = ifelse( (DMdur-as.numeric(ii)) >= 0, 166 DMdur-as.numeric(ii), NA ) ) 167pr.rates[,ia, ii ,"DM/Ins",] <- ci.pred( Ins.Dead, newdata = dnew ) 168pr.rates[,ia, ii ,"All" ,] <- ci.pred( All.Dead, newdata = dnew ) 169 } 170 } 171 172 173################################################### 174### code chunk number 16: mort-int 175################################################### 176par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1 ) 177plot( NA, xlim=c(40,82), ylim=c(5,300), bty="n", 178 log="y", xlab="Age", ylab="Mortality rate per 1000 PY" ) 179abline( v=seq(40,80,5), h=outer(1:9,10^(0:2),"*"), col=gray(0.8) ) 180for( aa in 4:7*10 ) for( ii in 1:4 ) 181 matshade( aa+as.numeric(dimnames(pr.rates)[[1]]), 182 cbind( pr.rates[,paste(aa),ii,"DM/Ins",], 183 pr.rates[,paste(aa),ii,"All" ,] )*1000, 184 type="l", lty=1, lwd=2, 185 col=c("red","limegreen") ) 186 187 188################################################### 189### code chunk number 17: Tr 190################################################### 191Tr <- list( "DM" = list( "Ins" = DM.Ins, 192 "Dead" = DM.Dead ), 193 "Ins" = list( "Dead(Ins)" = Ins.Dead ) ) 194 195 196################################################### 197### code chunk number 18: make-ini 198################################################### 199str( ini <- Si[NULL,1:9] ) 200 201 202################################################### 203### code chunk number 19: ini-fill 204################################################### 205ini[1:2,"lex.id"] <- 1:2 206ini[1:2,"lex.Cst"] <- "DM" 207ini[1:2,"Per"] <- 1995 208ini[1:2,"Age"] <- 60 209ini[1:2,"DMdur"] <- 5 210ini[1:2,"sex"] <- c("M","F") 211ini 212 213 214################################################### 215### code chunk number 20: simL 216################################################### 217set.seed( 52381764 ) 218Nsim <- 5000 219system.time( simL <- simLexis( Tr, 220 ini, 221 t.range = 12, 222 N = Nsim ) ) 223 224 225################################################### 226### code chunk number 21: sum-simL 227################################################### 228summary( simL, by="sex" ) 229 230 231################################################### 232### code chunk number 22: Tr.p-simP 233################################################### 234Tr.p <- list( "DM" = list( "Ins" = DM.Ins, 235 "Dead" = All.Dead ), 236 "Ins" = list( "Dead(Ins)" = All.Dead ) ) 237system.time( simP <- simLexis( Tr.p, 238 ini, 239 t.range = 12, 240 N = Nsim ) ) 241summary( simP, by="sex" ) 242 243 244################################################### 245### code chunk number 23: Cox-dur 246################################################### 247library( survival ) 248Cox.Dead <- coxph( Surv( DMdur, DMdur+lex.dur, 249 lex.Xst %in% c("Dead(Ins)","Dead")) ~ 250 Ns( Age-DMdur, knots=ad.kn ) + 251 I(lex.Cst=="Ins") + 252 I(Per-2000) + sex, 253 data = Si ) 254round( ci.exp( Cox.Dead ), 3 ) 255 256 257################################################### 258### code chunk number 24: TR.c 259################################################### 260Tr.c <- list( "DM" = list( "Ins" = Tr$DM$Ins, 261 "Dead" = Cox.Dead ), 262 "Ins" = list( "Dead(Ins)" = Cox.Dead ) ) 263system.time( simC <- simLexis( Tr.c, 264 ini, 265 t.range = 12, 266 N = Nsim ) ) 267summary( simC, by="sex" ) 268 269 270################################################### 271### code chunk number 25: nState 272################################################### 273system.time( 274nSt <- nState( subset(simL,sex=="M"), 275 at=seq(0,11,0.2), from=1995, time.scale="Per" ) ) 276nSt[1:10,] 277 278 279################################################### 280### code chunk number 26: pstate0 281################################################### 282pM <- pState( nSt, perm=c(1,2,4,3) ) 283head( pM ) 284par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) 285plot( pM ) 286plot( pM, border="black", col="transparent", lwd=3 ) 287text( rep(as.numeric(rownames(pM)[nrow(pM)-1]),ncol(pM)), 288 pM[nrow(pM),]-diff(c(0,pM[nrow(pM),]))/5, 289 colnames( pM ), adj=1 ) 290box( col="white", lwd=3 ) 291box() 292 293 294################################################### 295### code chunk number 27: pstatex 296################################################### 297clr <- c("limegreen","orange") 298# expand with a lighter version of the two chosen colors 299clx <- c( clr, rgb( t( col2rgb( clr[2:1] )*2 + rep(255,3) ) / 3, max=255 ) ) 300par( mfrow=c(1,2), las=1, mar=c(3,3,4,2), mgp=c(3,1,0)/1.6 ) 301# Men 302plot( pM, col=clx, xlab="Date of FU" ) 303lines( as.numeric(rownames(pM)), pM[,2], lwd=3 ) 304mtext( "60 year old male, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6) ) 305mtext( "Survival curve", side=3, line=1.5, adj=0 ) 306mtext( "DM, no insulin DM, Insulin", side=3, line=0.5, adj=0, col=clr[2] ) 307mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[1] ) 308axis( side=4 ) 309axis( side=4, at=1:19/20, labels=FALSE ) 310axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 ) 311# Women 312pF <- pState( nState( subset(simL,sex=="F"), 313 at=seq(0,11,0.2), 314 from=1995, 315 time.scale="Per" ), 316 perm=c(1,2,4,3) ) 317plot( pF, col=clx, xlab="Date of FU" ) 318lines( as.numeric(rownames(pF)), pF[,2], lwd=3 ) 319mtext( "60 year old female, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6) ) 320mtext( "Survival curve", side=3, line=1.5, adj=0 ) 321mtext( "DM, no insulin DM, Insulin", side=3, line=0.5, adj=0, col=clr[2] ) 322mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[1] ) 323axis( side=4 ) 324axis( side=4, at=1:19/20, labels=FALSE ) 325axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 ) 326 327 328################################################### 329### code chunk number 28: pstatey 330################################################### 331par( mfrow=c(1,2), las=1, mar=c(3,3,4,2), mgp=c(3,1,0)/1.6 ) 332# Men 333pM <- pState( nState( subset(simL,sex=="M"), 334 at=seq(0,11,0.2), 335 from=60, 336 time.scale="Age" ), 337 perm=c(1,2,4,3) ) 338plot( pM, col=clx, xlab="Age" ) 339lines( as.numeric(rownames(pM)), pM[,2], lwd=3 ) 340mtext( "60 year old male, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6) ) 341mtext( "Survival curve", side=3, line=1.5, adj=0 ) 342mtext( "DM, no insulin DM, Insulin", side=3, line=0.5, adj=0, col=clr[2] ) 343mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[1] ) 344axis( side=4 ) 345axis( side=4, at=1:19/20, labels=FALSE ) 346axis( side=4, at=1:19/20, labels=FALSE, tcl=-0.4 ) 347axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 ) 348# Women 349pF <- pState( nState( subset(simL,sex=="F"), 350 at=seq(0,11,0.2), 351 from=60, 352 time.scale="Age" ), 353 perm=c(1,2,4,3) ) 354plot( pF, col=clx, xlab="Age" ) 355lines( as.numeric(rownames(pF)), pF[,2], lwd=3 ) 356mtext( "60 year old female, diagnosed 1990, aged 55", side=3, line=2.5, adj=0, col=gray(0.6) ) 357mtext( "Survival curve", side=3, line=1.5, adj=0 ) 358mtext( "DM, no insulin DM, Insulin", side=3, line=0.5, adj=0, col=clr[2] ) 359mtext( "DM, no insulin", side=3, line=0.5, adj=0, col=clr[1] ) 360axis( side=4 ) 361axis( side=4, at=1:9/10, labels=FALSE ) 362axis( side=4, at=1:19/20, labels=FALSE, tcl=-0.4 ) 363axis( side=4, at=1:99/100, labels=FALSE, tcl=-0.3 ) 364 365 366################################################### 367### code chunk number 29: comp-0 368################################################### 369PrM <- pState( nState( subset(simP,sex=="M"), 370 at=seq(0,11,0.2), 371 from=60, 372 time.scale="Age" ), 373 perm=c(1,2,4,3) ) 374PrF <- pState( nState( subset(simP,sex=="F"), 375 at=seq(0,11,0.2), 376 from=60, 377 time.scale="Age" ), 378 perm=c(1,2,4,3) ) 379CoxM <- pState( nState( subset(simC,sex=="M"), 380 at=seq(0,11,0.2), 381 from=60, 382 time.scale="Age" ), 383 perm=c(1,2,4,3) ) 384CoxF <- pState( nState( subset(simC,sex=="F"), 385 at=seq(0,11,0.2), 386 from=60, 387 time.scale="Age" ), 388 perm=c(1,2,4,3) ) 389 390par( mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(3,1,0)/1.6 ) 391 plot( pM, border="black", col="transparent", lwd=3 ) 392lines( PrM, border="blue" , col="transparent", lwd=3 ) 393lines( CoxM, border="red" , col="transparent", lwd=3 ) 394text( 60.5, 0.05, "M" ) 395box( lwd=5, col="white" ) ; box( lwd=2, col="black" ) 396 397 plot( pF, border="black", col="transparent", lwd=3 ) 398lines( PrF, border="blue" , col="transparent", lwd=3 ) 399lines( CoxF, border="red" , col="transparent", lwd=3 ) 400text( 60.5, 0.05, "F" ) 401box( lwd=5, col="white" ) ; box( lwd=2, col="black" ) 402 403 404################################################### 405### code chunk number 30: CHANGE1 (eval = FALSE) 406################################################### 407## source( "../R/simLexis.R", keep.source=TRUE ) 408 409 410################################################### 411### code chunk number 31: CHANGE2 412################################################### 413simX <- Epi:::simX 414sim1 <- Epi:::sim1 415lint <- Epi:::lint 416get.next <- Epi:::get.next 417chop.lex <- Epi:::chop.lex 418 419 420################################################### 421### code chunk number 32: simLexis.rnw:972-975 422################################################### 423cbind( 424attr( ini, "time.scales" ), 425attr( ini, "time.since" ) ) 426 427 428################################################### 429### code chunk number 33: simLexis.rnw:1000-1001 430################################################### 431simLexis 432 433 434################################################### 435### code chunk number 34: simLexis.rnw:1018-1019 436################################################### 437simX 438 439 440################################################### 441### code chunk number 35: simLexis.rnw:1031-1032 442################################################### 443sim1 444 445 446################################################### 447### code chunk number 36: simLexis.rnw:1044-1045 448################################################### 449lint 450 451 452################################################### 453### code chunk number 37: simLexis.rnw:1055-1056 454################################################### 455get.next 456 457 458################################################### 459### code chunk number 38: simLexis.rnw:1065-1066 460################################################### 461chop.lex 462 463 464################################################### 465### code chunk number 39: simLexis.rnw:1083-1084 466################################################### 467nState 468 469 470################################################### 471### code chunk number 40: simLexis.rnw:1093-1094 472################################################### 473pState 474 475 476################################################### 477### code chunk number 41: simLexis.rnw:1098-1100 478################################################### 479plot.pState 480lines.pState 481 482 483