Case study: copula dependence between temperature and ozone

1. Pseudo-observations, the dependence stripped of the margins

The copula lives on the unit square, and the data enter it as pseudo-observations: each value is replaced by its rank divided by n + 1, which maps the marginal to a uniform and leaves only the dependence. gpum_copula() does this internally; here is the picture it works from.

u <- rank(x) / (length(x) + 1)
v <- rank(y) / (length(y) + 1)
plot(u, v, pch = 16, col = grDevices::adjustcolor("black", 0.5),
     xlab = "u = rank(temperature)", ylab = "v = rank(ozone)",
     main = "Pseudo-observations on the unit square")

The points look denser toward the upper-right corner than the lower-left, which would be upper-tail dependence. Whether that apparent asymmetry is real dependence or an artefact of the eye is what the formal comparison settles.

2. Fitting with automatic family selection

Rather than assume a family, family = "auto" fits all four candidates and selects the one preferred by predictive comparison (WAIC on the copula density, which is a proper density on the unit square). The choice of family is itself a hypothesis about tail behaviour, so letting the data choose is part of the answer.

fit <- gpum_copula(data.frame(temp = x, ozone = y), family = "auto",
                   n_iter = 12000, n_chains = 8, seed = 1, backend = "cpu")
fit
#> <gpum_copula>
#>   selected family : frank (by WAIC)
#>   theta           : 7.038 [5.507, 8.637]
#>   Kendall's tau   : 0.561 [0.487, 0.625]  (empirical 0.586)
#>   tail dependence : lower 0.000, upper 0.000
#>   diagnostic      : R-hat 1.000 (converged); empirical tau inside the 95% band (adequate fit)
#>   family comparison:
#>    family criterion  value  delta converged
#>     frank      WAIC -95.79  0.000      TRUE
#>  gaussian      WAIC -88.78  7.002      TRUE
#>    gumbel      WAIC -87.66  8.121      TRUE
#>   clayton      WAIC -60.69 35.094      TRUE

The printed object is the technical diagnostic: the selected family, the dependence parameter and Kendall’s tau with their credible intervals, the tail-dependence coefficients, and a verdict line that reports convergence (R-hat) and goodness of fit (whether the empirical Kendall tau falls inside the model’s posterior band). The comparison table ranks the four families by WAIC.

The data select Frank, not Gumbel. Frank is a symmetric, tail-independent copula: strong overall dependence (Kendall’s tau near 0.56) with no concentration in either tail. The upper-tail hint the scatter suggested did not survive the comparison; the Gumbel family, which would impose upper-tail dependence, ranks below Frank on WAIC. The fan in the raw scatter was heteroscedasticity of the ozone margin, higher variance at high temperature, not a coupling of the extremes. This is the discipline the copula enforces: the dependence structure is read after the marginal shapes are removed by the ranks, and the eye’s reading of the raw scatter, which mixes the two, is overruled.

3. Reading the dependence

Kendall’s tau is the family-free summary, comparable across copulas; the tail-dependence coefficients are the extreme-co-movement probabilities.

kendall_tau(fit)
#>      mean      q2.5    median     q97.5 
#> 0.5614716 0.4871316 0.5633997 0.6249870
tail_dependence(fit)
#> $lower
#>   mean   q2.5 median  q97.5 
#>      0      0      0      0 
#> 
#> $upper
#>   mean   q2.5 median  q97.5 
#>      0      0      0      0

Both tail-dependence coefficients are zero, the signature of the Frank family: temperature and ozone are strongly associated through the body of the distribution, but their extremes are asymptotically independent. Knowing that today is among the hottest on record does not, in the tail limit, make a record ozone reading more likely than the overall association already implies. That is a materially different risk statement from the upper-tail dependence the scatter suggested, and it is the copula, not the eye, that is entitled to make it.

4. Visual inspection

