Truncation in dcar_normal

[Mattia Falaschi]

Dear nimble users, 

I am working on model to describe the activity of an animal over the year.

I have a trigonometric model where I implemented a spatial random effect with the dcar_normal distribution. Here is an example:

model <- nimbleCode({
 
  # Trigonometric model
  for (i in 1:nrows){
    y[i] ~ dbin(p[i],Nrecords[i])
    logit(p[i]) <- b0 + b1*cos( ( doy[i]*2*pi/(366)  ) + ( -theta0 -  rho[cell[i]]) )
  }
 
  # CAR model
  rho[1:ncell] ~ dcar_normal(adj=adj[1:sumneigh],weights=weights[1:sumneigh],num=num_adj[1:ncell],tau=tau,zero_mean=1)
  # Variance of car model
  tau ~ T(dnorm(0,sd=0.5),0,)

  ## Priors of fixed effects
  b0 ~ T(dnorm(0,0.1),-10,10)
  b1 ~ T(dnorm(0,0.1),0,10)
  gamma0 ~ dunif(-pi, pi)
  theta0 <- pi*cos(gamma0)

})

However, since I have a trigonometric model, values at certain distances (2*pi) have the same meaning. The result is that the model fails to reach convergence and show multiple peaks at certain sites, e.g.:

So, my Idea was to truncate the dcar_normal distribution with T(dcar_normal(...)) but it is not possible, since it returns this error:

Rmodel <- nimbleModel(code = model, constants = constants, data = data, inits = inits)
Defining model
Error in processBoundsAndTruncation() :
  Cannot implement truncation for dcar_normal; 'p' and 'q' functions not available.

Do you have any workaround to truncate the dcar_normal distribution or other suggestions?

All the best,

Mattia 

[dirkdouw...@gmail.com]

Hi Mattia,

I have dealt with similar issues. You can use dconstraint, which is described on page 53 of the user manual https://r-nimble.org/manuals/NimbleUserManual.pdf. I am not exactly sure how you want to truncate it, but let's say you want rho[1] to be between 0 and 2pi, then just do,

 constraint_data ~ dconstraint(0<rho[1]  & rho[1]<2*pi),

where contraint_data is 1 and entered in as data.

I think that should work, but my understanding is that it truncates the joint prior distribution, which is not necessarily the same as truncating one of the marginal priors. So it will change the marginal priors of the other parameters as well. For instance, the marginal prior of tau might change to better reflect the constraint. This is probably not a big deal, but something to be aware of.

[Mattia Falaschi]

Hi Dirk and thank yuo for your suggestion.

I managed to use dconstraint, and it seems to be working properly.

The only issue is that I get this error message after the end of each chain:

running chain 1...
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
  [Error] Invalid model values were induced during the MCMC, which were caused by centering of the CAR (dcar_normal) process.  Centering takes place because of the argument zero_mean=1 to the dcar_normal distribution.  Results of this MCMC run may be invalid. This can be avoided by absorbing the mean into the CAR process, and omitting the zero_mean argument.

Have you ever encontered this? The results I obtain seem to be what expected (still no convergence, but this may be related to model complexity).

All the best, 

Mattia

[dirkdouw...@gmail.com]

Hi Mattia,

Someone else might know more than me, but my understanding is that zero_mean=1 will recenter the random effects at 0 after sampling. So, if you constrain them to be greater than 0, then obviously they cannot be recentered at 0 without violating that constraint. I would try to figure out with your model where you want rho to be centered. You can always do zero_mean=0, and then it will estimate the mean, but this may cause identifiability issues depending on your model. Or if you know rho should be centered at, say, a, then you could just replace rho with rho+a.

[Mattia Falaschi]

Thanks again!

I think that rho make sense to be centred on zero in my case (if more info are needed, I can better explain the model).

This is the model I am using:

model <- nimbleCode({

 
  for (i in 1:nrows){
   
    y[i] ~ dbin(p[i],Nrecords[i])
   
    logit(p[i]) <- b0 +

      b1*cos( ( doy[i]*2*pi/(366)  ) + ( -theta0 -  rho[cell[i]] ) ) +
      b2*cos( ( doy[i]*2*pi/(366/2)) + ( -theta1 - (rho[cell[i]])*2 ) )

   
  }
 
  # CAR model
  rho[1:ncell] ~ dcar_normal(adj=adj[1:sumneigh],weights=weights[1:sumneigh],num=num_adj[1:ncell],tau=tau,zero_mean=1)
  # Variance of car model
  tau ~ T(dnorm(0,sd=0.5),0,)

  for(i in 1:ncell){
   constraint_data[i] ~ dconstraint( (-pi/2) < rho[i]  & rho[i] < (pi/2) )
  }

 
 
  ## Priors of fixed effects
  b0 ~ T(dnorm(0,0.1),-10,10)
  b1 ~ T(dnorm(0,0.1),0,10)

  b2 ~ T(dnorm(0,0.1),0,10)

  gamma0 ~ dunif(-pi, pi)
  theta0 <- pi*cos(gamma0)

  gamma1 ~ dunif(-pi, pi)
  theta1 <- pi*cos(gamma1)

})

I highlited the dconstraint part, and as you can see, I am constraining rho between -pi/2 and pi/2, so I do not expect the issue of re-centering at zero.

Now, I am just concerned about whether these values outside the range are somehow influencing the estimates of parameters or not.

(anyway, the model now managed to reach convergence and the results are what I expexted, so thanks again!)

[Perry de Valpine]

Hi Mattia,

