Predict on output of lgcp from inlabru gives unexpected result

[Virginia Morera Pujol]

I am running an model in inlabru, using the function lgcp with a SPDE field. This is what my data looks like. These animals never fly over land, which is why I have removed the outer boundary to the east and inside the islands.The mesh is generated with this code:


boundary_in <- inla.sp2segment(b_in)
boundary_out <- inla.sp2segment(b_out)
mesh <- inla.mesh.2d(boundary = list(boundary_in, boundary_out), # these are SpatialPolygons (attached)
                     max.edge = c(70, 200),                      # max edge in the inner boundary is approx. the expected range/5

                       cutoff = 10)                                       # any value larger than this and the boundaries of the islands start to lose too much its shape

And the model is run with this:

# specify the 2d matérn field for the random effect
matern2d <- inla.spde2.pcmatern(mesh,
                                prior.sigma = c(0.1, 0.01),  #specified that sd is unlikely to be larger than 0.1
                                prior.range = c(20, 0.1))    # we don't want our model picking up spatial structures at ranges smaller than 20 km.


# specify the components
cmp0 <- coordinates ~ mySPDE(map = coordinates, model = matern2d) +
                      Intercept


# fit the model
fit0 <- lgcp(components = cmp0,
             data = ds.df2,     # spatial points dataframe (attached)
             samplers = bound)  # same polygon as the inner boundary (also attached)

            

Everything seems to run according to plan but when I try to predict the intensity with

lmbd0 <- predict(fit0, pixels(mesh, mask = T), ~exp(Intercept + mySPDE))

and try to plot it using the ggplot functions

ggplot() + gg(lmbd0["median"]) + coord_equal()

I get the error message

Error in if ((w[1] * sm + w[2] * cm + w[3] * dm + w[4]) < best$score) break : missing value where TRUE/FALSE needed


Trying to plot it outside ggplot with

