R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

Survival Functions

We are going to use different survival functions: - DeMoivre - Gompertz - Makeham

So first, lets create a vector cotaining values from 0 to 110, which will represent the age of the individuals.

age<-c(0:110)

Summary age

Just yo check the quantiles of age from 0 to 110:

summary(age)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0    27.5    55.0    55.0    82.5   110.0

Now we can create the survival functions. We will use the following formulas:

\[ \text{DeMoivre: } \; {}_t p_x = S_x(t) = \frac{\omega - x - t}{\omega - x} \]

\[ \text{Gompertz: } \; {}_t p_x = S_x(t) = e^{-\frac{B c^x (c^{t-1})}{\ln(c)}} \]

\[ \text{Makeham: } \; {}_t p_x = S_x(t) = e^{-A t - \frac{B c^x (c^{t-1})}{\ln(c)}} \]

DeMoivre:

Let’s start with DeMoivre survival function:

#``# Parameters
w<-110
# The function:
DeMoivre<-function(x,t){
  (w-x-t)/(w-x)
}
for (i in seq(from = 0, to = 90, by = 5)){
  surv<-DeMoivre(20,i)
  cat("The probability of surviving", i, "years from age 20 is", surv,"\n")}
## The probability of surviving 0 years from age 20 is 1 
## The probability of surviving 5 years from age 20 is 0.9444444 
## The probability of surviving 10 years from age 20 is 0.8888889 
## The probability of surviving 15 years from age 20 is 0.8333333 
## The probability of surviving 20 years from age 20 is 0.7777778 
## The probability of surviving 25 years from age 20 is 0.7222222 
## The probability of surviving 30 years from age 20 is 0.6666667 
## The probability of surviving 35 years from age 20 is 0.6111111 
## The probability of surviving 40 years from age 20 is 0.5555556 
## The probability of surviving 45 years from age 20 is 0.5 
## The probability of surviving 50 years from age 20 is 0.4444444 
## The probability of surviving 55 years from age 20 is 0.3888889 
## The probability of surviving 60 years from age 20 is 0.3333333 
## The probability of surviving 65 years from age 20 is 0.2777778 
## The probability of surviving 70 years from age 20 is 0.2222222 
## The probability of surviving 75 years from age 20 is 0.1666667 
## The probability of surviving 80 years from age 20 is 0.1111111 
## The probability of surviving 85 years from age 20 is 0.05555556 
## The probability of surviving 90 years from age 20 is 0

Using ggplot library we can use a better plot to visualize the survival probability since different ages and years:

result<-matrix(data=NA,nrow=length(age),ncol=w+1)
for (i in (seq_along(age)-1)){
    for (j in 0:(w-i)) {
      if (i<= w-j) {
        result[i+1,j+1]<-DeMoivre(i,j)
          } else {
        result[i+1,j+1]<-NA
        }
    }
}
result<-t(result)
result<-as.data.frame(result)
result$t<-(seq_along(0:110)-1)

library(ggplot2)
ggplot(result,aes(x=t))+
  geom_line(aes(y=V1, color="Age 0")) +
  geom_line(aes(y=V11, color="Age 10")) +
  geom_line(aes(y=V21, color="Age 20")) +
  geom_line(aes(y=V31, color="Age 30")) +
  geom_line(aes(y=V41, color="Age 40")) +
  geom_line(aes(y=V51, color="Age 50")) +
  geom_line(aes(y=V61, color="Age 60")) +
  geom_line(aes(y=V71, color="Age 70")) +
  geom_line(aes(y=V81, color="Age 80")) +
  geom_line(aes(y=V91, color="Age 90")) +
  labs(x = "Time (t)", y = "Survival Probability",
       title = "De Moivre's Law of Mortality",
       subtitle = "Probability of surviving t years") +
  scale_color_manual(values = rainbow(10)) +
  theme_minimal()
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 20 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 30 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 40 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 50 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 60 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 70 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 80 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 90 rows containing missing values or values outside the scale range
## (`geom_line()`).

Gompertz:

Now, counting with B=0.0000027 and c=1.124 for the survival probability, we’re going to estimate it th rough a plot to represent the probability since different ages and years:

At first we need to create the function:

# The function:
gompertz <- function(x, t, B, c) {
  exponent <- - (B * c^x * (c^t - 1)) / log(c)
  return(exp(exponent))
  }

To create the function we need to estimate specific probabilities:

for (i in seq(from = 0, to = 90, by = 5)){
  surv <- gompertz(20,i, 0.0000027, 1.124)
  cat("The probability of surviving", i, "years from age 20 is", surv,"\n")
}
## The probability of surviving 0 years from age 20 is 1 
## The probability of surviving 5 years from age 20 is 0.99981 
## The probability of surviving 10 years from age 20 is 0.9994693 
## The probability of surviving 15 years from age 20 is 0.9988583 
## The probability of surviving 20 years from age 20 is 0.9977631 
## The probability of surviving 25 years from age 20 is 0.9958012 
## The probability of surviving 30 years from age 20 is 0.9922913 
## The probability of surviving 35 years from age 20 is 0.9860252 
## The probability of surviving 40 years from age 20 is 0.9748827 
## The probability of surviving 45 years from age 20 is 0.9552072 
## The probability of surviving 50 years from age 20 is 0.9208987 
## The probability of surviving 55 years from age 20 is 0.8624068 
## The probability of surviving 60 years from age 20 is 0.7666245 
## The probability of surviving 65 years from age 20 is 0.6206611 
## The probability of surviving 70 years from age 20 is 0.4249039 
## The probability of surviving 75 years from age 20 is 0.2153065 
## The probability of surviving 80 years from age 20 is 0.06359116 
## The probability of surviving 85 years from age 20 is 0.007131258 
## The probability of surviving 90 years from age 20 is 0.0001407397

To visualize the results obtained in the last step:

result<-matrix(data=NA,nrow=length(age),ncol=w+1)
for (i in (seq_along(age)-1)){
    for (j in 0:(w-i)) {
      if (i<= w-j) {
        result[i+1,j+1]<-gompertz(i,j, 0.0000027, 1.124)
          } else {
        result[i+1,j+1]<-NA
        }
    }
}
result<-t(result)
result<-as.data.frame(result)
result$t<-(seq_along(0:110)-1)

library(ggplot2)
ggplot(result,aes(x=t))+
  geom_line(aes(y=V1, color="Age 0")) +
  geom_line(aes(y=V11, color="Age 10")) +
  geom_line(aes(y=V21, color="Age 20")) +
  geom_line(aes(y=V31, color="Age 30")) +
  geom_line(aes(y=V41, color="Age 40")) +
  geom_line(aes(y=V51, color="Age 50")) +
  geom_line(aes(y=V61, color="Age 60")) +
  geom_line(aes(y=V71, color="Age 70")) +
  geom_line(aes(y=V81, color="Age 80")) +
  geom_line(aes(y=V91, color="Age 90")) +
  labs(x = "Time (t)", y = "Survival Probability",
       title = "De Gompertz` Law of Mortality",
       subtitle = "Probability of surviving t years") +
  scale_color_manual(values = rainbow(10)) +
  theme_minimal()
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 20 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 30 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 40 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 50 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 60 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 70 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 80 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 90 rows containing missing values or values outside the scale range
## (`geom_line()`).

Makeham:

For our last step, we are going to analize the Makeham survival function!!!

Now, let’s apply the following values: A=0.00022, B=0.0000027 and c=1.124 for the survival probability function. We are going to estimate it through a plot to visually understand the relation between the probability of surviving with the passing of the years.

# The function:
makeham <- function(x, t, A, B, c) {
  exponent <- -(A*t)-(B*c^x*(c^t - 1))/log(c)
  return(exp(exponent))
}

We can create a function to estimate these specific probabilities:

for (i in seq(from = 0, to = 90, by = 5)){
  surv <- makeham(20,i, 0.00022, 0.0000027, 1.124)
  cat("The probability of surviving", i, "years from age 20 is", surv,"\n")
}
## The probability of surviving 0 years from age 20 is 1 
## The probability of surviving 5 years from age 20 is 0.9987108 
## The probability of surviving 10 years from age 20 is 0.9972729 
## The probability of surviving 15 years from age 20 is 0.9955675 
## The probability of surviving 20 years from age 20 is 0.9933826 
## The probability of surviving 25 years from age 20 is 0.9903394 
## The probability of surviving 30 years from age 20 is 0.9857637 
## The probability of surviving 35 years from age 20 is 0.978462 
## The probability of surviving 40 years from age 20 is 0.9663414 
## The probability of surviving 45 years from age 20 is 0.9457973 
## The probability of surviving 50 years from age 20 is 0.9108243 
## The probability of surviving 55 years from age 20 is 0.8520346 
## The probability of surviving 60 years from age 20 is 0.7565716 
## The probability of surviving 65 years from age 20 is 0.6118488 
## The probability of surviving 70 years from age 20 is 0.4184105 
## The probability of surviving 75 years from age 20 is 0.211783 
## The probability of surviving 80 years from age 20 is 0.06248174 
## The probability of surviving 85 years from age 20 is 0.006999142 
## The probability of surviving 90 years from age 20 is 0.0001379805