When zero_mean is TRUE (or 1), then after the special CAR updates that are done on each rho[i], the entire vector rho[1:ncell] is re-centered at 0. However, there is no way in that step to automatically respect dconstraint declarations on the rho[i]s. The warning message is likely indicating that when the re-centering is done, some of the dconstraints are violated, and the CAR sampler doesn't know what to do about it so it's simply telling you that it happened. You could do compiled_model$getLogProb() after the mcmc has run to see if it returns NaN or -Inf, which would indicate that the model has some invalid (impossible) values, and those would likely be from violations of the dconstraints.

I don't see a great way around this given your modeling need, so I would suggest using zero_mean=FALSE (or 0). You could then consider a couple of things after the MCMC has run. One thing would be to simply accept the results as the posterior, recognizing that the dconstraint on each element of rho[i] has effectively forced the rho[1:ncell] to be identifiable. However, that is a somewhat weird kind of identifiability. So another thing you could consider would be to post-process the posterior samples to center rho[1:ncell] on zero and move the mean of the b1*cos(...) and b2*sin(...) terms into b0. Evidently, you would then get some rho[i] values (after centering on 0) outside of the [-pi/2, pi/2] range, but perhaps those wouldn't turn out to be problematic.

HTH

Perry

--
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 visit https://groups.google.com/d/msgid/nimble-users/d3254e49-8043-491d-baa4-92d29de9cfacn%40googlegroups.com.

[Mattia Falaschi]

Hi Perry, and thank you for your reply.

As for the first part, I tried with compiled_model$getLogProb() and I get a value aroun -7000, so no NaN or -Inf.

I am not sure I completely got the meaning of the second part.

For now, I re-run the model with zero_mean=0 and it runs smoothly. I also added a derived parameter to estimate the mean of rho, by defining rho.mean <- mean(rho[1:ncell]) which in my trial was estimated to be around 0.1.

After that, I am not sue what you mean with "center rho[1:ncell] on zero and move the mean of the b1*cos(...) and b2*sin(...) terms into b0".

What I did, is to add rho.mean into the trigonometric functions b1*cos(...) and b2*cos(...) when predicting activity across the year, in order to have curve showing the average activity after accounting for the spatial random effect, but I am not sure this is what you intended.

Say "out" is the model output, I can do:

doys <- 1:366
ch <- rbind(out$chain1, out$chain2, out$chain3)
pred <- matrix(NA, nrow=length(doys),ncol=nrow(ch))
for(i in 1:nrow(ch)){
  pred[,i] <- plogis(ch[i,"b0"] + 

                     ch[i,"b1"]*cos((doys*2*pi/(366)  )+(-ch[i,"theta0"] - (ch[i,"rho.mean"]) ))+
                     ch[i,"b2"]*cos((doys*2*pi/(366/2))+(-ch[i,"theta1"] - (ch[i,"rho.mean"])*2 )))
}

Just for completeness, the predictions I get with the two methods (zero_mean=1 or zero_mean=0 with re-centering based on rho.mean) are similar but not identical, as you can see here:

[Perry de Valpine]

Hi Mattia,

Thanks for this. I see my suggestion was a bit hasty. I should have said theta0 instead of b0, and then I think we are on the same page.

Using a simpler case: if one has rho[1:n] following dcar with zero_mean=TRUE, then often one has also something like:

eta[i] <- intercept + <other linear terms> + rho[i]

# followed by use of a function of eta[i] (perhaps through a link function or other step) as a parameter of a data distribution.

In that case, since mean(rho[1:n]) will be held at 0, the intercept is also identifiable.

If one has the same setup but uses zero_mean=FALSE, then if one considers adding a constant c to intercept and subtracting the same constant c from all rho[i], one gets an equivalent model in terms of downstream likelihoods (although perhaps with different prior probabilities). For this reason, both the intercept and the rho[i] may look poorly mixed and wander randomly such that intercept + mean(rho[1:n]) is well identified because as the intercept wanders to higher values, mean(rho[1:n]) wanders to correspondingly lower values

In such a situation, it is actually reasonable to run the MCMC anyway and then post-process the samples to obtain the posterior of interest. That would be, for each i, setting intercept_prime = intercept + mean(rho[1:n]) and rho_prime[i] = rho[i] - mean(rho[1:n]), where the "prime" names are the identifiable quantities.

In your situation, I was trying to apply the same logic, but I looked at b0 as the relevant intercept, but now I realize it is theta0, and also theta1 which makes the problem also different in another way. Let's say you just had the theta0 part (as in the first code you sent). Then I think adding a constant c to theta0 and subtracting the same constant c from all rho[i] would give an equivalent model for downstream likelihoods (e.g. for y[i]). So you could post-process by setting theta0_prime = theta0 + mean(rho[1:n]) and rho_prime[i] = rho[i] - mean(rho[1:n]). 

Curiously, you would then also need to do similar steps for theta1 using 2*mean(rho[1:n]).

Also I guess what you really need is mean(rho[cell[1:n]]), i.e. the average over the nested rho indices from cell[], since I don't know if every cell is represented exactly once. That won't work directly in model code but you could do it in a nimbleFunction or after the MCMC (if you are monitoring rho).

If I follow that last code you sent, that might be doing essentially what I'm saying here. The code you sent that makes the figures seems to omit the actual rho[i] values, and I think you intended that for prediction purposes. If so, that matches treating theta0 + mean(rho[1:n]) as the equivalent theta0 you would obtain if the rho[1:n] were kept centered.

Perry

[Mattia Falaschi]

Thank you Perry, now it is all clear to me!

I think that in this specific case mean(rho[cell[1:n]]) and  mean(rho[1:ncell]) are equivalent, given that rho is defined as "rho[1:ncell] ~ dcar_normal(...)" . In the dataset I actually have more records for each "cell", but the dimension of the parameter is equal to the number of cells (ncell).

Cheers, 

Mattia