Create a gist now

Instantly share code, notes, and snippets.

@thenoviceoof /longevity.R Secret
Last active Sep 30, 2017

What would you like to do?
Produces a graph comparing Gompertz models, one fit to Dutch data and another with a RR=0.5 intervention
## 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment