Multivariate Analysis

Every distribution and model covered so far is univariate: each unit has a single event time, and units are treated as independent replicates. Multivariate survival analysis relaxes the independence: it models several correlated event-time series jointly. A pair of failure times on the two bearings of one shaft, the times to two related complications in one patient, or the lifetimes of two components sharing an environment are all naturally dependent, and pretending otherwise understates the joint risk.

How much can independence get wrong? Take two pumps in a redundant (parallel) pair, each with a 10% chance of failing within a year. If they fail independently, the chance that both fail — and the system goes down — is \(0.1 \times 0.1 = 1\%\). If they share a cause of early failure (the same batch of seals, the same contaminated fluid), the joint probability can be several times larger, even though each pump on its own is exactly as reliable as before. Nothing in the two marginal distributions reveals this; it lives entirely in the dependence between them.

The difficulty is that a joint distribution mixes two very different things: what each series looks like on its own, and how the series move together. The copula is the device that separates them. SurPyval’s surpyval.multivariate module implements bivariate copula models whose margins are ordinary SurPyval distributions; worked examples are on the Multivariate Modelling with SurPyval page and the API on Multivariate Modelling.

Copulas and Sklar’s theorem

The probability integral transform

The starting point is a fact every simulation relies on: if \(X\) has a continuous CDF \(F\), then \(U = F(X)\) is uniform on \([0, 1]\). Transforming each series by its own CDF therefore strips away its marginal shape — Weibull, LogNormal, whatever it was — and leaves a uniform variable. What remains between the transformed variables \(U_1 = F_1(X_1)\) and \(U_2 = F_2(X_2)\) is pure dependence: if the bearings tend to fail early together, \(U_1\) and \(U_2\) tend to be small together.

Definition and Sklar’s theorem

A copula \(C(u_1, u_2)\) is simply a joint distribution function whose margins are uniform on \([0, 1]\): \(C(u_1, u_2) = P(U_1 \leq u_1, U_2 \leq u_2)\). Its whole job is to encode dependence, stripped of any marginal shape. Sklar’s theorem [Sklar1959mv] says that any joint distribution can be written this way: given marginal CDFs \(F_1, F_2\), the joint CDF is

\[H(x_1, x_2) = C\big(F_1(x_1),\, F_2(x_2)\big).\]

For continuous margins the copula is unique, and conversely any copula combined with any margins through this formula gives a valid joint distribution. Because \(F_1\) and \(F_2\) map each series onto uniform margins, the copula \(C\) carries only the dependence structure. This is the key modelling freedom: the margins can be any survival distribution — a Weibull for one series, a LogNormal for the other — while the copula, chosen separately, governs how they are coupled. The marginal question (“how long does this component last?”) and the dependence question (“do the two fail together?”) are answered by different parts of the model.

Two copulas bracket all the others. Independence is \(C(u_1, u_2) = u_1 u_2\). Every copula lies between the Fréchet-Hoeffding bounds

\[\max(u_1 + u_2 - 1,\, 0) \;\leq\; C(u_1, u_2) \;\leq\; \min(u_1, u_2),\]

the lower bound describing perfect negative dependence (one series is a decreasing function of the other) and the upper bound perfect positive dependence (one is an increasing function of the other). A copula family with a parameter moves between independence and one or both of these extremes.

Everything else from the copula

Writing \(u_1 = F_1(x_1)\) and \(u_2 = F_2(x_2)\), every joint quantity used in survival analysis follows from \(C\) and the margins:

the joint CDF \(H(x_1, x_2) = P(X_1 \leq x_1, X_2 \leq x_2) = C(u_1, u_2)\) — both have failed by their respective times;
the joint survival function, by inclusion-exclusion,

