Why fuzzy survival analysis?

A survival record is often less precise than the analysis assumes. A patient “progressed during the second quarter” — the failure happened, but its date is only known to lie between two examinations. A subject is “lost to follow-up, last seen alive”: the censoring time is known, but an unobserved failure cannot be excluded. A date recalled as “around the middle of the year” is rounded.

Classical survival analysis forces such a record into a single number: the midpoint of the quarter, the censoring time as if the status were certain. Everything downstream then pretends the record was exact, and the resulting confidence intervals are too narrow.

Fuzzy survival analysis keeps the imprecision. A vague quantity is represented by a fuzzy number — a set of candidate values with a membership degree, or equivalently a family of nested alpha-cuts

\[F_\alpha = \bigl[\, f^-_\alpha,\; f^+_\alpha \,\bigr], \qquad 0 < \alpha \le 1,\]

with F_0 the support and F_1 the core. The parameter alpha reads as “how sure am I”: at alpha = 1 the most plausible values, at alpha = 0 everything still compatible with the information.

1 Representing vague information

Two shapes cover most of what a record looks like.

# the failure is entirely plausible anywhere in the second quarter
q2 <- fuzzy_number_interval(4, 6)
q2
#> <fuzzy_number> with 1 element
#> [1] "I[4, 6]"
alpha_cuts(q2, c(0, 0.5, 1))
#>   id level lower upper
#> 1  1   0.0     4     6
#> 2  1   0.5     4     6
#> 3  1   1.0     4     6

A flat core means the value is completely plausible over the whole interval — this is the natural model for interval-censored-by-examination data. Its consequence, discussed in section 7, is that every alpha-cut is identical.

# a rounded recollection: certainly 3 to 9, concentrated in the middle
recalled <- fuzzy_number_trapezoid(3, 5, 7, 9)
alpha_cuts(recalled, c(0, 0.25, 0.5, 0.75, 1))
#>   id level lower upper
#> 1  1  0.00   3.0   9.0
#> 2  1  0.25   3.5   8.5
#> 3  1  0.50   4.0   8.0
#> 4  1  0.75   4.5   7.5
#> 5  1  1.00   5.0   7.0

Here the alpha-cuts are strictly nested, so alpha really carries information.

plot(c(recalled, q2))

The membership function is recovered from the alpha-cuts, so the two representations never disagree:

round(membership(recalled, c(2, 3, 4, 5, 6, 7, 8, 9, 10)), 3)
#>      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,]    0    0  0.5    1    1    1  0.5    0    0

2 The data container

fuzzy_surv() is the analogue of survival::Surv(): it bundles fuzzy times with failure indicators.

tm <- fuzzy_number_crisp(c(2, 6, 9))
tm[2] <- fuzzy_number_interval(4, 6)   # "failed in the second quarter"
fuzzy_surv(tm, c(1, 1, 1))
#> fuzzy_surv: 3 observations
#>   event times: 2 crisp, 1 vague (mean width 0.667)
#>   status:      3 certain events, 0 certain censored, 0 uncertain

A status strictly between 0 and 1 is a degree of belief that the subject failed. Such a subject enters the risk set normally and contributes a fractional number of events, which is exactly the classical estimator when the degrees are 0 or 1.

fuzzy_surv(c(2, 5, 7, 9), c(1, 0.5, 1, 1))   # "lost to follow-up, maybe failed"
#> fuzzy_surv: 4 observations
#>   event times: 4 crisp, 0 vague (mean width 0)
#>   status:      3 certain events, 0 certain censored, 1 uncertain

3 The fuzzy Kaplan–Meier estimator

At level alpha, the data compatible with the record form the box

\[D_\alpha = \prod_i \bigl[\, \tilde t^-_i(\alpha),\; \tilde t^+_i(\alpha) \,\bigr].\]

Replacing every time by the lower end of its cut makes failures happen as early as possible, so that curve is a lower bound on survival; the upper end gives the upper bound. The alpha-cut of the fuzzy survival function is

\[\hat S_\alpha(t) = \bigl[\, \hat S^-_\alpha(t),\; \hat S^+_\alpha(t) \,\bigr].\]

data(fuzzy_dfs)
fit <- fkm(fuzzy_surv(ft, fs) ~ 1, data = fuzzy_dfs)
fit
#> Fuzzy Kaplan-Meier estimate
#>   n = 240, events = 81, uncertain = 40, alpha levels = 21
#>   mean width of the alpha-cuts at level 0: 0.675
summary(fit, times = c(6, 12, 24))[, c("time", "n.risk", "lower", "surv", "upper")]
#>   time n.risk  lower   surv  upper
#> 1    6    212 0.9781 0.9781 0.9781
#> 2   12    187 0.9360 0.9410 0.9460
#> 3   24    129 0.7668 0.7784 0.7900

