Wavelet significance and the red noise null

A scalogram is thousands of tests at once. Why a white noise null flags a quarter of the map, and how to test blobs rather than pixels in R.

Every scalogram has blobs. The transform is a filter bank, the data are noisy, and some cells will come out darker than others no matter what the series is. So the question that follows the map is: which of these blobs would surprise me if the population had no cycle at all?

That question has three separate traps in it, and the standard workflow steps on all three.

The theory, in one line

Torrence and Compo (1998) showed that if the series is Gaussian with background power spectrum \(P_k\), then at each scale the wavelet power is distributed as

\[\frac{|W_n(s)|^2}{\sigma^2} \Rightarrow \frac{1}{2} P_k \chi^2_2\]

so the 95% level at each period is var(x) * P_k * qchisq(0.95, 2) / 2. Two degrees of freedom, because the Morlet coefficient is complex and its real and imaginary parts each contribute one.

Everything hinges on \(P_k\): the spectrum of the background you are testing against. Get that wrong and the test does not degrade gracefully, it inverts.

The transform itself is the same twenty lines as in the previous post, folded away here. What is new is everything after it.

Trap one: the background is not white

Here is a population index with no cycle whatsoever. It is an AR(1) process: this year’s deviation is 0.7 of last year’s plus a shock. That is not an exotic model, it is what you get from any population with overlapping generations, carry-over of condition, or a slow-moving driver.

set.seed(285)
n <- 256
y <- as.numeric(arima.sim(list(ar = 0.7), n = n))
acf1 <- as.numeric(acf(y, lag.max = 1, plot = FALSE)$acf[2])
o <- cwt(y)
P <- Mod(o$W)^2
usable <- outer(o$period, o$coi, "<=")
c(true_phi = 0.7, sample_lag1_acf = round(acf1, 3),
  map_cells = length(P), usable_cells = sum(usable))

       true_phi sample_lag1_acf       map_cells    usable_cells 
          0.700           0.652       18176.000       10591.000 

Now test it two ways. The theoretical background spectrum of an AR(1) with parameter \(\alpha\) is a standard formula; setting \(\alpha = 0\) gives white noise, a flat line.

P_theor <- function(period, dt, alpha) {
  fr <- dt / period
  (1 - alpha^2) / (1 + alpha^2 - 2 * alpha * cos(2 * pi * fr))
}
sig95 <- function(o, x, alpha) var(x) * P_theor(o$period, o$dt, alpha) * qchisq(0.95, 2) / 2

lev_white <- sig95(o, y, 0)
lev_red <- sig95(o, y, acf1)
frac_white <- mean((P > lev_white)[usable])
frac_red <- mean((P > lev_red)[usable])
c(flagged_under_white_null = round(frac_white, 3),
  flagged_under_red_null = round(frac_red, 3))

flagged_under_white_null   flagged_under_red_null 
                   0.240                    0.046 

24% of the map is “significant at 95%” against a white background. There is no cycle in this series. A test that rejects 24% of the time when the null is true is not a 5% test, it is a machine for finding cycles in anything.

Against the red background the figure is 4.6%, which is what a 95% test is supposed to give.

Figure 1: Wavelet power for pure AR(1) noise with no cycle in it. The upper panel outlines everything the white noise null calls significant; the lower panel outlines what the red noise null calls significant. Same data, same transform, same 95% level.

One series proves nothing. Repeat it two hundred times.

mc <- function(B, n, phi, seed) {
  set.seed(seed)
  out <- matrix(NA_real_, B, 2)
  for (b in 1:B) {
    z <- as.numeric(arima.sim(list(ar = phi), n = n))
    oz <- cwt(z); Pz <- Mod(oz$W)^2
    ok <- outer(oz$period, oz$coi, "<=")
    a <- as.numeric(acf(z, lag.max = 1, plot = FALSE)$acf[2])
    out[b, ] <- c(mean((Pz > sig95(oz, z, 0))[ok]), mean((Pz > sig95(oz, z, a))[ok]))
  }
  colnames(out) <- c("white", "red")
  out
}
m <- mc(200, n, 0.7, 1285)
rate_white <- mean(m[, "white"])
rate_red <- mean(m[, "red"])
c(white_null_rejection_rate = round(rate_white, 3),
  red_null_rejection_rate = round(rate_red, 3),
  red_sd_across_series = round(sd(m[, "red"]), 3))

