### R code from vignette source 'fudat' ### Encoding: UTF-8 ################################################### ### code chunk number 1: fudat.rnw:2-6 ################################################### options( width=90, # prompt=" ", continue=" ", 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: fudat.rnw:44-47 ################################################### library( Epi ) data( DMlate ) str( DMlate ) ################################################### ### code chunk number 3: fudat.rnw:55-56 ################################################### head( DMlate ) ################################################### ### code chunk number 4: fudat.rnw:59-62 ################################################### with( DMlate, c(sum(!is.na(dodth)), sum(dox-dodm), sum(!is.na(dodth)) / sum(dox-dodm) * 1000 ) ) ################################################### ### code chunk number 5: fudat.rnw:278-284 ################################################### data( DMlate ) Ldm <- Lexis( entry = list( age=dodm-dobth ), exit = list( age=dox-dobth ), exit.status = factor( !is.na(dodth), labels=c("Alive","Dead") ), data = DMlate ) summary( Ldm ) ################################################### ### code chunk number 6: fudat.rnw:310-315 ################################################### Ldm <- Lexis( entry = list( age=dodm-dobth, per=dodm ), exit = list( age=dox-dobth ), exit.status = factor( !is.na(dodth), labels=c("Alive","Dead") ), data = DMlate ) ################################################### ### code chunk number 7: fudat.rnw:332-335 ################################################### oo <- options(digits=5) head( Ldm[,-(10:11)] ) options( oo ) ################################################### ### code chunk number 8: Ldm-default ################################################### plot( Ldm ) ################################################### ### code chunk number 9: Ldxtra ################################################### ypi <- 11 # years per inch pdf( "./graph/fudat-Ldm-xtra.pdf", width=25/ypi+1, height=100/ypi+1 ) par( mai=c(3,3,1,1)/4, mgp=c(3,1,0)/1.6 ) plot( Ldm, 2:1, lwd=1, grid=list(1990+0:5*5,1:19*5), lty.grid=1, col.grid=gray(0.7), xlim=c(1990,2015), xaxs="i", ylim=c( 0, 100), yaxs="i", las=1 ) points( Ldm, 2:1, pch=c(NA,16)[Ldm$lex.Xst], col="black", cex=0.6 ) dev.off() ################################################### ### code chunk number 10: fudat.rnw:466-467 ################################################### library( popEpi ) ################################################### ### code chunk number 11: fudat.rnw:469-477 ################################################### Xdm <- splitLexis( Ldm, "age", breaks=0:100 ) library( popEpi ) Mdm <- splitMulti( Ldm, age=0:100 ) summary( Ldm ) summary( Xdm ) summary( Mdm ) subset( Ldm, lex.id==151 )[,1:6] subset( Mdm, lex.id==151 )[,1:6] ################################################### ### code chunk number 12: fudat.rnw:512-518 ################################################### astab <- xtabs( cbind(D=(lex.Xst=="Dead"), Y=lex.dur) ~ floor(age) + sex, data=Mdm ) str( astab ) round( ftable( astab[51:55,,], row.vars=1 ),1 ) ################################################### ### code chunk number 13: fudat.rnw:542-548 ################################################### a1 <- 0:99 # the 1-year age-classes mold <- glm( astab[,"M","D"] ~ factor(a1) - 1, offset = log(astab[,"M","Y"]), family = poisson ) mm <- glm( astab[,"M",] ~ factor(a1) - 1, family = poisreg ) mf <- glm( astab[,"F",] ~ factor(a1) - 1, family = poisreg ) ################################################### ### code chunk number 14: fudat.rnw:581-585 ################################################### cbind( poisson = coef(mold), poisreg = coef(mm), astab[,"M",], empirical = log(astab[,"M","D"]/astab[,"M","Y"]) )[35:38,] ################################################### ### code chunk number 15: fudat.rnw:612-617 ################################################### A <- 0:99+0.5 # midpoints of age-classes sm <- glm( astab[,"M",] ~ Ns( A, knots=c(20,40,60,70,80,90) ), family=poisreg ) sf <- glm( astab[,"F",] ~ Ns( A, knots=c(20,40,60,70,80,90) ), family=poisreg ) ################################################### ### code chunk number 16: fudat.rnw:629-636 ################################################### nd <- data.frame( A=seq(10,99,0.5) ) pm <- ci.pred( sm, nd ) pf <- ci.pred( sf, nd ) matshade( nd$A, cbind( pm, pf )*1000, plot=TRUE, col=c("blue","red"), lwd=2, log="y", ylim=c(0.5,300), las=1, xlab="Age", ylab="DM mortality per 1000 PY" ) ################################################### ### code chunk number 17: smooth-rates ################################################### summary( Mdm, t=T ) xm <- glm( cbind(lex.Xst=="Dead",lex.dur) ~ Ns( age, knots=c(20,40,60,70,80,90) ), family=poisreg, data=subset(Mdm,sex=="M") ) xf <- glm( cbind(lex.Xst=="Dead",lex.dur) ~ Ns( age, knots=c(20,40,60,70,80,90) ), family=poisreg, data=subset(Mdm,sex=="F") ) nd <- data.frame( age=seq(10,99,0.5) ) xpm <- ci.pred( xm, nd ) xpf <- ci.pred( xf, nd ) matshade( nd$age, cbind( pm, pf, xpm, xpf )*1000, plot=TRUE, col=c("blue","red"), lwd=2, lty=c(1,1,2,2), log="y", ylim=c(0.5,400), las=1, xlab="Age", ylab="DM mortality per 1000 PY" ) matpoints( 0:99+0.5, cbind( pmax(astab[,"M","D"]/astab[,"M","Y"],0.00042), pmax(astab[,"F","D"]/astab[,"F","Y"],0.00045) )*1000, pch=16, col=c("blue","red"), cex=0.4, type="o", lty=1 ) ################################################### ### code chunk number 18: fudat.rnw:736-743 ################################################### names( Mdm ) Cdm <- cutLexis( data = Ldm, cut = Ldm$dooad, timescale = "per", new.state = "OAD", precursor.states = "Alive" ) subset( Cdm, lex.id %in% 35:38 )[,c(1:6,11)] ################################################### ### code chunk number 19: fudat.rnw:761-764 ################################################### CMdm <- cutLexis( data=Mdm, cut=Mdm$dooad, timescale="per", new.state="OAD", precursor.states="Alive" ) subset( CMdm, lex.id %in% 35:38 )[,c(1:6,11)] ################################################### ### code chunk number 20: fudat.rnw:809-827 ################################################### clin <- data.frame( lex.id = c(35,35,36,36,37), ex = c("one","two")[c(1,2,1,2,2)], per = c(2003.1,2005.2,2002.8,2005.4,2006.3), sbp = c(120,134,149,166,120), wt = c(90,93,107,103,75) ) clin Zdm <- subset( CMdm, lex.id %in% 35:38 )[,1:6] Kdm <- addCov.Lexis( Zdm, # Lexis data frame to be expanded clin, # dataframe with id, date and measures timescale = "per",# time scale variable exnam = "ex", # exam name variable addScales = TRUE ) od <- options(digits=4) Kdm[,c(timeScales(Kdm), "lex.id","lex.dur","lex.Cst","lex.Xst", "ex","sbp","wt")] timeSince(Kdm) options( od )