To visualize the results obtained in the last step:

result<-matrix(data=NA,nrow=length(age),ncol=w+1)
for (i in (seq_along(age)-1)){
    for (j in 0:(w-i)) {
      if (i<= w-j) {
        result[i+1,j+1]<-makeham(i,j, 0.00022, 0.0000027, 1.124)
          } else {
        result[i+1,j+1]<-NA
        }
    }
}
result<-t(result)
result<-as.data.frame(result)
result$t<-(seq_along(0:110)-1)

library(ggplot2)
ggplot(result,aes(x=t))+
  geom_line(aes(y=V1, color="Age 0")) +
  geom_line(aes(y=V11, color="Age 10")) +
  geom_line(aes(y=V21, color="Age 20")) +
  geom_line(aes(y=V31, color="Age 30")) +
  geom_line(aes(y=V41, color="Age 40")) +
  geom_line(aes(y=V51, color="Age 50")) +
  geom_line(aes(y=V61, color="Age 60")) +
  geom_line(aes(y=V71, color="Age 70")) +
  geom_line(aes(y=V81, color="Age 80")) +
  geom_line(aes(y=V91, color="Age 90")) +
  labs(x = "Time (t)", y = "Survival Probability",
       title = "De Makeham`s Law of Mortality",
       subtitle = "Probability of surviving t years") +
  scale_color_manual(values = rainbow(10)) +
  theme_minimal()
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 20 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 30 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 40 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 50 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 60 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 70 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 80 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 90 rows containing missing values or values outside the scale range
## (`geom_line()`).

Analysis:

x<- 0
A <- 0.00022
B <- 0.0000027
c <- 1.124

data <- data.frame(
  Tiempo = age,
  DeMoivre = sapply(age, function(t) DeMoivre(x, t)),
  Gompertz = sapply(age, function(t) gompertz(x, t, B, c)),
  Makeham = sapply(age, function(t) makeham(x, t, A, B, c))
)
data_long <- tidyr::pivot_longer(data, cols = -Tiempo, 
               names_to = "Model", 
               values_to = "Survival")


ggplot(data_long, aes(x = Tiempo, y = Survival, color = Model)) +
  geom_line(linewidth = 1.2) +
  labs(title = "Comparison of Survival Functions",
       subtitle = paste("DeMoivre (ω =", w, "), Gompertz (B =", B, ", c =", c, "), Makeham (A =", A, ", B =", B, ", c =", c, ")"),
       x = "Time (t)",
       y = "tP₀ - Survival Probability",
       color = "Model") +
  scale_color_manual(values = c("DeMoivre" = "turquoise3", 
                               "Gompertz" = "tomato3", 
                               "Makeham" = "khaki1")) +
  ylim(0, 1) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        plot.subtitle = element_text(hjust = 0.5),
        legend.position = "bottom")

Conclusion:

As we can observe, each graphic has its own values, even though some of them are very similar (such as Gompertz and Makeham). We can deduce that even if the formulas are not completely the same, they will give us back similar values since the possibilities to survive after n years is each time smaller or less than the last one. In this case this could be applied to us, students, since most of us are 20 years old. This graphics could be a great representation of how our probabilities of surviving after a certain period of time will be. In our last step, which we know wasn’t complete necessary since it wasn’t required, we wanted to do a graphic that included all three formulas so it would be easier to see how they behave in the same values for different ages with the same parameters for each one.

We know our last graphic is not the best one, but is the best we could do to include all three in just one, and as we already said, the Gompertz and Makeham are very similar, while the DeMoivre is more like a linear function. All these we can observe from each individual one: in the DeMoivre we can observe several lines as in a linear function, meanwhile the Gompertz and Makeham look very alike since the results obtained in each one are very similar (you could check each result in every section of the functions), they only differ a few decimals or even centecimals.

Now, applying our results to our case, the Gompertz and Makeham models show higher survival probabilities compared to DeMoivre for young adults. DeMoivre wouldn’t be adequate for our 20 year old person since it’s more pessimist for young people and too optimist for older ones. The other two are more realistic and consistent for the human behaviour.