spectrum matching in pomp

[Fuhan-Yang]

I'm playing with 'spect()' in pomp and comparing with 'spec.pgram()' in R. While I found the difference in the outputs using the same data.

The 'spec.pgram()' calculated the raw periodogram using discrete Fourier transform (fast = F) with no smoothing. The argument I used for 'spect()' is `kernel.width = 1'. I also exponential transformed the output of 'spect()' to make sure they have the same scales (By default, the spec from 'spect()' is log-transformed, right?) While I still could see the range and the pattern are a little bit different, could you tell me what I missed there?

Another question is that, I tried to calculate the metrics of 'spect()' using logLik(), while I got NA from all the models. Could you please tell me what could be the causes?

[kingaa]

In spect.R e.g., lines 299ff, you can see that pomp's spect() calls spect.pgram, with certain non-default options, on de-trended and possibly transformed data. Could this be the reason for the differences you see?

I am not sure what you mean by "calculate the metrics of spect() using loglik()". Perhaps you can explain what you are trying to accomplish and provide some code showing what you have done?

[Fuhan-Yang]

Sure! I'm interested in using logLik(spect()) in my data to calculate how similar the simulations are to the data in terms of power spectrum. I used to use spec.pgram(), thus I want to make sure they produce the same output for the same data at first.

Below is the code:

library(pomp)
library(dplyr)

po <- ricker()
y <- po@data[1,]

plot(y,type = 'l')

###############################
## get the full power spectrum #
###############################

#### spec.pgram ####
spec <- spec.pgram(y,spans = NULL, kernel = NULL, taper = 0, pad = 0,fast = F,plot = F)
cycle_spec <- data.frame(cycle = 1/spec$freq, spec = spec$spec)

#### spect in pomp ####
po |>
  spect(nsim = 100,
    params = c(r = exp(3.8), sigma = 0.3, phi = 10, c = 1, N_0 = 7),
    vars = 'y',
    kernel.width = 1,
    detrend = 'none') -> spec_po

datspec <- as.numeric(spec_po@datspec)

cycle_spec_po <- data.frame(cycle = 1/spec_po@freq, spec = datspec)

op <- par(mfrow = c(2,1),mar = c(3,3,2,2),mgp = c(2,1,0))
plot(cycle_spec$cycle,cycle_spec$spec,type = 'h',
  xlim = c(2,10),xlab = 'cycle',ylab = 'spec',main = 'spec.pgram')
plot(cycle_spec_po$cycle,exp(cycle_spec_po$spec),type = 'h',
  xlim = c(2,10),xlab = 'cycle',ylab = 'exp(spec)',main = 'spect')
par(op)

#####################################
## get the loglikelihood of spectrum #
#####################################

logLik(spec_po) # NA

And the output is:

The range and the pattern are different, and the logLik() returned NA. I wonder what caused the differences, and why NA was returned. Thank you!

[kingaa]

As I mentioned above, the difference in the plots is likely due to the fact that spect() uses different arguments in its call to spec.pgram() than you use in your own call. You can examine the code to see just what these differences are.

logLik() returns NA when you call it on the 'spectd_pomp' object created by your call to spect() because it is undefined for objects of that type. If you want an objective function for spectrum-matching, use spect_objfun to create one. For example:

library(pomp)

ricker() |>
  spect_objfun(
    nsim = 100,
    params = c(r = exp(3.8), sigma = 0.3, phi = 10, c = 1, N_0 = 7),
    vars = "y",
    kernel.width = 3,
    detrend = "none",
    est=c("r","sigma"),
    partrans=parameter_trans(log=c("r","sigma")),
    paramnames=c("r","sigma")
  ) -> spec_fn

optim(
  fn=spec_fn,
  par=c(3,log(0.1))
) -> fit

spec_fn(fit$par)

plot(spec_fn)
coef(spec_fn)
logLik(spec_fn)

See the manual for more about spectrum-matching.

[Fuhan-Yang]

Thank you so much for the reply. The difference has been resolved, while yet there is an issue in the likelihood. There is a model that always returned NA in logLik from the spect_objfun. When checking the simulated power spectrum, there are two potential issues.

The first one is that this model is affected by stochastic extinction. That means, most of the simulated spectrum (> 98%) returned -Inf because most of the simulations went extinct during the burn-in period (before the beginning of the data). Not sure if this affects likelihood calculation because the sample size is small (7 out of 500 simulations).

The second one is that the remaining 7 simulations always have a dominant semi-annual cycle. Not sure if this caused NA in the likelihood, such as zero covariance matrix. While the spectral density of semi-annual cycle varies across simulations but always center around a certain value.

Could you explain how is the likelihood calculated for the power spectrum in pomp? In addition, I would like to know your comments about how convincing the dynamics can explain the data if it is severely affected by stochastic extinction while the valid simulations fit the data well.

Thank you!

[kingaa]

It's important to understand that, although the spectrum-matching objective function value can be extracted using the logLik function, it is not a likelihood, mathematically speaking. Rather, it is a discrepancy measure, due to Reuman et al. (the reference is in the spect help page). You can look at the code to see just what it is: the spect_discrep function is defined here.

With respect to your question about extinction. I don't want to sound glib, and I'm not perfectly sure what you're asking, but it seems that if your model predicts extinction more than 98% of the time, and the data don't show extinction, then perhaps there's a better explanation out there somewhere.

[kingaa]

Of course, the other possibility is that you want to explain the data by saying we got really lucky (or unlucky).

[Fuhan-Yang]

My data seems to show that this disease is an endemic equilibrium that is under the weak effect of climate drivers. By saying weak, first, parameter estimates showed this. Also, the simulations very easily went extinct at the initial outbreak under the nonlinear effects of climates. That means, this data is one of the observations of a dynamic that are very likely to have already been gone extinct...From this perspective, maybe another mechanism is more feasible. However, the goodness-of-fit is good, the median from the survived simulations could capture the data trend well.

So I wonder what is the role of stochasticity. It seems like stochasticity could capture the trend that could not be explained by the explicit assumptions from dynamical models (the fit from stochastic model is better than the deterministic model), but the chosen model is still a little suspicious. Because we could conclude this trend is solely caused by stochasticity, it could be possible that stochasticity kind of mask the true and more complicated mechanism behind the data. Maybe a more complicated model could solve this. This is just an irrelevant thought, any comments are welcome!

[kingaa]

My data seems to show that this disease is an endemic equilibrium that is under the weak effect of climate drivers.

This is an interpretation. Are there other, equally plausible interpretations?

By saying weak, first, parameter estimates showed this. Also, the simulations very easily went extinct at the initial outbreak under the nonlinear effects of climates. That means, this data is one of the observations of a dynamic that are very likely to have already been gone extinct...

Only if your model is correct. However, would not a model under which extinction was less likely be on its face more plausible?

From this perspective, maybe another mechanism is more feasible. However, the goodness-of-fit is good, the median from the survived simulations could capture the data trend well.

What do you mean when you "goodness-of-fit" is good? If 98% of your simulations go extinct, but your data do not, is that a good fit?

Also, what does it mean to compare the median of simulations with the data? Remember that the fundamental assumption of statistics is that the data are one realization from some stochastic process. In particular, they are not the median of several realizations of any process.

So I wonder what is the role of stochasticity. It seems like stochasticity could capture the trend that could not be explained by the explicit assumptions from dynamical models (the fit from stochastic model is better than the deterministic model), but the chosen model is still a little suspicious.

These are indeed interesting questions. Would they still be interesting to you if you had a better model in hand?

Perhaps it would be worthwhile to think about the forces that prevent extinction in the actual process....

[Fuhan-Yang]

Thanks for the comments, this is helpful as a discussion point. Thank you so much!