set.seed(1215)trial_data <- cetux %>%  as_tibble() %>%  filter(treat == "Control") %>%  mutate(age = rnorm(n(), mean = 60, sd = 15),          age = cut(age,     breaks = c(0, 40, 60, Inf),     labels = c("Young", "Middle", "Elderly")  )) %>%  rename(event = d,         time = years) %>%  select(time, event, age) %>%  slice(sample(1:n(),80)) # take a random small sample to reduce runtime.# Reference data is the trial population.ref_data <- trial_data %>%  select(age)

mod_age_strata <- survextrap(Surv(time, event) ~ age,                              data = trial_data,                             fit_method = "opt",                             iter = 200, # low iterations to reduce runtime                             smooth_model = "random_walk") plot(mod_age_strata)

# Time vector to evaluate hazard and survival times.t_vec <- seq(from = 0.001, to = 40, length.out = 50)# Get the standardised survival and hazard rates.# takes a few minutes to run but not too long.test_haz <- hazard(mod_age_strata,                    standardise_to(ref_data),                   t = t_vec,                   summ_fns = list(mean = mean))test_surv <- survival(mod_age_strata,                       standardise_to(ref_data),                      t = t_vec,                      summ_fns = list(mean = mean))haz_stand1 <- test_haz %>%  mutate(method = "survextrap method") surv_stand1 <- test_surv %>%  mutate(method = "survextrap method") # To check# Get the hazards and survival rates for the ref_data without standardising.# takes a couple of minutes to run.test_haz2 <- hazard(mod_age_strata,                     newdata = ref_data,                    t = t_vec,                    summ_fns = list(mean=mean))test_surv2 <- survival(mod_age_strata,                        newdata = ref_data,                       t = t_vec,                       summ_fns = list(mean=mean))haz_stand2  <- test_haz2 %>%  rename(hazard = mean) %>%  bind_cols(test_surv2 %>% # add the survival rates.              select(mean) %>%              rename(survival = mean)) %>%  group_by(t) %>%  summarise(mean = sum(hazard * survival)/sum(survival)) %>% # weighting by survival times.  mutate(method = "manual") %>%  select(t, mean, method)surv_stand2 <- test_surv2 %>%  rename(survival = mean) %>%  group_by(t) %>%  summarise(mean = mean(survival)) %>%  mutate(method = "manual")

haz_stand1 %>%  bind_rows(haz_stand2) %>%  ggplot(aes(x = t, y = mean, colour = method, linetype = method))+  theme_classic()+  geom_line(linewidth = 1.2)+  scale_linetype_manual(values = c("solid", "11")) +  xlab("Time")+  scale_y_continuous("Hazard")

surv_stand1 %>%  bind_rows(surv_stand2) %>%  ggplot(aes(x = t, y = mean, colour = method, linetype = method))+  theme_classic()+  geom_line(linewidth = 1.2)+  scale_linetype_manual(values = c("solid", "11")) +  xlab("Time")+  scale_y_continuous("Survival", labels = scales::percent)