@@ -26,6 +26,10 @@ of a box subject to area limits:
2626#| label: setup
2727#| include: false
2828library(CVXR)
29+ library(ggplot2)
30+ library(survival)
31+ source("_check_status.R")
32+ theme_set(theme_bw())
2933```
3034
3135``` {r}
@@ -41,6 +45,147 @@ The objective `x * y * z` multiplies three variables --- not DCP at all --- yet
4145solves it to global optimality. See the
4246[ DGP Tutorial] ( https://cvxr.rbind.io/examples/dgp/tutorial.html ) .
4347
48+ ### A worked example: ordering survival curves
49+
50+ DGP is not only for engineering design; it appears in statistics wherever a model is
51+ built from products and ratios of positive quantities. Here is one such problem, from
52+ @LimKimWang :2009.
53+
54+ Suppose we estimate survival curves $S_1, \ldots, S_N$ for several populations and
55+ know on scientific grounds that they are ** stochastically ordered** : population $a$
56+ dies earlier than $b$, meaning $S_a(t) \le S_b(t)$ at every $t$ (written $S_a \preceq
57+ S_b$). The unconstrained nonparametric estimates --- the Kaplan--Meier curves --- often
58+ * cross* where they should not. @LimKimWang :2009 observed that imposing the ordering
59+ turns the estimation into a geometric program, which ` CVXR ` solves directly.
60+
61+ Take ** right-censored** data. For population $i$, at each distinct failure time
62+ $t_ {ij}$ let $n_ {ij}$ subjects be at risk and $d_ {ij}$ of them fail. Write the one-step
63+ conditional survival probability $p_ {ij} = S_i(t_ {ij}) / S_i(t_ {i,j-1})$ and $q_ {ij} =
64+ 1 - p_ {ij}$, so survival is the running product $S_i(t_ {ij}) = \prod_ {r \le j} p_ {ir}$
65+ (the Kaplan--Meier form). The likelihood is the ** monomial**
66+ $$
67+ \mathcal L = \prod_{i} \prod_{j} p_{ij}^{\,n_{ij} - d_{ij}}\, q_{ij}^{\,d_{ij}},
68+ $$
69+ which is log-log affine and so drops straight into a GP. Three ingredients complete the
70+ model:
71+
72+ - ** Sum-to-one** , relaxed from $p_ {ij} + q_ {ij} = 1$ to the posynomial inequality
73+ $p_ {ij} + q_ {ij} \le 1$ (a GP admits only monomial equalities). It is tight at the
74+ optimum, since both exponents are nonnegative.
75+ - ** Monotonicity** of each survival curve, which is * automatic* here: $S_i$ is a
76+ product of factors $p \le 1$, so it is nonincreasing with no extra constraint.
77+ - ** Ordering** $S_a \preceq S_b$: at each pooled failure time $t$, the monomial
78+ $\bigl(\prod_ {t_ {ar} \le t} p_ {ar}\bigr)\bigl(\prod_ {t_ {br} \le t} p_ {br}\bigr)^{-1}
79+ \le 1$.
80+
81+ Everything is a monomial except the sum-to-one posynomial, so ` is_dgp() ` holds. We use
82+ the oropharyngeal-carcinoma trial of Kalbfleisch and Prentice, grouping patients by
83+ tumor stage ` T ` and nodal stage ` N ` :
84+
85+ ``` {r}
86+ #| label: surv-data
87+ oropharynx <- readRDS("data/oropharynx.rds")
88+
89+ grp_counts <- function(T, N) { # discrete-hazard counts
90+ s <- subset(oropharynx, tstage == T & nstage == N)
91+ dt <- sort(unique(s$days[s$status == 1])) # distinct failure times
92+ data.frame(t = dt,
93+ d = sapply(dt, function(x) sum(s$days == x & s$status == 1)), # failures
94+ n = sapply(dt, function(x) sum(s$days >= x))) # at risk
95+ }
96+ ```
97+
98+ ` fit_rc() ` assembles the monomial likelihood, the relaxed sum-to-one constraints, and
99+ one ordering monomial per pooled failure time for each requested pair ` c(a, b) `
100+ (meaning $S_a \preceq S_b$):
101+
102+ ``` {r}
103+ #| label: surv-fit
104+ fit_rc <- function(groups, order = list()) {
105+ G <- length(groups)
106+ p <- lapply(groups, function(g) Variable(nrow(g), pos = TRUE))
107+ q <- lapply(groups, function(g) Variable(nrow(g), pos = TRUE))
108+ mono <- function(v, a) Reduce(`*`, lapply(which(a != 0), function(j) v[j]^a[j]))
109+ Scum <- function(i, ix) if (length(ix)) prod_entries(p[[i]][ix]) else 1 # empty prod = 1
110+
111+ objective <- Reduce(`*`, lapply(seq_len(G), function(i)
112+ mono(p[[i]], groups[[i]]$n - groups[[i]]$d) * mono(q[[i]], groups[[i]]$d)))
113+ constraints <- lapply(seq_len(G), function(i) p[[i]] + q[[i]] <= 1)
114+
115+ for (o in order) { # S_a(t) <= S_b(t) at every pooled failure time
116+ a <- o[1]; b <- o[2]
117+ for (tk in sort(unique(c(groups[[a]]$t, groups[[b]]$t)))) {
118+ ia <- which(groups[[a]]$t <= tk); ib <- which(groups[[b]]$t <= tk)
119+ constraints <- c(constraints, Scum(a, ia) * Scum(b, ib)^(-1) <= 1)
120+ }
121+ }
122+ prob <- Problem(Maximize(objective), constraints)
123+ val <- psolve(prob, gp = TRUE)
124+ check_solver_status(prob)
125+ list(is_dgp = is_dgp(prob), logLik = log(val),
126+ S = lapply(seq_len(G), function(i) cumprod(value(p[[i]])))) # S = running product
127+ }
128+ ```
129+
130+ With no ordering the GP must return the ordinary Kaplan--Meier estimator --- a
131+ reassuring check:
132+
133+ ``` {r}
134+ #| label: surv-km-check
135+ g <- grp_counts(3, 3)
136+ fit <- fit_rc(list(g))
137+ km <- survfit(Surv(days, status) ~ 1,
138+ data = subset(oropharynx, tstage == 3 & nstage == 3))
139+ max(abs(fit$S[[1]] - summary(km, times = g$t)$surv)) # ~ 0: GP == Kaplan-Meier
140+ ```
141+
142+ Now the ordering. Take the three T = 3 groups differing in nodal involvement,
143+ $(T,N) = (3,1), (3,2), (3,3)$, and impose $S_3 \preceq S_2 \preceq S_1$ (more nodal
144+ involvement, worse survival):
145+
146+ ``` {r}
147+ #| label: surv-fig
148+ #| fig-width: 8
149+ #| fig-height: 3.3
150+ g1 <- list("(3,1)" = grp_counts(3, 1),
151+ "(3,2)" = grp_counts(3, 2),
152+ "(3,3)" = grp_counts(3, 3))
153+
154+ rc_free <- fit_rc(g1)
155+ rc_order <- fit_rc(g1, list(c(3, 2), c(2, 1))) # S3 <= S2 <= S1
156+
157+ step_rc <- function(fit, groups, panel)
158+ do.call(rbind, lapply(seq_along(groups), function(i) {
159+ g <- groups[[i]]; S <- fit$S[[i]]
160+ data.frame(t = c(0, g$t), S = c(1, S),
161+ grp = paste0("(T,N)=", names(groups)[i]), panel = panel)
162+ }))
163+
164+ df <- rbind(step_rc(rc_free, g1, "(a) Kaplan-Meier (unconstrained)"),
165+ step_rc(rc_order, g1, "(b) ordered: S3 <= S2 <= S1"))
166+
167+ ggplot(df, aes(t, S, colour = grp)) + geom_step(linewidth = 0.6) +
168+ facet_wrap(~panel) + coord_cartesian(xlim = c(0, 1850), ylim = c(0, 1)) +
169+ labs(x = "time (days)", y = "survival probability", colour = NULL) +
170+ theme(legend.position = "inside", legend.position.inside = c(0.9, 0.8),
171+ legend.background = element_blank())
172+ ```
173+
174+ The unconstrained Kaplan--Meier curves cross --- $(3,2)$ sits above $(3,1)$ for much of
175+ the range, contradicting $S_2 \le S_1$ --- while the ordered fit removes the crossings,
176+ reproducing Figure 1 of @LimKimWang :2009. The ordering costs only a little likelihood:
177+
178+ ``` {r}
179+ #| label: surv-ll
180+ c(unconstrained = round(rc_free$logLik, 2), ordered = round(rc_order$logLik, 2))
181+ ```
182+
183+ The website's
184+ [ survival-ordering example] ( https://cvxr.rbind.io/examples/dgp/survival-ordering.html )
185+ carries this further --- partial (non-total) orders, a stronger ** uniform** ordering
186+ that constrains the hazard ratio, and the ** interval-censored** case, where the same GP
187+ recipe applies with a different parameterization.
188+
44189## Disciplined Quasiconvex Programming (DQCP)
45190
46191A ** quasiconvex** function has convex sub-level sets but need not be convex itself
0 commit comments