sum to zero constraints in GPs

123 views
Skip to first unread message

Garyfallos Konstantinoudis

unread,
Nov 7, 2023, 1:13:10 PM11/7/23
to nimble-users
Dear all, 

I am trying to model the effect of temperature on all-cause mortality using NIMBLE and Gaussian processes. I have been modifying this tutorial (https://r-nimble.org/nimbleExamples/gaussian_process.html). The code I am using is as follows:

expcov <- nimbleFunction(    
  run = function(dists = double(2), rho = double(0), sigma = double(0)) {
    returnType(double(2))
    n <- dim(dists)[1]
    result <- matrix(nrow = n, ncol = n, init = FALSE)
    sigma2 <- sigma*sigma
    for(i in 1:n)
      for(j in 1:n)
        result[i, j] <- sigma2*exp(-dists[i,j]/rho)
    return(result)
  })
cExpcov <- compileNimble(expcov)

# model
model_exposure <- nimbleCode(
  {
    for(k in 1:N){
     
      O[k] ~ dpois(lambda[k])
      lambda[k] <- exp(log(E[k]) + beta_0 + b[X[k]])
    }
   
    # The Gaussian process for temperature
   
    cov[1:Jb, 1:Jb] <- expcov(dists = dist_cov[1:Jb, 1:Jb],
                              rho = rho,
                              sigma = sigma)
   
    b[1:Jb] ~ dmnorm(zeros[1:Jb], cov = cov[1:Jb, 1:Jb])
   
    # sum-to-zero contraints
    sum.constr <- sum(b[1:Jb])
    one ~ dconstraint(sum.constr == 0)
   
    # set the priors of the fixed effects
    # intercept
    beta_0 ~ dnorm(0, sd = 10)
   
    # set the priors of the hyperparameters
    sigma ~ dunif(0, 100) 
    rho ~ dunif(0, 20)
  }
)

I have a couple of questions:
1. I would like to add sum to zero constraints on the temperature effect b but the moment I add the following lines:
     sum.constr <- sum(b[1:Jb])
    one ~ dconstraint(sum.constr == 0)
the b doesnt move. Am I approaching it in a wrong way?
2. I tried modifying the  expcov to make it a squared exponential covariance so it is kind of smoother but i get an error when taking the chol that the matrix is not positive definite. I also tried to use the absolute difference between the bins in the temperature domain and again is not pos def. Any ideas?

See attached the code and the data. 

Best wishes, 
Gary




findata.rds
nimbleRscript.R

Chris Paciorek

unread,
Nov 9, 2023, 8:59:49 PM11/9/23
to Garyfallos Konstantinoudis, nimble-users
Hi Gary,

Regarding the constraint, that won't work. The reason being that any standard sampler assigned to `b` will fail because any proposal that changes values of `b` arbitrarily will invalidate the constraint.
For it to work as written, you'd have to write a special sampler that proposes valid moves in the reduced dimension space implied by the constraint. Instead, for your particular model, I suggest moving `beta_0` so that it is the mean of the `b` vector. Then you won't need a constraint -- the reason you need it now is the possibility of trading off between the mean of `b` and the value of `beta_0`.

mn[1:Jb] <- beta_0 * ones[1:Jb]
b[1:Jb] ~ dmnorm(mn[1:Jb], cov = cov[1:Jb, 1:Jb])

Regarding the positive definiteness, that's a well-known issue with the squared exponential. Unless the locations are all somewhat far apart, the matrix, while mathematically positive definite, is not numerically positive definite because there is at least one eigenvalue that is very close to zero (mathematically) that because of floating point imprecision ends up being negative on the computer. One approach people use (often in machine learning) is to add a small epsilon on the diagonal of the covariance. Another standard approach in spatial stats would be to use the Matern covariance. You'd need to write a user-defined function to compute that but we do provide the Bessel function that is needed, so it should be reasonably straightforward. Happy to help if you try it and have questions. 

(Side note: if temperature is binned, I'm slightly surprised you have the problem with positive definiteness, but perhaps your bins are pretty narrow?)

-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 view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/05f41646-97e2-46de-b42e-f08a68bc1f2en%40googlegroups.com.

Garyfallos Konstantinoudis

unread,
Nov 10, 2023, 6:01:12 AM11/10/23
to nimble-users
Hi Chris, 

thank you very much for the reply. 

Regarding the constraints I will implement accordingly. 

I tried to write the Matern covariance in my code, see attached. Some things to validate:
- I couldnt find gamma function here: https://r-nimble.org/html_manual/cha-RCfunctions.html (but lgamma) but it seems that there is a gamma function built, so I tried gamma() (but also exp(lgamma()))
- I used the besselK() function, it seems that in the manual the argument is k, whereas i think is built with x
-  I fixed the diagonal to 1, as the besselK of zero is Inf

Now I have the following error when im running nimbleModel()

Defining model Error in getSymbolicParentNodesRecurse(code, constNames, indexNames, nimbleFunctionNames, : R function 'cMatcov' has arguments that cannot be evaluated in the code 'cMatcov(dists = dist_cov[1:100, 1:100], rho = rho, sigma = sigma, nu = 0.5)'. Either the function must be a nimbleFunction or values for the following inputs must be specified as constants in the model: dist_cov[1:100, 1:100],rho,sigma.

I tried a some things but I still get this error, do you have any ideas?

(see updated code attached)

Best wishes, 
Gary
nimbleRscript_matern.R

Chris Paciorek

unread,
Nov 12, 2023, 1:56:11 PM11/12/23
to Garyfallos Konstantinoudis, nimble-users
Hi Gary,

- Please use `lgamma` for the log of the gamma function, and exponentiate as needed. You may want to do various calculations on the log scale, exponentiate and then multiply by the result from the besselK call. Working on the log scale is often better numerically, though in this case you need the result on the original scale and the besselK is done on the original scale.
- Yes, the manual should say 'x' not 'k' as the argument of the besselK function.
- Yes, setting the correlation to one at distance zero makes sense.

The error you are getting is because you are compiling to `cMatCov` and then trying to use that in the model. Just use `matcov` in the model and nimble will deal with compiling `matcov`.
You should remove `cMatcov <- compileNimble(matcov)`.

-chris


Garyfallos Konstantinoudis

unread,
Nov 15, 2023, 12:30:48 PM11/15/23
to nimble-users
Hi Chris, 

thank you for this. I managed to make it work but in the end decided to go with the squared exponential covariance as it gives me the amount of smoothing I expect in the result. The problem is that it doesnt converge. I have tried a couple of things:
1. The exponential one as showed on the online tutorial, works and converges but the result is not smooth enough
2. I played with different priors for the squared exponential hyperparameters, but is seems that the hyperparameters converge (somewhat) but the chains of the actual random variables of the GP are really bad (bad mixing, high correlation). Do you have any ideas why is this the case?
3. I tried to use nimbleHMC, but I have an error: "failed to create library", when running the compileNimble, any ideas for that?

See my code attached. Thank you in advance!

Best wishes, 
Gary

nimbleRscript_matern.R

Daniel Turek

unread,
Nov 16, 2023, 10:42:55 AM11/16/23
to Garyfallos Konstantinoudis, nimble-users
Gary, I'll jump in RE nimbleHMC.  I tried running your script, noting that you load nimbleHMC library at the beginning, although nimbleHMC is never used, and I didn't encounter any compilation errors (excepting one error with the compiled function name "csqExpcov" vs. "cExpcov").

Are you able to create a reproducible script, demonstrating the compilation error you had when using nimbleHMC?

Thanks,
Daniel


Reply all
Reply to author
Forward
0 new messages