lower and upper are the bounds at the requested level; surv is the defuzzified curve. The numbers at risk are reported as usual.

4 Comparing two arms

fitg <- fkm(fuzzy_surv(ft, fs) ~ arm, data = fuzzy_dfs)
fitg
#> 
#> Group: Control 
#>   n = 114, alpha levels = 21
#> 
#> Group: Treated 
#>   n = 126, alpha levels = 21
summary(fitg, times = 12)[[1]][, c("time", "n.risk", "lower", "surv", "upper")]
#>   time n.risk  lower   surv  upper
#> 1   12     84 0.9411 0.9411 0.9411

The fuzzy counterpart of a survival difference plot shows where the difference itself is uncertain.

plot(fitg, type = "diff")

flogrank(fuzzy_surv(ft, fs) ~ arm, data = fuzzy_dfs)
#> Fuzzy log-rank test
#>   n = 240, groups: Control vs Treated
#>   observed / expected events in group 1: 49 / 65.66
#>   crisp chi-square (level 1): 9.3273
#>   fuzzy chi-square range:    [9.3273, 9.3273]
#>   crisp p-value: 0.002258

On crisp data the log-rank statistic reproduces survival::survdiff() exactly.

5 The fuzzy Cox model

The partial likelihood is refitted at the two ends of each alpha-cut, and the alpha-cut of a coefficient spans every value obtained at a level at or below it. Because the map from level to bounds is cumulative, the intervals are nested.

fct <- fcox(fuzzy_surv(ft, fs) ~ age + arm, data = fuzzy_dfs)
fct
#> Fuzzy Cox proportional hazards model
#> 
#>        coef         p fuzzy.lower fuzzy.upper
#> age -0.0465 7.914e-05     -0.0696     -0.0234
#> arm  0.5211  0.005565      0.1527      0.8895
#> 
#> 240 subjects, 121 failures
#> Alpha levels used: 0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1

The hazard ratio is exp(beta); the exponential is monotone, so exponentiating preserves the membership structure of the coefficient.

alpha_cuts(coef(fct, exponentiate = TRUE), c(0, 0.5, 1))
#>   id level     lower     upper
#> 1  1   0.0 0.9545862 0.9545862
#> 2  2   0.0 1.6838509 1.6838509
#> 3  1   0.5 0.9545862 0.9545862
#> 4  2   0.5 1.6838509 1.6838509
#> 5  1   1.0 0.9545862 0.9545862
#> 6  2   1.0 1.6838509 1.6838509

plot(fct)

plot(fct, type = "membership")

The p-value is itself fuzzy: it spans the smallest and the largest value attainable over the fuzzy coefficient, given the standard error of the crisp fit.

defuzz(fct$p.value, "centroid")
#> [1] 7.914476e-05 5.564754e-03

6 Reducing a fuzzy number to a number

A fuzzy number carries a whole family of alpha-cuts; a type-reduction rule produces a single number. Nine classical rules are implemented.

defuzz(recalled, "centroid")
#> [1] 6
defuzz(recalled, c("min", "centroid", "max"))
#>          [,1]
#> min         3
#> centroid    6
#> max         9

Every rule of a symmetric fuzzy number returns its centre, which is the easiest way to check them:

defuzz(fuzzy_number_triangle(2, 5, 8),
       c("midrange", "centroid", "gmirl", "sgmrl", "yager"))
#>          [,1]
#> midrange    5
#> centroid    5
#> gmirl       5
#> sgmrl       5
#> yager       5

The survival curve itself can be turned into fuzzy numbers and reduced: # 7 What the fuzzy bounds do and do not mean

Three caveats, stated plainly.

Interval data give level-invariant bounds. A flat-core fuzzy number has the same alpha-cut at every level, so the bounds do not shrink with alpha. The uncertainty is real but it is not graded — the package cannot invent a gradation the record does not contain. When the plausibility genuinely varies within the interval, use a trapezoid instead.

tn <- as_fuzzy_number(fuzzy_dfs$ft)
lo <- tn$support[, 1]; hi <- tn$support[, 2]; w <- hi - lo
graded <- fuzzy_number_trapezoid(lo, lo + 0.4 * w, lo + 0.6 * w, hi)
fitg2 <- fkm(fuzzy_surv(I(list(graded)), fuzzy_dfs$fs) ~ 1)
wid <- fitg2$upper[, match(24, fitg2$time)] - fitg2$lower[, match(24, fitg2$time)]
round(wid[c(1, 6, 11, 16, 21)], 5)
#> [1] 0.02322 0.00150 0.00114 0.00033 0.00011

The width of the cut of S(24) now shrinks as the level grows.

The bounds are evaluated at three points of the data box. At each level the Kaplan–Meier curve is computed on the lower end, the centre and the upper end of the data box, and the extremes are reported. Because the estimator is not monotone in the event times — moving one subject past another changes both risk sets — those three points do not provably bracket every admissible dataset. The intervals are a good outer approximation, not a rigorous envelope.