\[S(x_1, x_2) = P(X_1 > x_1,\, X_2 > x_2) = 1 - u_1 - u_2 + C(u_1, u_2);\]
the joint density \(f(x_1, x_2) = c(u_1, u_2)\, f_1(x_1)\, f_2(x_2)\), where \(c = \partial^2 C / \partial u_1 \partial u_2\) is the copula density and \(f_1, f_2\) the marginal densities;
the conditional distribution of one series given the exact value of the other, through the partial derivative of \(C\) (often called the h-function):

\[P(X_2 \leq x_2 \mid X_1 = x_1) = \frac{\partial C(u_1, u_2)}{\partial u_1}.\]

The h-function is also the key to simulation: draw \(u_1\) uniform, draw a second uniform \(w\), and solve \(\partial C/\partial u_1(u_1, u_2) = w\) for \(u_2\); then map back through the margins’ quantile functions, \(x_j = F_j^{-1}(u_j)\).

A note on orientation. SurPyval applies the copula to the marginal CDFs, as in Sklar’s theorem above. Some survival texts instead apply a copula to the marginal survival functions, \(S(x_1, x_2) = \tilde{C}(S_1(x_1), S_2(x_2))\) (a “survival copula”); for example, the shared gamma-frailty model is a Clayton copula applied in this second way. The same family name then describes a different model — the tail that is dependent is flipped — so when comparing with published results check which convention was used.

Families of dependence

Different copula families describe qualitatively different dependence, especially in the tails — whether two series tend to fail together at short lives (lower-tail dependence) or survive together to long lives (upper-tail dependence):

Copula	Parameter	Dependence
Independence	none	none (\(\tau = 0\))
Clayton	\(\theta > 0\)	lower-tail (joint early failure)
Gumbel	\(\theta \geq 1\)	upper-tail (joint long survival)
Frank	\(\theta \neq 0\)	symmetric, no tail
Gaussian	\(\rho \in (-1, 1)\)	symmetric, no tail
Joe	\(\theta \geq 1\)	upper-tail, stronger than Gumbel’s
AMH (Ali-Mikhail-Haq)	\(-1 \leq \theta \leq 1\)	weak only (\(-0.18 \leq \tau \leq 1/3\)), no tail
StudentT	\(\rho \in (-1, 1)\), \(\nu > 0\)	symmetric, both tails

The strength of dependence is summarised by rank measures that do not depend on the margins — Kendall’s \(\tau\) and Spearman’s \(\rho\) — and the tendency to fail (or survive) together in the extremes by the tail dependence coefficients. Choosing a family is largely a question of which tail behaviour matches the physics or the clinical reality.

Measuring dependence

The usual (Pearson) correlation is a poor summary for lifetimes: it depends on the margins, and it changes if you take logs of the times. The copula-based measures depend only on \(C\), so they are the same whatever the margins and whatever monotone transformation is applied to the times:

Kendall’s tau is the probability that two independent pairs are concordant (the unit that fails first on series 1 also fails first on series 2) minus the probability they are discordant: \(\tau = 4\,E[C(U_1, U_2)] - 1 \in [-1, 1]\).
Spearman’s rho is the ordinary correlation of the uniforms, \(\rho_S = \operatorname{corr}(U_1, U_2) = 12\int_0^1\!\int_0^1 C(u_1, u_2)\, du_1\, du_2 - 3\).
The tail-dependence coefficients measure clustering in the extremes:

\[\lambda_L = \lim_{q \to 0^+} P\big(U_1 \leq q \mid U_2 \leq q\big) = \lim_{q \to 0^+} \frac{C(q, q)}{q}, \qquad \lambda_U = \lim_{q \to 1^-} P\big(U_1 > q \mid U_2 > q\big).\]

\(\lambda_L > 0\) means that, given one unit is among the very earliest failures, there is a non-vanishing probability that its partner is too: early failures come in pairs. \(\lambda_U > 0\) is the same statement for the longest survivors.

Two copulas can share the same \(\tau\) and still differ greatly in the tails, which is exactly where the risk of joint failure is decided. Fit the family, not just the correlation.

The families in detail

The formulas below are standard; [Nelsen2006mv] derives them and many more.

