# This code available at:
# http://pauldickman.com/software/stnet/relsurv.R
# 
# R code for calculating the Pohar-Perme estimate of marginal relative survival 
# using the relsurv package applied to the colon data widely used in 
# Stata tutorials on my website.
# 
# The R code was written by Enoch Yi-Tung Chen, who has written a tutorial 
# (https://enochytchen.com/tutorials/relsurv/tied/) comparing the various 
# software implementations.
# 
# See http://pauldickman.com/software/stnet/compare_stns/ for comparable estimates
# in four Stata commands.
# 
# Enoch Yi-Tung Chen, Paul Dickman
# August 2022

library(haven)
library(tidyverse)
library(relsurv)

# Read the patient data (colon cancer)
colon <- read_dta("http://pauldickman.com/data/colon.dta")

# Restrict to localised
colon <- subset(colon, stage==1)

colon2 <- colon %>%
          mutate(sex     = as_factor(sex),      
                 status  = as.numeric(status),
                 stime   = as.numeric(exit-dx),  # survival time in days
                 )%>%
          select(id, sex, status, dx, age, stime)

str(colon2)
summary(colon2)

# ===============================================================================
# Convert Stata popmort file into an R ratetable for use in rs.surv()
popmort <- read_dta("http://pauldickman.com/data/popmort.dta")
colnames(popmort)[2] <- "year"
colnames(popmort)[3] <- "age"

popmort <- popmort %>% 
           mutate (sex  = as_factor(sex),
                   prob = as.numeric(prob),
                  )%>%
           select(sex, age, year, prob)

# Make the dataset wider instead of longer
popmort <- popmort %>% spread(year, prob)

# Split the popmort file by sex (1 = m; 2 = f)
popmort_m <- subset(popmort, sex == 1);
popmort_m <-as.matrix(popmort_m[-c(1:2)])
popmort_f <- subset(popmort, sex == 2);
popmort_f <-as.matrix(popmort_f[-c(1:2)])

popmort_rt<- transrate(popmort_m,popmort_f, yearlim=c(1951,2000), int.length=1)

# Confirm valid ratetable
is.ratetable(popmort_rt)

#===============================================================================
# Use rs.surv() to estimate non-parametric relative survival by time since diagnosis
# rs.surv(), the follow-up time must be specified in days.
# Vital status (status) is coded as follows
# 1 [Dead: cancer]
# 2 [Dead: other]
# 0 [Alive]
# 4 [Lost to follow-up]

rssurv <- rs.surv(Surv(time  =  stime, 
                     event = status == 1 | status == 2 ) ~ 1,
                rmap = list(age = age*365.241, year = dx),
                ratetable = popmort_rt,
                data = colon2,
                method = "pohar-perme")

summary(rssurv, time = c(0:10) * 365.241, scale = 365.241)