Smoother 1D CAR Process in NIMBLE

[Oliver Stoner]

Dear all,

I am currently working on a discrete-time series model where I would like a relatively smooth temporal process in the mean. The model is for epidemiological forecasting and involves the last few data being missing so that they can be predicted. I started using a penalised cubic spline set up with GAM, which was very efficient (good sampling after 1000 iterations of burn-in) but I feel it is inappropriate because, in the missing region at the end of the time series, the combination of the lack of data and the smoothing penalty meant that the spline just continued whatever trend it had before the missing data.

I have also tried using a multivariate normal process with a squared exponential or a matern covariance function of temporal distance, which produced an ideal process but was very inefficient using the default RW_block sampler, requiring 1 million samples of burn-in. I tried using an AF_slice sampler which worked better but it was not reliable, often reaching the maximum number of contractions. I don't really know why this happens despite adding a small amount to the diagonal of the covariance matrix to make it more numerically stable. Has anyone had this issue for gaussian processes using these samplers?

Sigma[1:N,1:N] <- exp(-(dist[1:N,1:N]/rho)^2)*sigma^2 + 0.0001*diag(N)

delta[1:N] ~ dmnorm(zero[1:N],cov=Sigma[1:N,1:N])

sigma ~ T(dnorm(0,1),0,)

rho ~ dlnorm(2,1)

My most recent attempt is to base the process on an improper CAR prior in one dimension, hoping to emulate a thin plate spline like the 'rw2' model in INLA. I tried setting order=2 in dcar_normal, as suggested in the manual, without changing any of the weights. However, this still results in quite a coarse process:

adj=c(2,as.numeric(rbind(1:(N-2),3:N)),N-1)

n_adj=c(1,rep(2,N-2),1)

delta[1:N] ~ dcar_normal(adj=adj[1:l_adj],num=n_adj[1:N],tau=nu,c=2,zero_mean=1)

sigma ~ T(dnorm(0,1),0,)

nu <- 1/sigma^2

To achieve a smoother CAR process, should I provide the model with a different adjacency structure and weights? Alternatively, does anyone have any suggestions for how to more efficiently implement a more appropriate temporal process in NIMBLE?

Thanks,

Oliver

[Chris Paciorek]

Hi Oliver, a couple options for the second-order smoothing:

1) if you look at the 'ice' model in the original BUGS examples (available in the classic-bugs/vol2/ice directory in the installed NIMBLE package) you'll see how to specify a 2nd-order random walk directly by conditioning on the past rather than using the CAR representation.  There's also some discussion of this example in the http://www.openbugs.net/Examples/Ice.html.

The representation above is equivalent to an IID model on the second differences in time. (Whereas a first-order model is equivalent to an IID model on first differences.)

2) As far as the CAR representation for a second order model, Section 6.3 of the classic Breslow and Clayton 1993 paper (https://www.jstor.org/stable/pdf/2290687.pdf)  discusses the CAR formulation relative to the conditioning on the past formulation. There is also info in this tutorial I wrote: https://www.stat.berkeley.edu/~paciorek/research/techVignettes/techVignette5.pdf.

So to implement that in NIMBLE: you need to have each time point have as neighbors the nearest two time points in both the future and the past (with fewer neighbors for points at or within one time point of the boundary). And the weights should be what is given in the Breslow and Clayton reference (or in the Rue & Held Gaussian Markov Random Fields book, Section 3.4, equation 3.40). So you'd have  this inputs for the dcar_normal specification, for an example of 6 time points:

adj = c(2, 3, 1, 3, 4, 1, 2, 4, 5, 2, 3, 5, 6, 3, 4, 6, 4, 5)

## weights are negative of the non-zero values show in the Breslow/Clayton and Rue/Held refs, with the diagonals ignored:

weights = c(2, -1, 2, 4, -1, -1, 4, 4, -1, -1, 4, 4, -1, -1, 4, 2, -1 , 2)

num = c(2, 3, 4, 4, 3, 2)

c = 2

zero_mean = 0

That said, I suspect the simple representation in #1 above will mix just as well or better and is simpler to implement and explain.

I'm somewhat surprised the RW_block does as poorly as you found. However squared exponential and Matern with large "nu" parameter (say bigger than 3-4) gives very smooth processes that can strongly constrain neighboring values when sampling (which could still be the case if the amount added to the diagonal is rather small). So that can slow mixing. It might also be that you need to carefully specify the initial proposal covariance matrix - e.g., making sure that the proposal variances are all reasonable.

-Chris

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To post to this group, send email to nimble...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/63efc770-022c-473c-8609-b4c5541a87dc%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Oliver Stoner]

Hi Chris,

 

Thanks for taking the time to provide such a detailed solution. I have tried out the first solution and it works well, providing a nice smooth process and it’s much faster than the squared exponential to run.

 

I’m also keen to try out the CAR solution, as I am hoping to at some point try out using it in a cyclic domain, to capture seasonal patterns. However, when I follow your instructions I get this error when running buildMCMC:

 