Independence — \(C(u_1, u_2) = u_1 u_2\). No parameter; the joint distribution is the product of the margins. It is the baseline against which the others are judged: if a dependent copula does not fit clearly better, the simpler independent model may do.

Clayton —

\[ \begin{align}\begin{aligned}C(u_1, u_2) = \big(u_1^{-\theta} + u_2^{-\theta} - 1\big)^{-1/\theta}, \qquad \theta > 0,\\\tau = \frac{\theta}{\theta + 2}, \qquad \lambda_L = 2^{-1/\theta}, \qquad \lambda_U = 0 .\end{aligned}\end{align} \]

\(\theta \to 0\) gives independence and \(\theta \to \infty\) perfect positive dependence. Dependence concentrates in the lower tail: the two series are most tightly linked among early failures. Use it for common-cause early failure — a shared manufacturing defect, a shared harsh start-up, a shared contamination event — where knowing one unit failed young is strong evidence its partner will too.

Gumbel (Gumbel-Hougaard) —

\[ \begin{align}\begin{aligned}C(u_1, u_2) = \exp\!\Big(-\big[(-\ln u_1)^{\theta} + (-\ln u_2)^{\theta}\big]^{1/\theta}\Big), \qquad \theta \geq 1,\\\tau = 1 - \frac{1}{\theta}, \qquad \lambda_L = 0, \qquad \lambda_U = 2 - 2^{1/\theta} .\end{aligned}\end{align} \]

\(\theta = 1\) is independence. Dependence concentrates in the upper tail: the series are most tightly linked among the longest lives. Use it when a shared benign factor (a gentle duty cycle, a robust batch) lets pairs survive to old age together. (The Gumbel copula is unrelated to the univariate Gumbel distribution, surpyval.Gumbel.)

Frank —

\[ \begin{align}\begin{aligned}C(u_1, u_2) = -\frac{1}{\theta}\ln\!\left(1 + \frac{(e^{-\theta u_1} - 1)(e^{-\theta u_2} - 1)}{e^{-\theta} - 1}\right), \qquad \theta \neq 0,\\\tau = 1 - \frac{4}{\theta}\big[1 - D_1(\theta)\big], \qquad \lambda_L = \lambda_U = 0,\end{aligned}\end{align} \]

where \(D_1(\theta) = \frac{1}{\theta}\int_0^{\theta} \frac{s}{e^s - 1}\,ds\) is the first Debye function. \(\theta > 0\) gives positive and \(\theta < 0\) negative dependence, symmetric in the two tails and with no tail dependence. Use it for moderate, “everywhere-alike” association, or whenever the dependence may be negative (one series tends to be long when the other is short, as when two failure modes compete for the same weakness).

Gaussian —

\[ \begin{align}\begin{aligned}C(u_1, u_2) = \Phi_2\big(\Phi^{-1}(u_1), \Phi^{-1}(u_2); \rho\big), \qquad -1 < \rho < 1,\\\tau = \frac{2}{\pi}\arcsin\rho, \qquad \rho_S = \frac{6}{\pi}\arcsin\frac{\rho}{2}, \qquad \lambda_L = \lambda_U = 0,\end{aligned}\end{align} \]

where \(\Phi\) is the standard normal CDF and \(\Phi_2(\cdot, \cdot; \rho)\) the bivariate normal CDF with correlation \(\rho\). It is the dependence of a bivariate normal distribution transplanted onto arbitrary margins, so \(\rho\) has the familiar meaning of a correlation between the normal scores \(\Phi^{-1}(u_j)\). Positive and negative dependence are both allowed. Its lack of tail dependence means that, however large \(\rho\), joint extreme events become asymptotically independent — a Gaussian copula can understate the risk of joint early failure when the true dependence is Clayton-like.

Joe —

