Special Supervised Tasks: Time Series and Survival

Healthcare data often include multiple visits per patient, changing exposures, incomplete follow-up, and seasonal demand. Ordinary random rows and ordinary numeric regression may be inappropriate.

1. Two different questions about time

Forecasting: How many fictional appointments will a clinic receive next month? Observations form a chronological sequence, often with seasonality and autocorrelation.
Survival analysis: When will a follow-up occur for a person, given some people have not yet had a follow-up by the end of observation? We must model censoring.

A patient having no observed return is not automatically the same as “never returns.”

Example: Five fictional patients are watched for 12 months. months is the time until a follow-up visit, and returned is 1 if the visit happened inside the 12 months.

toy <- data.frame(
  patient = c("A", "B", "C", "D", "E"),
  months = c(3, 12, 7, 12, 5),
  returned = c(1, 0, 1, 0, 1)
)
toy

Output

  patient months returned
1       A      3        1
2       B     12        0
3       C      7        1
4       D     12        0
5       E      5        1

Patients B and D were still waiting at month 12. Their true return time is unknown; we only know it is longer than 12 months. That is a survival question (time until an event, with some people not yet having had it). A forecasting question would instead be “how many visits will the clinic see in week 20?”

2. Time-series forecasting

A time series is indexed by date/time; successive observations often depend on one another. Choose a forecast horizon and evaluate on later periods, never by randomly mixing future periods into training.

set.seed(500)
week <- seq_len(104)
fake_counts <- 40 + 4 * sin(2 * pi * week / 13) +
  .15 * week + rnorm(length(week), sd = 3)
training_counts <- ts(fake_counts[1:80], frequency = 13)
test_counts <- fake_counts[81:104] # reserved, not used to fit
fit <- arima(training_counts, order = c(1, 0, 0))
forecast <- predict(fit, n.ahead = length(test_counts))$pred
c(test_mae = mean(abs(test_counts - forecast)),
  naive_last_value_mae = mean(abs(test_counts - fake_counts[80])))

Output

            test_mae naive_last_value_mae 
            7.170464             4.082597

ARIMA combines autoregressive terms, differencing and moving-average terms. The order = c(1, 0, 0) used here has one autoregressive term and no differencing or moving average. Other families: seasonal ARIMA, exponential smoothing/ETS, regression with seasonal features, Prophet, boosted trees with lagged variables, and neural sequence models.

The ARIMA model’s test MAE is 7.17 and the naive last-value forecast’s is 4.08, so the naive forecast is better on these 24 weeks. The reason is the upward trend. An ARIMA(1, 0, 0) forecast decays toward the training mean (about 46), which is below the trend, while the last training value (52.5) is close to the level of the test weeks (mean about 53.5). This holds for this simulated series only; it is the reason to always score a naive baseline next to any model.

This toy demonstration uses one fixed model selected before examining the late-period score. In a real study, choose lag features, hyperparameters, and transformations on earlier rolling-origin training windows.

3. Why lag features can leak

A model predicting visits in week 20 may use the count from week 19 (known beforehand) but not week 21. A random row split can put week 21 in training and week 20 in test, giving the model unrealistic future knowledge. Define a rolling window or expanding-window split matching the intended deployment.

Example: A lag feature is the previous week’s count, shifted down one row. Here are six typed weekly counts.

weekly <- data.frame(week = 1:6, visits = c(31, 35, 33, 38, 36, 40))
weekly$visits_last_week <- c(NA, weekly$visits[-6])
weekly

Output

  week visits visits_last_week
1    1     31               NA
2    2     35               31
3    3     33               35
4    4     38               33
5    5     36               38
6    6     40               36

Week 4 uses 33 from week 3, which is already known when week 4 starts. Week 1 has no earlier week, so its lag is NA. Using weekly$visits[-1] instead (the next week) would put the future into the row, which is the leak.

4. Survival analysis: time and event indicator

Censoring means the exact event time is not observed before follow-up ends (or is lost). A survival record needs both observed time and event status (1 event observed, 0 right-censored). Dropping censored rows produces a biased sample of people who experienced the event earlier.

if (!requireNamespace("survival", quietly = TRUE)) {
  stop("Install the recommended R package 'survival' to run this chapter.")
}
set.seed(500)
n <- 120
age <- 20 + (seq_len(n) * 7) %% 61
true_time <- rexp(n, rate = exp(-4.5 + 0.015 * age))
followup_end <- 55
observed_time <- pmin(true_time, followup_end)
event_observed <- as.integer(true_time <= followup_end)
surv_data <- data.frame(
  age = age, time = observed_time, event = event_observed
)
table(surv_data$event)

