# Read in some life table data
library( foreign )
lAUS <- read.dta( "../data/ltab-AUS.dta" )
lDK  <- read.dta( "../data/ltab-DK.dta" )
lFIN <- read.dta( "../data/ltab-FIN.dta" )

str( lAUS )
str( lDK )
str( lFIN )

# Function to plot the empirical mortality rates along with the estimated
# population mortality
plot.ltab <-
function( ltab, tabnam, tabper )
{
if( sum( !is.na( match( c("Age","MortM","MortF"), names( ltab ) ) ) ) < 3 )
stop( "Names 'Age','MortM' and 'MortF' must all be in the dataframe" )

# Regression of the log-mortality rates
mort <- with( ltab, c( MortM, MortF ) )
age <- with( ltab, c( Age, Age ) )
middle <- ( age > 40 & age < 85 )
sex <- factor( rep( c("m","f"), rep( nrow(ltab), 2 ) ) )
l.fin <- lm( log2( mort / 10^5 ) ~ - 1 + sex + age, subset=middle )
rm( mort, age, sex )
# Population mortality measured in deaths per year described as:
# log( lambda ) = \alpha_s + beta * a
# alpha_m, alpha_k, beta are the parameters from l.fin.
#
ap.m <- ci.lin( l.fin, subset="sexm" )[,1]
ap.f <- ci.lin( l.fin, subset="sexf" )[,1]
bp   <- ci.lin( l.fin, subset="age" )[,1]

# Plot of the empirical rates and the fitted lines and their formula.
plot( ltab$MortM ~ ltab$Age, log="y", type="n", col="blue",
      ylim=c(5,30000), xlab="Age", ylab="Mortality per 100,000 person years" )
abline( h=c(1:10,1:10*10,1:10*100,1:10*1000,1:10*10000),
        v=seq(0,100,5), col=gray(0.8) )
lines( ltab$MortM ~ ltab$Age, lwd=3, cex=0.7, col="blue" )
lines( ltab$MortF ~ ltab$Age, lwd=3, cex=0.7, col="red" )
abline( log10( 2 ) * ap.m + 5, log10( 2 ) * bp, col="blue" )
abline( log10( 2 ) * ap.f + 5, log10( 2 ) * bp, col="red" )
abline( v=c(40,85), col=gray(0.8), lwd=2 )
bred <- strwidth( expression( log[2]*"[mortality per "*10^5*" (40-85 years)]" ) )
la <- floor( ( 100 - bred ) / 5 ) * 5
rect( la, 1, 110, 20, border=gray( 0.8 ), col="white" )
bred <- strwidth( paste( tabnam, "life tables", tabper ) )
la <- ceiling( ( bred ) / 5 ) * 5
rect( -5, 20000, la, 50000, border=gray( 0.8 ), col="white" )
box()
text( -1, 30000, paste( tabnam, "life tables", tabper ), adj=c(0,0.5) )
text( cnr( 97, 15 ), expression( log[2]*"[mortality per "*10^5*" (40-85 years)]" ),
                     adj=1, font=2 )
text( cnr( 97, 10 ), paste(   "Men:", formatC( ap.m, format="f", digits=3 ),
                                 "+", formatC( bp  , format="f", digits=3 ),
                               "age" ), col="blue", font=2, adj=c(1,0.5) )
text( cnr( 97,  5 ), paste( "Women:", formatC( ap.f, format="f", digits=3 ),
                                 "+", formatC( bp  , format="f", digits=3 ),
                               "age" ), col= "red", font=2, adj=c(1,0.5) )
}

plt( "ltab-AUS", width=10 )
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.5 )
plot.ltab( subset(lAUS,Period=="1997-1999"), "Australian", "1997-1999" )
dev.off()

plt( "ltab-DK", width=10 )
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.5 )
plot.ltab( subset(lDK, Period=="1997-1998"), "Danish", "1997-1998" )
dev.off()

plt( "ltab-FI", width=10 )
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.5 )
plot.ltab( lFIN, "Finnish", "1998" )
dev.off()

# Get the Swedish mortality data subset them to 1997-1999 and make a lifetable
mSE <- read.dta( "../data/mort-SE.dta" )
subSE <- subset( mSE, P>1996 & P<2000 )
mort <- with( subSE, tapply( D, list(A,sex), sum ) ) /
        with( subSE, tapply( Y, list(A,sex), sum ) ) * 10^5
lSE <- as.data.frame( mort )
names( lSE ) <- c("MortM","MortF")
lSE$Age <- as.numeric( rownames( lSE ) )
str( lSE )

plt( "ltab-SE", width=10 )
par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.5 )
plot.ltab( lSE, "Swedish", "1997-99" )
dev.off()