---
title: "Fuzzy survival analysis with vague event times"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Fuzzy survival analysis with vague event times}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 7, fig.height = 4.5)
```

```{r}
library(fuzzySurvival)
library(survival)
```

# 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.

```{r}
# the failure is entirely plausible anywhere in the second quarter
q2 <- fuzzy_number_interval(4, 6)
q2
alpha_cuts(q2, c(0, 0.5, 1))
```

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.

```{r}
# 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))
```

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

```{r}
plot(c(recalled, q2))
```

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

```{r}
round(membership(recalled, c(2, 3, 4, 5, 6, 7, 8, 9, 10)), 3)
```

# 2 The data container

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

```{r}
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))
```

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.

```{r}
fuzzy_surv(c(2, 5, 7, 9), c(1, 0.5, 1, 1))   # "lost to follow-up, maybe failed"
```

# 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].$$

```{r}
data(fuzzy_dfs)
fit <- fkm(fuzzy_surv(ft, fs) ~ 1, data = fuzzy_dfs)
fit
summary(fit, times = c(6, 12, 24))[, c("time", "n.risk", "lower", "surv", "upper")]
```

`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

```{r}
fitg <- fkm(fuzzy_surv(ft, fs) ~ arm, data = fuzzy_dfs)
fitg
summary(fitg, times = 12)[[1]][, c("time", "n.risk", "lower", "surv", "upper")]
```

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

```{r}
plot(fitg, type = "diff")
```

```{r}
flogrank(fuzzy_surv(ft, fs) ~ arm, data = fuzzy_dfs)
```

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.

```{r}
fct <- fcox(fuzzy_surv(ft, fs) ~ age + arm, data = fuzzy_dfs)
fct
```

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

```{r}
alpha_cuts(coef(fct, exponentiate = TRUE), c(0, 0.5, 1))
```

```{r}
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.

```{r}
defuzz(fct$p.value, "centroid")
```

# 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.

```{r}
defuzz(recalled, "centroid")
defuzz(recalled, c("min", "centroid", "max"))
```

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

```{r}
defuzz(fuzzy_number_triangle(2, 5, 8),
       c("midrange", "centroid", "gmirl", "sgmrl", "yager"))
```

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.

```{r}
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)
```

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.

```{r}
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))
```

# 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.

```{r}
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:

```{r}
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)
```

# 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.

```{r}
s <- fkm_fuzzy_surv(fit, times = c(6, 12, 24, 36))
alpha_cuts(s, c(0, 0.5, 1))
defuzz(s, "centroid")
```
```{r}
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.

```{r}
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))))
```

**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`.

```{r}
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))
all(fitc$lower == fitc$upper)
```