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