LifeTable Sailing

LifeTableSailing.R

# R code for 4th of July jet lag blog post
# Tim Riffe
############################################

# library(devtools)
# install_github("timriffe/TR1/TR1/HMDHFDplus")
library(HMDHFDplus)

# it'll ask you to give your HMD password and username in the console (no quotes)
Mx <- readHMDweb("USA", "Mx_1x1")

mx2lx <- function(mx){
	# (lazy-man's approximation)
	# we're not sacrificing too much precision here. just pretend...
	c(1, exp(-cumsum(mx)))[1:length(mx)]
}
lx2med <- function(lx){
	# get median by using monotonic spline over l(x)
	a <- 0:length(lx)
	splinefun(a ~ c(lx, 0), method = "hyman")(.5)
}
mx2med <- function(mx){
	lx2med(mx2lx(mx))
}
#e(x)^\star = mean life assuming constant hazard implied by half life
lx2exstar <- function(lx){
	lx2med(lx) / log(2)
}
lx2mstar <- function(lx){
	log(2) / lx2med(lx)
}
mx2mstar <- function(mx){
	lx2mstar(mx2lx(mx))	
}

mx2mstarx <- function(mx){
	N   <- length(mx)
	out <- rep(NA, N)
	for (i in 1:N){
		out[i] <- mx2mstar(mx[i:N])
	}
	out
}
mx2medx <- function(mx){
	N   <- length(mx)
	out <- rep(NA, N)
	for (i in 1:N){
		out[i] <- mx2med(mx[i:N])
	}
	out
}
# now for plots

mx <- Mx$Male[Mx$Year == 1980]
a  <- 0:110

png("mx1980males.png")
plot(a,mx,type='l',log='y',xlab = "age", ylab = "rate (log scale)")
dev.off()

png("mx1980malesSeg.png")
plot(a, mx, type = 'l', log = 'y', xlab = "age", ylab = "rate (log scale)")
segments(a, mx2mstarx(mx), a + mx2medx(mx), mx2mstarx(mx))
dev.off()

png("mxstar1980males.png")
plot(a, mx, type = 'l', log = 'y', xlab = "age", ylab = "rate (log scale)")
lines(a, mx2mstarx(mx), col = "blue")
legend("bottomright", lty = 1, col = c("black", "blue"), legend = c("m(x)", "m*(x)"))
dev.off()

png("sail1980.png")
plot(a, mx, type = 'n', xlab = "", ylab = "", log = 'y')
polygon(c(a, rev(a + mx2medx(mx))),
		c(mx2mstarx(mx), rev(mx2mstarx(mx))),
		col = "#CCCCAA")
dev.off()

#install.packages("animation")
library(animation)
saveGIF(for (yr in sort(unique(Mx$Year))){
			plot(a, mx, type = 'n', xlab = "", ylab = "", log = 'y', ylim = c(5e-3, 1.5), axes = FALSE)
			mx <- Mx$Male[Mx$Year == yr]
			polygon(c(a, rev(a + mx2medx(mx))), 
					c(mx2mstarx(mx), rev(mx2mstarx(mx))),
					col = "#CCCCAA50")
			text(15, 1.2, yr, cex = 4)
		}, movie.name = "USmalesSailing.gif", img.name = "Rplot", convert = "convert", 
		cmd.fun = system, clean = TRUE)


#saveHTML(for (yr in sort(unique(Mx$Year))){
#			plot(a, mx, type = 'n',xlab = "", ylab = "", log = 'y', ylim = c(5e-3, 1.5), axes = FALSE)
#			mx <- Mx$Male[Mx$Year == yr]
#			polygon(c(a, rev(a + mx2medx(mx))),
#					c(mx2mstarx(mx), rev(mx2mstarx(mx))),
#					col = "#CCCCAA50")
#			text(15, 1.2, yr, cex = 4)
#		}, img.name = "norm_plot", single.opts = "utf8: false", autoplay = FALSE, 
#		interval = 0.05, imgdir = "norm_dir", htmlfile = "USmalesSailing.html", 
#		ani.height = 500,verbose=FALSE, ani.width = 500, title = "US men's lifetable sailing",navigator=FALSE)
#