Wavelet coherence and phase

Do two ecological series cycle together, and which one leads? Build wavelet coherence by hand in R and learn why the answer needs smoothing.

You have two series: hare counts and lynx counts, host and parasite, rainfall and recruitment. They both wobble. The questions are whether they wobble together, at which timescale, during which years, and who moves first.

The cross-correlation function answers a blunter version of this: one number per lag, averaged over the whole record, assuming the relationship never changed. Wavelet coherence is the version that keeps the time axis. It is also the version with the most ways to fool you, and the first of them is built into the algebra.

The transform is the same one from the first post in this series, folded away here.

The cross-wavelet transform, and why its modulus is empty

Multiply one transform by the conjugate of the other and you get the cross-wavelet transform, \(W_{xy} = W_x W_y^*\). Its argument is the phase difference between the two series at that time and scale, which is genuinely useful. Its modulus is called cross-wavelet power, gets plotted constantly, and contains no information about the relationship whatsoever:

set.seed(286)
n <- 300
x1 <- as.numeric(arima.sim(list(ar = 0.8), n))
x2 <- as.numeric(arima.sim(list(ar = 0.8), n))   # completely independent
o1 <- cwt(x1); o2 <- cwt(x2)
Wxy <- o1$W * Conj(o2$W)
lhs <- Mod(Wxy)^2
rhs <- Mod(o1$W)^2 * Mod(o2$W)^2
max_rel_diff <- max(abs(lhs - rhs) / rhs)
c(max_relative_difference = signif(max_rel_diff, 3))

max_relative_difference 
               9.61e-16 

\(|W_{xy}|^2 = |W_x|^2 |W_y|^2\), to 9.6e-16. That is not an empirical finding, it is an identity: the modulus of a product is the product of the moduli. The cross-power map is the two power maps multiplied together. If both species happen to have a strong cycle at period 8, the cross-power map lights up at period 8 whether or not the species have ever met.

usable <- outer(o1$period, o1$coi, "<=")
lhs_ok <- replace(lhs, !usable, -Inf)
i_max <- which(lhs_ok == max(lhs_ok), arr.ind = TRUE)[1, ]
c(largest_cross_power = round(max(lhs_ok), 1),
  at_period = round(o1$period[i_max[1]], 2), at_year = as.numeric(i_max[2]))

largest_cross_power           at_period             at_year 
            4529.00               46.75               66.00 

The largest cross-power on the map of two unrelated series is 4529, at period 46.8, sitting where both series happened to have a large excursion at the same time. Maraun and Kurths (2004) put this plainly: the wavelet cross spectrum is not suitable for significance testing the relationship between two processes. Use coherence.

Coherence, and why it needs smoothing

Coherence is cross-power normalised by the two individual powers, which is exactly the structure of a correlation coefficient:

\[R^2 = \frac{|S(W_{xy}/s)|^2}{S(|W_x|^2/s) \cdot S(|W_y|^2/s)}\]

The \(S\) is a smoothing operator, and without it the whole thing collapses. Watch:

R2_raw <- Mod(Wxy)^2 / (Mod(o1$W)^2 * Mod(o2$W)^2)
c(max_deviation_from_one = signif(max(abs(R2_raw - 1)), 3),
  min_coherence = min(R2_raw), max_coherence = max(R2_raw))

max_deviation_from_one          min_coherence          max_coherence 
              8.88e-16               1.00e+00               1.00e+00 

Unsmoothed coherence is identically 1, everywhere, for any two series, to 8.9e-16. It is the previous identity divided by itself. Two independent AR(1) processes are perfectly coherent at every time and every scale, which tells you the quantity is measuring nothing.

Smoothing is not a cosmetic step, then. It is what creates the statistic. Following Torrence and Webster (1999), smooth along time with a Gaussian whose width grows with the scale, then along scale with a short boxcar:

smooth_wt <- function(M, scale, dt, dj) {
  n <- ncol(M)
  npad <- 2^ceiling(log2(n))
  kk <- seq_len(floor(npad / 2)) * (2 * pi / npad)
  k2 <- c(0, kk, -rev(kk[seq_len(floor((npad - 1) / 2))]))^2
  S <- matrix(0 + 0i, nrow(M), n)
  for (a in seq_len(nrow(M))) {                       # gaussian in time, width proportional to scale
    z <- c(M[a, ], rep(0, npad - n))
    S[a, ] <- (fft(exp(-0.5 * (scale[a] / dt)^2 * k2) * fft(z), inverse = TRUE) / npad)[seq_len(n)]
  }
  m <- 0.6 / (dj * 2)                                 # boxcar of 0.6 decorrelation length (Morlet)
  w <- c(m %% 1, rep(1, 2 * round(m) - 1), m %% 1); w <- w / sum(w)
  h <- (length(w) - 1) / 2
  out <- matrix(0 + 0i, nrow(M), n)
  for (a in seq_len(nrow(M))) {
    idx <- pmin(pmax(seq(a - h, a + h), 1), nrow(M))
    out[a, ] <- colSums(w * S[idx, , drop = FALSE])
  }
  out
}

