Survival

R Markdown

This is an R Markdown document. The idea of this document is to serve as an introduction to R, but also to show how to plot some survival functions.

The firs thing we will do is to download and install R and RStudio (Now Posit). You can find the download links at R Project and Posit.

Once you have installed R and RStudio, we can start working. You can open RStudio and create a new R Markdown document by going to File -> New File -> R Markdown...

You should do this, so it will be easiear to share your results. You can also use the R console, but R Markdown is more powerful and allows you to create documents that include both code and text.

Now, let’s start with some basic R code.

Changing working directory: To do this, you have to specify R the location on your computer of the documents or the file where you will save all your work.

To do this you go to Session -> Set Working Directory -> Choose Directory... and select the folder where you want to save your work.Now we can star working.

setwd("D:/Curso 2025/R")

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)

Just to check what we input

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:

• DeMoivre: \(_tp_x=S_x(t) = \frac{\omega -x-t}{\omega -x}\)
• Gompertz: \(_tp_x=S_x(t) = e^{-\frac{Bc^x(c^t-1)}{\ln(c)}}\)
• Makeham: \(_tpx=S_x(t) = e^{-At-\frac{Bc^x(c^t-1)}{\ln(c)}}\)

Now we can create the survival functions in R.Lets start with DeMoivre.

#``# Parameters
omega<-110
# The function:
DeMoivre<-function(x,t){
  (omega-x-t)/(omega-x)
}

Now, we can give any value of x and t and we will find the probability of surviving. Lets estimate the probability of a person age 20 surviving 0,5,10,20 years.

DeMoivre(20,0)

## [1] 1

DeMoivre(20,5)

## [1] 0.9444444

DeMoivre(20,10)

## [1] 0.8888889

DeMoivre(20,20)

## [1] 0.7777778

We can create a function to estimate specific probabilities:

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

We can also have a plot, lets plot the probability of surviving one year for any age in our vector of ages.

result<-numeric(length(age))
for (i in seq_along(age)){
  result[i]<-DeMoivre(i,1)
}

plot(age, result, type="l", col="blue", lwd=2,
     xlab="Age", ylab="Survival Probability",
     main="De Moivre's Law of Mortality",
     sub="Probability of surviving 1 year")

And we can change the age, now lets change our function to surviving 5 years for any age.

result<-numeric(length(age))
for (i in seq_along(age)){
  result[i]<-DeMoivre(i,5)
}

plot(age, result, type="l", col="blue", lwd=2,
     xlab="Age", ylab="Survival Probability",
     main="De Moivre's Law of Mortality",
     sub="Probability of surviving 5 years")

As you can see we have a problem… negative probabilities? We have to add a restriction.

result<-numeric(length(age))
for (i in seq_along(age)){
  if (i<=105){
    result[i]<-DeMoivre(i,5)
  } else {
    result[i]<-NA
  }
}

plot(age, result, type="l", col="blue", lwd=2,
     xlab="Age", ylab="Survival Probability",
     main="De Moivre's Law of Mortality",
     sub="Probability of surviving 5 years")

Now lets plot for different ages and years. We will also use ggplot library to get a nicer plot.

result<-matrix(data=NA,nrow=length(age),ncol=omega+1)
for (i in (seq_along(age)-1)){
    for (j in 0:(omega-i)) {
      if (i<= omega-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)

## Warning: package 'ggplot2' was built under R version 4.5.1

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()`).

ACTIVITY:

Plot Gompertz and Makeham functions using ggplot and compare the results with DeMoivre. You just need the final plot (the one in colors). Your activity should include DeMoivre, Gompertz and Makeham function plots; however, grading will be based only on Gompertz and Makeham.

Use the following values: \(A=0.00022, B=0.0000027, c=1.124\))

For this activity you can work in pairs of teams of 3. All of you have to submit the activity individually so it will be reflected on the system.

Grading: - 50% Correct code - 30% Nice plots - 20% Discussion