calculate Schoen's del for the entire HMD

Schoens_del.R

# Tim Riffe
# HT Andrew Noymer
# Link to original Schoen article (Demography, 1970):
# http://link.springer.com/article/10.2307/2060150


#library(devtools)
#install_github("timriffe/TR1/TR1/HMDHFDplus")
library(HMDHFDplus)
library(data.table)

# a nify del function:
Schoen_del <- function(mx1,mx2=rep(1,length(mx1))){
	stopifnot(length(mx1)==length(mx2))
	mx1[mx1 == 0] <- NA
	mx2[mx2 == 0] <- NA
	prod(mx1 / mx2, na.rm = TRUE) ^ (1 / length(mx1))
}

# all HMD country codes:
XYZ <- getHMDcountries()

# you need to define two character strings, us and pw as your
# HMD username and password. If you use these functions, just run 
# it, enter the username/password in the console (no quotes) and press enter
# and the value will be assigned to the variable. This way it doesn't get
# saved in your script ....


# us <- userInput()
# pw <- userInput()

# now grab all the lifetables and mesh together..
# grab them all
HMDL <- do.call(rbind,lapply(XYZ, function(x, us, pw){
		    cat(x,"\n")
			Males        <- readHMDweb(x,"mltper_1x1",username=us,password=pw)
			Females      <- readHMDweb(x,"fltper_1x1",username=us,password=pw)
			Males$Sex    <- "m"
			Females$Sex  <- "f"
            CTRY         <- rbind(Females, Males)
			CTRY$PopName <- x
			CTRY		
		}, us = us, pw = pw))

# make it a data.table for slick operations
HMDL <- data.table(HMDL)

# tapply() would also do the trick, but doesn't come out so nicely organized.
del     <- HMDL[, Schoen_del(mx), by = list(PopName, Year, Sex)]
del     <- as.data.frame(del)

# take care of the Belgium missing years...
del$V1[del$V1 == 1] <- NA

# plot everything
png("HMDschoensdel.png", width=800,height=800)
plot(NULL, 
		type = 'n', 
		xlim = range(del$Year), 
		ylim = c(10,max(1/del$V1, na.rm=TRUE)), 
        log = 'y',
		ylab = "del")
for (ctry in XYZ){
	lines(del$Year[del$PopName == ctry & del$Sex == "f"], 
			1/del$V1[del$PopName == ctry & del$Sex == "f"], 
			col = "#FF000050")
	lines(del$Year[del$PopName == ctry & del$Sex == "m"], 
			1/del$V1[del$PopName == ctry & del$Sex == "m"], 
			col = "#0000FF50")
}
legend("topleft",lty=1,col=c("#FF000050","#0000FF50"),legend=c("females","males"))
dev.off()