problem with pbeta() and random effects

[Bob Douma]

Dear all,

I am experiencing problems when I use pbeta in combination with random effects. 

I want to estimate the mean of two random variables with random intercepts and the correlation between the two random variables; one variable is normally distributed and the other is beta distributed. The correlation between the two is estimated through a gaussian copula (a copula allows to separate the margins from the dependence). This requires the use of pbeta. 

When I run this model without random effects it works fine, and estimates are as expected (see nonnormal_beta.R and nonnormal_beta.cpp in attach). However, when I add random effects to the means, the estimates for the correlation parameter and the standard deviation of the normal are far off (see second half of nonnormal_beta.R and  and mixed_model_two_random_int_nonnormal.cpp). Strangely enough it works fine when I fit a correlated gamma-normal with random effects. Is there a problem with pbeta?

Any suggestions how I can solve this? Any help is appreciated

cheers, Bob 

p.s. I am quite new to TMB, so my coding style is not very compact

[Hans Skaug]

Hi, 

Interesting model. I doubt that there is a problem with pbeta. The two most likely explanations are:

i) Your random effect model is not identifiable, i.e. not even an infinite amount of data would determine the parameters.

ii) The Laplace approximation is inaccurate for your model.

My prior guess is on i), although you say that it works when you replace the beta distribution with a gamma.

My recommendation here is to use the "map" argument (as I see you are already using) to nail down exactly which

parameter is causing the problem.

ii) can be ruled out, or eliminated, by using a more accurate integration method than the Laplace approximation:

http://kaskr.github.io/adcomp/namespacegauss__kronrod.html#ad2a3bebe64e675f09612aa6211ef1998

but it takes some work to set it up the first time.

Hans

[Cole Monnahan - NOAA Federal]

It's also worth integrating with MCMC via tmbstan which can provide some insight into issues like non-identifiability. 

[Bob Douma]

Thanks a lot Hans and Cole for your suggestions. I will try them and report back to you.

Bob 

[Bob Douma]

Dear Hans,

I have tried your suggestions. Below my findings:

i) When I constrain the model through the map argument it appears that all parameters to the beta distribution need to be constrained and the correlation parameter need to be set to zero in order to let the model run.

When I run stan the model converges to values that seem decent. When I integrate over the random effects I also get decent estimates (I programmed this in R, not in TMB because I am more fluent in R - so I am sure that I do not make mistakes - see attach Mixed_effectmodel_integrate_beta.R). So I think we can conclude that the model is identifiable. 

ii) When I toggle laplace = TRUE in tmbstan it does not work anymore. Interestingly, the TMB works if I run the optimiser without specifying a gradient. Also in this case I get decent estimates. These estimates are very similar as the ones found by stan, though the random effects estimates are further off than I expected.  See attach for all the files (including the correlated gamma-normal with random effects)

What I do not understand is that the nlimb works when you do not specify the gradient. Does that support the conclusion that the laplace approximation is not accurate enough? Or does that mean that something goes wrong in the automatic differentiation?

Again thanks for your help. And I am happy with the outcome that the model is identified :). My ultimate goal is to fit two correlated mixed effect models with random intercepts (and slopes possibly) and this is a nice step ahead.

Bob 

[Cole Monnahan - NOAA Federal]

Hi Bob, that's quite interesting. I haven't looked at or thought about your model much, but this seems to suggest something is awry with the LA. 

From the tmbstan output can you confirm that all parameters pass checks (Rhat, ESS)? Also, you can do pairs(fit, pars=...) and specify a subset of parameters. I'd recommend playing with that to see if there are any obviously unusual geometries going on. 

What do you mean the laplce=TRUE option doesn't work? It crashes? Or takes too long to run? Or runs but has not converged?

optim will use finite difference if no gradient fn is supplied. What is more surprising is that it works to optimize.

--
To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to a topic in the Google Groups "TMB Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/tmb-users/NSqL2xVQJmw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to tmb-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/49d546f8-9f2d-4207-8160-5fb188f64b77n%40googlegroups.com.

[Bob Douma]

Hi Cole,

When I use the argument laplace=TRUE in the tmbstan function I get the following warning after some time:

[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                        

[2] "  No acceptably small step size could be found. Perhaps the posterior is not continuous?"

[1] "error occurred during calling the sampler; sampling not done

When I use algorithm HMC in tmbstan it is still at the first iteration after 12 hours. 

For the tmbstan (without laplace=TRUE) the parameters pass checks (Rhat is in all cases equal to 1). Lowest ESS is 1056, but on average 2445 with 3000 iterations, 1500 warmup, thin=1 and 4 chains). Traceplot also looks fine (though maybe a bit of autocorrelation in the beta estimates). See attach. 

Bob   

[Hans Skaug]

Hi Bob,