Output

 0  1 
42 78

The table counts the simulated people who had the event (1) and those still event-free when follow-up ended (0, right-censored).

All 120 times were simulated; none represent actual patient encounters.

5. Kaplan–Meier and Cox proportional hazards

Kaplan–Meier estimates the survival curve (probability of not having had the event by each time). Cox regression models how predictor differences relate to the hazard under a proportional hazards assumption.

surv_obj <- survival::Surv(surv_data$time, surv_data$event)
km <- survival::survfit(surv_obj ~ 1, data = surv_data)
plot(km, xlab = "Simulated follow-up time",
     ylab = "Estimated event-free probability")

Kaplan-Meier curve of estimated event-free probability against simulated follow-up time, as a solid step line with dashed confidence bands. The curve starts at 1 and falls steadily as follow-up time increases, for made-up data.

cox <- survival::coxph(surv_obj ~ age, data = surv_data)
summary(cox)

Output

Call:
survival::coxph(formula = surv_obj ~ age, data = surv_data)

  n= 120, number of events= 78 

        coef exp(coef) se(coef)   z Pr(>|z|)  
age 0.015804  1.015929 0.006322 2.5   0.0124 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    exp(coef) exp(-coef) lower .95 upper .95
age     1.016     0.9843     1.003     1.029

Concordance= 0.591  (se = 0.033 )
Likelihood ratio test= 6.29  on 1 df,   p=0.01
Wald test            = 6.25  on 1 df,   p=0.01
Score (logrank) test = 6.35  on 1 df,   p=0.01

The summary reports 78 events among the 120 people. The age coefficient is 0.0158, which is a hazard ratio of 1.016 per year of age (95% interval 1.003 to 1.029, p = 0.012). The Kaplan-Meier plot falls from 1 as follow-up time passes, and the dashed bands show the uncertainty around that estimate.

A hazard ratio is not a risk ratio or a prediction probability. Check proportional hazards, competing events, time-varying effects, and censoring mechanisms before interpretation. With repeated patients or center effects, account for dependence appropriately.

6. Advanced survival families

Parametric survival: exponential, Weibull, log-normal; explicit time-distribution assumptions.
Accelerated failure time models: describe multiplicative effects on event time.
Flexible / competing risk models: account for multiple mutually exclusive event types.
Random survival forests and gradient boosted survival: nonlinear survival prediction with censoring.
Recurrent event / multi-state models: repeated transitions and events.

R’s survival, flexsurv, cmprsk, randomForestSRC, ranger, and tidymodels ecosystem support different portions; packages and censoring support vary.

Example: A parametric model assumes a time distribution. This fits a Weibull accelerated failure time model to the simulated data with survreg() from the survival package.

weib <- survival::survreg(surv_obj ~ age, data = surv_data, dist = "weibull")
round(coef(weib), 3)

Output

(Intercept)         age 
      4.797      -0.018

exp(coef(weib)["age"])

Output

      age 
0.9823779

The age coefficient is -0.018, so exp() of it is about 0.98. In an accelerated failure time model, that is the multiplier on the event time for one more year of age: older people have slightly shorter times to the event, about 2% shorter per year. This points the same way as the Cox result in section 5 (more hazard with age), under a different assumption about the time distribution.

Specify index date, time horizon, minimum follow-up, administrative censoring, patient grouping, and temporal validation. For survival, use censored-aware measures such as concordance, time-dependent AUC, calibration, and integrated Brier score with appropriate censoring adjustments. Ordinary classification at an arbitrary time point can be biased if follow-up is incomplete.

Example: Concordance is the share of comparable pairs of people for whom the model gives the higher risk to the person with the earlier event. 0.5 is a coin flip and 1 is perfect ordering. It handles censored people correctly.

survival::concordance(cox)

Output

Call:
concordance.coxph(object = cox)

n= 120 
Concordance= 0.5912 se= 0.03317
concordant discordant     tied.x     tied.y    tied.xy 
      3686       2541         52          0          0

The concordance is 0.59 (standard error 0.033). Of the 6,227 comparable pairs without ties in age, the Cox model ranked 3,686 correctly. Pairs tied in age count as half, so the printed 0.591 is (3686 + 52/2) / (3686 + 2541 + 52), slightly different from 3686 / 6227 = 0.592. That is only a little better than a coin flip, so age alone orders these simulated people weakly.

Try it

Explain why a row-wise random split may leak future information in forecasting.

Show solution

