Produces a graph comparing Gompertz models, one fit to Dutch data and another with a RR=0.5 intervention

longevity.R

## Licensed under CC-0, "No Rights Reserved"
## (https://creativecommons.org/publicdomain/zero/1.0/legalcode).
## Following on Gwern's licensing.

library(ggplot2)
library(reshape2)

## Years to calculate.
y = 0:120

## The next 2 lines are from Gwern's longevity page
## (https://www.gwern.net/Longevity) (licensed under CC-0
## https://www.gwern.net/About#license),
## with parameters fit on Dutch population data from
## "Mortality hazard rates and life expectancy", Cramer and Kaas (2013).
S <- function(t, RR=1) { exp(-((RR*0.000016443)/log(1.1124) * (1.1124^t - 1))) }
df = data.frame(year=y, normal=sapply(y, S), half=sapply(y, S, 0.5))

## Transform the frame to include columns names as values.
ndf = melt(df, id.vars=c("year"))

plot = ggplot(ndf, aes(x=year, y=value, color=variable)) +
    geom_point(shape=1) +
    ggtitle("Impact of RR on Gompertz-model longevity") +
    xlab("Years") +
    ylab("Survival curve")

ggsave(filename="plot.png", width=7, height=4.5, plot)