Thanks for doing experiments. Your findings seem to indicate that the Laplace approximation is failing. 

This led me to test:

> obj$fn(obj$par)

iter: 1  mgc: 4.327788e-09 

[1] 1083.761

attr(,"logarithm")

[1] TRUE

> obj$gr(obj$par)

iter: 1  mgc: 4.327788e-09 

outer mgc:  NaN 

[1] NaN NaN NaN NaN NaN NaN NaN

So the gradient of the Laplace approximation is NA in all components, but the Laplace approximation itself can be evaluated.

I have not come across such a situation, and I do not see anything in your code that would indicate a problem.

I do not know how to debug this, but I suggest that you simplify the code as much as possible. For instance,

only have a single y,z pair, check if all fixed effects can be removed (without loosing the NAs)

In conclusion: the way forward is to investigate what makes obj$gr produce NA.

Hans

[Bob Douma]

Hi Hans,

Thanks for spending time on this problem. Ok. I will do some more work to understand what is causing the problem. With a single y,z pair, you mean having one observation for y and z?

Bob 

[Bob Douma]

Hi Hans,

I tried to simplify the problem as much as possible. I did so by a little trick: call pbeta(y,s1,s2), subsequently qbeta(), and afterwards dbeta().  Whether it breaks depends on the data used. Proportions that are close to one or zero lead to an NA evaluation in the gradient. See attach for two subsets of the data (dataset1: values between 0.25 and 0.75 and dataset2 with values above 0.75 and below 0.25).

In fact if you only fit a model to two datapoints whose value are above 0.75 it breaks, and it works when it is below.

So I think the problem is caused by pbeta, because by integrating over the random effect the pbeta values become zero or one. I guess this is a floating point precision problem? This could probably be solved by adding a log.p=T argument to pbeta similar as in R. 

Cheers, Bob 

On Friday, May 28, 2021 at 3:53:08 PM UTC+2 Hans Skaug wrote:

[Hans Skaug]

Hi,

Thanks for simplifying your code. Your trick is clever,   because it proves that the must

be a problem with either pbeta or qbeta. Your two lines:

    p(0) = pbeta(x1[i], s1, s2);

    q(0) = qbeta(p(0),s1,s2); 

means that q(0)=x1[i]. Hence the "likelihood contribution"      

ans -= dbeta(q(0),s1 , s2, true);

does not depend on the random effect u. This again means that the Laplace

approximation is exact (due the normal prior on u). This is of course

assuming infinite numerical precission.

What goes wrong in finite precission with the AD is not clear to me.

You are differenting a code for q(0) which, as a computer program depends on u,

but mathematically does not.  The solution must lie in the underlying

numerical procedure procedure for pbeta(), which I have not looked into.

It is hard to say if how difficult it is to fix.

Sorry for not  being able to come up with a more constructive answer.

Hans

[Bob Douma]

Hi Hans,

Thanks for the response. It is good that we have nailed down the problem to something more specific. 

It is ok that you cannot give a more constructive answer right now. I can imagine the pbeta is not a top-priority issue to solve either. For now, I know what I need to know: that a correlated mixed-effect model is identifiable and can be fitted. 

Let me know if you or others have worked on a solution - do you want me to address this as an issue on the github page?  And thanks for developing such a nice interface to AD. 

Cheers, Bob 

[Hans Skaug]

Hi,

Yes, you can file an issue, not so much under the expectation that someone will look into it soon,

but rather to warn other users that there may be a problem with using pbeta in combination with the Laplace approximation.

(I am not aware that this has been tested before). For the record,

if you look into the github-folder \adcomp\TMB\inst\include\tiny_ad\beta you will

see that TMB uses the same code to evaluate pbeta as R does. 

As for alternative ways of doing this, I do not see an easy fix. It would be nice

to have "copula library" in TMB. The copula density is a composite function (involving pbeta in your case),

for which one can calculate 1,2,3nd order derivatives by hand, and which can be feeded into TMB via "atomic functions".

(However, this requires understanding of how atomic functions work in TMB). This

may be more computationally effient and more numerically robust. However,

this is not to say that the issue with pbeta is not important.

Hans

 

On Monday, May 31, 2021 at 8:43:28 PM UTC+2 bob....@gmail.com wrote:

Hi Hans,

[Bob Douma]

Hi Hans,

I submitted my finding as an issue to TMB to warn future users. I fully agree that it would be nice to have a copula library. I am not sure if I can contribute to this function since I do not understand how the atomic functions work. Maybe I have some time later to dive into it. 

Bob 

--
To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to a topic in the Google Groups "TMB Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/tmb-users/NSqL2xVQJmw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to tmb-users+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/f01f2274-57b6-4f8d-9b7a-e604b1a4909en%40googlegroups.com.