Passing indices to excursions.inla

[Former user]

Hi,

This is a message that I posted in INLA discussion group. Finn suggested to left an issue here.

I would like to use the excursions.inla function, but the available documentation is scarce.

This is my model formula for estimating PM10 daily concentrations in Italy, for the month of January.

as.formula(lpm10~Intercept+dust+aod550.s+log.pbl00.s+log.pbl12.s+sp.s+t2m.s+tp.s+ptp.s+q_dem.s+i_surface.s+d_a1.s-1)->myformula

#spde with autoregressive AR(1) structure: update(myformula,.~.+f(id_centralina,model="iid")+f(i,model=spde,group = i.group,control.group = list(model="ar1",hyper=list(theta=theta_hyper))))->myformula

#spde definition inla.spde2.pcmatern(mesh=mesh,alpha=2,constr=FALSE,prior.range = c(150,0.8),prior.sigma = c(0.8,0.2))->spde

#n_giorni: the number of days in January, 1..31 inla.spde.make.index(name="i",n.spde=spde$n.spde,n.group = n_giorni)->iset

inla(myformula, data=inla.stack.data(mystack,spde=spde), family ="gaussian", verbose=TRUE, control.compute = list(openmp.strategy="pardiso.parallel",cpo=TRUE,waic=TRUE,dic=TRUE,config=TRUE), control.fixed = list(prec.intercept = 0.001, prec=1,mean.intercept=0), control.predictor =list(A=inla.stack.A(mystack),compute=TRUE) )->inla.out

Now, I have my inla.out object and I would like to find the excursion sets for one specific day of January, similarly to Cameletti's Piemonte example.

I supposed that the the "ind" argument is what I need (am I wrong?)

I tried something like this:

##if I want the excursions sets for the 1st of January or 2 for the second which(iset$i.group==1)->indici / which(iset$i.group==2)->indici

#This works excursions.inla(inla.out,stack = mystack,tag="training",method = "QC",type=">",u=log(50),ind = indici)->exc.out

#But if I do: which(iset$i.group==29)->indici #29 of January

I get:

Error in seq.default(xmin, xmax, len = n) : 'from' must be a finite number In addition: Warning messages: 1: In max(marginal$y) : no non-missing arguments to max; returning -Inf 2: In min(x, na.rm = na.rm) : no non-missing arguments to min; returning Inf 3: In max(x, na.rm = na.rm) : no non-missing arguments to max; returning -Inf 4: In min(m$x) : no non-missing arguments to min; returning Inf 5: In max(m$x) : no non-missing arguments to max; returning -Inf

Finally: if I run the same command without the "ind" argument , the commmand runs and terminates with no error, but then I am not able to retrieve the results for one specific day.

Any suggestion?

Regards Guido

[David Bolin]

Hi Guido,

The ind-argument is a bit tricky with the INLA interface, and I agree that there should be more documentation about it. I believe that the problem you are facing is that excursions.inla already extracts one part of the model through the "tag" argument, and because of this, it is assumed that "ind" contains indices within that component. So my guess is that which(iset$i.group==29) creates indices which are larger than the number of elements in the component that is extracted from the "training" part of the stack.

If you could create a reproducable example, or send me the inla object and the code you use I can take a look at it.