# Copyright (c) Longevitas Ltd.  All rights reserved.

# Model of left-truncation and right-censorship
frame()
par(mar=c(0,0,0,0))
plot.window(xlim=c(0, 9), ylim=c(0, 9))
#rect(0, 8, 1, 8.5, col="grey", border=NA)
rect(1, 8, 7, 8.5, col="lightgrey", border=NA)
lines(x=c(0, 9), y=c(8, 8), lwd=2, col="darkblue")
lines(x=c(1, 1), y=c(7.5, 8.5), lwd=2, col="darkblue")
lines(x=c(7, 7), y=c(7.5, 8.5), lwd=2, col="darkblue")
text(1, 7, expression(y[i]), cex=1.5)
text(7, 7, expression(y[i]+t[i]), cex=1.5)
text(7, 6, expression(d[i]==0), cex=1.5)
text(8.5, 8, "X", cex=2, col="red")

#rect(0, 3, 1, 3.5, col="grey", border=NA)
rect(1, 3, 5, 3.5, col="lightgrey", border=NA)
lines(x=c(0, 9), y=c(3, 3), lwd=2, col="darkblue")
lines(x=c(1, 1), y=c(2.5, 3.5), lwd=2, col="darkblue")
text(5, 3, "X", cex=2, col="red")
text(1, 2, expression(y[i]), cex=1.5)
text(5, 2, expression(y[i]+t[i]), cex=1.5)
text(5, 1, expression(d[i]==1), cex=1.5)

# R source for drawing the five mortality laws in laws.pdf
ages<-seq(50,120)
alpha<--12.0
beta<-0.09
epsilon<--10.0
rho<--1.0

# Now calculate qx according to the five laws
Cox<-rep(exp(epsilon), length(ages))
Gompertz<-exp(alpha+beta*ages)
Makeham<-exp(epsilon)+Gompertz
Perks<-Gompertz/(1.0+Gompertz)
Beard<-Gompertz/(1+exp(rho+alpha+beta*ages))

par(mar=c(4.5, 4, 0, 1))
par(mfrow=c(1,2), oma=c(0,0,0,0))
plot(x=ages, y=Gompertz, type="n", xlab="Age", ylab="Force of mortality", font.lab=2, frame=FALSE, axes=FALSE)
lines(x=ages, y=Cox, lwd=2, lty=1)
lines(x=ages, y=Gompertz, lwd=2, lty=2)
lines(x=ages, y=Makeham, lwd=2, lty=3)
lines(x=ages, y=Perks, lwd=2, lty=4)
lines(x=ages, y=Beard, lwd=2, lty=5)
axis(1, font=2)
axis(2, las=1, font=2)

names<-c(expression(paste("Gompertz (1825) ", q[x] == exp(alpha + beta * x))), expression(paste("Perks (1932) ", q[x] == frac(exp(alpha + beta * x), 1+exp(alpha + beta * x)))), expression(paste("Extreme value:  ", q[x] == 1-exp(-exp(alpha + beta * x)))), expression(paste("Probit-Gompertz: ", q[x] == Phi(alpha + beta * x))), "Reversed extreme value: ", expression(paste(q[x] == exp(frac(-1, exp(alpha + beta * x))))))
lines<-c(1, 2, 3, 4, 5, NA)
legend(60, 0.65, names, lty=lines, lwd=2, bty="n", cex=0.85, y.intersp=0)

# Now do same on log scale
plot(x=ages, y=log(Gompertz), type="n", xlab="Age", ylab="log(mortality rate)", font.lab=2, frame=FALSE, axes=FALSE)
lines(x=ages, y=log(Cox), lwd=2, lty=1)
lines(x=ages, y=log(Gompertz), lwd=2, lty=2)
lines(x=ages, y=log(Makeham), lwd=2, lty=3)
lines(x=ages, y=log(Perks), lwd=2, lty=4)
lines(x=ages, y=log(Beard), lwd=2, lty=5)
axis(1, font=2)
axis(2, las=1, font=2)

# Bootstrap plots
d<-read.csv(file="c:\\research\\baj\\laws\\Bootstrapping weighted models.csv", header=TRUE)
par(mar=c(4, 5, 0, 0))
plot(d$ActualLives, d$ExpectedLives, pch=20, frame=FALSE, axes=FALSE, xlab="Actual deaths", ylab="Expected deaths")
axis(1, font=2, las=1)
axis(2, font=2, las=1)

