### R code from vignette source 'ms-steno2.rnw' ################################################### ### code chunk number 1: ms-steno2.rnw:2-5 ################################################### 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: ms-steno2.rnw:14-23 ################################################### library(survival) library(Epi) library(popEpi) # popEpi::splitMulti returns a data.frame rather than a data.table options("popEpi.datatable" = FALSE) library(tidyverse) setwd("c:/bendix/teach/AdvCoh/courses/Melb.2023/pracs") getwd() clear() ################################################### ### code chunk number 3: ms-steno2.rnw:35-40 ################################################### data(steno2) steno2 <- cal.yr(steno2) steno2 <- transform(steno2, doEnd = pmin(doEnd, doDth, na.rm = TRUE)) str(steno2) ################################################### ### code chunk number 4: ms-steno2.rnw:49-59 ################################################### 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, scale.R = 100) ################################################### ### code chunk number 5: ms-steno2.rnw:93-99 ################################################### data(st2alb) cut2 <- cal.yr(st2alb) names(cut2) names(cut2) <- c("lex.id", "cut", "new.state") str(cut2) head(cut2) ################################################### ### code chunk number 6: ms-steno2.rnw:103-105 ################################################### addmargins(table(table(cut2$lex.id))) cut2$lex.id %>% table %>% table %>% addmargins ################################################### ### code chunk number 7: boxL3 ################################################### L3 <- rcutLexis(L2, cut2) summary(L3) boxes(L3, boxpos = TRUE, cex = 0.8) ################################################### ### code chunk number 8: ms-steno2.rnw:140-144 ################################################### (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 9: ms-steno2.rnw:154-160 ################################################### 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 10: ms-steno2.rnw:167-170 ################################################### L4 <- rcutLexis(L3, xcut) L4 <- Relevel(L4, c("Norm","Mic","Mac","D(CVD)","D(oth)")) summary(L4) ################################################### ### code chunk number 11: b4 ################################################### clr <- c("forestgreen","orange","red","blue",gray(0.3)) boxes(L4, boxpos = list(x = c(20,20,20,80,80), y = c(10,50,90,75,25)), show.BE = "nz", scale.R = 100, digits.R = 2, cex = 0.9, pos.arr = 0.3, col.bg = clr, col.border = clr, col.txt = c("white","black")[c(1,2,1,1,1)]) ################################################### ### code chunk number 12: ms-steno2.rnw:212-215 ################################################### S4 <- splitMulti(L4, tfi = seq(0, 25, 1/12)) summary(L4) summary(S4) ################################################### ### code chunk number 13: ms-steno2.rnw:224-227 ################################################### ma <- glm.Lexis(S4, ~ Ns(tfi, knots = seq( 0, 20, 5)) + Ns(age, knots = seq(50, 80, 10)) + lex.Cst) ################################################### ### code chunk number 14: b4m ################################################### clr <- c("forestgreen","orange","red",gray(0.3)) summary(Relevel(L4, list("Dead" = 4:5), first = FALSE)) boxes(Relevel(L4, list("Dead" = 4:5), first = FALSE), boxpos = list(x = c(20,20,20,80), y = c(10,50,90,50)), show.BE = "nz", scale.R = 100, digits.R = 2, cex = 0.9, pos.arr = 0.3, col.bg = clr, col.border = clr, col.txt = c("white","black")[c(1,2,1,1)]) ################################################### ### code chunk number 15: ms-steno2.rnw:262-270 (eval = FALSE) ################################################### ## m1 <- 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:273-281 (eval = FALSE) ################################################### ## m2 <- 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:294-295 ################################################### round(ci.exp(ma), 2) ################################################### ### code chunk number 18: ms-steno2.rnw:307-308 ################################################### Wald(ma, subset = "lex.Cst") ################################################### ### code chunk number 19: ms-steno2.rnw:316-329 ################################################### # other causes of death 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) # # CVD death 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:335-337 ################################################### Wald(mo, subset = "Cst") Wald(mC, subset = "Cst") ################################################### ### code chunk number 21: ms-steno2.rnw:346-347 ################################################### summary(L2$age) ################################################### ### code chunk number 22: ms-steno2.rnw:355-357 ################################################### 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, 19, 0.5)), ain = c(45, 55, 65))[-1,], age = ain + tfi, lex.Cst = "Mic") prf[ 1: 5,] 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.02,20), plot = TRUE, xlab = "Age at follow-up (years)", ylab = "Mortality rate per 100 PY") abline(v = c(45, 55, 65), lty = 3) ################################################### ### code chunk number 24: mort3 ################################################### par(mfrow=c(1,3)) for(st in c("Norm","Mic","Mac")) { matshade(prf$age, pmin(pmax( cbind(ci.pred(mo, transform(prf, lex.Cst = st)), ci.pred(mC, transform(prf, lex.Cst = st))) * 100, 0.05), 60), 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:449-454 ################################################### 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:460-466 ################################################### 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:474-479 ################################################### 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:497-511 ################################################### 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("Mac","Mic"), to = c("Mic","Norm")) 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:547-558 ################################################### 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 30: ms-steno2.rnw:611-615 ################################################### ini <- L2[,c("per", "age", "tfi")] ini <- rbind(transform(ini, lex.Cst = factor("Mic"), allo = factor("Int")), transform(ini, lex.Cst = factor("Mic"), allo = factor("Conv"))) str(ini) ################################################### ### code chunk number 31: ms-steno2.rnw:625-632 ################################################### set.seed(1952) system.time( Sorg <- simLexis(Tr = Tr, # models for each transition init = ini, # cohort of straters N = 100, # how many copies of each person in ini t.range = 21, # how long should we simulate before censoring n.int = 200))# how many intervals for evaluating rates ################################################### ### code chunk number 32: ms-steno2.rnw:636-639 ################################################### Sorg <- Relevel(Sorg, c("Norm", "Mic", "Mac", "D(CVD)", "D(oth)")) summary(Sorg, by = "allo") subset(Sorg, lex.id %in% 28:32) ################################################### ### code chunk number 33: ms-steno2.rnw:644-645 ################################################### addmargins(table(table(Sorg$lex.id))) ################################################### ### code chunk number 34: ms-steno2.rnw:652-659 ################################################### system.time( Nst <- nState(Sorg, at = seq(0, 20, 0.2), from = 0, time.scale = "tfi")) str(Nst) head(Nst) ################################################### ### code chunk number 35: ms-steno2.rnw:664-675 ################################################### Nint <- nState(subset(Sorg, allo == "Int"), at = seq(0, 20.2, 0.1), from = 0, time.scale = "tfi") Nconv<- nState(subset(Sorg, allo == "Conv"), at = seq(0, 20.2, 0.1), from = 0, time.scale = "tfi") cbind( head(Nint), NA, head(Nconv)) ################################################### ### code chunk number 36: ms-steno2.rnw:682-686 ################################################### Pint <- pState(Nint ) Pconv <- pState(Nconv) str(Pint) head(Pint) ################################################### ### code chunk number 37: ms-steno2.rnw:694-698 ################################################### clr <- c("forestgreen", "orange", "red", "blue", gray(0.4)) par(mfrow = c(1,2), mar=c(3,3,2,2)) plot(Pint , col = clr, xlim = c(0, 20)) plot(Pconv, col = clr, xlim = c(20, 0)) ################################################### ### code chunk number 38: pStates ################################################### clr <- c("forestgreen", "orange", "red", "blue", gray(0.4)) par(mfrow = c(1,2), mar=c(3,3,2,2)) plot(Pint, col = clr, xlim = c(0, 20)) # the survival curve lines(as.numeric(rownames(Pint)), Pint[,"Mac"], lwd = 3, col = "black") lines(as.numeric(rownames(Pint)), Pint[,"Mac"], lwd = 1, col = "white") text(rownames(Pint)[150], Pint[150,] - diff(c(0, Pint[150,]))/2, colnames(Pint), col = "white", cex = 0.8) plot(Pconv, col = clr, xlim = c(20, 0)) # the survival curve lines(as.numeric(rownames(Pconv)), Pconv[,"Mac"], lwd = 3, col = "black") lines(as.numeric(rownames(Pconv)), Pconv[,"Mac"], lwd = 1, col = "white") text(rownames(Pconv)[150], Pconv[150,] - diff(c(0, Pconv[150,]))/2, colnames(Pint), col = "white", cex = 0.8) mtext(c("Intensive care","Conventional care"), side = 3, at = c(1,3)/4, outer = TRUE, line = -2) ################################################### ### code chunk number 39: ms-steno2.rnw:758-769 ################################################### 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 40: ms-steno2.rnw:783-791 ################################################### system.time( Sdef <- simLexis(Tr = Tr, init = ini, N = 1000, t.range = 21, n.int = 200)) summary(Sdef, by = "allo") subset(Sdef, lex.id < 5) ################################################### ### code chunk number 41: ms-steno2.rnw:799-805 ################################################### 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 42: ms-steno2.rnw:809-849 ################################################### P45c <- pState(nState(subset(Sdef, ain == 45 & allo == "Conv"), at = seq(0, 21, 0.1), from = 0, time.scale = "tfi")) P45i <- pState(nState(subset(Sdef, ain == 45 & allo == "Int"), at = seq(0, 21, 0.1), from = 0, time.scale = "tfi")) P50c <- pState(nState(subset(Sdef, ain == 55 & allo == "Conv"), at = seq(0, 21, 0.1), from = 0, time.scale = "tfi")) P50i <- pState(nState(subset(Sdef, ain == 55 & allo == "Int"), at = seq(0, 21, 0.1), from = 0, time.scale = "tfi")) P55c <- pState(nState(subset(Sdef, ain == 55 & allo == "Conv"), at = seq(0, 21, 0.1), from = 0, time.scale = "tfi")) P55i <- pState(nState(subset(Sdef, ain == 55 & allo == "Int"), at = seq(0, 21, 0.1), from = 0, time.scale = "tfi")) P60c <- pState(nState(subset(Sdef, ain == 55 & allo == "Conv"), at = seq(0, 21, 0.1), from = 0, time.scale = "tfi")) P60i <- pState(nState(subset(Sdef, ain == 55 & allo == "Int"), at = seq(0, 21, 0.1), from = 0, time.scale = "tfi")) P65c <- pState(nState(subset(Sdef, ain == 65 & allo == "Conv"), at = seq(0, 21, 0.1), from = 0, time.scale = "tfi")) P65i <- pState(nState(subset(Sdef, ain == 65 & allo == "Int"), at = seq(0, 21, 0.1), from = 0, time.scale = "tfi")) ################################################### ### code chunk number 43: 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 44: ms-steno2.rnw:887-893 ################################################### (ain <- seq(45, 65, 5)) (alo <- levels(S4$allo)) pdef <- NArray(c(list(ain = ain, allo = alo), dimnames(P45i))) str(pdef) ################################################### ### code chunk number 45: ms-steno2.rnw:900-909 ################################################### for(aa in ain) for(gg in alo) 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 46: panel3 ################################################### ain <- seq(45, 65, 10) par(mfrow = c(length(ain),2), mar = c(1,2,1,2)/5, oma = c(2,4,2,3), mgp = c(3,1,0) / 1.6) for(aa in ain) { mat2pol(pdef[paste(aa),"Int" ,,], col = clr, xlim = c(0,20), xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n") axis(side = 2) axis(side = 4, labels = NA) mat2pol(pdef[paste(aa),"Conv",,], col = clr, xlim = c(20,0), xaxs = "i", yaxs = "i", xaxt = "n", yaxt = "n") axis(side = 2, labels = NA) axis(side = 4) } 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 = 2) ################################################### ### code chunk number 47: ms-steno2.rnw:969-972 ################################################### AaJepi <- AaJ.Lexis(L4, timeScale = "tfi") class(AaJepi) AaJepi ################################################### ### code chunk number 48: ms-steno2.rnw:977-979 ################################################### AaJepi$transitions summary(L4, simp = FALSE) ################################################### ### code chunk number 49: ms-steno2.rnw:985-990 ################################################### names(AaJepi) AaJepi$states head(AaJepi$pstate) head(AaJepi$lower) head(AaJepi$upper) ################################################### ### code chunk number 50: AaJ ################################################### par(mfrow = c(1,1)) mat2pol(AaJepi$pstate, perm = c(1:3,5,4), x = AaJepi$time, col = clr) lines(AaJepi$time, apply(AaJepi$pstate[,1:3], 1, sum), lwd = 5) ################################################### ### code chunk number 51: ms-steno2.rnw:1002-1004 ################################################### AaJallo <- AaJ.Lexis(L4, ~ allo, timeScale = "tfi") AaJallo ################################################### ### code chunk number 52: ms-steno2.rnw:1010-1014 ################################################### AaJallo$states AaJallo$strata wh <- rep(substr(names(AaJallo$strata), 6, 9), AaJallo$strata) table(wh) ################################################### ### code chunk number 53: AaJstates ################################################### par(mfrow = c(1,2), mar=c(3,3,2,0), oma = c(0,0,0,3), las = 1, bty = "n") mat2pol(AaJallo$pstate[wh=="Int",], x = AaJallo$time[wh=="Int"], col = clr, xlim = c(0,21), xaxs = "i", yaxs = "i") lines(AaJallo$time[wh=="Int"], apply(AaJallo$pstate[,1:3], 1, sum)[wh=="Int"], lwd = 4) mat2pol(AaJallo$pstate[wh=="Conv",], x = AaJallo$time[wh=="Conv"], col = clr, xlim = c(21,0), xaxs = "i", yaxs = "i") lines(AaJallo$time[wh=="Conv"], apply(AaJallo$pstate[,1:3], 1, sum)[wh=="Conv"], lwd = 4) axis(side = 4) mtext(c("Int","Conv"), side = 3, at = c(1,3)/4, outer = TRUE, line = -2) ################################################### ### code chunk number 54: ms-steno2.rnw:1080-1085 ################################################### tLive <- xtabs(lex.dur ~ ain + allo + lex.Cst, data = Sdef) / nid(Sdef) * 10 str(mtLive <- addmargins(tLive, 3)) round(ftable(mtLive , col.vars = c(3,2)), 1) round(ftable(mtLive[,,-(4:5)], col.vars = c(3,2)), 1) ################################################### ### code chunk number 55: ms-steno2.rnw:1089-1090 ################################################### round((mtLive[,"Int",-(4:5)] - mtLive[,"Conv",-(4:5)]), 1) ################################################### ### code chunk number 56: ms-steno2.rnw:1100-1106 ################################################### tDead <- xtabs((20 - tfi - lex.dur) ~ ain + allo + lex.Xst, data = subset(Sdef, lex.Xst %in% c("D(oth)", "D(CVD)"))) / nid(Sdef) * 10 str(mtDead <- addmargins(tDead[,,4:5], 3)) round(ftable(mtDead , col.vars = c(3,2)), 1) round(ftable(mtDead[,,-(4:5)], col.vars = c(3,2)), 1) ################################################### ### code chunk number 57: ms-steno2.rnw:1111-1112 ################################################### round((mtDead[,"Int",] - mtDead[,"Conv",]), 1) ################################################### ### code chunk number 58: ms-steno2.rnw:1128-1135 ################################################### data(st2clin) str(st2clin) st2clin <- cal.yr(st2clin) names(st2clin) names(st2clin)[1:2] <- c("lex.id","per") summary(st2clin) addmargins(table(table(st2clin$lex.id))) ################################################### ### code chunk number 59: ms-steno2.rnw:1141-1150 ################################################### S5 <- addCov.Lexis(S4, st2clin, "per") tt <- table(st2clin$lex.id) (who <- names(tt[tt == 3])[1]) subset(st2clin, lex.id == who) subset(S5, lex.id == who, select = c(lex.id,per,tfi,tfc,exnam,a1c,chol,crea)) timeScales(S5) timeSince(S5) ################################################### ### code chunk number 60: ms-steno2.rnw:1169-1182 ################################################### 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) ################################################### ### code chunk number 61: ms-steno2.rnw:1186-1188 ################################################### round(ci.exp(detc, subset = -(1:8), pval = T), 3) round(ci.exp(impc, subset = -(1:8), pval = T), 3) ################################################### ### code chunk number 62: ms-steno2.rnw:1191-1203 ################################################### 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, subset = -(1:8), pval = T), 3) round(ci.exp(mCc, subset = -(1:8), pval = T), 3) ################################################### ### code chunk number 63: ms-steno2.rnw:1249-1254 ################################################### 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 64: ms-steno2.rnw:1260-1262 ################################################### subset(S4 , lex.id == 102)[,1:8] subset(St4, lex.id == 102)[,1:9] ################################################### ### code chunk number 65: ms-steno2.rnw:1271-1272 ################################################### cbind(with(subset(St4, grepl("D", lex.Tr)), table(lex.Tr))) ################################################### ### code chunk number 66: ms-steno2.rnw:1276-1283 ################################################### 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, data = subset(St4, grepl("D", lex.Tr))) round(ci.exp(stD)[,1,drop=F],3)