# ========================================================== 
# A.4.2 Lexis Diagrams
# R code to create Figure A.3
#
# December 11, 2018
# ==========================================================


events <- NULL
events$id   <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
events$bday <- c(1942, 1904, 1965, 1910, 1923, 1952, 1938, 1921, 1925, 1949)
events$age.entry <- c(39.461, 78.689, 15.540, 70.286, 60.011, 28.753, 42.705, 59.398, 55.680, 33.667)
events$age.ps  <- c(18, 22, 15, 50, 52, 12, 37, 42, 48, 15)
events$age.psa <- c(33, 31, 15, 50, 52, 19, 42, 42, 48, 28)
events$age.last.contact <- c(74.853, 86.746, 52.107, 80.821, 74.229, 65.227, 63.184, 83.077, 59.469, 45.990)
events$death.status <- c(0, 1, 0, 1, 1, 0, 1, 1, 1, 1)


library(Epi)

followup <- Lexis(entry=list(cal = age.entry + bday, age = age.entry), 
                  exit=list(cal = age.last.contact + bday), 
                  exit.status = death.status, id=id, 
                  data=events)
print(followup[,c("cal","age","lex.dur","lex.Cst","lex.Xst","lex.id")])

events$age.exit <- ifelse(events$age.entry <= events$age.ps, events$age.ps + 0.5, events$age.entry)
psc.to.entry <- Lexis(entry=list(cal = age.ps + bday, age = age.ps), 
                      exit=list(cal = age.exit + bday), 
                      exit.status = 1, id=id, data=events)

events$age.exit <- ifelse(events$age.ps == events$age.psa, events$age.psa + 0.5, events$age.psa)
psc.to.psa <- Lexis(entry=list(cal = age.ps + bday, age = age.ps), 
                    exit=list(cal = age.exit + bday), 
                    exit.status = 1, id=id, data=events)


pdf("figureA3.pdf", width=10, height=6) 
par(mai=c(0.8,0.8,0.3,0.3))

plot(0, 0, type="n", axes=F, xlim=c(1920, 2020), ylim=c(0,100), xlab="", ylab="")
axis(side=1, at=seq(1920,2020,by=10), label=T)
axis(side=2, at=seq(0,100,by=20), label=T, las=2, adj=1)
mtext("CALENDAR  YEAR", side=1, line=2.6)
mtext("AGE  OF  INDIVIDUAL", side=2, line=2.8)

lines(psc.to.psa, lty=1, lwd=2)
lines(psc.to.entry, lty=2, lwd=1)
points(psc.to.entry, pch=1, cex=1)
lines(followup, lty=2, lwd=1)
points(subset(followup, status(followup) %in% 1), pch=19, cex=1)
box()

dev.off()