Error in CAR_normal_checkAdjWeightsNum(adj, weights, num) :

  weights argument to dcar_normal() should only contain positive values

 

Is there something I’m missing?

 

Thanks,

 

Oliver

[Chris Paciorek]

That looks like an overzealous check on our part. I'll file a bug
report and we'll look into it for our next release. If we fix it
soon, we can let you know how to build the fixed version from our
Github repository.

> To achieve a smoother CAR process, should I provide the model with a different adjacency structure and weights? Alternatively, does anyone have any suggestions for how to more efficiently implement a more appropriate temporal process in NIMBLE?
>
>
>
> Thanks,
>
>
>
> Oliver
>
> --
> You received this message because you are subscribed to the Google Groups "nimble-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
> To post to this group, send email to nimble...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/63efc770-022c-473c-8609-b4c5541a87dc%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to the Google Groups "nimble-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
> To post to this group, send email to nimble...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/002601d447be%240b2f13e0%24218d3ba0%24%40gmail.com.

[Oliver Stoner]

Sounds good, thanks! 

[Daniel Turek]

Oliver, apologies for the over-sight not allowing negative weights in the dcar_normal distribution.  This has now been corrected on our development branch of NIMBLE, which will now correctly allow negative values in the "weights" vector argument to dcar_normal, and also negative values in the "C" vector argument (specifying  normalized weights) to the dcar_proper distribution.

Until our next updated release of NIMBLE, you can install the development branch from GitHub (which includes these changes) using:

> remove.packages('nimble')

> library(devtools)

> install_github('nimble-dev/nimble', ref = 'devel', subdir = 'packages/nimble')

> library(nimble)

Please let me know if you have a chance to try this out, and whether or not it solves your problem.

Cheers,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CANP6ef%3Doqu2sHdbpCP3rvNQB93ktbCZ2nqoxuZFyRjVu50s%3D1g%40mail.gmail.com.

[Oliver Stoner]

Hi all,

Thanks for fixing it so quickly. It does seem to have fixed the problem. The two approaches seem to produce very similar results, but conceptually I need to think more about the implications of conditioning only on the past versus conditioning on both the past and the future. I think conditioning only on the past might be easier to justify in terms of forecasting disease counts into the future, though I can't say for certain.

Both are slightly more efficient than the squared exponential model, though they don't mix well compared to the first order random walk, at least for this model.

Thanks again for the advice and the quick fix. I'm now going to try to replace my cyclic cubic spline which captures the seasonal patterns with a cyclic second-order CAR.

Oliver

[Chris Paciorek]

Just so I'm clear here, I wanted to say that the posterior inference on the latent process conditions on both the future and the past even when representing the model as conditioning on the past. It's just a matter of writing the model in a different representation.

And a side note that in most cases when one increases the smoothness of a process, mixing gets worse because posterior dependence increases.

Chris

--

You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To post to this group, send email to nimble...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/fd658e4a-419f-442b-8476-caba04d64bd7%40googlegroups.com.

[Oliver Stoner]

Hi all,

I have a follow up question about this. My data is one count for each week over a few years i.e. t=1,2,...,N.

Currently in the mean of my model (lambda_t) I have a combination of a first-order random walk over the whole series f(t) and a cyclic second order CAR g(w_t) where w_t=1,2,...,52 representing the time of year.

log(lambda_t) = f(t) + g(w_t)

The latter is designed to capture seasonal structure which is quite evident in this data. By doing this I get a nice smooth seasonal effect from the second order CAR and the first-order random walk does a good job of capturing non-seasonal temporal variation.

However, the model is for forecasting beyond the range of the data and I would really prefer to have a smoother non-seasonal effect where the trend is extrapolated into the future. Thanks to your advice in September I can change the first order random walk in t to a second order random walk. However, when I do this, regardless of what initial values I choose I always end up in a situation where the second-order random walk in t absorbs all of the seasonal variation and the second order CAR in w does nothing. I think this might have something to do with prior choice for the RW1 and CAR models. I have tried changing the second order CAR seasonal effect to a penalised cyclic cubic spline, with the same result.

Have any of you experienced this kind of problem before? Does anybody have a suggestion for the priors or an alternative formulation so that I can have a smooth but distinct temporal and seasonal component?

Best wishes,

Oliver

[Chris Paciorek]

Hmm, interesting. You might be able to write the model in terms of
basis functions and orthogonalize the 2nd-order terms relative to the
cyclic terms. There's work by John Hughes and by Murali Haran among
others on orthogonalization in CAR-type models (in the context of
spatial confounding), and possibly something in Simon Wood's book as
well.

Chris

> --
> You received this message because you are subscribed to the Google Groups "nimble-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
> To post to this group, send email to nimble...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/8dec7164-d4df-4e1f-9b20-ca18275ff5d7%40googlegroups.com.