# Lives-weighted analysis
d<-read.csv(file="c:\\research\\baj\\laws\\Lives- v. Amounts weighting.csv", header=TRUE)
plotchar<-20
par(mar=c(0, 3, 2, 0), mfrow=c(1, 2))
plot(d$LivesWeightedAmounts, pch=plotchar, frame=FALSE, axes=FALSE, xlab="", ylab="Ratio of expected v. actual", ylim=range(0.75, 1.6), col="grey", type="n")
lines(d$LivesWeightedLives, col="grey")
points(d$LivesWeightedAmounts, pch=plotchar, col="black", cex=0.6)
axis(2, font=2, las=1)
title("(a) Lives-weighted model fit")
plot(d$AmountsWeightedAmounts, pch=plotchar, frame=FALSE, axes=FALSE, xlab="", ylab="", ylim=range(0.75, 1.6), col="grey", type="n")
lines(d$AmountsWeightedLives, col="grey")
points(d$AmountsWeightedAmounts, pch=plotchar, col="black", cex=0.6)
axis(2, font=2, las=1)
title("(b) Amounts-weighted model fit")
legend(1000, 1.55, pch=c(NA, plotchar), col=c("grey", "black"), lty=c(1, NA), legend=c("Ratio by lives", "Ratio by pension amount"), bty="n")

# Plots of models in Table 6
ages<-seq(60, 100)
rAge<-data.frame(X=c(60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100), Events=c(358,449,440,504,543,979,1042,1079,1208,1218,1343,1414,1509,1644,1730,1879,1984,2019,2067,2207,2176,2211,2119,1938,1903,1734,1605,1390,1213,1055,996,867,715,552,438,369,254,180,122,82,1), Exposures=c(73831.376731,77760.254447,75742.920653,71966.273796,69138.266094,104605.282005,101322.694692,95590.648040,89857.985503,84005.496840,79461.079133,74729.316472,70054.301688,66010.102860,62196.092731,59268.633746,54778.716778,50950.725570,46980.816431,43498.219980,40037.247638,34744.596137,29725.450318,25294.938364,21332.615226,17692.047556,14150.656323,11101.899638,9365.200914,7823.414084,6404.780607,4996.599332,3861.602540,2950.861435,2241.410447,1677.122034,1249.148130,883.223322,643.826907,469.680919,0.813748), Expected=c(425.200131,499.274624,542.870874,576.356510,618.411421,1045.791044,1130.491425,1191.275804,1250.845062,1306.053779,1379.800053,1449.631281,1517.709751,1597.622730,1681.092467,1789.136358,1847.145199,1919.037226,1976.268493,2044.179645,2100.448968,2035.515914,1944.862266,1848.478322,1741.205142,1612.577286,1439.888080,1262.547837,1189.763370,1110.130863,1014.601451,883.910836,763.040919,651.216948,552.394462,461.638727,383.969616,303.168890,246.947120,201.129788,0.369319), DevianceResidual=c(-3.350963,-2.289421,-4.566902,-3.080532,-3.097476,-2.087951,-2.667396,-3.306186,-1.218448,-2.464683,-0.995149,-0.939716,-0.223783,1.154747,1.187119,2.107104,3.146113,2.262507,2.025638,3.554929,1.638745,3.835590,3.891818,2.065717,3.819556,2.986904,4.271829,3.528989,0.671488,-1.668643,-0.585780,-0.570630,-1.757897,-3.993550,-5.051617,-4.469431,-7.072398,-7.658840,-8.822520,-9.545348,0.854885))
muGompertz<-exp(-11.8512 + 0.110611 * ages)
muMakeham<-exp(-31.7401) + exp(-11.8512 + 0.110611 * ages)
muPerks<-exp(-12.4393 + 0.119024 * ages) / (1.0 + exp(-12.4393 + 0.119024 * ages))
muBeard<-exp(-13.4434 + 0.133545 * ages) / (1.0 + exp(-13.4434 + 1.00715 + 0.133545 * ages))
muMakehamPerks<-(exp(-30.0494) + exp(-12.4393 + 0.119024 * ages)) / (1.0 + exp(-12.4393 + 0.119024 * ages))
muMakehamBeard<-(exp(-6.7183) + exp(-14.2129 + 0.143227 * ages)) / (1.0 + exp(-14.2129 + 1.1407 + 0.143227 * ages))
par(mar=c(3, 4, 0, 0), mfrow=c(1, 2))
plot(x=ages, y=log(muGompertz), ylim=range(-5.5, -1), type="n", xlab="Age", ylab="log(force of mortality)", font.lab=2, frame=FALSE, axes=FALSE)
points(x=rAge$X, y=log(rAge$Events / rAge$Exposure), pch=19, lwd=2, col="red")
lines(x=ages, y=log(muGompertz), lwd=2, lty=1, col="blue")
axis(1, font=2)
axis(2, las=1, font=2)
legend(60, -1, pch=c(NA, 19), lwd=c(2, 2), col=c("blue", "red"), lty=c(1, NA), legend=c("Fitted Gompertz law", "Crude observed force of mortality"), bty="n")

