Using nimbleEcology to get MLE

[gesta...@gmail.com]

I am interested in using nimbleEcology to fit Nmixture models by MLE, and began playing around with a toy example to try to learn the ropes. A few questions cropped up for me.

1.      Does it matter if the model code includes priors, even though they are not used (I think). In the occupancy example at (https://cran.r-project.org/web/packages/nimbleEcology/vignettes/Introduction_to_nimbleEcology.html it looks like priors are included in the model.

2.      In my toy example, when I return coefficients for p (in the call to optim), or when I return the derived detection probabilities, everything looks OK. But when I try both, the parameter estimates for the coefficients are off, and the hessian has a bunch of zeroes. How can I return both? I fiddled around with the getDependencies function, but couldn’t get it to work the way I expected. I have the same question about the site-level abundances – what if I want to return those?

3.      Is it possible to include random effects, e.g., a random site effect for lambda or p?

 

I suspect the problem is just my lack of understanding about getDependencies or maybe even my inexperience with optim.

Code for the toy example is attached.

 

Thanks,

Glenn

[Benjamin R Goldstein]

Hi Glenn--

Thanks for the questions and glad to see your interest!

1. You should not include priors in the model code when doing MLE. They'll contribute to the log likelihood that's being optimized, which is not what you want. In the vignette we're writing models for MCMC, which we might want to make more explicit.

2. The problem I'm seeing in your code is that you're including "p" in your CalcLogLik function, and providing values for p separately from the values of b1 and b0, which determine p. The only parameters in your model are b0, b1, and lambda, so try just finding maximum likelihood estimates of those. When I fix that on my end, I get sensible estimates and SEs from your first example. In your version that breaks, the optimizer thinks that each "p" is a separate parameter to be estimated. If you want to estimate values of "p" at each site under the model you've written, you can calculate them using the data and the MLE estimates of b0 and b1, but you can't do derived quantities in an MLE setting the way you can in MCMC (and note that the error propagation for those quantities becomes more complicated too).

3. At present, the main way to include random effects when estimating a NIMBLE model with MLE is using the new Laplace approximation feature part of the NIMBLE 1.0.0 release. I'll admit I haven't yet used this new feature so I don't have advice handy, but the NIMBLE User Manual is the best place to start for that (Chapter 16.2).

Hope this helps,

Ben

--
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/1e081886-2260-4bfd-897d-b903b6c20261n%40googlegroups.com.

[gesta...@gmail.com]

OK, thanks Ben. I think that all makes sense. If I understand you correctly, I am already almost there, and all I need to do (ignoring for now the random effects bit!) is delete from the model the lines specifying priors, and then do the MLE fit that I already specify in my third try (where I specify only lambda, b0, and b1). Correct?

Something like:

nc3 <- nimbleCode({
for(j in 1:noccs){
logit(p[j]) <- b0 + b1*X[j]
}
for(i in 1:nSites){
y[i,1:noccs] ~ dNmixture_v(lambda, prob = p[1:noccs],
Nmin = -1, Nmax = -1, len = noccs)
}
})
nim <- nimbleModel(code=nc3,name='nim',constants=Consts,data=Data,inits=Inits,calculate=FALSE)
Cnim <- compileNimble(nim)

# fit by MLE
CalcLogLik <- nimbleFunction(
  setup = function(model, nodes)
    calcNodes <- model$getDependencies(nodes, self = FALSE),
  run = function(v = double(1)) {
    values(model, nodes) <<- v
    return(model$calculate(calcNodes))
    returnType(double(0))
  }
)

LogLik <- CalcLogLik(nim, c("b0","b1","lambda"))
CLogLik <- compileNimble(LogLik, project = nim)
out2 <- optim(c(0.5,0.1,2), CLogLik$run, control = list(fnscale = -1),hessian = TRUE)

[Benjamin R Goldstein]

What you've got looks good. One additional thing--the dynamic bounds (Nmin = -1 and Nmax = -1) should only be used in the MCMC context as they can introduce bias in the MLE; the bounds of the truncated infinite sum shouldn't change during MLE. Instead, you should set Nmin = 0 and set Nmax = K where K is some value much higher than the possible maximum abundance at a site. (This is the upper bound of the truncated infinite sum used in likelihood calculation.)

Ben

--
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/03dd4c45-94c2-4678-ba19-4ab627adf25an%40googlegroups.com.

[Wei Zhang]

Just a quick follow up on Ben's point 3 in the first email. MLE now can be done very easily using the Laplace functionality nimble 1.0.0+; see buildLaplace. You only need the model object for using Laplace. If there are no random effects in the model, Laplace does not involve any approximations; otherwise marginal likelihood of top-level model parameters is obtained via the Laplace approximation. Gradient of the loglikelihood will be obtained via automatic differentiation and used for optimization which should help for complicated models. Also standard errors or covariance matrix of the ML estimators can be returned for uncertainty measurement. 

Best wishes

Wei

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CAEbSdx1Zc1R%2BHaviyF0ZS%2BGq4D41EkXHAOODFLXqA-%3DpCiA4Dg%40mail.gmail.com.

[Perry de Valpine]

However, just to add to the discussion: we haven't yet released an update to nimbleEcology with support for AD in the marginal distributions provided in nimbleEcology. We will be doing that.

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CAL5Jjd%3DpTKaz00VQ8naMkxpvB1Q%3D%3Da32H9dav3pSN1G%3D_HhqkQ%40mail.gmail.com.

[Glenn Stauffer]

Thanks all!

I will take a look at chapter 16, and anticipate the support for nimbleEcology distributions.

Ben, I liked the idea of the dynamic bounds, but didn't think about implications for MLE. I've always before used nimble for mcmc.

Glenn