###################################################
### chunk number 1: 
###################################################
library( Epi )


###################################################
### chunk number 2: 
###################################################
st <- read.table( "../data/pm-dk.txt", header=T, as.is=T,
                 na.strings="." )
st
str( st )


###################################################
### chunk number 3: 
###################################################
# Change the character variables with dates to fractional calendar
# years
for( i in 2:5 ) st[,i] <- cal.yr( st[,i], format="%d/%m/%Y" )
st$exit[nrow(st)] <- cal.yr(Sys.Date())
# Attach the data for those still alive
st$fail <- !is.na(st$death)
st[!st$fail,"death"] <- 2009
st


###################################################
### chunk number 4: 
###################################################
attach( st )


###################################################
### chunk number 5: 
###################################################
# Lexis object
L <- Lexis( entry = list(per=birth), 
             exit = list(per=death,age=death-birth),
             exit.status = fail, 
             data = st )
# Plot Lexis diagram
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, xaxt="n" ) # Omit x-labels
plot( L, xlim=c(1945,2010), ylim=c(32,88), lwd=3, las=1,grid=0:20*5, col="black",
xlab="Calendar time", ylab="Age" )
points( L, pch=c(NA,16)[L$lex.Xst+1] )
#put names of the prime ministers on the plot
text( death, death-birth, Name, adj=c(1.05,-0.05), cex=0.7 )
par( xaxt="s" )
axis( side=1, at=seq(1950,2010,10) ) # x-labels at nice places


###################################################
### chunk number 6: 
###################################################
# New Lexis object describing periods in an office 
# and lines added to a picture
in_office <- c( rep(FALSE,nrow(st)-1), TRUE )
st <- cbind( st, in_office )
Lo <- Lexis( entry = list(per=entry), 
              exit = list(per=exit,age=exit-birth),
              exit.status = in_office, 
              data = st )
lines( Lo, lwd=3, las=1, col="red" )
# the same may be plotted using command segments
box()
segments( birth, 0, death, death-birth, lwd=2 )
segments( entry, entry-birth, exit, exit-birth, lwd=4, col="red" )


###################################################
### chunk number 7: 
###################################################
abline( h=50 )


###################################################
### chunk number 8: 
###################################################
age_entry <- Lo$age
age_exit <- Lo$age+Lo$lex.dur
n_birthday<- sum( ( age_entry<50) & ( age_exit>50 ) )
n_birthday


###################################################
### chunk number 9: 
###################################################
abline( v=cal.yr( "2/10/1972", format="%d/%m/%Y" ) )


###################################################
### chunk number 10: 
###################################################
alive <- (L$death >=2004)
n_alive <- sum(alive)
n_alive
#Anker Jorgensen - 1 person has got 2 lex.id's
levels( as.factor( subset( L$Name,alive==T ) ) )


###################################################
### chunk number 11: 
###################################################
# New lexis object - since entry to the office to the death
Ln <- Lexis( entry=list(per=entry), exit=list(per=death,age=death-entry),
exit.status=fail, data=st )
ny <- 2008-1945
n_alive <- vector( "numeric", ny )
for (i in 1:ny)
{
alive <- ( (Ln$death >=(1944+i)) & (Ln$entry<=(1944+i)) )
n_alive[i] <- nlevels( as.factor( subset( Ln$Name, alive==T ) ) )
}


###################################################
### chunk number 12: MAXALIVE
###################################################
plot( n_alive~seq(1945,(1945+ny-1),1), type="l", xlab="Calendar year", 
ylab="Number of former prime ministers alive" )


###################################################
### chunk number 13: 
###################################################
rect( 1970, 50, 1990, 70, lwd=2, border="green",col="lightgreen"  )	


###################################################
### chunk number 14: 
###################################################
polygon( c(1955,2005,2005,1965,1955), c(30,80,70,30,30), lwd=2,
     border="blue", col="lightblue" )
# Now draw the Lexis diagram again on top of the shaded areas


###################################################
### chunk number 15: 
###################################################
# Prime-minister years lived by persons
#  aged 50 to 70 in the period 1 January 1970 through 1 January 1990.
x1 <- splitLexis ( Lo ,breaks=c(0,50,70,100), time.scale="age" )
x2 <- splitLexis ( x1, breaks=c(1900,1970,1990,2010), time.scale="per" )
summary( x2 )
tapply( status(x2,"exit")==1, list( timeBand(x2,"age","left"),
                                    timeBand(x2,"per","left") ), sum )
tapply( dur(x2),  list( timeBand(x2,"age","left"),
                        timeBand(x2,"per","left") ), sum )
# Computing the person-years in the 1925-35 cohort
x3 <- subset( Lo , birth>1925 & birth<=1935 )
summary( x3 )
dur( x3 )
# Computing person years in the intersection
x4 <- subset( x2 , birth>1925 & birth<=1935 )
summary( x4 )
dur( x4 )