plot(lmbd0["median])
plot(bound, border = "white", add = T)

gets me what you can see in the attached image "intensity.jpeg", a completely flat surface with intensity = 0 and one single "pixel" of intensity 5e+232

and indeed,  summary(lmbd0$median) gives me this:

      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
 0.000e+00  0.000e+00  0.000e+00 4.847e+258  0.000e+00 5.387e+262 

I have tried changing the values for the SPDE priors for range and sd, and their limits, but I get the same problem.  I suspected it has something to do with my points laying too close tot the boundary, so I tried making another mesh with an inner boundary that had a buffer of 2 km around what you see, but the problem persisted... 

Any suggestions will be much appreciated! 

Virginia. 

[Finn Lindgren]

Meant to send the below reply to the list!

The gmap plot require the mesh to be constructed with
inla.mesh.2d(..., crs=CRS(proj4string(b_in)))
so that it knows how to convert it to longitude/latitude for plotting.

Finn

---------- Forwarded message ----------
From: Finn Lindgren <finn.l...@gmail.com>
Date: 8 May 2018 at 20:10
Subject: Re: [r-inla] Predict on output of lgcp from inlabru gives
unexpected result
To: Virginia Morera Pujol <morera....@gmail.com>


Hi,

I believe in your suspicion that it runs into problems due to the
geometry being tricky with so many points close to the boundary, but
I'll keep investigating why it fails so spectacularly.
Also, could you send me the example code for the 2km expanded buffer
as well? 2km however is a tiny buffer, since cutoff is 10 km. A 20km
buffer would be more reasonable (but see below about other model
options).

Predicting on the log-scale shows the problem more clearly:
lmbd1 <- predict(fit0, pixels(mesh), ~(Intercept + mySPDE))
ggplot() + gg(lmbd1["median"]) + coord_equal()
As you noted, there's a single hotspot near the coast where the model
breaks down. In the past, such hotspots were the cause of an invalid
numerical integration scheme in the lgcp likelihood approximation, but
the current method shouldn't behave like this.
Predicting log(abs(...)), just to see if there is any pattern hiding
underneath, does show a pattern that connects reasonably well to the
actual point pattern, so it seems it's just the one hotspot that cases
some numerical issue.

Since the current CRAN version, we have fixed a few minor bugs,
available in the development version; see
https://github.com/fbachl/inlabru/tree/devel for installation
instructions,
but since I was running the development version, this clearly did not
help in this case.

You should take a look at Haakon's "barrier" SPDE models, that deal
much better with this kind of situation; several recent threads
discuss those.
As far as I'm aware of, noone has used them in combination with
inlabru yet, but I'm keen to get those working, so if you try them and
run into problems, please let me know!
In principle, they should work seamlessly with inlabru::lgcp(), but
the inner workings of the inlabru&inla code may complain a bit...

Btw, you might find this code helpful (at least I did, since it helped
me verify that your coordinate system makes sense!)

inlabru::gmap(ds.df2) + gm(ds.df2,mapping=aes_(~x,~y)) + gm(mesh)

Resulting plot is attached.

Finn

On 8 May 2018 at 19:40, Finn Lindgren <finn.l...@gmail.com> wrote:
> Hi,
> turning on verbose output with options=list(verbose=TRUE) reveals that
> inla() is struggling with the optimisation due to NaN or Inf
> appearing; I'm not sure of the precise reason for that.
>
> Finn

>> --
>> 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 post to this group, send email to
>> r-inla-disc...@googlegroups.com.
>> Visit this group at https://groups.google.com/group/r-inla-discussion-group.
>> For more options, visit https://groups.google.com/d/optout.
>
>
>
> --
> Finn Lindgren
> email: finn.l...@gmail.com



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


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

[Finn Lindgren]

Hi,

I think I've identified the issue:
Since the mesh building simplifies the boundary, there is a mismatch
between the actual mesh boundary and the specified "samplers" domain,
which is the original, high-resolution boundary shape.
In principle, this shouldn't be an issue as long as the mismatch is
small, but it seems it results in some badly scaled integration
weights near the boundary; e.g., a too small integration weight in the
vicinity of an observed point should result in precisely the behaviour
we see, where the intensity is dragged towards infinity; it simply
isn't penalised in the approximate likelihood like it would be in the
true likelihood.

There are some related issues for calculating the integration weights
that interacts with this, making it difficult to debug, so a proper
solution may need to wait until I have more time to clean up the code
and make it more robust to edge cases like this (our original inlabru
package developer has moved on to work in a bank...).

The "barrier" models can have meshes that extend inland, which should
eliminate the problem even with the current inlabru integration code,
so I would suggest trying that, in particular since it removes the
edge effects that may also confuse your model, as it has an isolated
small "lake" in the southeast corner!

I'd still be interested in seeing the case with an extended domain you
mentioned, so I have more test cases to use for debugging.

Finn

On 8 May 2018 at 19:05, Virginia Morera Pujol <morera....@gmail.com> wrote:
>

[Finn Lindgren]

On 9 May 2018 at 00:52, Finn Lindgren <finn.l...@gmail.com> wrote:
> The "barrier" models can have meshes that extend inland, which should
> eliminate the problem even with the current inlabru integration code,
> so I would suggest trying that, in particular since it removes the
> edge effects that may also confuse your model, as it has an isolated
> small "lake" in the southeast corner!

Small caveat: those models should have _less_ issues of this type, but
currently, any model where the sampling domain doesn't align with the
mesh edges is liable to have this problem.
Alternative workaround is therefore to _first_ simplify the boundary
_before_ constructing the mesh, instead of letting the mesh builder do
it; that way, you have more control over the simplification, and can
supply the simplified boundary as the sampler parameter, which should
solve the problem.
inla.simplify.curve() may, or may not, be helpful for this (it doesn't
do exactly what we want in this case; it doesn't eliminate close
points)

Finn

[Virginia Morera Pujol]

Hello!

Thanks for all the suggestions, it seems I have a lot to do now!

The code I used to make the 2km buffer around the boundaries (which gives the same result) is:

b_in_2buff <- rgeos::gBuffer(b_in, width = 2)
b_out_2buff <- rgeos::gBuffer(b_out, width = 2)
boundary_in <- inla.sp2segment(b_in_2buff)
boundary_out <- inla.sp2segment(b_out_2buff)
mesh <- inla.mesh.2d(boundary = list(boundary_in, boundary_out), max.edge = c(70, 200), cutoff = 10)

I tried running it again with a buffer of 20 as you suggested, but I lose the islands completely (see attached image) and inla crashes (see other attached image) giving errors related to NAs and Infs as well:

GMRFLib version 3.0-0-snapshot, has recived error no [12] 
Reason : The Newton-Raphson optimizer did not converge 
Message : Condition `lambda < 1.0 / lambda_lim' is not TRUE 
Function : GMRFLib_init_GMRF_approximation_store__intern 
File : approx-inference.c Line : 2890 
RCSId : file: approx-inference.c hgid: f83f00640846 date: Tue May 08 10:16:23 2018 +0300

(All this is run with the developement version of inlabru now, and R version 3.5.0)

I suspect this might be because the boundary sort of "bends on itself" with this buffer, causing weird geometries.

I'll start trying to apply your suggestions and keep you posted of developments.

Thank you again!

Virginia Morera

PhD Student

Department of Animal Biology

University of Barcelona

Aquest correu electrònic i els annexos poden contenir informació confidencial o protegida legalment i està adreçat exclusivament a la persona o entitat destinatària. Si no sou  el destinatari final o la persona encarregada de rebre’l, no esteu autoritzat a llegir-lo, retenir-lo, modificar-lo, distribuir-lo, copiar-lo ni a revelar-ne el contingut. Si heu rebut aquest correu electrònic per error, us preguem que n’informeu al remitent i que elimineu del sistema el missatge i el material annex que pugui contenir. Gràcies per la vostra col·laboració.

Este correo electrónico y sus anexos pueden contener información confidencial o legalmente protegida y está exclusivamente dirigido a la persona o entidad destinataria. Si usted no es el destinatario final o la persona encargada de recibirlo, no está autorizado a leerlo, retenerlo, modificarlo, distribuirlo, copiarlo ni a revelar su contenido. Si ha recibido este mensaje electrónico por error, le rogamos que informe al remitente y elimine del sistema el mensaje y el material anexo que pueda contener. Gracias por su colaboración.

This email message and any documents attached to it may contain confidential or legally protected material and are intended solely for the use of the individual or organization to whom they are addressed. We remind you that if you are not the intended recipient of this email message or the person responsible for processing it, then you are not authorized to read, save, modify, send, copy or disclose any of its contents. If you have received this email message by mistake, we kindly ask you to inform the sender of this and to eliminate both the message and any attachments it carries from your account.Thank you for your collaboration.

[Finn Lindgren]

In the figure, there are some issues where the blue (inner boundary) runs close but not identical to the outer mesh boundary. What I meant was to only extend the outer boundary, not the inner one.

Finn

<inla_error.JPG>

<mesh_20buff.jpeg>

[Virginia Morera Pujol]

Ah.

I just did that, and the result is exactly the same as without the buffer (see image).

to build the mesh:

b_out_buff <- rgeos::gBuffer(b_out, width = 20)
boundary_in <- inla.sp2segment(b_in)
boundary_out <- inla.sp2segment(b_out_buff)

mesh <- inla.mesh.2d(boundary = list(boundary_in, boundary_out), max.edge = c(70, 200), cutoff = 10)

Virginia Morera

PhD Student

Department of Animal Biology

University of Barcelona

Aquest correu electrònic i els annexos poden contenir informació confidencial o protegida legalment i està adreçat exclusivament a la persona o entitat destinatària. Si no sou  el destinatari final o la persona encarregada de rebre’l, no esteu autoritzat a llegir-lo, retenir-lo, modificar-lo, distribuir-lo, copiar-lo ni a revelar-ne el contingut. Si heu rebut aquest correu electrònic per error, us preguem que n’informeu al remitent i que elimineu del sistema el missatge i el material annex que pugui contenir. Gràcies per la vostra col·laboració.

Este correo electrónico y sus anexos pueden contener información confidencial o legalmente protegida y está exclusivamente dirigido a la persona o entidad destinataria. Si usted no es el destinatario final o la persona encargada de recibirlo, no está autorizado a leerlo, retenerlo, modificarlo, distribuirlo, copiarlo ni a revelar su contenido. Si ha recibido este mensaje electrónico por error, le rogamos que informe al remitente y elimine del sistema el mensaje y el material anexo que pueda contener. Gracias por su colaboración.

This email message and any documents attached to it may contain confidential or legally protected material and are intended solely for the use of the individual or organization to whom they are addressed. We remind you that if you are not the intended recipient of this email message or the person responsible for processing it, then you are not authorized to read, save, modify, send, copy or disclose any of its contents. If you have received this email message by mistake, we kindly ask you to inform the sender of this and to eliminate both the message and any attachments it carries from your account.Thank you for your collaboration.

[Virginia Morera Pujol]

I have been trying Haakon's code for the Barrier model (script and data attached) and everything seems to run smoothly until I get to actually construct the model, using the new mesh I've created and the barrier.triangles object. When running

barrier.model = inla.barrier.pcmatern(mesh2, barrier.triangles = barrier.triangles,
                                      prior.range = c(6, 0.5),
                                      prior.sigma = c(3, 0.01))

I get the error message

Error in solve.default(t(Ts) %*% Ts) :
system is computationally singular: reciprocal condition number = 2.89737e-21

I am pretty sure it has to do with the polygon used for my mesh having self intersections (see code attached) but I don't know how to fix that. Some posts suggest running gBuffer with width = 0, but this just gets me the outer polygon, and I lose all the islands and shapes inside.

Virginia Morera

PhD Student

Department of Animal Biology

University of Barcelona

Aquest correu electrònic i els annexos poden contenir informació confidencial o protegida legalment i està adreçat exclusivament a la persona o entitat destinatària. Si no sou  el destinatari final o la persona encarregada de rebre’l, no esteu autoritzat a llegir-lo, retenir-lo, modificar-lo, distribuir-lo, copiar-lo ni a revelar-ne el contingut. Si heu rebut aquest correu electrònic per error, us preguem que n’informeu al remitent i que elimineu del sistema el missatge i el material annex que pugui contenir. Gràcies per la vostra col·laboració.

Este correo electrónico y sus anexos pueden contener información confidencial o legalmente protegida y está exclusivamente dirigido a la persona o entidad destinataria. Si usted no es el destinatario final o la persona encargada de recibirlo, no está autorizado a leerlo, retenerlo, modificarlo, distribuirlo, copiarlo ni a revelar su contenido. Si ha recibido este mensaje electrónico por error, le rogamos que informe al remitente y elimine del sistema el mensaje y el material anexo que pueda contener. Gracias por su colaboración.

This email message and any documents attached to it may contain confidential or legally protected material and are intended solely for the use of the individual or organization to whom they are addressed. We remind you that if you are not the intended recipient of this email message or the person responsible for processing it, then you are not authorized to read, save, modify, send, copy or disclose any of its contents. If you have received this email message by mistake, we kindly ask you to inform the sender of this and to eliminate both the message and any attachments it carries from your account.Thank you for your collaboration.

[Haakon Bakka]

As an aside: With the newest version of INLA, you can run (assuming you constructed the mesh from a polygon poly.water):

water.tri = inla.over_sp_mesh(poly.water, y = mesh, type = "centroid", ignore.CRS = T)

poly.water2 = inla.barrier.polygon(mesh, barrier.triangles = barrier.tri)

The polygon you are actually using, when you are modeling the data is poly.water2. You may want to use this for plotting and determining in-out points, rather than using your original polygon. I do.

If you do the mesh for a stationary model or a barrier model, the problems you have encountered probably will not lead to any serious issues, but you may want to resolve them anyway.

Additional tips can be found at

https://haakonbakka.bitbucket.io/btopic104.html

and

https://haakonbakka.bitbucket.io/btopic127.html

Best,

Haakon

--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@googlegroups.com.

[Haakon Bakka]

And for your new question, I will run your code and check.

- Haakon

--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@googlegroups.com.

[Finn Lindgren]

I believe that generated polygon is what could/should be used a samplers polygon in inlabru as well.

Finn

[Haakon Bakka]

Small error in the sample code. Should have been

water.tri = inla.over_sp_mesh(poly.water, y = mesh, type = "centroid", ignore.CRS = T)

poly.water2 = inla.barrier.polygon(mesh, barrier.triangles = water.tri)

To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsubscr...@googlegroups.com.

[Virginia Morera Pujol]

Yep, tried that. the error persists:

Error in gUnaryUnion(mesh.polys) :
  TopologyException: Input geom 0 is invalid: Self-intersection at or near point 7103.8250221099997 1307.92857018 at 7103.8250221099997 1307.92857018

Virginia Morera

PhD Student

Department of Animal Biology

University of Barcelona

Aquest correu electrònic i els annexos poden contenir informació confidencial o protegida legalment i està adreçat exclusivament a la persona o entitat destinatària. Si no sou  el destinatari final o la persona encarregada de rebre’l, no esteu autoritzat a llegir-lo, retenir-lo, modificar-lo, distribuir-lo, copiar-lo ni a revelar-ne el contingut. Si heu rebut aquest correu electrònic per error, us preguem que n’informeu al remitent i que elimineu del sistema el missatge i el material annex que pugui contenir. Gràcies per la vostra col·laboració.

Este correo electrónico y sus anexos pueden contener información confidencial o legalmente protegida y está exclusivamente dirigido a la persona o entidad destinataria. Si usted no es el destinatario final o la persona encargada de recibirlo, no está autorizado a leerlo, retenerlo, modificarlo, distribuirlo, copiarlo ni a revelar su contenido. Si ha recibido este mensaje electrónico por error, le rogamos que informe al remitente y elimine del sistema el mensaje y el material anexo que pueda contener. Gracias por su colaboración.

This email message and any documents attached to it may contain confidential or legally protected material and are intended solely for the use of the individual or organization to whom they are addressed. We remind you that if you are not the intended recipient of this email message or the person responsible for processing it, then you are not authorized to read, save, modify, send, copy or disclose any of its contents. If you have received this email message by mistake, we kindly ask you to inform the sender of this and to eliminate both the message and any attachments it carries from your account.Thank you for your collaboration.

[Haakon Bakka]

Hi Virginia,

I actually did not get your error when I ran the code. I was also able to plot the poly.barrier.

I think it is the self-intersection that is the problem (also causing the next error). Could you try some of the tricks at https://haakonbakka.bitbucket.io/btopic127.html ?

Here, we use gBuffer and gIntersection and gSimplify. Maybe if you have a more complex polygon you can make sure that one is more well-defined before you make the mesh?

PS: If you plot plot(mesh2), you see that it looks quite bad. Could you try to get a mesh that looks a bit nicer (as plot(mesh) does)?

Best,

Haakon

--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@googlegroups.com.

[Finn Lindgren]

gSimplify uses a version of the same algorithm as inla.siplify.curve that I mentioned previously, and hence is liable to the same issue: it tries to preserve non-linear/curved features while removing points that fall on straight lines, but what we would like is a method that _smooths_ the curve. I’m not aware of a good function for that.

Finn

[Virginia Morera Pujol]

Ok,

I did manage to run the rest of the code to set the barrier model by increasing the tolerance in gSimplify to 10 (which makes sense since 10 is the value for the cutoff in the mesh anyway).

I redid the mesh with the same edge parameters as the mesh1 (see attached).

However, I am not sure how to translate Haakon's code from INLA into inlabru. I tried this:

barrier.model = inla.barrier.pcmatern(mesh2, barrier.triangles = barrier.triangles,
                                      prior.range = c(6, 0.5),
                                      prior.sigma = c(3, 0.01))

# specify the components
cmp0 <- coordinates ~ myBARRIER(map = coordinates, model = barrier.model) + Intercept

# fit the model
fit0 <- lgcp(components = cmp0, # "model formula"
             data = ds.df2,
             samplers = new_sp,#
             options=list(verbose=TRUE))

And it failed pretty spectacularly, fatal error included (see attached image). This is the text output I got.

 

Warning message:
package "methods" in options("defaultPackages") was not found 
During startup - Warning messages:
1: package 'datasets' in options("defaultPackages") was not found 
2: package 'utils' in options("defaultPackages") was not found 
3: package 'grDevices' in options("defaultPackages") was not found 
4: package 'graphics' in options("defaultPackages") was not found 
5: package 'stats' in options("defaultPackages") was not found 
6: package 'methods' in options("defaultPackages") was not found 
Error in sort.int(x, na.last = na.last, decreasing = decreasing, ...) : 
  there is no .Internal function 'sorted_fpass'
Error in (function (data, model, stackmaker, n = 10, result = NULL, family,  : 
  INLA returned message: The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
  If this does not help, please contact the developers at <he...@r-inla.org>.
In addition: Warning message:
In inla.inlaprogram.has.crashed() :
  The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
  If this does not help, please contact the developers at <he...@r-inla.org>.

So I figure I probably did not specify the model in a way inlabru understands...

Virginia Morera

PhD Student

Department of Animal Biology

University of Barcelona

Aquest correu electrònic i els annexos poden contenir informació confidencial o protegida legalment i està adreçat exclusivament a la persona o entitat destinatària. Si no sou  el destinatari final o la persona encarregada de rebre’l, no esteu autoritzat a llegir-lo, retenir-lo, modificar-lo, distribuir-lo, copiar-lo ni a revelar-ne el contingut. Si heu rebut aquest correu electrònic per error, us preguem que n’informeu al remitent i que elimineu del sistema el missatge i el material annex que pugui contenir. Gràcies per la vostra col·laboració.

Este correo electrónico y sus anexos pueden contener información confidencial o legalmente protegida y está exclusivamente dirigido a la persona o entidad destinataria. Si usted no es el destinatario final o la persona encargada de recibirlo, no está autorizado a leerlo, retenerlo, modificarlo, distribuirlo, copiarlo ni a revelar su contenido. Si ha recibido este mensaje electrónico por error, le rogamos que informe al remitente y elimine del sistema el mensaje y el material anexo que pueda contener. Gracias por su colaboración.

This email message and any documents attached to it may contain confidential or legally protected material and are intended solely for the use of the individual or organization to whom they are addressed. We remind you that if you are not the intended recipient of this email message or the person responsible for processing it, then you are not authorized to read, save, modify, send, copy or disclose any of its contents. If you have received this email message by mistake, we kindly ask you to inform the sender of this and to eliminate both the message and any attachments it carries from your account.Thank you for your collaboration.

To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsubscr...@googlegroups.com.

[Finn Lindgren]

Hi, that’s indeed a spectacular failure message, which isn’t only due to inlabru, but likely also issues with R 3.5.0...  but I do see some things that inlabru will definitely need: for regular spde models, it fins the mesh in the spde object, but since the barrier model is a different kind of model object that doesn’t work. You have to supply it manually to lgcp(), I believe as domain=mesh.

Can you send me the updated code so I can check on R 3.4.2?

Finn

<inla_fatal_error.JPG>

<meshes.jpeg>

[Virginia Morera Pujol]

New code and data file attached.

Running it with domain = mesh2 in the lgcp call gives me the error

Error in ipoints(NULL, domain[[nm]], name = nm) : object 'ips' not found

Virginia Morera

PhD Student

Department of Animal Biology

University of Barcelona

Aquest correu electrònic i els annexos poden contenir informació confidencial o protegida legalment i està adreçat exclusivament a la persona o entitat destinatària. Si no sou  el destinatari final o la persona encarregada de rebre’l, no esteu autoritzat a llegir-lo, retenir-lo, modificar-lo, distribuir-lo, copiar-lo ni a revelar-ne el contingut. Si heu rebut aquest correu electrònic per error, us preguem que n’informeu al remitent i que elimineu del sistema el missatge i el material annex que pugui contenir. Gràcies per la vostra col·laboració.

Este correo electrónico y sus anexos pueden contener información confidencial o legalmente protegida y está exclusivamente dirigido a la persona o entidad destinataria. Si usted no es el destinatario final o la persona encargada de recibirlo, no está autorizado a leerlo, retenerlo, modificarlo, distribuirlo, copiarlo ni a revelar su contenido. Si ha recibido este mensaje electrónico por error, le rogamos que informe al remitente y elimine del sistema el mensaje y el material anexo que pueda contener. Gracias por su colaboración.

This email message and any documents attached to it may contain confidential or legally protected material and are intended solely for the use of the individual or organization to whom they are addressed. We remind you that if you are not the intended recipient of this email message or the person responsible for processing it, then you are not authorized to read, save, modify, send, copy or disclose any of its contents. If you have received this email message by mistake, we kindly ask you to inform the sender of this and to eliminate both the message and any attachments it carries from your account.Thank you for your collaboration.

[Finn Lindgren]

Thanks!

the lgcp parameter should be
domain = list(coordinates = mesh2)
instead of
domain = mesh2

With that correction, it runs with no errors on R 3.4.2; the missing
"methods" etc warning I belive is entirely unrelated, and is issue
with either 3.5.0 or how you run the code; the command line Rscript
method used to have this problem, where it didn't load the core R
packages, to which methods and the other listed packages in your
warning message belong. I thought they'd fixed that problem though.
How did you run the code when that warning came up? Can you run
something much more simple and see if you get the same warning (i.e.
something that doesn't involve INLA or inlabru at all).

Haakon, how should one interpret the Theta1 and Theta2 parameters for
the barrier model output? Since it's not stored as an inla.spde
object, inla.spde.result won't work, and for regular pcmatern models
inla converts the output to range and sigma.

The summary output indicates that it managed to make an estimation,
but there are clearly issues with it.
I wonder if the covering square that shows up in new_sp is causing
problems by confusing the samplers setup, see attached figure of
new_sp and the estimated myBARRIER model component.

The inlabru predict() function doesn't manage to do a spatial
prediction (just comes out as constant), so I think there is still
some code in inlabru that only does the right thing for ordinary spde
model objects; the main issue tends to be how and where it needs to
know about the mesh, which isn't present in the barrier.model object
(at least not where inlabru would expect it to be).

Finn

>>>> an email to r-inla-discussion...@googlegroups.com.

>>>> To post to this group, send email to

>>>> r-inla-disc...@googlegroups.com.

>>>> Visit this group at
>>>> https://groups.google.com/group/r-inla-discussion-group.
>>>> For more options, visit https://groups.google.com/d/optout.
>>>
>>>
>>
>> <inla_fatal_error.JPG>
>>
>> <meshes.jpeg>
>
>

[Virginia Morera Pujol]

I feel like we're slowly (but surely) getting there!

The issue was indeed with R version 3.5.0. I did the following tests

- R Version 3.5.0 worked perfectly with other not-inla stuff, and also with very simple inlabru models like the ones you sent us to test before the workshop (attached script)

- R Version 3.4.4 works perfectly to run the barrier model with the script you have sent me, and I manage to replicate your results.

I'll try and construct a mesh without that square and see what changes.

Virginia Morera

PhD Student

Department of Animal Biology

University of Barcelona

Aquest correu electrònic i els annexos poden contenir informació confidencial o protegida legalment i està adreçat exclusivament a la persona o entitat destinatària. Si no sou  el destinatari final o la persona encarregada de rebre’l, no esteu autoritzat a llegir-lo, retenir-lo, modificar-lo, distribuir-lo, copiar-lo ni a revelar-ne el contingut. Si heu rebut aquest correu electrònic per error, us preguem que n’informeu al remitent i que elimineu del sistema el missatge i el material annex que pugui contenir. Gràcies per la vostra col·laboració.

Este correo electrónico y sus anexos pueden contener información confidencial o legalmente protegida y está exclusivamente dirigido a la persona o entidad destinataria. Si usted no es el destinatario final o la persona encargada de recibirlo, no está autorizado a leerlo, retenerlo, modificarlo, distribuirlo, copiarlo ni a revelar su contenido. Si ha recibido este mensaje electrónico por error, le rogamos que informe al remitente y elimine del sistema el mensaje y el material anexo que pueda contener. Gracias por su colaboración.

This email message and any documents attached to it may contain confidential or legally protected material and are intended solely for the use of the individual or organization to whom they are addressed. We remind you that if you are not the intended recipient of this email message or the person responsible for processing it, then you are not authorized to read, save, modify, send, copy or disclose any of its contents. If you have received this email message by mistake, we kindly ask you to inform the sender of this and to eliminate both the message and any attachments it carries from your account.Thank you for your collaboration.

>>>> an email to r-inla-discussion-group+unsub...@googlegroups.com.

>>>> To post to this group, send email to

>>>> r-inla-discussion-group@googlegroups.com.

>>>> Visit this group at
>>>> https://groups.google.com/group/r-inla-discussion-group.
>>>> For more options, visit https://groups.google.com/d/optout.
>>>
>>>
>>
>> <inla_fatal_error.JPG>
>>
>> <meshes.jpeg>
>
>

[Finn Lindgren]

Note: the square wasn't in the mesh, just in the new_sp object.
Finn

>> >>>> an email to r-inla-discussion...@googlegroups.com.

>> >>>> To post to this group, send email to

>> >>>> r-inla-disc...@googlegroups.com.

[Finn Lindgren]

Some more clues:

The new_sp object seems to have two components; one rectangle, and one
set of polygons that indicate what is _land_. When used as "samplers"
information, this means that lgcp will try to define the lgcp
likelihood across anywhere that isn't land, which includes part of the
mesh domain that lies outside the study domain. You'll need to
construct a samplers polygon that covers precisely the study domain;
it should be possible to get that as the difference between the
rectangle and the land polygons (that are hiding in new_sp) somehow.

You also need to supply mesh=mesh2 to the myBARRIER component definition.

Unfortunately, there is some code that then expect this to be a
regular inla.sdpe2 model object, and fails because of that. I'll keep
tracing that issue.

Finn

[Virginia Morera Pujol]

Quick note:

I constructed a polygon without the outside square to use as boundary and it does indeed change the results.

See attachments:

boundary.Rdata contains the new polygon without the rectangle

mesh_barrier_model_v2.R contains the bit of code to construct such polygon, but it's just a gDifference

run_barrier_model_v2.R has the code to actually run the model with the updated boundary for "samplers"

Virginia Morera

PhD Student

Department of Animal Biology

University of Barcelona

Aquest correu electrònic i els annexos poden contenir informació confidencial o protegida legalment i està adreçat exclusivament a la persona o entitat destinatària. Si no sou  el destinatari final o la persona encarregada de rebre’l, no esteu autoritzat a llegir-lo, retenir-lo, modificar-lo, distribuir-lo, copiar-lo ni a revelar-ne el contingut. Si heu rebut aquest correu electrònic per error, us preguem que n’informeu al remitent i que elimineu del sistema el missatge i el material annex que pugui contenir. Gràcies per la vostra col·laboració.

Este correo electrónico y sus anexos pueden contener información confidencial o legalmente protegida y está exclusivamente dirigido a la persona o entidad destinataria. Si usted no es el destinatario final o la persona encargada de recibirlo, no está autorizado a leerlo, retenerlo, modificarlo, distribuirlo, copiarlo ni a revelar su contenido. Si ha recibido este mensaje electrónico por error, le rogamos que informe al remitente y elimine del sistema el mensaje y el material anexo que pueda contener. Gracias por su colaboración.

This email message and any documents attached to it may contain confidential or legally protected material and are intended solely for the use of the individual or organization to whom they are addressed. We remind you that if you are not the intended recipient of this email message or the person responsible for processing it, then you are not authorized to read, save, modify, send, copy or disclose any of its contents. If you have received this email message by mistake, we kindly ask you to inform the sender of this and to eliminate both the message and any attachments it carries from your account.Thank you for your collaboration.

>> >>>> an email to r-inla-discussion-group+unsub...@googlegroups.com.

>> >>>> To post to this group, send email to

>> >>>> r-inla-discussion-group@googlegroups.com.

[Finn Lindgren]

Ah, just what was needed. :-)
There will still be some issues though; even if the code might run
without mesh=mesh2, I don't think it will compute the right thing
anyway, and it fails with it there due to incorrect assumptions
elsewhere in the code.
I'll have a look after dinner...

Finn

>> >> >>>> an email to r-inla-discussion...@googlegroups.com.

>> >> >>>> To post to this group, send email to

>> >> >>>> r-inla-disc...@googlegroups.com.

[Finn Lindgren]

I've been tracing parts of the inlabru code, and there are a couple of
issues with incorrect assumptions that makes the model not do quite
the right thing, including not allowing the predict() method to do its
job properly, and handling the coordinate systems in an unforseen way
(it wants straight lines to be great circles, which deforms the
vertical line of the "bound" polygon when constructing the
likelihood).

The now too coarse domain boundary leaves some of the points outside
of the allowed domain, and it appears the code doesn't check for
that.. At first I though this is what cased the two artefacts in the
results, but that turns out not to be the case, because those effects
are still there when I remove the stray points.

Raw output plot, ignoring the Intercept:
ggplot() + gg(mesh2, color =
exp(fit0$summary.random$myBARRIER$`0.5quant`)) + gg(bound) +
coord_equal() ## check the rectangle is not there

This shows more of what actually happens within the study domain (with
predict(), we would be able to mask it out more easily):
ggplot() + gg(mesh2, color = exp(pmin(0,
fit0$summary.random$myBARRIER$`0.5quant`))) + gg(bound) +
coord_equal() ## check the rectangle is not there
It still doesn't look like a good fit to the data, but we're
definitely getting closer;
the barrier is clearly doing its job of decoupling at least point on
the mainland from the model on the ocean. The large spatial estimated
correlation range makes it hard to see what effect the islands have;
now that we know what we're looking for in terms of boundary curves
etc, you could start increasing the detail on the bound again.

I'll look into updating the parts of the inlabru code that aren't
quite doing the right thing, but I can't make any promises regarding
when; I've already spent more time on this than I really had time for
this week.

Finn

[Finn Lindgren]

On 9 May 2018 at 21:16, Finn Lindgren <finn.l...@gmail.com> wrote:
> (it wants straight lines to be great circles, which deforms the
> vertical line of the "bound" polygon when constructing the
> likelihood).

This issue can be worked around by doing the _opposite_ of what
gSimpify does: instead of removing points that lie along straight
lines, _add_ points at regular intervals, so that you have more
control over the true shape of the polygon; a straight line in your
coordinate system is _not_ a straight line on the sphere (a great
circle) and _not_ a straight line in other projections used internally
by inlabru to construct an integration scheme in square km units.
Another possible option is to remove the crs information on all the
objects. I believe that inlabru will then simply use the raw
coordinate system, without trying to be clever.

I've attached my latest testing code.

Finn

[Finn Lindgren]

One more test attached; swapping back to a regular pcmatern model, but
ignoring the sampling domain both when constructing the model and when
analysing the results.
It also ignores all CRS information. I also changed the prior for
the range; 50% chance of having a range smaller than 6 km was very
unreasonable; the estimate is 700 km now.

This is a sanity check that shows that the mesh is sensible, and that
the original hotspot artefact is not present in this model, and that
the spde models should in principle give a reasonable result.
In fact, you could use this approach while waiting for me to fix the
inlabru bugs (and also possible barrier model bugs!). You just need to
be careful with the study domain boundary and mesh construction;
you'll get some "pull towards zero" effect near the coast, but as a
temporary result it's perhaps still useful for you.

Finn

[Haakon Bakka]

Hi,

Reply to: Haakon, how should one interpret the Theta1 and Theta2 parameters for

the barrier model output? 

The first is log(sigma) the second is log(range), so use inla.tmarginal to transform them with exponentials.

Best,

Haakon

>>>> an email to r-inla-discussion-group+unsub...@googlegroups.com.

>>>> To post to this group, send email to

>>>> r-inla-discussion-group@googlegroups.com.

>>>> Visit this group at
>>>> https://groups.google.com/group/r-inla-discussion-group.
>>>> For more options, visit https://groups.google.com/d/optout.
>>>
>>>
>>
>> <inla_fatal_error.JPG>
>>
>> <meshes.jpeg>
>
>

[Virginia Morera Pujol]

Hi Finn,

I increased gradually the complexity of the boundary (decreasing the cutoff value and the tolerance in the gsimplify function) and tested the output from the barrier model. The three images attached go from very simple boundary (cutoff = 10, the one we were using) to the more complex one (cutoff = 1). The effect seems to concentrate towards the point near the boundary were the initial model failed. 

I will keep working with the normal SPDE model using this boundary (cutoff = 10) and ignoring the sampling domain as you suggested, as I still have a long way to go until I can run the definitive models... But please keep me posted of developments re using barrier models with inlabru or fixing the instabilities near the boundary. 

Thank you, both, for your invaluable help!

Virginia Morera

PhD Student

Department of Animal Biology

University of Barcelona

Aquest correu electrònic i els annexos poden contenir informació confidencial o protegida legalment i està adreçat exclusivament a la persona o entitat destinatària. Si no sou  el destinatari final o la persona encarregada de rebre’l, no esteu autoritzat a llegir-lo, retenir-lo, modificar-lo, distribuir-lo, copiar-lo ni a revelar-ne el contingut. Si heu rebut aquest correu electrònic per error, us preguem que n’informeu al remitent i que elimineu del sistema el missatge i el material annex que pugui contenir. Gràcies per la vostra col·laboració.

Este correo electrónico y sus anexos pueden contener información confidencial o legalmente protegida y está exclusivamente dirigido a la persona o entidad destinataria. Si usted no es el destinatario final o la persona encargada de recibirlo, no está autorizado a leerlo, retenerlo, modificarlo, distribuirlo, copiarlo ni a revelar su contenido. Si ha recibido este mensaje electrónico por error, le rogamos que informe al remitente y elimine del sistema el mensaje y el material anexo que pueda contener. Gracias por su colaboración.

This email message and any documents attached to it may contain confidential or legally protected material and are intended solely for the use of the individual or organization to whom they are addressed. We remind you that if you are not the intended recipient of this email message or the person responsible for processing it, then you are not authorized to read, save, modify, send, copy or disclose any of its contents. If you have received this email message by mistake, we kindly ask you to inform the sender of this and to eliminate both the message and any attachments it carries from your account.Thank you for your collaboration.

>>>> an email to r-inla-discussion-group+unsubscr...@googlegroups.com.

>>>> To post to this group, send email to
>>>> r-inla-discussion-group@googlegroups.com.
>>>> Visit this group at
>>>> https://groups.google.com/group/r-inla-discussion-group.
>>>> For more options, visit https://groups.google.com/d/optout.
>>>
>>>
>>
>> <inla_fatal_error.JPG>
>>
>> <meshes.jpeg>
>
>

[Virginia Morera Pujol]

Quick update:

I have been trying to add covariates to the model using exactly the same data, mesh, and pcmatern model for the spatial field, and using 1d matérn models for the covariates (since their are not "linear" to the intensity) and for some of the covariates I run into "Variable evaluates to NAN or INF" messages again, and the predicted intensity looks suspiciously "concentrated" around the same spot that was causing trouble before. I have attached images of the covariate itself (with meshpoints overlaid) and the predicted intensity (and sd) but I can send code and data if you think it's worth a look.

Virginia Morera

PhD Student

Department of Animal Biology

University of Barcelona

Aquest correu electrònic i els annexos poden contenir informació confidencial o protegida legalment i està adreçat exclusivament a la persona o entitat destinatària. Si no sou  el destinatari final o la persona encarregada de rebre’l, no esteu autoritzat a llegir-lo, retenir-lo, modificar-lo, distribuir-lo, copiar-lo ni a revelar-ne el contingut. Si heu rebut aquest correu electrònic per error, us preguem que n’informeu al remitent i que elimineu del sistema el missatge i el material annex que pugui contenir. Gràcies per la vostra col·laboració.

Este correo electrónico y sus anexos pueden contener información confidencial o legalmente protegida y está exclusivamente dirigido a la persona o entidad destinataria. Si usted no es el destinatario final o la persona encargada de recibirlo, no está autorizado a leerlo, retenerlo, modificarlo, distribuirlo, copiarlo ni a revelar su contenido. Si ha recibido este mensaje electrónico por error, le rogamos que informe al remitente y elimine del sistema el mensaje y el material anexo que pueda contener. Gracias por su colaboración.

This email message and any documents attached to it may contain confidential or legally protected material and are intended solely for the use of the individual or organization to whom they are addressed. We remind you that if you are not the intended recipient of this email message or the person responsible for processing it, then you are not authorized to read, save, modify, send, copy or disclose any of its contents. If you have received this email message by mistake, we kindly ask you to inform the sender of this and to eliminate both the message and any attachments it carries from your account.Thank you for your collaboration.