failing to fit random walk along year by location

[Thierry Onkelinx]

Dear all,

I'm trying to fit a dataset with multiple non-zero count data time series. About 250 locations, sampled once a month (6 winter months) during 30 years. About 34% of the design has a non-zero count. Every location has at least 5 years (8 observations) with non-zero counts, the median location as 18 years (49 observations) with non-zero counts.

The full model has a binomial component to fit the probability of non-zero counts. The idea is to fit a type 0 zero inflated distribution with a variable point mass of zero's.

The model (below the text) runs fine when I include a global trend (cyear), global seasonal effect (imonth), average location effect (location) and a change in season effect by year (imonth2).

I'd like to allow for a changing trend by location (cyear2). That model fails to fit. I get a bunch of "GMRFLib_2order_approx: rescue NAN/INF values in logl for idx=46787". After some time the model is killed.

Any suggestions?

Best regards,

Thierry

inla(

  count ~ 0 + intercept +
    f(
      cyear, model = "rw1", scale.model = TRUE,
      hyper = list(theta = list(prior = "pc.prec", param = c(1, 0.01)))
    ) +
    f(
      cyear2, model = "rw1", constr = TRUE, replicate = ilocation,
      hyper = list(theta = list(prior = "pc.prec", param = c(1, 0.01)))
    ) +
    f(
      imonth, model = "rw1", constr = TRUE,
      hyper = list(theta = list(prior = "pc.prec", param = c(0.5, 0.01)))
    ) +
    f(
      imonth2, model = "rw1", constr = TRUE, replicate = cyear,
      hyper = list(theta = list(prior = "pc.prec", param = c(0.5, 0.01)))
    ) +
    f(
      location, model = "iid", constr = TRUE,
      hyper = list(theta = list(prior = "pc.prec", param = c(1, 0.01)))
    ),
  data = truncated_count, family = "zeroinflatednbinomial0",
  control.compute = list(config = TRUE, waic = TRUE),
  control.predictor = list(link = 1),
  control.family = list(
    list(hyper = list(theta2 = list(initial = -20, fixed = TRUE)))
  )
)

ir. Thierry Onkelinx

Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry....@inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////

[Helpdesk (Haavard Rue)]

Somtimes these zero-inflated models can be a little tricky to get
started, is the likelihood can be not log-concave. you can force this
adding

control.inla=list(cmin=0)

if this run goes well, can try to remove it with a rerun like

r=inla(...,control.inla=list(cmin=0))
r$.args$control.inla$cmin = -Inf
r=inla.rerun(r)


btw, there is a new implementatino of the zeroinflated binomial and
poisson, which allow for more options, see

inla.doc("0poisson")

nbinom, can be achived adding an iid random effect

if this does not help, please share so I can rerun it here

thx
H



On Fri, 2023-05-19 at 15:28 +0200, 'Thierry Onkelinx' via R-inla

> --
> You received this message because you are subscribed to the Google
> Groups "R-inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to r-inla-discussion...@googlegroups.com.
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/CAJuCY5zMFa7RStcR1qH8S_O9Czvne_PTFjs_k66ryOjgfjmmKQ%40mail.gmail.com
> .

--
Håvard Rue
he...@r-inla.org

[Thierry Onkelinx]

Dear Havard,

control.inla=list(cmin=0) doesn't solve the problem. The issue is probably with the data from some locations. Some locations have overall moderate counts (tens of individuals) in combination with one or two observations with high counts (thousands of individuals).

Some subsets fit without issues, other subsets fit with some messages "GMRFLib_2order_approx: rescue NAN/INF values in logl for idx". The subsets with a large portion of the data yield lots of "GMRFLib_2order_approx: rescue NAN/INF values in logl for idx" and fail to fit.

How can I capture the "GMRFLib_2order_approx" messages from the output? They are not present in the logfile.

The 0poisson is not an option. I use a more complex structure for the binomial part of the model. A zero trunctaced poisson and negative binomial distribution would be a nice thing to have.

Best regards,

ir. Thierry Onkelinx

Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry....@inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////

Op za 20 mei 2023 om 10:28 schreef Helpdesk (Haavard Rue) <he...@r-inla.org>:

[Helpdesk (Haavard Rue)]

thx for the info.

zero truncated poisson you can get from zeroinflatedpoisson0 setting 'p'
to a fixed low value.

same for the nbinomial

[Sylvan Benaksas]

I just wanted to follow up on this, has there be any implementation or is this still the best advice for a truncated poisson and negative binomial,

Here is the code I am trying for zero truncated nbinomial,

inla(family = c("binomial", "zeroinflatednbinomial0"), #zeroinflatednbinomial0
                control.family = list(
                  list(),   
                  list(
                    hyper = list(
                      prob = list(
                        fixed = TRUE,
                        initial = -20)))), ...)

Many thanks

[Finn Lindgren]

From INLA 23.10.19-1, there is an "nzpoisson" for the zero-truncated Poisson; see
  https://inlabru-org.github.io/inlabru/articles/zip_zap_models.html

for an example in a hurdle model context.
I'm not aware of a Negative Binomial version (I don't see one in the c-code)

Finn



Finn

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/4083b2dc-4037-4934-a189-ab3bef5f6d1bn%40googlegroups.com.

--

Finn Lindgren
email: finn.l...@gmail.com

[Helpdesk (Haavard Rue)]

There is no corresponding one for the nbinomial, but 'nbinomial = poisson +
iid', approximately, so... If important it could be added, or added by user
using 'cloglik'

Håvard Rue
he...@r-inla.org