Vagueness is not bias. Coarsening failures to examination dates biases the estimate; widening intervals reflects the imprecision but does not undo it. On the example data the gap between the fuzzy curve and the curve of the unobserved exact times is much larger than the width of the fuzzy band.

truefit <- survfit(Surv(true_time, true_status) ~ 1, data = fuzzy_dfs)
tt2 <- summary(truefit)$time
truth <- vapply(tt2, function(t) summary(truefit, times = t)$surv, numeric(1))
gap <- abs(defuzz(fit, "centroid", times = fit$time) -
             fuzzySurvival:::km_eval(tt2, truth, fit$time))
data.frame(fuzziness = max(fit$upper[1, ] - fit$lower[1, ]),
           coarsening_bias = max(gap))
#>    fuzziness coarsening_bias
#> 1 0.05669086       0.5283677

8 Validating your own setting

sim_fuzzy_surv() generates exact times, then coarsens them exactly as a record does, so the fuzzy estimators can be checked against a known truth.

set.seed(7)
d <- sim_fuzzy_surv(n = 300, beta = c(age = 0.35, arm = -0.55),
                    prop_vague = 0.7, prop_uncertain = 0.2)
fit_sim <- fkm(fuzzy_surv(ft, fs) ~ 1, data = d)
plot(fit_sim)

The curve built on the midpoints is one of the representative datasets, so the level 0 cut must contain it:

cut_sim <- as_fuzzy_number(d$ft)$acut(0)
mid <- signif((cut_sim$lower[, 1] + cut_sim$upper[, 1]) / 2, 12)
obs_km <- fuzzySurvival:::km_fractional(mid, d$fs)
rv <- fuzzySurvival:::km_eval(obs_km$time, obs_km$surv, fit_sim$time)
all(rv >= fit_sim$lower[1, ] - 1e-8) && all(rv <= fit_sim$upper[1, ] + 1e-8)
#> [1] TRUE

9 Reference

fuzzy_surv() builds the response; fkm(), fcox() and flogrank() estimate; fkm_fuzzy_surv() turns a fitted curve into fuzzy numbers; defuzz() reduces them; alpha_cuts(), membership(), support(), core() and fuzzy_spread() inspect them. sim_fuzzy_surv() generates test problems.

s <- fkm_fuzzy_surv(fit, times = c(6, 12, 24, 36))
alpha_cuts(s, c(0, 0.5, 1))
#>    id level     lower     upper
#> 1   1   0.0 0.9780630 0.9780630
#> 2   2   0.0 0.9359943 0.9459864
#> 3   3   0.0 0.7667985 0.7900214
#> 4   4   0.0 0.6062923 0.6257605
#> 5   1   0.5 0.9780630 0.9780630
#> 6   2   0.5 0.9359943 0.9459864
#> 7   3   0.5 0.7667985 0.7900214
#> 8   4   0.5 0.6062923 0.6257605
#> 9   1   1.0 0.9780630 0.9780630
#> 10  2   1.0 0.9359943 0.9459864
#> 11  3   1.0 0.7667985 0.7900214
#> 12  4   1.0 0.6062923 0.6257605
defuzz(s, "centroid")
#> [1] 0.9780630 0.9409903 0.7784100 0.6160264

plot(fit)

The guaranteed properties

Two properties are enforced by construction and covered by the test suite.

Nested alpha-cuts. The lower bound never decreases and the upper bound never increases as alpha grows.

L <- fit$lower; U <- fit$upper
c(nested_lower = all(vapply(seq_len(ncol(L)),
                           function(k) all(diff(L[, k]) >= -1e-9), logical(1))),
  nested_upper = all(vapply(seq_len(ncol(U)),
                           function(k) all(diff(U[, k]) <= 1e-9), logical(1))))
#> nested_lower nested_upper 
#>         TRUE         TRUE

Exact reduction on crisp data. When every event time is known, the fuzzy estimator is the classical Kaplan–Meier estimator, which makes the implementation checkable against survival.

set.seed(42)
n <- 80
tt <- rexp(n, 0.05) * 10
st <- rbinom(n, 1, 0.7)
tt[!st] <- tt[!st] + runif(sum(!st), 0, 5)

fitc <- fkm(fuzzy_surv(tt, st) ~ 1)
ref <- vapply(fitc$time, function(t) summary(survfit(Surv(tt, st) ~ 1),
                                            times = t)$surv, numeric(1))
ref[is.na(ref)] <- 1
max(abs(defuzz(fitc, "midrange", level = 1, times = fitc$time) - ref))
#> [1] 0
all(fitc$lower == fitc$upper)
#> [1] TRUE

Fuzzy survival analysis with vague event times