Proposal for adjusted lexis surface made up of equilateral Lexis triangles

EqLexis.R

# the function:

Rates must be an all-numeric MATRIX with 4 columns: "Age", "Year", "Cohort" and a 4th column, such as "ASFR",
col = a vector of colors produced by your favorite color ramp (there is a default inside the function). must be 1 shorter than the number of breaks
breaks = the cut points for determinig colors- must be 1 longer than the number of colors
... additional arguments to pass to plot()

function no works fine with pre-logged mortality data. 

Legend (from fields package) now goes to right margin automatically, and you can pass on arguments to make it friendly for logged data. example to come soon. see ?fields:::image.plot for an idea of how to pass special arguments to the legend.

Also new: multiple adjusment parameters for the different axis labels and tick lables, each is a 2 element vector indicating c(x shift,yshift). For most HFD fertility data labels automatically get pretty close to where they need to be, but for larger surfaces you'll need to tinker so that it looks good. 

You can now also change the reference line spacing to be every 10 years or every 5 years (or any other multiple). 5 or 10 are the only ones that make sense. If you want something else then specify guides = FALSE and draw them yourself using segments(). some pointers for hand-drawing triangle grid lines:

the y coordinates need to be scaled by sin(pi/3) (or sqrt(3)/2, whichever is easier to remember)
each x coordinate of the period lines needs to be shifted left by .5 times it's age position. e.g. a 1960 period line going from age 12 to age 55 would have x coordinates 1960-.5*12 and 1960-.5*55, and the y positions would be sin(pi/3)*12 and sin(pi/3)*55

cohort lines are shifted by the same amount but to the right, e.g. the 1960 cohort from age 12 to 55 would be: 1960+.5*12, 1960+.5*55 and the same y coordinates.

It may take fiddling to get right..., hence the labor into making them draw automatically.