Answer: With a random split, some test weeks come before training weeks, so the model is trained on information that would not exist yet at prediction time.

Step by step

Forecasting means predicting weeks that have not happened yet, so training data must all lie before the test data.
Split 30 weeks at random and check, for each test week, whether a later week is in the training set.

set.seed(1)
test_weeks <- sort(sample(1:30, 6))
train_weeks <- setdiff(1:30, test_weeks)
test_weeks

Output

[1]  1  2  4  7 23 25

sapply(test_weeks, function(w) any(train_weeks > w))

Output

[1] TRUE TRUE TRUE TRUE TRUE TRUE

3. Every TRUE is a test week with at least one later week in training. A model could learn from those later weeks (and from lag features built on them), which looks like skill but is not available in real use. A chronological split (weeks 1 to 24 train, 25 to 30 test) has no such week.

Recalculate the MAE for a naive forecast predicting the last training value for every test week.

Show solution

Answer: The naive last-value forecast has a test MAE of 4.08, which is lower (better) than the ARIMA model’s 7.17 on these 24 weeks.

Step by step

The last training value is week 80, fake_counts[80]. Predict that same number for each of the 24 test weeks.
MAE is the mean absolute difference between the test counts and the forecast.

mean(abs(test_counts - fake_counts[80]))

Output

[1] 4.082597

3. The output is the naive MAE that section 2 also printed. Compare it with the model’s test_mae there: the smaller number is the better forecast on these 24 weeks.

Increase followup_end and observe how many simulated events are censored.

Show solution

Answer: Longer follow-up means fewer censored people: 75 of 120 at followup_end = 20, 42 at 55 (the chapter’s value), 15 at 100, and 1 at 200.

Step by step

Reuse the simulated true_time from section 4. Anyone with true_time greater than followup_end is censored.
Count the censored people and the events for several follow-up lengths.

ends <- c(20, 55, 100, 200)
data.frame(followup_end = ends,
           censored = sapply(ends, function(e) sum(true_time > e)),
           events = sapply(ends, function(e) sum(true_time <= e)))

Output

  followup_end censored events
1           20       75     45
2           55       42     78
3          100       15    105
4          200        1    119

3. Read across the rows: censored plus events always equals 120, and the longer we watch, the more events we see. The row for 55 matches the chapter’s table(surv_data$event).

Why is “not returned by study end” different from “will never return”?

Show solution

Answer: “Not returned by study end” means the return time is unknown and at least as long as the follow-up. “Will never return” is a claim about the whole future, which the data cannot show.

Step by step

Follow-up stops on a fixed date, so a person with no return yet may return the day after.
Treating those people as non-returners (or dropping them) changes the average. Compare the mean time using only people with an event against all rows.

with(surv_data, c(events_only = mean(time[event == 1]),
                  all_rows = mean(time),
                  censored_people = sum(event == 0)))

Output

    events_only        all_rows censored_people 
       19.22794        31.74816        42.00000

3. The mean time among only the 78 people with an event is 19.2, against 31.7 over all 120 rows. The 42 censored people were left out of the first mean, and they are the ones who took longest, so dropping them makes returns look much quicker than they are. Survival methods keep them and use the fact that they were event-free up to their last time.

Explain what a Cox hazard ratio does not mean.

Show solution

Answer: It is not a risk ratio and not a predicted probability. It is a ratio of instantaneous event rates (hazards) for a one-unit difference in the predictor, assuming the hazards stay proportional over time. In this model the hazard ratio is 1.016 per year of age (interval 1.003 to 1.029).

Step by step

Get the hazard ratio per year of age, with its confidence interval, from the fitted Cox model.

round(exp(cbind(HR = coef(cox), confint(cox))), 4)

Output

        HR  2.5 % 97.5 %
age 1.0159 1.0034 1.0286

2. Read it: for each extra year of age, the hazard (the event rate at a given moment among people still event-free) is multiplied by the HR value. That is a rate comparison, not “the chance of an event by month 55” and not “X% more people get the event”. 3. To get a probability you need the survival curve (section 5), which also depends on time. The hazard ratio is also only trustworthy if the proportional hazards assumption is reasonable.

Checkpoint

You can distinguish forecasting from survival, understand time-ordered evaluation and right censoring, and identify ARIMA, Cox, Kaplan–Meier, survival forests, and more specialized methods.

1. Two different questions about time

2. Time-series forecasting

3. Why lag features can leak

4. Survival analysis: time and event indicator

5. Kaplan–Meier and Cox proportional hazards

6. Advanced survival families

7. Validation for time-related tasks

Try it

Checkpoint