The plot() method is the visual companion of the printed diagnostic: the fitted density contours over the pseudo-observations, the posteriors of the parameter and of Kendall’s tau against the empirical value, and the observed-against-generated check on the unit square.

plot(fit)

The Frank contours are symmetric about the anti-diagonal, no corner favoured, the posterior of tau sits on the empirical mark, and the generated blue cloud overlays the observed black one. That overlay is the confirmation that a symmetric, tail-independent copula reproduces the observed dependence; had the upper tail been genuinely coupled, the generated cloud would fall short of the observed points in the upper-right corner.

5. Against the classical estimator

The copula package fits the same families by maximum pseudo-likelihood. The Bernstein-von Mises theorem says the Bayesian posterior mean and the maximum-likelihood estimate must agree for a regular model at this sample size; the check confirms the fit is not an artefact of the sampler.

if (requireNamespace("copula", quietly = TRUE)) {
  cop <- switch(fit$family,
    gumbel = copula::gumbelCopula(), clayton = copula::claytonCopula(),
    frank = copula::frankCopula(), gaussian = copula::normalCopula())
  ml <- copula::fitCopula(cop, cbind(u, v), method = "mpl")
  cat(sprintf("family %s: posterior mean %.3f vs pseudo-MLE %.3f\n",
              fit$family, fit$param["mean"], as.numeric(coef(ml))))
  cat(sprintf("Kendall tau: posterior %.3f vs empirical %.3f\n",
              fit$tau["mean"], fit$tau_empirical))
}
#> family frank: posterior mean 7.038 vs pseudo-MLE 7.131
#> Kendall tau: posterior 0.561 vs empirical 0.586

6. Forcing the eye’s hypothesis, and seeing it rejected

The scatter argued for upper-tail dependence, which is the Gumbel family. Fitting Gumbel by force shows what imposing that hypothesis would claim, and the comparison of step 2 already recorded that it fits worse. Read the two side by side: Gumbel would assert a positive upper tail-dependence coefficient the data do not support, at a worse WAIC, while agreeing with Frank on the overall strength.

fit_gumbel <- gpum_copula(data.frame(temp = x, ozone = y),
                          family = "gumbel", n_iter = 8000, n_chains = 8,
                          seed = 1, backend = "cpu")
rbind(
  selected_frank = c(tau = fit$tau["mean"],
                     tail_upper = fit$tail_upper["mean"]),
  forced_gumbel = c(tau = fit_gumbel$tau["mean"],
                    tail_upper = fit_gumbel$tail_upper["mean"])
)
#>                 tau.mean tail_upper.mean
#> selected_frank 0.5614716       0.0000000
#> forced_gumbel  0.5295618       0.6140167

The two families agree on the overall strength of dependence and disagree on where it lives: Frank places none in the tails, Gumbel would place a substantial share in the upper tail. Choosing the family is choosing a hypothesis about the extremes, and the copula workflow makes that choice an inference the data adjudicate rather than an assumption the analyst imports.

7. Where this stops

The workflow here is bivariate. Three or more variables need the vine decomposition of a later release, which factors a high-dimensional dependence into a cascade of bivariate copulas, each of the kind fit here. The marginal distributions were handled non-parametrically through the ranks; the companion release on marginal auto-selection fits an explicit distribution to each column so the two halves of Sklar’s theorem, margins and copula, are both modelled. And a caveat of method: heavy ties in either variable blur the pseudo-observations and the copula’s identifiability (Genest and Neslehova 2007); the ozone and temperature readings here are near-continuous, so the ranks are clean.

References

• Genest, C., Neslehova, J. (2007). A primer on copulas for count data. ASTIN Bulletin, 37(2), 475-515.
• Nelsen, R. B. (2006). An Introduction to Copulas, 2nd ed. Springer.
• Sklar, A. (1959). Fonctions de repartition a n dimensions et leurs marges. Publications de l’Institut de Statistique de l’Universite de Paris, 8, 229-231.