### R code from vignette source 'ms-steno2.rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: ms-steno2.rnw:2-7 ################################################### options( width=90, # prompt=" ", continue=" ", # Absence of prompts makes it easier for # students to copy from the final pdf document 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: ms-steno2.rnw:36-43 ################################################### library(survival) library(Epi) library(popEpi) # popEpi::splitMulti returns a data.frame rather than a data.table options("popEpi.datatable" = FALSE) library(tidyverse) clear() ################################################### ### code chunk number 3: ms-steno2.rnw:48-55 ################################################### nround <- function(df, dec = 2) { wh.num <- sapply(df, is.numeric) df[,wh.num] <- round(df[,wh.num], dec) print(df) } ################################################### ### code chunk number 4: ms-steno2.rnw:67-73 ################################################### data(steno2) steno2 <- transform(cal.yr(steno2), doEnd = ifelse(!is.na(doDth), doDth, doEnd)) str(steno2) ################################################### ### code chunk number 5: ms-steno2.rnw:82-92 ################################################### L2 <- Lexis(entry = list(per = doBase, age = doBase - doBth, tfi = 0), exit = list(per = doEnd), exit.status = factor(deathCVD + !is.na(doDth), labels=c("Mic","D(oth)","D(CVD)")), id = id, data = steno2) summary(L2, t = TRUE) boxes(L2, boxpos = TRUE, show.BE = TRUE) ################################################### ### code chunk number 6: ms-steno2.rnw:126-132 ################################################### data(st2alb) cut2 <- rename(cal.yr(st2alb), lex.id = id, cut = doTr, new.state = state) str(cut2) ################################################### ### code chunk number 7: ms-steno2.rnw:135-136 ################################################### with(cut2, addmargins(table(table(lex.id)))) ################################################### ### code chunk number 8: boxL3 ################################################### L3 <- rcutLexis(L2, cut2) summary(L3) boxes(L3, boxpos = TRUE, cex = 0.8) ################################################### ### code chunk number 9: ms-steno2.rnw:166-170 ################################################### (jump <- subset(L3, (lex.Cst == "Norm" & lex.Xst == "Mac") | (lex.Xst == "Norm" & lex.Cst == "Mac"))[, c("lex.id", "per", "lex.dur","lex.Cst", "lex.Xst")]) ################################################### ### code chunk number 10: ms-steno2.rnw:177-183 ################################################### set.seed(1952) xcut <- select(transform(jump, cut = per + lex.dur * runif(per, 0.1, 0.9), new.state = "Mic"), c(lex.id, cut, new.state)) xcut ################################################### ### code chunk number 11: ms-steno2.rnw:190-192 ################################################### L4 <- rcutLexis(L3, xcut) summary(L4) ################################################### ### code chunk number 12: b4 ################################################### boxes(L4, boxpos = list(x = c(20,20,20,80,80), y = c(50,90,10,75,25)), show.BE = "nz", scale.R = 100, cex = 0.9) ################################################### ### code chunk number 13: ms-steno2.rnw:231-234 ################################################### S4 <- splitMulti(L4, tfi = seq(0, 25, 1/2)) summary(L4) summary(S4) ################################################### ### code chunk number 14: ms-steno2.rnw:243-247 ################################################### ma <- glm.Lexis(S4, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst) ci.exp(ma) ################################################### ### code chunk number 15: ms-steno2.rnw:259-266 (eval = FALSE) ################################################### ## ma <- glm(cbind(lex.Xst %in% c("D(oth)", "D(CVD)") & lex.Cst != lex.Xst, ## lex.dur) ## ~ Ns(tfi, knots = seq( 0, 20, 5)) + ## Ns(age, knots = seq(50, 80, 10)) + ## lex.Cst, ## family = poisreg, ## data = subset(S4, lex.Cst %in% c("Norm","Mic","Mac"))) ################################################### ### code chunk number 16: ms-steno2.rnw:269-276 (eval = FALSE) ################################################### ## ma <- glm((lex.Xst %in% c("D(oth)", "D(CVD)") & lex.Cst != lex.Xst) ## ~ Ns(tfi, knots = seq( 0, 20, 5)) + ## Ns(age, knots = seq(50, 80, 10)) + ## lex.Cst, ## offset = log(lex.dur), ## family = poisson, ## data = subset(S4, lex.Cst %in% c("Norm","Mic","Mac"))) ################################################### ### code chunk number 17: ms-steno2.rnw:284-285 ################################################### round(ci.exp(ma), 2) ################################################### ### code chunk number 18: ms-steno2.rnw:291-292 ################################################### Wald(ma, subset = "lex.Cst") ################################################### ### code chunk number 19: ms-steno2.rnw:300-310 ################################################### mo <- glm.Lexis(S4, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst, to = "D(oth)") round(ci.exp(mo), 3) mC <- glm.Lexis(S4, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst, to = "D(CVD)") round(ci.exp(mC), 3) ################################################### ### code chunk number 20: ms-steno2.rnw:316-318 ################################################### Wald(mo, subset = "Cst") Wald(mC, subset = "Cst") ################################################### ### code chunk number 21: ms-steno2.rnw:325-326 ################################################### summary(L2$age) ################################################### ### code chunk number 22: ms-steno2.rnw:332-334 ################################################### expand.grid(tfi = c(NA, seq(0, 20, 5)), ain = c(45, 55, 65)) ################################################### ### code chunk number 23: mort1 ################################################### prf <- transform(expand.grid(tfi = c(NA, seq(0, 20, 0.5)), ain = c(45, 55, 65))[-1,], age = ain + tfi, lex.Cst = "Mic") head(prf) prf[40:44,] matshade(prf$age, cbind(ci.pred(mo, prf), ci.pred(mC, prf)) * 100, lwd = 3, col = c("black","blue"), log = "y", ylim = c(0.01,50), plot = TRUE) ################################################### ### code chunk number 24: mort3 ################################################### par(mfrow=c(1,3)) for(st in c("Norm","Mic","Mac")) { matshade(prf$age, cbind(ci.pred(mo, transform(prf, lex.Cst = st)), ci.pred(mC, transform(prf, lex.Cst = st))) * 100, lwd = 3, col = c("black","blue"), log = "y", ylim = c(0.1,50), plot = TRUE) text(60, 50, st, adj = 0) } ################################################### ### code chunk number 25: ms-steno2.rnw:403-408 ################################################### mix <- glm.Lexis(S4, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst * allo, to = "D(oth)") round(ci.exp(mix), 3) ################################################### ### code chunk number 26: ms-steno2.rnw:414-420 ################################################### mox <- glm.Lexis(S4, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst / allo, to = "D(oth)") round(ci.exp(mox), 3) c(deviance(mox), deviance(mix)) ################################################### ### code chunk number 27: ms-steno2.rnw:428-433 ################################################### mCx <- glm.Lexis(S4, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst / allo, to = "D(CVD)") round(ci.exp(mCx), 3) ################################################### ### code chunk number 28: ms-steno2.rnw:441-455 ################################################### det <- glm.Lexis(S4, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst / allo, from = c("Norm","Mic"), to = c("Mic","Mac")) imp <- glm.Lexis(S4, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst / allo, from = c("Mic","Mac"), to = c("Norm","Mic")) round(ci.exp(det), 3) round(ci.exp(imp), 3) round(ci.exp(det, subset="al"), 1) round(ci.exp(imp, subset="al"), 1) ################################################### ### code chunk number 29: ms-steno2.rnw:504-508 ################################################### ini <- L2[,c("per", "age", "tfi")] ini <- rbind(transform(ini, lex.Cst = "Mic", allo = "Int"), transform(ini, lex.Cst = "Mic", allo = "Conv")) str(ini) ################################################### ### code chunk number 30: ms-steno2.rnw:516-527 ################################################### Tr <- list(Norm = list(Mic = det, "D(oth)" = mox, "D(CVD)" = mCx), Mic = list(Mac = det, Norm = imp, "D(oth)" = mox, "D(CVD)" = mCx), Mac = list(Mic = imp, "D(oth)" = mox, "D(CVD)" = mCx)) lapply(Tr, names) ################################################### ### code chunk number 31: ms-steno2.rnw:540-547 ################################################### set.seed(1952) system.time( Sorg <- simLexis(Tr = Tr, # models for each transition init = ini, # cohort of straters N = 10, # how many copies of each person in ini t.range = 20, # how long should we simulate before censoring n.int = 200))# how many intervals for evaluating rates ################################################### ### code chunk number 32: ms-steno2.rnw:551-554 ################################################### Sorg <- Relevel(Sorg,c("Norm", "Mic", "Mac", "D(CVD)", "D(oth)")) summary(Sorg, t = T) nround(subset(Sorg, lex.id %in% 28:32), 2) ################################################### ### code chunk number 33: ms-steno2.rnw:558-559 ################################################### addmargins(table(table(Sorg$lex.id))) ################################################### ### code chunk number 34: ms-steno2.rnw:566-573 ################################################### system.time( Nst <- nState(Sorg, at = seq(0, 20, 0.1), from = 0, time.scale = "tfi")) str(Nst) head(Nst) ################################################### ### code chunk number 35: ms-steno2.rnw:578-588 ################################################### Nint <- nState(subset(Sorg, allo == "Int"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi") Nconv<- nState(subset(Sorg, allo == "Conv"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi") head(Nint) head(Nconv) ################################################### ### code chunk number 36: ms-steno2.rnw:596-600 ################################################### Pint <- pState(Nint ) Pconv <- pState(Nconv) str(Pint) head(Pint) ################################################### ### code chunk number 37: pStates ################################################### par(mfrow = c(1,2), mar=c(3,3,2,2)) clr <- c("forestgreen", "orange", "red", "blue", gray(0.4)) plot(Pint, col = clr, xlim = c(0, 20)) # the survival curve lines(as.numeric(rownames(Pint)), Pint[,"Mac"], lwd = 4, col = "white") text(rownames(Pint)[150], Pint[150,] - diff(c(0, Pint[150,]))/2, colnames(Pint), col = "white") plot(Pconv, col = clr, xlim = c(20, 0)) # the survival curve lines(as.numeric(rownames(Pconv)), Pconv[,"Mac"], lwd = 4) text(rownames(Pconv)[150], Pconv[150,] - diff(c(0, Pconv[150,]))/2, colnames(Pconv), col = "white") mtext(c("Int","Conv"), side = 3, at = c(1,3)/4, outer = TRUE, line = -2) ################################################### ### code chunk number 38: ms-steno2.rnw:655-666 ################################################### ini <- S4[1:10,c("lex.id", "per", "age", "tfi", "lex.Cst", "allo")] ini[,"lex.id"] <- 1:10 ini[,"per"] <- 1993 # not used but it is a time scale in S4 ini[,"age"] <- ini[,"ain"] <- rep(seq(45,65,5), 2) ini[,"tfi"] <- 0 ini[,"lex.Cst"] <- factor("Mic", levels = c("Norm","Mic","Mac","D(CVD)","D(oth)")) ini[,"allo"] <- factor(rep(c("Int","Conv"), each = 5)) ini str(ini) ################################################### ### code chunk number 39: ms-steno2.rnw:678-687 ################################################### system.time( Sdef <- simLexis(Tr = Tr, init = ini, N = 100, t.range = 20, n.int = 200)) # str(Sdef) summary(Sdef) nround(head(Sdef)) ################################################### ### code chunk number 40: ms-steno2.rnw:695-701 ################################################### P45i <- nState(subset(Sdef, ain == 45 & allo == "Int"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi") head(P45i) head(pState(P45i)) ################################################### ### code chunk number 41: ms-steno2.rnw:705-745 ################################################### P45c <- pState(nState(subset(Sdef, ain == 45 & allo == "Conv"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi")) P45i <- pState(nState(subset(Sdef, ain == 45 & allo == "Int"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi")) P50c <- pState(nState(subset(Sdef, ain == 55 & allo == "Conv"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi")) P50i <- pState(nState(subset(Sdef, ain == 55 & allo == "Int"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi")) P55c <- pState(nState(subset(Sdef, ain == 55 & allo == "Conv"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi")) P55i <- pState(nState(subset(Sdef, ain == 55 & allo == "Int"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi")) P60c <- pState(nState(subset(Sdef, ain == 55 & allo == "Conv"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi")) P60i <- pState(nState(subset(Sdef, ain == 55 & allo == "Int"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi")) P65c <- pState(nState(subset(Sdef, ain == 65 & allo == "Conv"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi")) P65i <- pState(nState(subset(Sdef, ain == 65 & allo == "Int"), at = seq(0, 20, 0.1), from = 0, time.scale = "tfi")) ################################################### ### code chunk number 42: panel5 ################################################### par(mfrow = c(5,2), mar = c(1,1,0,0), oma = c(3,3,1,0), mgp=c(3,1,0)/1.6) plot(P45i, col = clr, xlim = c(0,20)) plot(P45c, col = clr, xlim = c(20,0)) plot(P50i, col = clr, xlim = c(0,20)) plot(P50c, col = clr, xlim = c(20,0)) plot(P55i, col = clr, xlim = c(0,20)) plot(P55c, col = clr, xlim = c(20,0)) plot(P60i, col = clr, xlim = c(0,20)) plot(P60c, col = clr, xlim = c(20,0)) plot(P65i, col = clr, xlim = c(0,20)) plot(P65c, col = clr, xlim = c(20,0)) mtext(c("Int","Conv"), side = 3, at = c(1,3)/4, outer = TRUE, line = 0) mtext(paste(seq(45,65,5)), side = 2, at = (5:1*2-1)/10, outer = TRUE, line = 0) ################################################### ### code chunk number 43: ms-steno2.rnw:781-787 ################################################### (ain <- seq(45, 65, 5)) (all <- levels(S4$allo)) pdef <- NArray(c(list(ain = ain, allo = all), dimnames(P45i))) str(pdef) ################################################### ### code chunk number 44: ms-steno2.rnw:793-802 ################################################### for(aa in ain) for(gg in all) pdef[paste(aa), gg, ,] <- nState(subset(Sdef, ain == aa & allo == gg), at = as.numeric(dimnames(pdef)[["when"]]), from = 0, time.scale = "tfi") pdef <- sweep(pdef, 1:2, pdef[,,1,"Mic"], "/") str(pdef) ################################################### ### code chunk number 45: panel3 ################################################### ain <- seq(45, 65, 10) par(mfrow = c(length(ain),2), mar = c(3,3,1,1), oma = c(0,2,1,0), mgp = c(3,1,0) / 1.6) for(aa in ain) { mat2pol(pdef[paste(aa),"Int" ,,], col = clr, xlim = c(0,20)) mat2pol(pdef[paste(aa),"Conv",,], col = clr, xlim = c(20,0)) } mtext(c("Int","Conv"), side = 3, at = c(1,3)/4, outer = TRUE, line = 0) mtext(ain, side = 2, at = (length(ain):1 * 2 - 1) / (length(ain) * 2), outer = TRUE, line = 0) ################################################### ### code chunk number 46: ms-steno2.rnw:829-831 ################################################### str(pdef) ftable(pdef[,,1:5,], row.vars = c(1,3)) ################################################### ### code chunk number 47: ms-steno2.rnw:845-852 ################################################### mid <- function(x) x[-1] - diff(x) / 2 str(pdef) pmid <- apply(pdef, c(1,2,4), mid) str(pmid) pyr <- apply(pmid, 2:4, sum) * 0.1 str(pyr) round(ftable(pyr, col.vars = 3:2), 1) ################################################### ### code chunk number 48: ms-steno2.rnw:865-866 ################################################### round(pyr[,"Int",] - pyr[,"Conv",], 1) ################################################### ### code chunk number 49: ms-steno2.rnw:883-888 ################################################### AaJ <- survfit(Surv(tfi, tfi + lex.dur, lex.Xst) ~ 1, id = lex.id, data = L4) AaJ$transitions summary(L4) ################################################### ### code chunk number 50: ms-steno2.rnw:892-896 (eval = FALSE) ################################################### ## survfit(Surv(tfi, tfi + lex.dur, lex.Xst) ~ 1, ## id = lex.id, ## istate = lex.Cst, ## data = L4) ################################################### ### code chunk number 51: ms-steno2.rnw:905-913 ################################################### R4 <- sortLexis(L4) last <- rev(!duplicated(rev(R4$lex.id))) R4$lex.Xst <- ifelse(last & R4$lex.Cst == R4$lex.Xst, "cens", as.character(R4$lex.Xst)) R4 <- Relevel(factorize(R4), "cens") summary(L4) summary(R4) ################################################### ### code chunk number 52: ms-steno2.rnw:935-939 ################################################### AaJ <- survfit(Surv(tfi, tfi + lex.dur, lex.Xst) ~ 1, id = lex.id, istate = lex.Cst, data = R4) ################################################### ### code chunk number 53: ms-steno2.rnw:943-945 ################################################### AaJ$transitions[,c(6,1:5)] summary(R4) ################################################### ### code chunk number 54: ms-steno2.rnw:950-955 ################################################### names(AaJ) AaJ$states head(AaJ$pstate) head(AaJ$lower) head(AaJ$upper) ################################################### ### code chunk number 55: ms-steno2.rnw:958-961 ################################################### mat2pol(AaJ$pstate, perm = c(2,1,3,5,4), x = AaJ$time, col = clr) lines(AaJ$time, apply(AaJ$pstate[,1:3], 1, sum), lwd = 5) ################################################### ### code chunk number 56: ms-steno2.rnw:966-970 ################################################### AaJ <- survfit(Surv(tfi, tfi + lex.dur, lex.Xst) ~ allo, id = lex.id, istate = lex.Cst, data = R4) ################################################### ### code chunk number 57: ms-steno2.rnw:975-979 ################################################### AaJ$states AaJ$strata wh <- rep(substr(names(AaJ$strata), 6, 9), AaJ$strata) table(wh) ################################################### ### code chunk number 58: AaJstates ################################################### par(mfrow = c(1,2), mar=c(3,3,2,2)) mat2pol(AaJ$pstate[wh=="Int",], perm = c(2,1,3:5), x = AaJ$time[wh=="Int"], col = clr, xlim = c(0,21), xaxs = "i", yaxs = "i") lines(AaJ$time[wh=="Int"], apply(AaJ$pstate[,1:3], 1, sum)[wh=="Int"], lwd = 4) mat2pol(AaJ$pstate[wh=="Conv",], perm = c(2,1,3:5), x = AaJ$time[wh=="Conv"], col = clr, xlim = c(21,0), xaxs = "i", yaxs = "i") lines(AaJ$time[wh=="Conv"], apply(AaJ$pstate[,1:3], 1, sum)[wh=="Conv"], lwd = 4) mtext(c("Int","Conv"), side = 3, at = c(1,3)/4, outer = TRUE, line = -2) ################################################### ### code chunk number 59: ms-steno2.rnw:1037-1044 ################################################### data(st2clin) str(st2clin) st2clin <- rename(cal.yr(st2clin), lex.id = id, per = doV) summary(st2clin) addmargins(table(table(st2clin$lex.id))) ################################################### ### code chunk number 60: ms-steno2.rnw:1050-1059 ################################################### S5 <- addCov.Lexis(S4, st2clin, "per") tt <- table(st2clin$lex.id) (who <- names(tt[tt == 3])[1]) subset(st2clin, lex.id == who) nround(subset(S5, lex.id == who, select = c(lex.id,per,tfi,tfc,exnam,a1c,chol,crea))) timeScales(S5) timeSince(S5) ################################################### ### code chunk number 61: ms-steno2.rnw:1078-1104 ################################################### detc <- glm.Lexis(S5, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst / allo + a1c + chol + log2(crea), from = c("Norm","Mic"), to = c("Mic","Mac")) impc <- glm.Lexis(S5, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst / allo + a1c + chol + log2(crea), to = c("Norm","Mic"), from = c("Mic","Mac")) round(ci.exp(detc), 3) round(ci.exp(impc), 3) moc <- glm.Lexis(S5, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst / allo + a1c + chol + log2(crea), to = "D(oth)") mCc <- glm.Lexis(S5, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst / allo + a1c + chol + log2(crea), to = "D(CVD)") round(ci.exp(moc), 3) round(ci.exp(mCc), 3) ################################################### ### code chunk number 62: ms-steno2.rnw:1149-1154 ################################################### St4 <- stack(S4) c(nrow(S4), nrow(St4)) table(S4$lex.Cst) table(St4$lex.Tr, St4$lex.Cst) ftable(St4$lex.Tr, St4$lex.Xst, St4$lex.Fail, col.vars = 2:3) ################################################### ### code chunk number 63: ms-steno2.rnw:1159-1162 ################################################### table(tt <- table(S4$lex.id)) nround(subset(S4 , lex.id == 102)[,1:8], 1) nround(subset(St4, lex.id == 102)[,1:9], 1) ################################################### ### code chunk number 64: ms-steno2.rnw:1170-1171 ################################################### cbind(with(subset(St4, grepl("D", lex.Tr)), table(lex.Tr))) ################################################### ### code chunk number 65: ms-steno2.rnw:1174-1182 ################################################### stD <- glm(cbind(lex.Fail, lex.dur) ~ Ns(tfi, knots = seq( 0, 20, 5)) * lex.Tr + Ns(age, knots = seq(50, 80, 10)) * lex.Tr + lex.Tr / allo + sex, family = poisreg, offset = log(lex.dur), data = subset(St4, grepl("D", lex.Tr))) round(ci.exp(stD)[,1,drop=F],3)