Last active
May 1, 2016 18:54
-
-
Save timriffe/51f360fdce29314a7574cf9adcd4668d to your computer and use it in GitHub Desktop.
dunno what to call it. It's a lexis surface, but using prospective age instead of age, and age itself is the z value. Rather convoluted.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Author: tim | |
# this code appears in a blog post here: | |
# https://sites.google.com/site/timriffepersonal/DemogBlog/aprospectiveagelexissurfaceconvolutedspielen | |
############################################################################### | |
# contains 1) a prospective age calculator | |
# 2) a quantile age calculator | |
# 3) a prospective age surface, but with age swapped for e(x). concoluted. | |
# note: If you don't know what a prospective age is, check out Sanderson & Scherbov (2007): | |
# http://www.demographic-research.org/volumes/vol16/2/ | |
library(HMDHFDplus) | |
# define username and password in console | |
ltf <- readHMDweb("USA","fltper_1x1",username=us,password=pw) | |
# a prospective age calculator, given e(x) | |
ex2age <- function(ex,Age = (1:length(ex)-1), ex.at = seq(50,5,by=-5)){ | |
# splines make it easier to interpolate to exact values. | |
out <- splinefun(Age~ex)(ex.at) | |
names(out) <- ex.at | |
out | |
} | |
# a quantile age calculator, given l(x) | |
lxquantile2age <- function(lx, Age = (1:length(lx)-1), p = seq(.9,.1,by=-.1), closeout = 2){ | |
# includes closing out the lifetable at open age + 2, whatever. Only OK | |
# if your lifetable goes to 110+ (could also assume exponential or something) | |
# same use of spline, but we force monotonic. | |
out <- splinefun(c(Age,max(Age)+closeout)~c(lx,0),method="monoH.FC")(p) | |
names(out) <- p | |
out | |
} | |
Age <- 0:110 | |
yr1 <- 1950 ; yr2 <- 2010 | |
# a look at selected period prospective ages | |
ex1 <- ltf$ex[ltf$Year == yr1] | |
ex2 <- ltf$ex[ltf$Year == yr2] | |
prosp <- cbind(ex = seq(50,5,by=-5), | |
age1 = round(ex2age(ex1),1), | |
age2 = round(ex2age(ex2),1)) | |
head(prosp) | |
# ex age1 age2 | |
#50 50 24.1 32.3 | |
#45 45 29.4 37.5 | |
#40 40 34.8 42.8 | |
#35 35 40.3 48.2 | |
#30 30 45.9 53.8 | |
#25 25 51.9 59.6 | |
# or a look at period quantile ages | |
lx1 <- ltf$lx[ltf$Year == yr1] / 1e5 | |
lx2 <- ltf$lx[ltf$Year == yr2] / 1e5 | |
quant <- cbind(prob = seq(.9,.1,by=-.1), | |
age1 = round(lxquantile2age(lx1),1), | |
age2 = round(lxquantile2age(lx2),1)) | |
head(quant) | |
# prob age1 age2 | |
#0.9 0.9 48.1 62.5 | |
#0.8 0.8 60.8 72.1 | |
#0.7 0.7 67.5 77.6 | |
#0.6 0.6 72.2 81.6 | |
#0.5 0.5 75.8 84.7 | |
#0.4 0.4 79.0 87.5 | |
# you didn't think I'd not try to make a surface out | |
# of this stuff, did you? | |
# [period assumption sinner here] | |
library(reshape2) | |
# an age-period matrix of e(x) | |
EX <- acast(ltf,Age~Year,value.var = "ex") | |
# a big matrix of prospective ages, sort of: | |
# we're using thresholds rather than comparing | |
# with a standard, but in the end such comparisons | |
# will be possible. | |
Sand <- apply(EX,2,ex2age,ex.at=5:50) | |
# define a color ramp, because colors | |
colRamp <- colorRampPalette(RColorBrewer::brewer.pal(9,"Blues"),space="Lab") | |
# a custom cohort line function, | |
# not necessarily robust or user-friendly. | |
# made for the sake of this single plot... | |
LexLines <- function(EX,N=5,...){ | |
Sand <- apply(EX,2,ex2age,ex.at=5:50) | |
# make it hmmm. | |
ages <- as.integer(rownames(EX)) | |
aN <- ages[ages %% N == 0] | |
aN <- aN[aN < max(Sand)] | |
aN <- aN[aN > min(Sand)] | |
years <- as.integer(colnames(EX)) | |
yN <- years[years %% N == 0] | |
EXN <- EX[as.character(aN),as.character(yN)] | |
for (m in 1:(nrow(EXN)-1)){ | |
for (n in 1:(ncol(EXN)-1)){ | |
segments(yN[n],EXN[m,n],yN[n+1],EXN[m+1,n+1], xpd=FALSE, ...) | |
} | |
} | |
# left side: | |
offl <- min(yN) - min(years) | |
if (offl > 0){ | |
EXleft <- EX[as.character(aN+offl),1] | |
for (m in 1:(nrow(EXN)-1)){ | |
segments(min(years),EXleft[m],yN[1],EXN[m+1,1], xpd=FALSE, ...) | |
} | |
} | |
# right side | |
offr <- max(years) - max(yN) | |
if(offr > 0){ | |
EXright <- EX[as.character(aN-offl),ncol(EX)] | |
for (m in 1:(nrow(EXN)-1)){ | |
segments(max(yN),EXN[m,ncol(EXN)], max(years), EXright[m+1], xpd=FALSE, ...) | |
} | |
} | |
cohL <- min(yN) - aN | |
for (i in 1:(length(cohL)-1)){ | |
rise <- EXN[i+1,1] - EXN[i,1] | |
run <- N | |
slope <- rise / run | |
angle <- atan(slope) * 180 / pi | |
text(min(yN), EXN[i,1]+1, cohL[i], col = "yellow", cex = .8, srt = angle) | |
} | |
cohR <- max(yN) - aN | |
for (i in 1:(length(cohR)-1)){ | |
rise <- EXN[i+1,ncol(EXN)] - EXN[i,ncol(EXN)] | |
run <- N | |
slope <- rise / run | |
angle <- atan(slope) * 180 / pi | |
text(max(yN), EXN[i,ncol(EXN)]+1, cohR[i], col = "yellow", cex = .8, srt = angle) | |
} | |
} | |
# the complicated plot, will describe in blog. | |
filled.contour(as.integer(colnames(Sand)),as.integer(rownames(Sand)),t(Sand), | |
asp=1,color=colRamp, xlab = "Year", | |
ylab = "e(x)",key.title = title("Age"), | |
# but quick tip: since filled.contour() spits back a device with crazy coords, | |
# if you want to plot more in the figure, you need to do it by passing arguments | |
# to plot.axes (non-intuitively). It's just like panel.first = list(bla bla) | |
plot.axes = { LexLines(EX,N=5,col="#FFFF0080"); | |
contour(as.integer(colnames(Sand)),as.integer(rownames(Sand)),t(Sand), | |
drawlabels = TRUE, axes = FALSE, | |
frame.plot = FALSE, add = TRUE, | |
labcex = 1); | |
axis(1); axis(2); | |
# grids are helpful in this case to match prospective ages... | |
abline(h=seq(5,50,by=5),col="#FFFFFF30"); | |
abline(v=seq(1935,2015,by=5),col="#FFFFFF30")}) | |
# one could do this for quantiles as well, maybe some other time. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment