# R code for graphing 2D
# As used in "The Importance of Year of Birth in Two-Dimensional Mortality Data", October 2005.

# Hacked version of filled.contour() which drops the key and is therefore better able to handle missing values.
newcontour<-function (x = seq(0, 1, len = nrow(z)), y = seq(0, 1, len = ncol(z)), 
    z, xlim = range(x, finite = TRUE), ylim = range(y, finite = TRUE), 
    zlim = range(z, finite = TRUE), levels = pretty(zlim, nlevels), 
    nlevels = 20, color.palette = cm.colors, col = color.palette(length(levels) - 
        1), plot.title, plot.axes, key.title, key.axes, asp = NA, 
    xaxs = "i", yaxs = "i", las = 1, axes = TRUE, frame.plot = axes, 
    ...) 
{
    if (missing(z)) {
        if (!missing(x)) {
            if (is.list(x)) {
                z <- x$z
                y <- x$y
                x <- x$x
            }
            else {
                z <- x
                x <- seq(0, 1, len = nrow(z))
            }
        }
        else stop("no 'z' matrix specified")
    }
    else if (is.list(x)) {
        y <- x$y
        x <- x$x
    }
    if (any(diff(x) <= 0) || any(diff(y) <= 0)) 
        stop("increasing 'x' and 'y' values expected")
    plot.new()
    plot.window(xlim, ylim, "", xaxs = xaxs, yaxs = yaxs, asp = asp)
    if (!is.matrix(z) || nrow(z) <= 1 || ncol(z) <= 1) 
        stop("no proper 'z' matrix specified")
    if (!is.double(z)) 
        storage.mode(z) <- "double"
    .Internal(filledcontour(as.double(x), as.double(y), z, as.double(levels), 
        col = col))
    if (missing(plot.axes)) {
        if (axes) {
            title(main = "", xlab = "", ylab = "")
            Axis(x, side = 1)
            Axis(y, side = 2)
        }
    }
    else plot.axes
    if (frame.plot) 
        box()
    if (missing(plot.title)) 
        title(...)
    else plot.title
    invisible()
}

# Read in the data, age x year of observation.
r<-read.csv("f:\\Research\\International\\EW_females_ap_YoBvYoO.csv", header=TRUE)
EWfemalerates<-cbind(r$X1962,r$X1963,r$X1964,r$X1965,r$X1966,r$X1967,r$X1968,r$X1969,r$X1970,r$X1971,r$X1972,r$X1973,r$X1974,r$X1975,r$X1976,r$X1977,r$X1978,r$X1979,r$X1980,r$X1981,r$X1982,r$X1983,r$X1984,r$X1985,r$X1986,r$X1987,r$X1988,r$X1989,r$X1990,r$X1991,r$X1992,r$X1993,r$X1994,r$X1995,r$X1996,r$X1997,r$X1998,r$X1999,r$X2000,r$X2001,r$X2002,r$X2003)
r<-read.csv("f:\\Research\\International\\EW_males_ap_YoBvYoO.csv", header=TRUE)
EWmalerates<-cbind(r$X1962,r$X1963,r$X1964,r$X1965,r$X1966,r$X1967,r$X1968,r$X1969,r$X1970,r$X1971,r$X1972,r$X1973,r$X1974,r$X1975,r$X1976,r$X1977,r$X1978,r$X1979,r$X1980,r$X1981,r$X1982,r$X1983,r$X1984,r$X1985,r$X1986,r$X1987,r$X1988,r$X1989,r$X1990,r$X1991,r$X1992,r$X1993,r$X1994,r$X1995,r$X1996,r$X1997,r$X1998,r$X1999,r$X2000,r$X2001,r$X2002,r$X2003)
r<-read.csv("f:\\Research\\International\\DE_males_ac_YoBvYoO.csv", header=TRUE)
DEmalerates<-cbind(r$X1962,r$X1963,r$X1964,r$X1965,r$X1966,r$X1967,r$X1968,r$X1969,r$X1970,r$X1971,r$X1972,r$X1973,r$X1974,r$X1975,r$X1976,r$X1977,r$X1978,r$X1979,r$X1980,r$X1981,r$X1982,r$X1983,r$X1984,r$X1985,r$X1986,r$X1987,r$X1988,r$X1989,r$X1990,r$X1991,r$X1992,r$X1993,r$X1994,r$X1995,r$X1996,r$X1997,r$X1998,r$X1999,r$X2000,r$X2001,r$X2002,r$X2003)
r<-read.csv("f:\\Research\\International\\DE_females_ac_YoBvYoO.csv", header=TRUE)
DEfemalerates<-cbind(r$X1962,r$X1963,r$X1964,r$X1965,r$X1966,r$X1967,r$X1968,r$X1969,r$X1970,r$X1971,r$X1972,r$X1973,r$X1974,r$X1975,r$X1976,r$X1977,r$X1978,r$X1979,r$X1980,r$X1981,r$X1982,r$X1983,r$X1984,r$X1985,r$X1986,r$X1987,r$X1988,r$X1989,r$X1990,r$X1991,r$X1992,r$X1993,r$X1994,r$X1995,r$X1996,r$X1997,r$X1998,r$X1999,r$X2000,r$X2001,r$X2002,r$X2003)

# Bind together the columns in the data file.  Bit repetitive this, but you only have to do it once.

# Create vectors for the years of observation and years of birth you want to graph.
years<-c(seq(1962,2003))
births<-c(seq(1900,1960))

# Set the text size (0.8 = 80% of normal)
textsize<-0.8

# Set plotting margins around the graph
par(mar=c(4, 4, 2, 2), mfrow=c(2,2))

# Plot 2d image
levels<-c(-10.0, 0.00, 0.01, 0.02, 0.03, 0.04)
labels<-c("0%", "1%", "2%", "3%", "4%")
newcontour(births, years, EWmalerates, levels=levels, col=c(gray(0:6 / 6)), axes=FALSE, font.lab=2, cex.lab=1.2, ylab="Year of observation", xlab="Year of birth", frame=FALSE, lwd=2)
axis(1, font=2)
axis(2, font=2)
axis(3, font=2)
# This draws a dashed line at age 65
lines(x=c(1800, 2000), y=c(1865, 2065), lty=2, lwd=2)
text(x=1945, y=1970, "E&W, males", font=2)

newcontour(births, years, EWfemalerates, levels=levels, col=c(gray(0:6 / 6)), axes=FALSE, font.lab=2, cex.lab=1.2, ylab="Year of observation", xlab="Year of birth", frame=FALSE, lwd=2)
axis(1, font=2)
axis(2, font=2)
axis(3, font=2)
axis(4, font=2)
# This draws a dashed line at age 65
lines(x=c(1800, 2000), y=c(1865, 2065), lty=2, lwd=2)
text(x=1945, y=1970, "E&W, females", font=2)

par(mar=c(4, 4, 0, 2))
newcontour(births, years, DEmalerates, levels=levels, col=c(gray(0:6 / 6)), axes=FALSE, font.lab=2, cex.lab=1.2, ylab="Year of observation", xlab="Year of birth", frame=FALSE, lwd=2)
axis(1, font=2)
axis(2, font=2)
# This draws a dashed line at age 65
lines(x=c(1800, 2000), y=c(1865, 2065), lty=2, lwd=2)
text(x=1945, y=1970, "Germany, males", font=2)

newcontour(births, years, DEfemalerates, levels=levels, col=c(gray(0:6 / 6)), axes=FALSE, font.lab=2, cex.lab=1.2, ylab="Year of observation", xlab="Year of birth", frame=FALSE, lwd=2)
axis(1, font=2)
axis(2, font=2)
axis(4, font=2)
# This draws a dashed line at age 65
lines(x=c(1800, 2000), y=c(1865, 2065), lty=2, lwd=2)
legend(x=1940, y=1980, pch=rep(19, 6), col=c(gray(0:6 / 6)), legend=c("<0%", "0-1%", "1-2%", "2-3%", "3-4%", ">4%"), bty="n")

# You now have an R graph.  You can maximise the size of the window and pick "Save as->PDF" to
# save the graph as print-quality graphics.

par(mar=c(4.5, 4, 0, 0))
d<-read.csv("c:\\Research\\International\\cohort_v_splines.csv", header=TRUE)
plot(x=d$Age, y=d$Actual, pch=20, frame=FALSE, xlab="Age", ylab="Annual improvement", axes=FALSE, font=2, ylim=range(0.0, 0.05), lwd=2)
axis(1, las=1, font=2)
axis(2, las=1, font=2, labels=c("0%", "1%", "2%", "3%", "4%", "5%"), at=c(0.0, 0.01, 0.02, 0.03, 0.04, 0.05))
points(x=d$Age, y=d$Medium, pch=21, lwd=2)
lines(x=d$Age, y=d$P.spline, lwd=2)
legend(x=50, y=0.02, pch=c(20, 21, NA), lty=c(NA, NA, 1), lwd=c(2, 2, 2), legend=c("Smoothed actual improvements", "Medium-cohort projection", "P-spline projection"), bty="n")