wcoh <- function(o1, o2) {
  s <- o1$scale
  Sxy <- smooth_wt(o1$W * Conj(o2$W) / s, s, o1$dt, o1$dj)
  Sxx <- Re(smooth_wt((Mod(o1$W)^2 + 0i) / s, s, o1$dt, o1$dj))
  Syy <- Re(smooth_wt((Mod(o2$W)^2 + 0i) / s, s, o1$dt, o1$dj))
  list(R2 = Mod(Sxy)^2 / (Sxx * Syy), phase = Arg(Sxy))
}

Two details that are easy to get wrong. The division by s inside every smoothed term is the scale rectification of Liu et al. (2007), and it has to be inside the smoothing, not outside. And the boxcar width 0.6 / dj is a property of the Morlet wavelet’s decorrelation length in scale, so changing w0 or dj changes it.

cc <- wcoh(o1, o2)
c(mean_coherence_independent_pair = round(mean(cc$R2[usable]), 3),
  max_coherence = round(max(cc$R2[usable]), 3),
  in_range = all(cc$R2 >= 0 & cc$R2 <= 1 + 1e-9))

mean_coherence_independent_pair                   max_coherence 
                          0.355                           0.966 
                       in_range 
                          1.000 

Now the same two independent series give a mean coherence of 0.355 rather than 1. The quantity exists.

A pair that really is coupled

Prey drives predator with a quarter-cycle delay: the classic signature, and a good test of whether the machinery recovers a lag we planted.

set.seed(1286)
tt <- 1:n
per <- 10
true_lag <- 2.5                                       # a quarter of the period
prey <- sin(2 * pi * tt / per) + rnorm(n, 0, 0.4)
pred <- sin(2 * pi * (tt - true_lag) / per) + rnorm(n, 0, 0.4)
q1 <- cwt(prey); q2 <- cwt(pred)
cq <- wcoh(q1, q2)
band <- which.min(abs(q1$period - per))
sel <- q1$period[band] <= q1$coi
coh_band <- mean(cq$R2[band, sel])
ph <- mean(cq$phase[band, sel])
lag_hat <- ph / (2 * pi) * q1$period[band]
c(band_period = round(q1$period[band], 2), mean_coherence = round(coh_band, 3),
  mean_phase_rad = round(ph, 4), implied_lag = round(lag_hat, 3), true_lag = true_lag)

   band_period mean_coherence mean_phase_rad    implied_lag       true_lag 
        10.170          0.947          1.544          2.500          2.500 

Coherence 0.947 at the period-10.17 band, and the phase implies a lag of 2.50 steps against a planted 2.5. The phase converts to a lag by \(\text{lag} = \phi \cdot \text{period} / 2\pi\), and it does not have to be an integer. Compare what the cross-correlation gets:

cf <- ccf(prey, pred, lag.max = 12, plot = FALSE)
c(ccf_peak_lag = cf$lag[which.max(cf$acf)], ccf_peak_value = round(max(cf$acf), 3))

  ccf_peak_lag ccf_peak_value 
        -3.000          0.714 

The cross-correlation lands on -3, because it can only report whole steps. (R’s sign convention for ccf(a, b) is that a negative lag means a leads b, which catches everyone once.) The wavelet gets 2.50, and it gets it separately for every year, which is the part the cross-correlation cannot do at all.

Figure 1: The prey and predator series (top) and their wavelet coherence (bottom). Arrows show the phase difference at the period-10 band: pointing right and down means the prey leads the predator by a quarter of a cycle. The pale wedges are the cone of influence.

What counts as high coherence

0.947 sounds impressive. Is 0.6? Simulate the null: two independent AR(1) series, two hundred times, and look at the whole distribution of coherence you get from nothing.

mc_coh <- function(B, n, phi, seed) {
  set.seed(seed)
  out <- numeric(0)
  for (b in 1:B) {
    z1 <- as.numeric(arima.sim(list(ar = phi), n))
    z2 <- as.numeric(arima.sim(list(ar = phi), n))
    p1 <- cwt(z1); p2 <- cwt(z2)
    r <- wcoh(p1, p2)$R2
    out <- c(out, r[outer(p1$period, p1$coi, "<=")])
  }
  out
}
nullR2 <- mc_coh(200, 200, 0.7, 2286)
c(null_mean = round(mean(nullR2), 3), null_median = round(median(nullR2), 3),
  null_95th = round(quantile(nullR2, 0.95), 3), null_99th = round(quantile(nullR2, 0.99), 3),
  share_of_null_above_0.5 = round(mean(nullR2 > 0.5), 3))

              null_mean             null_median           null_95th.95% 
                  0.337                   0.309                   0.742 
          null_99th.99% share_of_null_above_0.5 
                  0.859                   0.250 

Two series with no connection at all produce a median coherence of 0.31, and 25% of the null map exceeds 0.5. The 95th percentile is 0.74. So a patch at coherence 0.7 is not evidence of anything; it is the twentieth percentile of noise. This is why coherence figures need a simulated level, not a rule of thumb, and it is why Maraun and Kurths (2004) concluded that the widely reported ENSO to North Atlantic Oscillation coherence is an artefact for most of the twentieth century.

Figure 2: Coherence values from 200 pairs of independent red noise series, restricted to the usable part of the map. High coherence is the default output of the method, not a discovery.

Three ways the phase lies

The wrap. Phase is an angle, so it identifies the lag only up to a whole period. Plant a lag of 7 in a period-10 cycle:

set.seed(2287)
big_lag <- 7
w1 <- sin(2 * pi * tt / per) + rnorm(n, 0, 0.4)
w2 <- sin(2 * pi * (tt - big_lag) / per) + rnorm(n, 0, 0.4)
v1 <- cwt(w1); v2 <- cwt(w2)
cw <- wcoh(v1, v2)
selw <- v1$period[band] <= v1$coi
phw <- mean(cw$phase[band, selw])
c(true_lag = big_lag, mean_phase_rad = round(phw, 3),
  apparent_lag = round(phw / (2 * pi) * v1$period[band], 3),
  coherence = round(mean(cw$R2[band, selw]), 3))

      true_lag mean_phase_rad   apparent_lag      coherence 
         7.000         -1.901         -3.079          0.976 

The coherence is 0.98, so the relationship is unmistakable, and the phase says the lag is -3.08. Negative. The follower appears to lead. A lag of 7 and a lag of -3 are the same angle in a period-10 cycle, and no amount of data will separate them. Phase gives you the lag modulo the period, and choosing between the candidates is ecology, not statistics.

The common driver. Two species, no interaction of any kind, both responding to the same climate index:

set.seed(3286)
driver <- sin(2 * pi * tt / 12)
sp_a <- 1.2 * driver + rnorm(n, 0, 0.5)
sp_b <- 1.2 * driver + rnorm(n, 0, 0.5)   # sp_b never sees sp_a
d1 <- cwt(sp_a); d2 <- cwt(sp_b)
cd <- wcoh(d1, d2)
bd <- which.min(abs(d1$period - 12))
seld <- d1$period[bd] <= d1$coi
c(coherence_at_period_12 = round(mean(cd$R2[bd, seld]), 3),
  mean_phase_rad = round(mean(cd$phase[bd, seld]), 4))

coherence_at_period_12         mean_phase_rad 
                0.9700                -0.0539 

Coherence 0.970, phase -0.054 radians: in lockstep, indistinguishable from a tight mutualism. Coherence measures association in the time-frequency plane. It has no more to say about mechanism than a correlation coefficient does, and the arrows make it look like it does.

The edge. Coherence is inflated near the ends of the record for the same reason power is attenuated there, and the smoothing spreads the contamination further than the cone suggests:

c(mean_coherence_usable = round(mean(cc$R2[usable]), 3),
  mean_coherence_in_the_cone = round(mean(cc$R2[!usable]), 3))

     mean_coherence_usable mean_coherence_in_the_cone 
                     0.355                      0.483 

On the independent pair, the mean coherence inside the cone is 0.483 against 0.355 in the usable region. The corners of a coherence map are the least trustworthy part and the most colourful.

Honest limits

Everything above uses a pointwise level, and coherence maps have the same patchiness problem as power maps, so the areawise thinking applies here too and matters more, because the smoothing makes the patches larger by construction. The AR(1) surrogate null assumes each series is separately AR(1), which is a modelling choice about the background and not a fact about voles. And coherence, high and beautifully arrowed, still cannot tell a shared driver from an interaction: for that you need an experiment, or a method built for the question.

Where to go next

Four posts in, you can build the map, test it, and pair it with a second series. The last one in the series is about the ways the whole apparatus can be wrong in ways the map does not show: a bias that decides which cycle looks bigger, a single event that impersonates a cycle, and a resolution knob you chose without noticing.

Related tutorials

• Wavelet analysis of population cycles
• Wavelet significance and the red noise null
• Checking a wavelet analysis
• Convergent cross mapping

References

Grinsted, Moore and Jevrejeva 2004 Nonlinear Processes in Geophysics 11(5/6):561-566 (10.5194/npg-11-561-2004)

Torrence and Webster 1999 Journal of Climate 12(8):2679-2690

Maraun and Kurths 2004 Nonlinear Processes in Geophysics 11(4):505-514 (10.5194/npg-11-505-2004)

Liu, Liang and Weisberg 2007 Journal of Atmospheric and Oceanic Technology 24(12):2093-2102 (10.1175/2007JTECHO511.1)

Cazelles, Chavez, Berteaux, Menard, Vik, Jenouvrier and Stenseth 2008 Oecologia 156(2):287-304 (10.1007/s00442-008-0993-2)

Grenfell, Bjornstad and Kappey 2001 Nature 414(6865):716-723 (10.1038/414716a)