\[ \begin{align}\begin{aligned}C(u_1, u_2) = 1 - \big(\bar u_1^{\theta} + \bar u_2^{\theta} - \bar u_1^{\theta} \bar u_2^{\theta}\big)^{1/\theta}, \qquad \bar u = 1 - u, \qquad \theta \geq 1,\\\tau = 1 - \frac{2}{\theta}\, \frac{\psi(2/\theta + 1) - \psi(2)}{2/\theta - 1}, \qquad \lambda_L = 0, \qquad \lambda_U = 2 - 2^{1/\theta},\end{aligned}\end{align} \]

with \(\psi\) the digamma function (at \(\theta = 2\) the quotient is \(\psi'(2)\)). Like the Gumbel it links the long lives, but more strongly for the same Kendall’s tau: at \(\tau = 0.5\) the Joe copula has \(\lambda_U = 0.73\) (\(\theta = 2.86\)), the Gumbel 0.59. The parameterisation is that of R’s copula::joeCopula and of VineCopula’s family 6.

Ali-Mikhail-Haq (AMH) —

\[ \begin{align}\begin{aligned}C(u_1, u_2) = \frac{u_1 u_2}{1 - \theta(1 - u_1)(1 - u_2)}, \qquad -1 \leq \theta \leq 1,\\\tau = 1 - \frac{2\{\theta + (1 - \theta)^2 \ln(1 - \theta)\}} {3\theta^2}, \qquad \lambda_L = \lambda_U = 0 \ (\theta < 1).\end{aligned}\end{align} \]

Spearman’s rho has a closed form too, through the dilogarithm. The family reaches only weak dependence of either sign, \(-0.1817 \leq \tau \leq 1/3\), and is a cheap closed-form model for it (R’s copula::amhCopula). Data more strongly dependent than that drive the fit to the bound \(\theta = \pm 1\) (a valid copula; \(\theta = 1\) is the Clayton copula with \(\theta = 1\), with \(\lambda_L = 1/2\)), which it returns, as Clayton returns independence for negatively dependent data: compare the empirical Kendall’s tau with that range first.

Student-t —

\[ \begin{align}\begin{aligned}C(u_1, u_2) = T_{2,\nu}\big(T_\nu^{-1}(u_1), T_\nu^{-1}(u_2); \rho\big), \qquad -1 < \rho < 1, \quad \nu > 0,\\\tau = \frac{2}{\pi}\arcsin\rho, \qquad \lambda_L = \lambda_U = 2\,T_{\nu + 1}\!\left(-\sqrt{\frac{(\nu + 1) (1 - \rho)}{1 + \rho}}\right),\end{aligned}\end{align} \]

where \(T_\nu\) is the t CDF with \(\nu\) degrees of freedom and \(T_{2,\nu}\) the bivariate t CDF. It is the Gaussian copula with heavier joint tails: the same Kendall’s tau for a given \(\rho\), but extreme events — the earliest failures and the longest survivals — come together, the more so the smaller \(\nu\). As \(\nu \to \infty\) it becomes the Gaussian copula; data with no tail dependence send the fitted \(\nu\) there, and the fit then warns (“No finite maximum”) and recommends the Gaussian copula. The parameters are R’s copula::tCopula(param = rho, df = nu).

Rotations. Turning a copula round gives a new one: the 180-degree rotation, \(C_{180}(u_1, u_2) = u_1 + u_2 - 1 + C(1 - u_1, 1 - u_2)\), is the survival copula, the copula of \((1 - U_1, 1 - U_2)\), with the lower and upper tail dependence swapped; the 90- and 270-degree rotations, \(u_2 - C(1 - u_1, u_2)\) and \(u_1 - C(u_1, 1 - u_2)\), negate Kendall’s tau and Spearman’s rho and put the family’s tail in a corner where one series is short and the other long. SurPyval rotates the Clayton, Gumbel and Joe copulas with the rotation option of fit and from_params (R’s VineCopula convention, with the parameter kept in its own range as in pyvinecopulib); a survival Clayton links long lives, a survival Gumbel or Joe early failures.

In SurPyval, Kendall’s tau is computed in closed form for every family (Frank’s through a numerically evaluated Debye integral), and so are Spearman’s rho for the Gaussian, Independence, Frank and AMH copulas (Frank’s is \(1 - 12\{D_1(\theta) - D_2(\theta)\}/\theta\), with \(D_k\) the Debye functions) and the tail-dependence coefficients. Spearman’s rho for the Clayton, Gumbel and Joe copulas is the integral above, taken by Gauss-Legendre quadrature (accurate to about \(10^{-11}\)), and the Student-t’s is \(12\,E[U_1 U_2] - 3\) integrated through its conditional distribution (to about \(10^{-13}\)). (Before version 0.22 it was estimated from 50,000 simulated pairs, up to \(5 \times 10^{-3}\) off.)

Choosing a family

Start from the physics: is there a mechanism that makes early failures cluster (Clayton), long lives cluster (Gumbel), or a diffuse association with no special tail (Frank, Gaussian)?
Joint extremes in both tails point to the Student-t; joint long lives more strongly than a Gumbel allows point to the Joe.
Only Frank, Gaussian, Student-t and (weakly) AMH can express negative dependence. Clayton (\(\theta > 0\)), Gumbel and Joe (\(\theta \geq 1\)) are positive-only in SurPyval; fitted to negatively dependent data they are pushed to their independence boundary, so check the sign of the empirical Kendall’s tau first.
Compare fitted families by their log-likelihood or AIC (the fitted model’s log_likelihood and aic(), which use the full censored and truncated joint likelihood below) and by plotting simulated samples against the data. The how-to page shows both.

Estimation

Fitting a copula model means estimating both the marginal parameters and the copula parameter. Two strategies trade off robustness against efficiency:

IFM (Inference Functions for Margins) fits each margin independently and then estimates the copula parameter with the margins held fixed. It is fast, a poor margin cannot spoil the other one, and it is the usual default.
MLE optimises the copula parameter jointly with all marginal parameters. It is more efficient when the model is well specified, at a higher computational cost, and it is the one that stays correct when the truncation or censoring of one series depends on the other (see IFM or MLE? below).

Both maximise a likelihood, so the likelihood comes first.

The likelihood with censored data

Because the margins are ordinary survival distributions, the joint likelihood inherits survival analysis’s treatment of incomplete data: each series of a joint observation can be observed, right, left or interval censored, whatever the status of the other series, using the same codes as the univariate models (c of 0, 1, -1, 2). One idea covers all of them.

Every censored row is a rectangle. What we know about a row is that \(X_1\) lies in some set \(I_1\) and \(X_2\) in some set \(I_2\), and its likelihood is the probability of that. For a censored series the set is an interval: \((x, \infty)\) if right censored, \((-\infty, x]\) if left censored, \((x_l, x_r]\) if interval censored. On the copula scale an interval \((a_j, b_j]\) of \(X_j\) becomes \((F_j(a_j), F_j(b_j)]\) of \(U_j\), with \(F_j(-\infty) = 0\) and \(F_j(\infty) = 1\). The probability of a rectangle then follows from the joint CDF by inclusion-exclusion, as it does for any bivariate distribution:

\[P(a_1 < X_1 \leq b_1,\; a_2 < X_2 \leq b_2) = C(B_1, B_2) - C(A_1, B_2) - C(B_1, A_2) + C(A_1, A_2),\]

with \(A_j = F_j(a_j)\) and \(B_j = F_j(b_j)\). Every copula has \(C(0, \cdot) = C(\cdot, 0) = 0\), \(C(u, 1) = u\) and \(C(1, v) = v\), so an infinite end of an interval simply drops terms or turns \(C\) into a margin. For example, both series right censored gives \(1 - u_1 - u_2 + C(u_1, u_2)\) with \(u_j = F_j(x_j)\), the joint survival function.

An exact observation is a very thin rectangle. If \(X_1\) is observed at \(x_1\), shrink its interval to \((x_1, x_1 + \mathrm{d}x_1]\). The difference \(C(F_1(x_1 + \mathrm{d}x_1), \cdot) - C(F_1(x_1), \cdot)\) becomes \(\partial C/\partial u_1 \cdot f_1(x_1)\,\mathrm{d}x_1\): the difference in that argument turns into a derivative, multiplied by the marginal density. (The \(\mathrm{d}x_1\) is the same for every value of the parameters, so it is dropped, exactly as a univariate likelihood uses the density for an observed failure.)

So each series applies one operation to its argument of \(C\):

What is known about \(X_j\)	Censoring code `c`	Operation on argument \(j\) of \(C\)
observed exactly at \(x\)	`0`	differentiate at \(u = F_j(x)\), and multiply by \(f_j(x)\)
right censored, \(X_j > x\)	`1`	\(C(\cdot, 1) - C(\cdot, F_j(x))\)
left censored, \(X_j \leq x\)	`-1`	evaluate at \(F_j(x)\)
interval censored, \(x_l < X_j \leq x_r\)	`2`	\(C(\cdot, F_j(x_r)) - C(\cdot, F_j(x_l))\)

Applying the operations for both series gives the row’s likelihood; with four codes per series there are sixteen combinations, all built from the four functions \(C\), \(\partial C/\partial u_1\), \(\partial C/\partial u_2\) and the copula density \(c\). Some examples, with \(u_j = F_j(x_j)\):

\[\begin{split}\text{both observed:}\quad & c(u_1, u_2)\, f_1(x_1)\, f_2(x_2) \\ \text{1 observed, 2 right censored:}\quad & f_1(x_1)\Big[1 - \frac{\partial C}{\partial u_1}(u_1, u_2)\Big] \\ \text{1 observed, 2 left censored:}\quad & f_1(x_1)\,\frac{\partial C}{\partial u_1}(u_1, u_2) \\ \text{both right censored:}\quad & 1 - u_1 - u_2 + C(u_1, u_2) \\ \text{both left censored:}\quad & C(u_1, u_2) \\ \text{1 left censored, 2 interval censored:}\quad & C\big(u_1, F_2(x_{r,2})\big) - C\big(u_1, F_2(x_{l,2})\big)\end{split}\]

The second and third lines read naturally: the density of seeing series 1 fail at \(x_1\), times the conditional probability that series 2 had not (or had) failed by \(x_2\) given that. This is the h-function at work, and it is how a censored partner still carries information about the dependence: a row where series 1 failed early and series 2 was still running at a late censoring time is evidence against strong lower-tail dependence.

Counts n weight each row’s log-likelihood, so a row standing for five identical units counts five times.

For the likelihood to be valid the censoring must not depend on the unobserved lifetimes. For a single series this is the usual independent-censoring assumption. Jointly it allows more: series 2 may be censored at a time that depends on what happened to series 1 (for example, the unit is retired a fixed time after series 1 fails), because \(x_1\) is part of the observed row. As explained below, IFM does not allow this.

Truncation: conditioning on being seen

Truncation is different from censoring. Censoring hides part of an observation that is in the sample; truncation decides which rows are in the sample at all. If series \(j\) could only be seen inside the window \((t_{l,j}, t_{r,j})\), a row appears in the data only if both series landed in their windows. Each row’s likelihood (from the rectangle rule above) is therefore divided by the probability of that event, which is the copula mass of the truncation rectangle, by the same inclusion-exclusion:

\[P(\text{row observable}) = C(r_1, r_2) - C(l_1, r_2) - C(r_1, l_2) + C(l_1, l_2), \qquad l_j = F_j(t_{l,j}),\; r_j = F_j(t_{r,j}).\]

In SurPyval the window is given per row and per series (t of shape (N, 2, 2)), with \(\pm\infty\) for “no limit” (\(l_j = 0\), \(r_j = 1\)).

A common special case is truncation of one series only, say a burn-in of length \(b\): rows are seen only if \(X_1 > b\). The divisor is then \(C(1, 1) - C(F_1(b), 1) = 1 - F_1(b)\), a marginal probability that does not involve the copula parameter. It is tempting to conclude that the truncation can be dealt with inside margin 1 alone. It cannot, because the selection also changes which values of \(X_2\) are seen. The density of \(X_2\) among the rows that pass the burn-in is

\[f_2(x_2)\,\frac{P(X_1 > b \mid X_2 = x_2)}{P(X_1 > b)} = f_2(x_2)\, \frac{1 - \dfrac{\partial C}{\partial u_2}\big(F_1(b), F_2(x_2)\big)} {1 - F_1(b)} .\]

Under independence the fraction is one and nothing changes. Under positive dependence a long \(X_2\) makes passing the burn-in more likely, so the observed \(X_2\) values are shifted towards long lives, although series 2 was never truncated itself. The joint likelihood accounts for this automatically, through the copula density \(c(u_1, u_2)\) that pairs each \(x_2\) with its \(x_1\). A fit of margin 2 alone cannot.

How SurPyval fits

Every family has the same fit(x, c=None, n=None, t=None, margins=None, how="IFM", xl=None, xr=None, init=None); margins (one per series) is required and exactly two series are supported. With how="IFM" [JoeXu1996mv]:

each margin is fitted by the usual univariate maximum likelihood to its own series, honouring that series’ censoring codes (interval-censored entries through xl/xr), the row counts n and that series’ own truncation window; margins passed as already-fitted models are used as they are (under how="MLE" they supply the starting values and are re-estimated with the configuration they were fitted with: the same offset, limited-failure or zero-inflated option, and any fixed parameters kept at their values). A margin can also be non-parametric (the class surpyval.KaplanMeier, say, or a fitted non-parametric model), which gives the semi-parametric estimator of [GenestGhoudiRivest1995mv]: the margin’s step CDF, rescaled by \(N/(N+1)\) so the largest values stay inside the unit square, replaces \(F_j\), and the margin contributes no density term (a step function has none, and it would not depend on \(\theta\)). The copula parameter then rests on no assumption about the margins’ shapes; the likelihood is a pseudo-likelihood, so its value and the criteria compare copula families with the same margins only, and k counts the copula parameters (plus any parametric margin’s);
with the margins fixed, the copula parameter is chosen to maximise the joint log-likelihood above (censoring operations and truncation divisor included). The search is a Nelder-Mead search on an unconstrained transformation of the parameter (for the Gaussian copula, \(\rho = \tanh(\cdot)\)). It starts from the value that matches the empirical Kendall’s tau of the rows where both series are observed (or from near-independence if fewer than three such rows exist), unless a starting value is passed as init.

With how="MLE" the IFM solution is the starting point for a joint Nelder-Mead search over the copula parameter and every margin parameter, maximising the full likelihood. This is slower (about ten times the IFM time in the examples on the how-to page) and, like any joint fit, it lets a misspecified margin pull on the copula parameter and on the other margin.

IFM or MLE?

IFM’s first stage fits each series as if it were the only one. That is correct when the data of each series, looked at on its own, is an honest sample of that margin: every series censored independently of its own lifetime, and no rows selected by what happened to the other series. The second stage then needs only the correct margins, so IFM is consistent, and it loses little efficiency when the dependence is moderate [JoeXu1996mv]. Two common designs break the first stage:

Truncation of a partner series. As shown above, a burn-in on series 1 shifts the observed values of series 2. Margin 2 is fitted to the shifted sample as if it were the population. The copula stage then has to explain the data with a distorted margin, and it does so with the wrong amount of dependence (too little, in the how-to page’s example). The margin of the truncated series itself is fine: its own window is exactly the selection it is subject to.
Censoring that depends on the partner. If series 2 is censored at a time set by series 1 (a unit retired soon after its first part fails, a patient whose follow-up for one complication stops after the other), the censoring of series 2 is related to \(X_2\) through the dependence. Viewed one series at a time this is informative censoring, and margin 2 is biased. The joint likelihood is valid, because the censoring depends only on the observed \(x_1\).

how="MLE" maximises the joint likelihood, in which the margins, the copula, the censoring operations and the truncation divisor all appear together, so it is right in both cases. It costs more time and is more exposed to a badly chosen margin. A practical check: fit both. If they agree, the cheaper IFM fit is fine; if they differ, look for truncation or censoring of one series that is driven by the other. The how-to page shows both failures of IFM and their MLE correction.

Some further points worth knowing:

The copula parameter is estimated on the scale \(u_j = F_j(x_j)\), so a poorly chosen margin distorts it. Check the margins with the univariate tools first (see Parametric SurPyval Modelling).
The fitted model reports the point estimates, without standard errors. Its log_likelihood (neg_ll()) is the full joint log-likelihood above at those estimates, and aic()/bic() count as parameters the copula’s and those of every margin the fit estimated (a margin passed already fitted is not re-estimated under IFM, so it does not count); bic() uses the number of rows in which at least one series failed (was not right-censored), or of all rows when none did – the sample-size rule of every SurPyval model (see Comparing models: information criteria). A from_params model has no likelihood and raises a ValueError.
Margin probabilities are kept a tiny distance (\(10^{-10}\)) inside \((0, 1)\) to keep the Archimedean formulas finite.
The Frank, Clayton, Gumbel and Joe formulas are evaluated in log space, from terms that neither cancel nor overflow, so they stay accurate however strong the dependence (Frank’s and the Student-t’s samplers also invert their h-functions in closed form). The Gaussian and Student-t formulas are written in \(1 - |\rho|\) and in the difference of the two quantiles, so they are accurate for any \(|\rho| < 1\) (the Gaussian CDF is scipy’s, Genz’s algorithm, within \(10^{-14}\) of a 40-digit integration up to \(\rho = 1 - 2^{-52}\)); \(\rho = \pm 1\), the comonotone or countermonotone copula, is not a member of either family and from_params refuses it. A fit whose \(\rho\) runs to \(\pm 1\) warns that the likelihood has no finite maximum.
The Student-t copula’s CDF (needed for rows censored in both series and for truncation) is the integral of its closed-form h-function, taken by tanh-sinh quadrature: within \(10^{-11}\) of Genz’s exact algorithm (R’s mvtnorm, integer \(\nu\)) and \(10^{-15}\) of a 30-digit integration at non-integer \(\nu\). scipy’s multivariate_t.cdf is a randomised quasi-Monte Carlo integration (about \(10^{-4}\) off, different on every call), which an optimiser cannot use.
Only bivariate models are supported; more than two series raise a NotImplementedError.

For worked examples — fitting a copula, checking each censoring pattern’s likelihood contribution, IFM against MLE under truncation and dependent censoring, querying the joint distribution and dependence measures, simulating correlated lifetimes and defining a new copula family — see the Multivariate Modelling with SurPyval page.

References

[Sklar1959mv]

Sklar, A. (1959). Fonctions de répartition à n dimensions et leurs marges. Publications de l’Institut de Statistique de l’Université de Paris, 8, 229–231.

[JoeXu1996mv] (1,2)

Joe, H., & Xu, J. J. (1996). The estimation method of inference functions for margins for multivariate models. Technical Report 166, Department of Statistics, University of British Columbia.

[Nelsen2006mv]

Nelsen, R. B. (2006). An Introduction to Copulas (2nd ed.). Springer. The standard reference for the families, dependence measures and tail dependence used on this page.

[GenestGhoudiRivest1995mv]

Genest, C., Ghoudi, K. and Rivest, L.-P. (1995). A semiparametric estimation procedure of dependence parameters in multivariate families of distributions. Biometrika, 82(3), 543-552.