white_null_rejection_rate   red_null_rejection_rate      red_sd_across_series 
                    0.201                     0.051                     0.015 

Over 200 cycle-free AR(1) series, the white null flags 20.1% of the map on average and the red null flags 5.1%. The red-noise null is not a refinement. It is the difference between a test and a rubber stamp.

Why the gap is so large is visible in the levels themselves:

i4 <- which.min(abs(o$period - 4)); i32 <- which.min(abs(o$period - 32))
c(red_level_at_period_4 = round(lev_red[i4], 2),
  red_level_at_period_32 = round(lev_red[i32], 2),
  level_ratio = round(lev_red[i32] / lev_red[i4], 1))

 red_level_at_period_4 red_level_at_period_32            level_ratio 
                  2.10                  19.79                   9.40 

An AR(1) background puts 9.4 times more expected power at period 32 than at period 4. The white null sets one flat bar across all periods, so the whole long-period end of the map clears it automatically. Long-period blobs are exactly the ones people get excited about.

Figure 2: The 95% levels under the two backgrounds, against the observed time-averaged power of the cycle-free series. The white level is flat by construction; the red level rises with period and tracks what the data actually do.

Trap two: the edge eats the power

The cone of influence is not only a caveat about long periods. It attenuates power near the ends of the record, which means the same cycle is harder to detect at year 5 than at year 128, and a test that ignores the cone is less powerful exactly where records begin and end.

s16 <- 20 + 4 * sin(2 * pi * (1:n) / 16)   # a cycle of constant amplitude throughout
o16 <- cwt(s16); P16 <- Mod(o16$W)^2
b16 <- which.min(abs(o16$period - 16))
c(power_at_year_128 = round(P16[b16, 128], 2),
  power_at_year_20 = round(P16[b16, 20], 2),
  power_at_year_5 = round(P16[b16, 5], 2),
  year_5_as_share_of_mid = round(P16[b16, 5] / P16[b16, 128], 3),
  cone_reaches_period_16_at_year = ceiling(16 / (fourier_factor(6) / sqrt(2))))

             power_at_year_128               power_at_year_20 
                       209.390                        174.330 
               power_at_year_5         year_5_as_share_of_mid 
                        87.460                          0.418 
cone_reaches_period_16_at_year 
                        22.000 

The sine has constant amplitude for all 256 years. Its measured period-16 power at year 5 is 41.8% of its mid-series value, purely because the wavelet is hanging over the edge into padding. Never read the first or last 22 years of a period-16 band as evidence that a cycle started or stopped.

Trap three: pointwise is not areawise

Here is the deepest one, and it survives fixing the other two. Look again at the red-null panel above. The test is calibrated, 4.6% of cells flagged, exactly as advertised. But those cells are not scattered like confetti. Wavelet power is smooth in both directions, because neighbouring scales use overlapping wavelets and neighbouring times use overlapping windows. False positives therefore arrive in clumps.

Label them and count. Four-connected components, no package needed:

lab_patches <- function(M) {
  lab <- matrix(0L, nrow(M), ncol(M)); cur <- 0L
  for (i in seq_len(nrow(M))) for (j in seq_len(ncol(M))) {
    if (!M[i, j] || lab[i, j]) next
    cur <- cur + 1L
    stack <- list(c(i, j)); lab[i, j] <- cur
    while (length(stack)) {
      p <- stack[[length(stack)]]; stack[[length(stack)]] <- NULL
      for (d in list(c(-1, 0), c(1, 0), c(0, -1), c(0, 1))) {
        a <- p[1] + d[1]; b <- p[2] + d[2]
        if (a >= 1 && a <= nrow(M) && b >= 1 && b <= ncol(M) && M[a, b] && !lab[a, b]) {
          lab[a, b] <- cur; stack[[length(stack) + 1]] <- c(a, b)
        }
      }
    }
  }
  lab
}
flagged <- (P > lev_red) & usable
lab <- lab_patches(flagged)
sizes <- sort(tabulate(lab[lab > 0]), decreasing = TRUE)
c(flagged_cells = sum(flagged), patches = length(sizes),
  largest_patch = sizes[1], second = sizes[2], third = sizes[3])

