### R code from vignette source 'bcMS' ### Encoding: UTF-8 ################################################### ### code chunk number 1: bcMS.rnw:4-7 ################################################### options( width=90, SweaveHooks=list( fig=function() par(mar=c(3,3,1,1),mgp=c(3,1,0)/1.6,las=1,bty="n") ) ) ################################################### ### code chunk number 2: bcMS.rnw:42-85 (eval = FALSE) ################################################### ## library( foreign ) ## bc <- read.dta( "http://fmwww.bc.edu/repec/bocode/m/multistate_example.dta" ) ## names( bc ) ## head( bc ) ## ## ## We will only record dates of relapse and metastasis where they ## ## actually occur; we use the os as the exit date relying on the ## ## factor osi as the death indicator; so we define the times of ## ## relapse (without metastasis), metastasis, exit and death --- ## ## counting the follow-up in years: ## ## bc <- transform( bc, tor = ifelse( rfi==0 , NA, rf )/12, ## tom = ifelse( mfi=="no" , NA, mf )/12, ## tox = os/12, ## tod = ifelse( osi=="alive", NA, os )/12, ## xst = factor( osi, labels=c("Alive","Dead") ), ## hormon = factor(hormon), ## grade = factor(grade) ) ## ## ## correct typo in level name ## levels( bc$size )[2] <- ">20-50 mm" ## ## ## If metastatis at relapse, then metastasis, ## ## if death and metastasis at same time, metastatis 1 week before death ## bc <- transform(bc, tor = ifelse(!is.na(tor) & !is.na(tom) & tor==tom, ## NA, ## tor ), ## tom = ifelse(!is.na(tom) & !is.na(tod) & tom==tod, ## tom-1/52, ## tom), ## pr.tr = pr_1) ## ## ## Select relevant variables for inclusion ## bc <- bc[,c("pid","year","age","meno", ## "size","grade","nodes","pr","pr.tr","er","hormon","chemo", ## "tor","tom","tod","tox","xst","rf","rfi","mf","mfi","os","osi")] ## ## Check ## subset( bc, !is.na(tor) & !is.na(tom) & !is.na(tod) )[1:5,] ## subset( bc, !is.na(tor) & is.na(tom) )[1:5,] ## subset( bc, is.na(tor) & !is.na(tom) )[1:5,] ## BrCa <- bc[,1:17] ## head( BrCa ) ## save( BrCa, file="BC.rda" ) # not included in the Epi package ################################################### ### code chunk number 3: read ################################################### library(Epi) library(popEpi) data(BrCa) str(BrCa) ################################################### ### code chunk number 4: defLex ################################################### set.seed( 1952 ) Lbc <- Lexis( entry = list( tfd = 0, A = age + runif(nrow(BrCa)), P = year + runif(nrow(BrCa)) ), exit = list( tfd = tox ), exit.status = xst, id = pid, data = BrCa ) summary( Lbc ) names( Lbc ) ################################################### ### code chunk number 5: cut ################################################### Rbc <- cutLexis(Lbc, cut = pmin(Lbc$tor, Lbc$tom, na.rm=TRUE), timescale = "tfd", new.state = "Rel", split.states = TRUE, new.scale = "tfr") summary(Rbc, timeScale = TRUE) ################################################### ### code chunk number 6: box-r ################################################### boxes(Rbc, boxpos = list(x = c(15,15,85,85), y = c(85,15,85,15)), show.BE = TRUE, scale.R = 100, cex = 1.1) ################################################### ### code chunk number 7: bcMS.rnw:159-162 ################################################### system.time( Sbc <- splitLexis(Rbc, breaks = seq(0, 100, 1/12), "tfd")) summary(Sbc) ################################################### ### code chunk number 8: bcMS.rnw:166-169 ################################################### system.time( Sbc <- splitMulti(Rbc, tfd = seq(0, 100, 1/12))) summary(Sbc) ################################################### ### code chunk number 9: bcMS.rnw:177-183 ################################################### Stbc <- stack(Sbc) round( ftable(xtabs( cbind('No. records'=lex.id/lex.id, lex.Fail, lex.dur) ~ lex.Tr + lex.Xst, data=Stbc), row.vars=c(3,1)), 1) ################################################### ### code chunk number 10: bcMS.rnw:219-229 ################################################### ( kd.ad <- with( subset( Sbc, lex.Cst=="Alive" & lex.Xst=="Dead"), quantile( tfd+lex.dur, probs=(1:4-0.5)/4) ) ) ( kd.ar <- with( subset( Sbc, lex.Cst=="Alive" & lex.Xst=="Rel"), quantile( tfd+lex.dur, probs=(1:4-0.5)/4) ) ) ( kd.rd <- with( subset( Sbc, lex.Cst=="Rel" & lex.Xst=="Dead(Rel)"), quantile( tfd+lex.dur, probs=(1:4-0.5)/4) ) ) ( kr.rd <- with( subset( Sbc, lex.Cst=="Rel" & lex.Xst=="Dead(Rel)"), quantile( tfr+lex.dur, probs=(1:4-0.5)/4) ) ) ( ka.rd <- with( subset( Sbc, lex.Cst=="Rel" & lex.Xst=="Dead(Rel)"), quantile( tfd-tfr, probs=(1:4-0.5)/4) ) ) ################################################### ### code chunk number 11: oldglm ################################################### m.ad <- glm(cbind(lex.Xst=="Dead", lex.dur) ~ Ns(tfd, knots = kd.ad), family = poisreg, data = subset(Sbc, lex.Cst=="Alive")) m.ar <- glm(cbind(lex.Xst=="Rel", lex.dur) ~ Ns(tfd, knots = kd.ar), family = poisreg, data = subset(Sbc, lex.Cst=="Alive")) m.rd <- glm(cbind(lex.Xst=="Dead(Rel)", lex.dur) ~ Ns(tfd, knots = kd.rd), family = poisreg, data = subset(Sbc, lex.Cst=="Rel")) x.rd <- update(m.rd, . ~ . + Ns(tfr, knots=kr.rd)) r.rd <- update(x.rd, . ~ . - Ns(tfd, knots=kd.rd)) anova( m.rd, x.rd, r.rd, test="Chisq" ) ################################################### ### code chunk number 12: newglm ################################################### M.ad <- glm.Lexis(Sbc, ~ Ns(tfd, knots=kd.ad), to = "Dead") round( cbind( ci.exp(m.ad), ci.exp(M.ad) ), 4 ) attr( M.ad, 'Lexis' ) ################################################### ### code chunk number 13: rate-pr ################################################### nd <- data.frame(tfd = seq(0, 15, 0.1)) ad.rate <- ci.pred(m.ad, nd) ar.rate <- ci.pred(m.ar, nd) rd.rate <- ci.pred(m.rd, nd) ################################################### ### code chunk number 14: pr-pl ################################################### clr <- rainbow(3) ; yl <- c(0.3,200) matshade(nd$tfd, cbind(ad.rate, ar.rate, rd.rate) * 100, plot = TRUE, type = "l", lty = 1, lwd = 3, col = clr, las = 1, log = "y", xlab = "Time since diagnosis (years)", ylim = yl, ylab = "Rate per 100 PY" ) text(par("usr")[2]*0.95, (10^par("usr"))[3]*1.4^(1:3), c("Alive -> Dead","Alive -> Relapse","Relapse -> Dead"), col = clr, adj = 1, font = 2 ) matshade(nd$tfd, ci.ratio(rd.rate, ad.rate), lty = 1, lwd = c(3,1,1), col = gray(0.6)) abline(h = 1, col = gray(0.6)) ################################################### ### code chunk number 15: bcMS.rnw:321-328 ################################################### M.rd <- glm.Lexis(Sbc, ~ Ns(tfd , knots = kd.rd), to = "Dead(Rel)") X.rd <- glm.Lexis(Sbc, ~ Ns(tfd , knots = kd.rd) + Ns( tfr, knots = kr.rd), to = "Dead(Rel)") XX.rd <- glm.Lexis(Sbc, ~ Ns(tfd , knots = kd.rd) + Ns( tfr, knots = kr.rd) + Ns(tfd-tfr, knots = ka.rd), to = "Dead(Rel)") anova( M.rd, X.rd, XX.rd, test = "Chisq" ) ################################################### ### code chunk number 16: bcMS.rnw:338-341 ################################################### x.rd <- update( m.rd, . ~ . + Ns( tfr, knots=kr.rd) ) xx.rd <- update( x.rd, . ~ . + Ns( tfd-tfr, knots=ka.rd) ) anova( m.rd, x.rd, xx.rd, test="Chisq" ) ################################################### ### code chunk number 17: bcMS.rnw:364-365 ################################################### with( subset(Sbc,lex.Xst=="Dead(Rel)"), quantile(tfr+lex.dur,37:39/40) ) ################################################### ### code chunk number 18: rd-int ################################################### nd <- data.frame(expand.grid(tfd = c(NA, seq(0, 15, 0.1)), tad = c(0, 0.5, 1, 2, 3, 5, 8))) nd <- subset(transform(nd, tfr = tfd - tad ), (tfr>=0 & tfr<7) | is.na(tfr) ) head( nd ) matshade(nd$tfd, cbind(ci.pred( x.rd, nd), ci.pred(xx.rd, nd)) * 100, plot = TRUE, type = "l", lty = c("11", "solid"), lend = "butt", lwd = 3, col = gray(c(5, 0) / 10), alpha = c(0, 0.07), las = 1, log = "y", xlab = "Time since diagnosis (years)", ylim = c(5,100), ylab = "Mortality rate per 100 PY" ) matshade(tt <- seq(0,15,0.1), ci.pred( m.rd, data.frame(tfd = tt) )*100, lwd = 3, lty = 1, col = clr[3] ) ################################################### ### code chunk number 19: bcMS.rnw:403-408 ################################################### c.ar <- update( m.ar, . ~ . + age + size + nodes + pr.tr + hormon) c.ad <- update( m.ad, . ~ . + age + size + nodes + pr.tr + hormon) c.rd <- update( m.rd, . ~ . + age + size + nodes + pr.tr + hormon) cx.rd <- update( x.rd, . ~ . + age + size + nodes + pr.tr + hormon) cxx.rd<- update(xx.rd, . ~ . + age + size + nodes + pr.tr + hormon) ################################################### ### code chunk number 20: bcMS.rnw:413-420 ################################################### nd <- transform(nd, age = 54, size = "<=20 mm", nodes = 1, pr.tr = 3, hormon = "no") head( nd ) ################################################### ### code chunk number 21: rd-int-cov ################################################### matshade(nd$tfd, cbind(ci.pred( cx.rd, nd), ci.pred(cxx.rd, nd)) * 100, plot = TRUE, type = "l", lty = c("11", "solid"), lend = "butt", lwd = 3, col = gray(c(5, 0) / 10), alpha = c(0, 0.07), las = 1, log = "y", xlab = "Time since diagnosis (years)", ylim = c(5,100), ylab = "Mortality rate per 100 PY" ) nt <- data.frame(tfd = seq(0, 15, 0.1)) nt <- transform(nt, age = 54, size = "<=20 mm", nodes = 1, pr.tr = 3, hormon = "no") matshade(nt$tfd, ci.pred(c.rd, nt) * 100, lwd = 3, lty = 1, col = clr[3] ) ################################################### ### code chunk number 22: int-test (eval = FALSE) ################################################### ## int.test <- NArray( list( model=c("c.ar","c.ad","c.rd","cx.rd"), ## var=c("age","size","nodes","pr.tr","hormon"), ## what=c("d.f.","Dev","P") ) ) ## str( int.test ) ## ## int.test[1,1,]<-as.numeric(anova( c.ar,update( c.ar,.~.+log(tfd+0.5):age ),test="Chisq")[2,3:5]) ## int.test[1,2,]<-as.numeric(anova( c.ar,update( c.ar,.~.+log(tfd+0.5):size ),test="Chisq")[2,3:5]) ## int.test[1,3,]<-as.numeric(anova( c.ar,update( c.ar,.~.+log(tfd+0.5):nodes ),test="Chisq")[2,3:5]) ## int.test[1,4,]<-as.numeric(anova( c.ar,update( c.ar,.~.+log(tfd+0.5):pr.tr ),test="Chisq")[2,3:5]) ## int.test[1,5,]<-as.numeric(anova( c.ar,update( c.ar,.~.+log(tfd+0.5):hormon),test="Chisq")[2,3:5]) ## ## int.test[2,1,]<-as.numeric(anova( c.ad,update( c.ad,.~.+log(tfd+0.5):age ),test="Chisq")[2,3:5]) ## int.test[2,2,]<-as.numeric(anova( c.ad,update( c.ad,.~.+log(tfd+0.5):size ),test="Chisq")[2,3:5]) ## int.test[2,3,]<-as.numeric(anova( c.ad,update( c.ad,.~.+log(tfd+0.5):nodes ),test="Chisq")[2,3:5]) ## int.test[2,4,]<-as.numeric(anova( c.ad,update( c.ad,.~.+log(tfd+0.5):pr.tr ),test="Chisq")[2,3:5]) ## int.test[2,5,]<-as.numeric(anova( c.ad,update( c.ad,.~.+log(tfd+0.5):hormon),test="Chisq")[2,3:5]) ## ## int.test[3,1,]<-as.numeric(anova( c.rd,update( c.rd,.~.+log(tfd+0.5):age ),test="Chisq")[2,3:5]) ## int.test[3,2,]<-as.numeric(anova( c.rd,update( c.rd,.~.+log(tfd+0.5):size ),test="Chisq")[2,3:5]) ## int.test[3,3,]<-as.numeric(anova( c.rd,update( c.rd,.~.+log(tfd+0.5):nodes ),test="Chisq")[2,3:5]) ## int.test[3,4,]<-as.numeric(anova( c.rd,update( c.rd,.~.+log(tfd+0.5):pr.tr ),test="Chisq")[2,3:5]) ## int.test[3,5,]<-as.numeric(anova( c.rd,update( c.rd,.~.+log(tfd+0.5):hormon),test="Chisq")[2,3:5]) ## ## int.test[4,1,]<-as.numeric(anova(cx.rd,update(cx.rd,.~.+log(tfd+0.5):age ),test="Chisq")[2,3:5]) ## int.test[4,2,]<-as.numeric(anova(cx.rd,update(cx.rd,.~.+log(tfd+0.5):size ),test="Chisq")[2,3:5]) ## int.test[4,3,]<-as.numeric(anova(cx.rd,update(cx.rd,.~.+log(tfd+0.5):nodes ),test="Chisq")[2,3:5]) ## int.test[4,4,]<-as.numeric(anova(cx.rd,update(cx.rd,.~.+log(tfd+0.5):pr.tr ),test="Chisq")[2,3:5]) ## int.test[4,5,]<-as.numeric(anova(cx.rd,update(cx.rd,.~.+log(tfd+0.5):hormon),test="Chisq")[2,3:5]) ## save( int.test, file="int-test.Rda") ################################################### ### code chunk number 23: bcMS.rnw:499-503 ################################################### load( file="/home/bendix/teach/AdvCoh/00/examples/bcMS/int-test.Rda") round( int.test[,,2], 2 ) round( int.test[,,3], 4 ) round( int.test[,,1], 0 ) ################################################### ### code chunk number 24: bcMS.rnw:528-538 ################################################### i.ar <- update( c.ar, . ~ . + log(tfd+0.5):size + log(tfd+0.5):pr.tr + log(tfd+0.5):hormon) i.ad <- c.ad i.rd <- update( c.rd, . ~ . + log(tfd+0.5):pr.tr) ix.rd <- update(cx.rd, . ~ . + log(tfd+0.5):pr.tr) round(ci.lin(i.ad ), 4) round(ci.lin(i.ar ), 4) round(ci.lin(i.rd ), 4) round(ci.lin(ix.rd), 4) ################################################### ### code chunk number 25: bcMS.rnw:553-571 ################################################### nd.size <- data.frame(tfd = rep( c(NA,seq(0,15,0.1)), 3 ), age = 45, size = rep( levels(Lbc$size), each=152 ), nodes = 5, pr.tr = 3, hormon = levels(Lbc$hormon)[1] ) nd.pr <- data.frame(tfd = rep( c(NA,seq(0,15,0.1)), 6 ), age = 45, size = levels(Lbc$size)[2], nodes = 5, pr.tr = rep( 0:5, each=152 ), hormon = levels(Lbc$hormon)[1] ) nd.hormon <- data.frame(tfd = rep( c(NA,seq(0,15,0.1)), 2 ), age = 45, size = levels(Lbc$size)[2], nodes = 5, pr.tr = 3, hormon = rep( levels(Lbc$hormon), each=152 ) ) ################################################### ### code chunk number 26: int-size ################################################### clr <- rainbow(3) ; yl <- c(0.03,60) ad.c.rate <- ci.pred(c.ad, nd.size) ad.i.rate <- ci.pred(i.ad, nd.size) ar.c.rate <- ci.pred(c.ar, nd.size) ar.i.rate <- ci.pred(i.ar, nd.size) rd.c.rate <- ci.pred(c.rd, nd.size) rd.i.rate <- ci.pred(i.rd, nd.size) matplot(nd.size$tfd, cbind(ad.c.rate, ad.i.rate, ar.c.rate, ar.i.rate, rd.c.rate, rd.i.rate) * 100, type = "l", lty = rep(c("22", "solid"),each = 3), lwd = c(2,0,0), col = rep(clr, each = 6), las = 1, lend = "butt", log = "y", xlab = "Time since diagnosis (years)", ylim = yl, ylab = "Rate per 100 PY" ) text( par("usr")[2]*0.95, (10^par("usr"))[3]*1.4^(1:3), c("A->D","A->R","R->D"), col = clr, adj = 1, font = 2 ) ################################################### ### code chunk number 27: int-pr ################################################### ad.c.rate <- ci.pred(c.ad, nd.pr) ad.i.rate <- ci.pred(i.ad, nd.pr) ar.c.rate <- ci.pred(c.ar, nd.pr) ar.i.rate <- ci.pred(i.ar, nd.pr) rd.c.rate <- ci.pred(c.rd, nd.pr) rd.i.rate <- ci.pred(i.rd, nd.pr) matplot(nd.pr$tfd, cbind(ad.c.rate, ad.i.rate, ar.c.rate, ar.i.rate, rd.c.rate, rd.i.rate) * 100, type = "l", lty = rep(c("22", "solid"),each = 3), lwd = c(2,0,0), col = rep(clr, each = 6), las = 1, lend = "butt", log = "y", xlab = "Time since diagnosis (years)", ylim = yl, ylab = "Rate per 100 PY" ) text( par("usr")[2]*0.95, (10^par("usr"))[3]*1.4^(1:3), c("A->D","A->R","R->D"), col = clr, adj = 1, font = 2 ) ################################################### ### code chunk number 28: int-hormon ################################################### ad.c.rate <- ci.pred(c.ad, nd.hormon) ad.i.rate <- ci.pred(i.ad, nd.hormon) ar.c.rate <- ci.pred(c.ar, nd.hormon) ar.i.rate <- ci.pred(i.ar, nd.hormon) rd.c.rate <- ci.pred(c.rd, nd.hormon) rd.i.rate <- ci.pred(i.rd, nd.hormon) matplot(nd.hormon$tfd, cbind(ad.c.rate, ad.i.rate, ar.c.rate, ar.i.rate, rd.c.rate, rd.i.rate) * 100, type = "l", lty = rep(c("22","solid"),each = 3), lwd = c(2,0,0), col = rep(clr, each = 6), las = 1, lend = "butt", log = "y", ylim = yl, ylab = "Rate per 100 PY", xlab = "Time since diagnosis (years)") text(par("usr")[2]*0.95, (10^par("usr"))[3]*1.4^(1:3), c("A->D","A->R","R->D"), col = clr, adj = 1, font = 2 ) ################################################### ### code chunk number 29: bcMS.rnw:699-715 ################################################### names( Rbc ) Lini <- Rbc[NULL,c("tfd","A","P","tfr", "lex.Cst","lex.Xst","lex.dur","lex.id", "age","size","nodes","pr.tr","hormon")] pr.nodes <- c(0,1,5,10) npr <- nlevels(Rbc$size) * length(pr.nodes) Lini[1:npr,"tfd"] <- 0 Lini[1:npr,"tfr"] <- NA Lini[1:npr,"lex.Cst"] <- "Alive" Lini[1:npr,"age"] <- 54 Lini[1:npr,"size"] <- rep(levels(Rbc$size), length(pr.nodes)) Lini[1:npr,"nodes"] <- rep(pr.nodes, each=nlevels(Rbc$size)) Lini[1:npr,"pr.tr"] <- 3 Lini[1:npr,"hormon"] <- "no" Lini str( Lini ) ################################################### ### code chunk number 30: bcMS.rnw:728-736 ################################################### TR <- list( Alive = list( Dead = i.ad, Rel = i.ar ), Rel = list( "Dead(Rel)" = i.rd ) ) TRx <- list( Alive = list( Dead = i.ad, Rel = i.ar ), Rel = list( "Dead(Rel)" = ix.rd ) ) lapply( TR, names ) lapply( TR, lapply, class ) ################################################### ### code chunk number 31: simx ################################################### system.time( xx <- simLexis( Tr=TR , init=Lini, N=200, t.range=16 ) ) ################################################### ### code chunk number 32: simL (eval = FALSE) ################################################### ## # not evaluated, run interactively before final compilation ## set.seed( 1952 ) ## x0 <- simLexis(Tr = TR , init = Lini, N = 2000, t.range = 16) ## x1 <- simLexis(Tr = TR , init = Lini, N = 2000, t.range = 16) ## x2 <- simLexis(Tr = TR , init = Lini, N = 2000, t.range = 16) ## x3 <- simLexis(Tr = TR , init = Lini, N = 2000, t.range = 16) ## x4 <- simLexis(Tr = TR , init = Lini, N = 2000, t.range = 16) ## sL <- rbind(x0, transform(x1, lex.id = lex.id+25000 ), ## transform(x2, lex.id = lex.id+50000 ), ## transform(x3, lex.id = lex.id+75000 ), ## transform(x4, lex.id = lex.id+100000)) ## summary( sL ) ## s0 <- simLexis(Tr = TRx, init = Lini, N = 2000, t.range = 16) ## s1 <- simLexis(Tr = TRx, init = Lini, N = 2000, t.range = 16) ## s2 <- simLexis(Tr = TRx, init = Lini, N = 2000, t.range = 16) ## s3 <- simLexis(Tr = TRx, init = Lini, N = 2000, t.range = 16) ## s4 <- simLexis(Tr = TRx, init = Lini, N = 2000, t.range = 16) ## sLx<- rbind(s0, transform(s1, lex.id = lex.id+25000 ), ## transform(s2, lex.id = lex.id+50000 ), ## transform(s3, lex.id = lex.id+75000 ), ## transform(s4, lex.id = lex.id+100000)) ## summary( sLx ) ## save(sL, sLx, ## file = "/home/bendix/teach/AdvCoh/00/examples/bcMS/sL.Rda") ################################################### ### code chunk number 33: loadsim ################################################### load(file = "/home/bendix/teach/AdvCoh/00/examples/bcMS/sL.Rda") summary(sLx) ################################################### ### code chunk number 34: nStstr ################################################### nn <- nState( sLx[1:1000,], at=seq(0,16,0.1), from=0, time.scale="tfd" ) pp <- pState( nn, perm=c(1,2,4,3) ) str( pp ) ################################################### ### code chunk number 35: defNArr ################################################### (tt <- with(sLx, table(nodes, size))) prX <- prA <- NArray(c(dimnames(tt), dimnames(pp))) str(prA) ################################################### ### code chunk number 36: fillNArr (eval = FALSE) ################################################### ## for( nn in dimnames(prA)[[1]] ) ## for( ss in dimnames(prA)[[2]] ) ## { ## prA[nn,ss,,] <- pState( nState( subset( sL , nodes==as.numeric(nn) & ## size==ss ), ## at = seq(0,16,0.1), ## from = 0, ## time.scale = "tfd" ), ## perm = c(1,2,4,3) ) ## prX[nn,ss,,] <- pState( nState( subset( sLx, nodes==as.numeric(nn) & ## size==ss ), ## at = seq(0,16,0.1), ## from = 0, ## time.scale = "tfd" ), ## perm = c(1,2,4,3) ) ## } ## save(prA, prX, file = "pr.Rda") ################################################### ### code chunk number 37: states ################################################### load(file = "/home/bendix/teach/AdvCoh/00/examples/bcMS/pr.Rda") clr <- c("forestgreen","maroon") clr <- cbind(clr, adjustcolor(clr[2:1], 0.5)) par(mfcol = c(3,4), mar = c(1,1.5,1,1), mgp = c(3,1,0)/1.6, oma = c(2,2,2,2), las = 1, bty="n") nnn <- dimnames(prA)[[1]] sss <- dimnames(prA)[[2]] for( nn in nnn ) for( ss in sss ) { plot.pState(prX[nn,ss,,], col = clr, xlim = c(0,15), ylab = "", xlab = "" ) lines(as.numeric(dimnames(prX)[[3]]), prX[nn,ss,, 2], lwd = 3, lty = 1, col = "black" ) matlines(as.numeric(dimnames(prA)[[3]]), prA[nn,ss,,1:3], lwd = 1, lty = 1, col = "white" ) axis(side = 2, at = 0:10/10, labels = NA, tcl = -0.4) axis(side = 4, at = 0:10/10, labels = NA, tcl = -0.4) axis(side = 2, at = 0:50/50, labels = NA, tcl = -0.2) axis(side = 4, at = 0:50/50, labels = NA, tcl = -0.2) } mtext(paste( "Nodes =",nnn), side = 3, at = (1:4*2-1)/8, outer = TRUE, line = 0, cex = 0.66, las = 0 ) mtext(paste( "Size" ,sss), side = 4, at = (3:1*2-1)/6, outer = TRUE, line = 0, cex = 0.66, las = 0 ) mtext("Time since diagnosis (years)", side = 1, outer = TRUE, line = 1, cex = 0.66, las = 0 ) mtext("Probability", side = 2, outer = TRUE, line = 1, cex = 0.66, las = 0 ) ################################################### ### code chunk number 38: sojourntimes ################################################### intgr <- function(y, x) (sum(y[-1]) + sum(y[-length(y)])) / 2 * x cA <- apply(prA[,,1:151,1:3], c(1,2,4), intgr, 0.1) cA[,,3] <- cA[,,2] - cA[,,1] dimnames( cA )[[3]] <- c("noRel","Total","Rel") cA <- cA[,,c(1,3,2)] round(ftable(cA, row.vars = c(3,2)), 2) cX <- apply( prX[,,1:151,1:3], c(1,2,4), intgr, 0.1) cX[,,3] <- cX[,,2] - cX[,,1] dimnames( cX )[[3]] <- c("noRel","Total","Rel") cX <- cX[,,c(1,3,2)] round(ftable(cX, row.vars = c(3,2)), 2) ################################################### ### code chunk number 39: sjt-comp ################################################### round(ftable( cX - cA , row.vars = c(3,2) ), 2) round(ftable((cX - cA) / cA*100, row.vars = c(3,2) ), 1) ################################################### ### code chunk number 40: mcut ################################################### mbc <- mcutLexis(Lbc, timescale = "tfd", wh = c("tor", "tom"), new.states = c("Rel", "Met"), seq.states = TRUE, new.scales = c("tfr", "tfm")) summary( mbc, timeScale = TRUE ) mbc <- Relevel(mbc, list("Alive", "Rel", Met = c("Met","Rel-Met"), "Dead")) summary(mbc) # or: mbc <- Relevel(mbc, list(1, 3, Met = 4:5, 2)) print(subset(mbc, lex.id %in% (1328+0:2))[,1:10], row.names = FALSE, digits = 4) ################################################### ### code chunk number 41: sub-death ################################################### xbc <- transform(mbc, lex.Xst = factor(ifelse(lex.Cst != "Alive" & lex.Xst == "Dead", paste("D(", lex.Cst, ")", sep=""), as.character(lex.Xst)))) xbc <- factorize( xbc ) levels(xbc) xbc <- Relevel(xbc, c(1:4,6,5)) levels(xbc) summary(xbc) ################################################### ### code chunk number 42: box-rmx ################################################### boxes(xbc, boxpos = list(x = c(15,40,15,85,85,85), y = c(85,50,15,85,50,15)), show.BE = "nz", scale.R = 100, cex = 1.1) ################################################### ### code chunk number 43: bcMS.rnw:1051-1052 (eval = FALSE) ################################################### ## Tr <- lapply( Tr, lapply, get ) ################################################### ### code chunk number 44: bcMS.rnw:1059-1060 (eval = FALSE) ################################################### ## bootTr <- lapply( Tr, lapply, function(x) paste("BOOT",x,sep="") ) ################################################### ### code chunk number 45: bcMS.rnw:1064-1070 (eval = FALSE) ################################################### ## unique.models <- unique( unlist( Tr ) ) ## for( m in unique.models ) ## { ## assign( paste("BOOT",m,sep=""), ## update( get(m), data=boot.Lexis(data) ) ) ## }