plot(x=ages, y=log(muGompertz), ylim=range(-5.5, -1), type="n", xlab="Age", ylab="log(force of mortality)", font.lab=2, frame=FALSE, axes=FALSE)
points(x=rAge$X, y=log(rAge$Events / rAge$Exposure), pch=19, lwd=2, col="red")
lines(x=ages, y=log(muMakehamBeard), lwd=2, lty=1, col="blue")
axis(1, font=2)
axis(2, las=1, font=2)
legend(60, -1, pch=c(NA, 19), lwd=c(2, 2), col=c("blue", "red"), lty=c(1, NA), legend=c("Fitted Makeham-Beard law", "Crude observed force of mortality"), bty="n")

# Draw the reliability theory illustration
frame()
par(mar=c(0,0,0,0))
plot.window(xlim=c(0, 15), ylim=c(0, 8))
drawelement<-function(x, y)
{
  lines(x=c(x, x+1), y=c(y, y), lwd=2, col="darkblue")
  rect(x+1, y-0.5, x+2, y+0.5, lwd=2)
  lines(x=c(x+2, x+3), y=c(y, y), lwd=2, col="darkblue")
}
drawblock<-function(x, y)
{
  drawelement(x+1, y+3.0)
  drawelement(x+1, y+1.0)
  drawelement(x+1, y-1.0)
  drawelement(x+1, y-3.0)
  lines(x=c(x, x+1), y=c(y, y), lwd=2, col="darkblue")
  lines(x=c(x+4, x+5), y=c(y, y), lwd=2, col="darkblue")
  lines(x=c(x+1, x+1), y=c(y-3.0, y+3.0), lwd=2, col="darkblue")
  lines(x=c(x+4, x+4), y=c(y-3.0, y+3.0), lwd=2, col="darkblue")
}
drawblock(0, 4)
drawblock(5, 4)
drawblock(10, 4)
text(2.5, 3, "X", cex=2, col="red")
text(2.5, 5, "X", cex=2, col="red")
text(2.5, 7, "X", cex=2, col="red")
text(7.5, 5, "X", cex=2, col="red")
text(12.5, 1, "X", cex=2, col="red")
text(12.5, 3, "X", cex=2, col="red")
text(12.5, 5, "X", cex=2, col="red")
text(12.5, 7, "X", cex=2, col="red")
text(2.5, 0, "Block A")
text(7.5, 0, "Block B")
text(12.5, 0, "Block C")

# Draw the reliability plot
d<-read.csv(file="c:\\research\\baj\\laws\\reliability.csv", header=TRUE)
par(mar=c(3,4,0,0))
muGompertz<-log(d$Mu[31])+(d$Age-50.0)*(log(d$Mu[51]) - log(d$Mu[31]))/(70.0 - 50.0)
plot(x=d$Age, y=log(d$Mu), ylim=range(-5.5, -1), xlim=range(50, 100), type="n", xlab="Age", ylab="log(force of mortality)", font.lab=2, frame=FALSE, axes=FALSE)
lines(x=d$Age, y=muGompertz, lwd=2, lty=1, col="grey")
lines(x=d$Age, y=log(d$Mu), lwd=2, lty=1, col="darkblue")
axis(1, font=2)
axis(2, las=1, font=2)
q<-read.csv(file="c:\\Consulting\\survival\\PCA00Base.csv", header=TRUE)
mu<--log(1-q$Males)
points(x=q$Age+0.5, y=log(mu), pch=20, lwd=2, col="red")
legend(50, -1, lwd=c(2, 2, 2), lty=c(1, 1, NA), pch=c(NA, NA, 20), col=c("darkblue", "grey", "red"), legend=c("Force of mortality according to reliability model", "Gompertz extrapolation of force of mortality between 50 and 70", "Force of mortality under PCA00Base (males)"), bty="n")

# Plot hazard functions for elements and blocks
ages<-seq(0, 100, 0.25)
lambda<-0.15
block<-function(lambda, n, ages)
{
  lambda * n * exp(-lambda * ages) * (1-exp(-lambda * ages))^(n-1) / (1- (1 - exp(-lambda * ages)) ^ n)
}
block(lambda, 10, ages)
par(mar=c(4.5,4,0,0))
plot(x=ages, y=block(lambda, 10, ages), type="n", xlab="Age", ylab="hazard", font.lab=2, frame=FALSE, axes=FALSE)
lines(x=ages, y=rep(lambda, length(ages)), lwd=2, lty=1, col="grey")
lines(x=ages, y=block(lambda, 100, ages), lwd=2, lty=3, col="black")
lines(x=ages, y=block(lambda, 1000, ages), lwd=2, lty=1, col="darkblue")
axis(1, font=2)
axis(2, las=1, font=2)
names<-c(expression(paste("Hazard for element, ", lambda==0.15)), expression(paste("Hazard for block, ", lambda==0.15, " , ", n==100)), expression(paste("Hazard for block, ", lambda==0.15, " , ", n==1000)))
legend(60, 0.1, lwd=c(2, 2, 2), lty=c(1, 3, 1), col=c("grey", "black", "darkblue"), legend=names, bty="n")

# Copyright (c) Longevitas Ltd.  All rights reserved.