flagged_cells       patches largest_patch        second         third 
          492            13           231            62            39 

492 false-positive cells, and they form only 13 patches. The largest is 231 cells: a solid, contiguous, eye-catching red island in a series with no cycle. Print that figure and the discussion writes itself.

So the honest test is not “is this cell above the level” but “is this patch bigger than the patches noise makes”. Simulate the null again, and this time record the largest patch each time.

mc_patch <- function(B, n, phi, seed) {
  set.seed(seed); big <- integer(B)
  for (b in 1:B) {
    z <- as.numeric(arima.sim(list(ar = phi), n = n))
    oz <- cwt(z); Pz <- Mod(oz$W)^2
    ok <- outer(oz$period, oz$coi, "<=")
    a <- as.numeric(acf(z, lag.max = 1, plot = FALSE)$acf[2])
    L <- lab_patches((Pz > sig95(oz, z, a)) & ok)
    tb <- tabulate(L[L > 0])
    big[b] <- if (length(tb)) max(tb) else 0L
  }
  big
}
null_big <- mc_patch(200, n, 0.7, 2285)
c(null_largest_patch_median = median(null_big),
  null_largest_patch_95th = round(quantile(null_big, 0.95), 0),
  null_largest_patch_max = max(null_big))

  null_largest_patch_median null_largest_patch_95th.95% 
                      187.5                       432.0 
     null_largest_patch_max 
                      583.0 

Under the correct, calibrated red-noise null the largest patch on the map is typically 187.5 cells, and one map in twenty has a patch of 432 cells or more. Our cycle-free series produced 231, comfortably inside the ordinary range. Pointwise significance told us 231 cells were remarkable; the areawise view says they are Tuesday.

What does pass

The test has to be able to say yes, or it is just pessimism. Put a real transient cycle into the same red noise: a period-10 oscillation running only from year 80 to year 140.

set.seed(3285)
env <- ifelse((1:n) >= 80 & (1:n) <= 140, 1, 0)
amp <- 3.0
z <- as.numeric(arima.sim(list(ar = 0.7), n = n)) + amp * env * sin(2 * pi * (1:n) / 10)
oz <- cwt(z); Pz <- Mod(oz$W)^2
okz <- outer(oz$period, oz$coi, "<=")
az <- as.numeric(acf(z, lag.max = 1, plot = FALSE)$acf[2])
Lz <- lab_patches((Pz > sig95(oz, z, az)) & okz)
tz <- sort(tabulate(Lz[Lz > 0]), decreasing = TRUE)
p_area <- mean(null_big >= tz[1])
proc_sd <- sqrt(1 / (1 - 0.7^2))
c(cycle_amplitude = amp, noise_process_sd = round(proc_sd, 2),
  flagged_cells = sum(Lz > 0), largest_patch = tz[1],
  areawise_p = round(p_area, 3))

 cycle_amplitude noise_process_sd    flagged_cells    largest_patch 
           3.000            1.400          560.000          513.000 
      areawise_p 
           0.015 

The transient produces a largest patch of 513 cells against a null 95th percentile of 432, for an areawise p-value of 0.015. Note the flagged-cell count, 560, is barely different from the cycle-free series’ 492. The number of significant cells carried no information at all. Their arrangement carried all of it.

Figure 3: Distribution of the largest contiguous patch of pointwise-significant cells across 200 cycle-free red noise series, with the 95th percentile marked. The transient cycle’s patch sits beyond it; the cycle-free series’ patch sits inside the pile.

How big does a cycle have to be? Sweep the amplitude against the same null:

sweep_one <- function(amp, seed = 3285) {
  set.seed(seed)
  zz <- as.numeric(arima.sim(list(ar = 0.7), n = n)) + amp * env * sin(2 * pi * (1:n) / 10)
  oo <- cwt(zz); PP <- Mod(oo$W)^2
  ok <- outer(oo$period, oo$coi, "<=")
  a <- as.numeric(acf(zz, lag.max = 1, plot = FALSE)$acf[2])
  LL <- lab_patches((PP > sig95(oo, zz, a)) & ok)
  tb <- tabulate(LL[LL > 0])
  bg <- if (length(tb)) max(tb) else 0L
  c(amplitude = amp, amp_over_noise_sd = round(amp / sqrt(2) / proc_sd, 2),
    largest_patch = bg, areawise_p = round(mean(null_big >= bg), 3))
}
t(sapply(c(0, 1.0, 1.6, 2.2, 3.0, 4.0), sweep_one))

     amplitude amp_over_noise_sd largest_patch areawise_p
[1,]       0.0              0.00           172      0.545
[2,]       1.0              0.50           172      0.545
[3,]       1.6              0.81           325      0.135
[4,]       2.2              1.11           431      0.055
[5,]       3.0              1.51           513      0.015
[6,]       4.0              2.02           583      0.005

The amp_over_noise_sd column is the cycle’s own standard deviation divided by the background’s. The areawise test crosses 0.05 at about one to one, and the cycle had six full oscillations in which to make its case. Below that it is invisible: at amplitude 1.0 the largest patch on the map is not even the cycle’s, it is the same noise patch the cycle-free run produced. Note also that amplitude 2.2 lands at p = 0.055 while 3.0 lands at 0.015. With 200 null replicates the p-value near the boundary carries visible Monte Carlo error of its own, so a result at 0.05 is a result at 0.05, not a finding.

A check that came back clean

Red noise is itself a model, and ecological time series are often redder than AR(1) can manage: population records show 1/f-like spectra with more long-period variance than a single autoregressive parameter allows (Halley 1996; Vasseur and Yodzis 2004). If the AR(1) null is misspecified in that direction, does it break?

pink <- function(n, beta = 1) {
  m <- 2^ceiling(log2(n)) * 2
  z <- fft(rnorm(m))
  f <- c(1, seq_len(m / 2), rev(seq_len(m / 2 - 1)))
  x <- Re(fft(z / f^(beta / 2), inverse = TRUE))[seq_len(n)]
  as.numeric(scale(x))
}
set.seed(4285)
pk_rate <- replicate(200, {
  z <- pink(n)
  oz <- cwt(z); Pz <- Mod(oz$W)^2
  ok <- outer(oz$period, oz$coi, "<=")
  mean((Pz > sig95(oz, z, as.numeric(acf(z, lag.max = 1, plot = FALSE)$acf[2])))[ok])
})
c(pink_noise_rejection_rate = round(mean(pk_rate), 3), nominal = 0.05)

pink_noise_rejection_rate                   nominal 
                    0.034                     0.050 

It does not. Fitting an AR(1) null to 1/f noise flags 3.4% of the map against a nominal 5%: conservative, not the other way. This one is worth reporting precisely because it is a negative result. The AR(1) null’s weakness is not the colour of the noise. It is that people leave it set to white.

Honest limits

The lag-1 autocorrelation is estimated from the same series that may contain the cycle, and a cycle inflates it, which lifts the null and costs power. The areawise test above uses one fixed patch definition (four-connected cells above a pointwise 95% level) and a purpose-built simulation; the published versions define patch geometry more carefully (Maraun et al. 2007). And none of this touches the deeper point that a “significant” patch means only that the background model does not explain it, which is a long way short of a cycle.

Where to go next

One series at a time is the easy case. The next post asks whether two species cycle together and which one leads, where the null is even less forgiving and the answer to “is this coherence high” turns out to be “compared to what”.

Related tutorials

• Wavelet analysis of population cycles
• Checking a wavelet analysis
• Spectral analysis of population cycles
• The false discovery rate and Benjamini-Hochberg

References

Torrence and Compo 1998 Bulletin of the American Meteorological Society 79(1):61-78

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

Maraun, Kurths and Holschneider 2007 Physical Review E 75:016707 (10.1103/PhysRevE.75.016707)

Cazelles, Cazelles and Chavez 2014 Journal of the Royal Society Interface 11(91):20130585 (10.1098/rsif.2013.0585)

Halley 1996 Trends in Ecology and Evolution 11(1):33-37 (10.1016/0169-5347(96)81067-6)

Vasseur and Yodzis 2004 Ecology 85(4):1146-1152 (10.1890/02-3122)