EqLexis <- function(Rates, 
		value = "ASFR", 
		col, 
		breaks, 
		guides = TRUE, 
		guide.int = 5,
		legend = TRUE, 
		lab.breaks = NULL, 
		legend.args = NULL, 
		axis.args = NULL, 
		a.lab=NULL,
		c.lab=NULL,
		y.lab=NULL,
		a.lab.adj=c(0,0),
		c.lab.adj=c(0,0),
		y.lab.adj=c(0,0),
		a.ax.lab.adj=c(0,0),
		c.ax.lab.adj=c(0,0),
		y.ax.lab.adj=c(0,0),
		...){
	# try to coerce to numeric matrix
	if (is.data.frame(Rates)){
		Rates <- as.matrix(Rates[,colnames(Rates) %in% c("Age","Year","Cohort",value)])
	}
	
	# define polygon plot function:
	EQLexTriTiltLeft <- function(x,brks=breaks,cols=col,val=value){
		coli <- cols[(brks-x[val]) >= 0][1]
		ys <- sin(pi/3)
		# lower
		if (diff(c(x["Age"],x["Year"]))==x["Cohort"]){
			xcoord <- c(x["Year"]-.5*x["Age"],
					x["Year"]+1-.5*x["Age"],
					x["Year"]-.5*x["Age"]+.5)
			ycoord <- c(x["Age"]*ys,
					x["Age"]*ys,
					(x["Age"]+1)*ys)
		} else { #upper
			xcoord <- c(x["Year"]-.5*x["Age"],
					(x["Year"]+.5)-.5*x["Age"],
					(x["Year"]-.5)-.5*x["Age"])
			ycoord <- c(x["Age"]*ys,
					(x["Age"]+1)*ys,
					(x["Age"]+1)*ys)
		}
		polygon(xcoord,ycoord,col=coli,border="transparent")
	}	
	
	# default breaks and colors: 
	# if specifying different breaks or colors, remember, as with other
	# image()-like function, 
	if (missing(breaks)) {
		breaks <- approx(x=range(pretty(range(Rates[,value]))),n=201)$y
	}
	if (missing(col)) {
		colramp <- grDevices:::colorRampPalette(c("white","blue","green","yellow","orange"),space = "rgb")
		col <- colramp((length(breaks)-1))
	}
	# allow plotting in margins
	par(xpd=TRUE)
	
	# choose decent x and y ranges
	
	# ysc used to scale y everywhere = sin(pi/3), remember the unit circle!
	ysc <- sqrt(3)/2
	xlim <- c(min(Rates[,"Year"])-2-(diff(range(Rates[,"Age"]))*.5),max(Rates[,"Year"])+1)
	
	ages <- unique(Rates[,"Age"])
	ages <- c(ages,max(ages)+1)
	yl <- range(ages)
	ylim <- ysc*yl
	
	plot(NULL,type="n",ylim=ylim,xlim=xlim,asp=1,axes=FALSE,xlab="",ylab="", ...)
	
	# EQLexTriTiltLeft is the main polygon function
	apply(Rates,1,EQLexTriTiltLeft,cols=col,brks=breaks)
	
	# age references for lines
	agerefs <- sort(ages[ages%%guide.int==0])
	
	# year lines
	yrs <- c(unique(Rates[,"Year"]),max(Rates[,"Year"])+1)
	yrrefs <- sort(yrs[yrs %% guide.int == 0])
	x1 <- yrrefs - yl[1] * .5
	x2 <- yrrefs - yl[2] * .5
	if (guides){
		segments(x1 + .25, ylim[1] - ysc * .5, x2, ylim[2], col = "#00000030")
	} else {
		segments(x1 + .25, ylim[1] - ysc * .5, x1 - .5, ylim[1] + ysc, col = "black")
	}
	text(x1 + 2 + y.ax.lab.adj[1], ylim[1] - 1 + y.ax.lab.adj[2], yrrefs, pos = 1, srt = -60)
	
	# cohort lines to draw 
	cohs <- unique(Rates[,"Cohort"])
	cohrefs <- sort(cohs[cohs%%guide.int==0])
	# on left we want them to stop on the earliest year
	brange 	<- c(min(Rates[,"Year"]),min(Rates[,"Age"]),max(Rates[,"Year"])+1.5,max(Rates[,"Age"])+1.5)
	# shift x position left by .5 for each age
	x1 	<- pmax(cohrefs+brange[2],brange[1])
	x2	<- pmin(cohrefs+brange[4],brange[3])
	# rescale y positions using 'ysc'
	y1 	<- pmax(brange[2],x1-cohrefs)
	y2 	<- pmin(brange[4],x2-cohrefs)
	x1 	<- x1-.5*y1
	x2	<- x2-.5*y2
	y1 	<- y1*ysc
	y2 	<- y2*ysc
	if (guides){
		segments(x1,y1,x2,y2,col="#00000030")
	} else {
		segments(x2-1,y2-ysc*2,x2,y2,col="black")
	}
	text(x2 + 1 + c.ax.lab.adj[1], y2 + 1 + c.ax.lab.adj[2], cohrefs, srt = 60, pos = 3)
	
	# drawing age refs needed some year info
	if (guides) {
		segments(min(yrs)-agerefs*.5-.5,agerefs*ysc,max(xlim)-agerefs*.5,agerefs*ysc,col="#00000030")
	} else {
		segments(min(yrs)-agerefs*.5-.5,agerefs*ysc,min(yrs)-agerefs*.5+.5,agerefs*ysc,col="black")
	}
	
	text(min(yrs) - agerefs * .5 - .5 + a.ax.lab.adj[1], agerefs * ysc + a.ax.lab.adj[2], labels = agerefs, pos = 2)
	
	text(min(yrs) - 10 - agerefs[4] * .5 + a.lab.adj[1], agerefs[3] + a.lab.adj[2], labels = a.lab, cex = 1.5)
	text(max(yrrefs) - agerefs[2] * .5 + c.lab.adj[1], agerefs[length(agerefs) - 2] + c.lab.adj[2], labels = c.lab, srt = 60, cex = 1.5, pos = 3)
	text(median(yrrefs)+y.lab.adj[1],0+y.lab.adj[2],y.lab,cex=1.5,srt=-60)
	
	# legend:
	if (legend){
		fields:::image.plot(z = brks, lab.breaks = lab.breaks, col = cols, legend.args = legend.args, axis.args = axis.args, legend.only = TRUE)
	}
	# FINISH bounding box
	# give parallelogram box
	polygon(c(min(Rates[,"Year"])-.5*(max(Rates[,"Age"])+1),
					min(Rates[,"Year"])-.5*min(Rates[,"Age"]),
					max(Rates[,"Year"])+1-.5*min(Rates[,"Age"]),
					max(Rates[,"Year"])+1-.5*(max(Rates[,"Age"])+1)),
			ysc*c(max(Rates[,"Age"])+1,min(Rates[,"Age"]),min(Rates[,"Age"]),max(Rates[,"Age"])+1),
			
	)

}




# test data comes from the HFD, where "GBR_SCOasfrTR.txt" are age-period-cohort fertility rates from Scotland. you'll have to sign up at the HFD to get the data, but it's free. This code will work with the data taken as downloaded

Rates <- read.table("GBR_SCOasfrTR.txt", header = TRUE, skip = 2, na.strings = ".", as.is = TRUE)
# cut out open ages (you can leave them if you just make the column numeric)
Rates <- Rates[Rates$Age != "12-" & Rates$Age != "55+", ]
cnames <- colnames(Rates)
Years <- unique(Rates$Year)
Rates <- as.matrix(Rates)
Rates <- matrix(as.numeric(Rates),ncol=4)
colnames(Rates) <- cnames
brks <- approx(x=range(pretty(range(Rates[,"ASFR"]))),n=201)$y

# plotting code- will produce example image from blog post
dev.new(height=7,width=11)
EqLexis(Rates=Rates,
		main="Scotland Fertility, 1945 - 2009,\nCohorts 1890-1995, HMD",
		sub="Proposal for a non-distorted rendering of APC data")

# any questions, write me: tim.riffe@gmail.com