Re: [r-inla] Error when fitting no nugget SPDE

[Finn Lindgren]

Hi Joshua,

There's been a small misunderstanding about the mesh; the comment about the data locations was that to avoid some numerical issues with large precision nugget, they should be located at mesh nodes, but _not_ that only such mesh nodes were allowed.  On the contrary, you almost certainly need more mesh nodes to get a good approximation to the continuous model. Also, having only the observations as mesh nodes can _also_ lead to bad numerics due to badly shaped triangulation.

Secondly, in principle, if the nugget, as seen by the discretised spde model, really is small, shouldn't a "free" estimation pick that up?  The numerics issue related to mesh nodes in relation to data locations is related to is; the sum of discretised model and nugget may be a better approximation to the continuous model that the discretised model on its own.

Finally; the prior specifications often confuse me; double-check the PDF documentation for setting the nugget value; some model parameters are set in logscale, others not.

Finn

On 7 Apr 2014, at 22:59, Joshua French <joshua...@gmail.com> wrote:

Hello,

I am working on fitting a SPDE model with (essentially) no nugget.  Based on discussion in a previous thread, this can be accomplished by specifying a large, fixed precision for the nugget.  It was also noted that this could be numerically unstable unless the mesh was only included the observed data locations.  I have created a model to do this (with the mesh created using inla.delaunay).  Unfortunately, INLA is still crashing with the following error:

Note: method with signature ‘TsparseMatrix#Matrix#missing#replValue’ chosen for function ‘[<-’,
 target signature ‘dgTMatrix#ngCMatrix#missing#numeric’.
 "Matrix#nsparseMatrix#missing#replValue" would also be valid


    GMRFLib version 3.0-0-snapshot, has recived error no [3]
    Reason    : Matrix is (numerical) singular
    Function  : GMRFLib_comp_posdef_inverse
    File      : lapack-interface.c
    Line      : 71
    RCSId     : file: lapack-interface.c  hgid: eac68d10cc0f  date: Mon Feb 10 14:18:28 2014 +0100

/bin/sh: line 1: 34776 Abort trap: 6           '/Library/Frameworks/R.framework/Versions/3.0/Resources/library/INLA/bin/mac/64bit/inla' -b -s -v '/var/folders/x7/cm9jq3nn6mxdn31r9sw3sbv80000gn/T//Rtmp0pohci/file87d26c27d5f2/Model.ini' > /var/folders/x7/cm9jq3nn6mxdn31r9sw3sbv80000gn/T//Rtmp0pohci/file87d26c27d5f2/Logfile.txt
Error 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 help; please contact the developers at <he...@r-inla.org>.

I've tried scaling the coords/pcoords, changing the spde model (including alpha), but the problem persists.  The model runs if I change the control.family parameter to control.family=list(hyper=list(theta=list(initial=20, fixed=TRUE)) in the inla function (see p. 70 of the spde tutorial) which apparently also fixes the nugget precision to be a large value, but my post model checks on the predicted values are very bad (around 7% instead of 95% coverage), so I'm not suspicious about whether this is correct.

Does anyone here have suggestions for how I might correct this problem, fit the model another way, etc.?  Reproducible code for the problem is attached and provided below. 

library(INLA)

# (approximately) center coordinates
y = c(-1.523756,-2.22315,6.306653,-3.370748,3.157801,2.340856,7.686473,0.2584694,7.636864,-1.733082,2.893784,8.303231,2.698877,3.329094,6.857603,9.002041,-0.9022448,-0.05520936,2.502314,5.874786,4.692085,8.73891,4.202348,7.754542,2.94552,-0.886113,7.841755,8.476143,3.605204,2.705603,3.696793,8.136205,0.07849022,10.65243,5.662148,11.3402,-2.672724,8.99015,10.64247,-2.060972,-0.4470968,10.00641,8.573999,7.037847,7.308873,3.456908,6.926109,11.85081,0.3283561,2.568085,0.5851343,-0.783349,2.437771,2.174987,5.163221,0.09572735,-4.131473,-1.019985,8.216652,-0.9120677,1.386178,4.804644,9.977686,8.664915,4.104794,4.001969,5.571503,1.422291,-1.016264,5.529922,2.66573,0.3342805,4.807335,9.998714,5.448492,8.101185,8.026452,7.247551,8.773828,-2.183535,5.357466,5.746667,7.746981,11.13623,-4.302845,6.744897,7.680085,2.163964,-0.7716895,6.438806,6.720176,-2.824282,2.329994,2.190169,10.25891,3.372272,4.533381,-1.655705,-3.852587,7.3127)

set.seed(100)
coords <- matrix(runif(2 * 100), ncol = 2);
coords1 = coords[,1]; coords2 = coords[,2]

g = seq(0.0125, 0.9875, length = 40)
pgrid = expand.grid(g, g)
pcoords = as.matrix(pgrid)
pcoords1 = pgrid[,1]; pcoords2 = pgrid[,2]

# create mesh with observed data locations only
mesh = inla.delaunay(coords)

# create A matrices for observed and prediction coordinates
A = inla.spde.make.A(mesh, loc = coords)
Ap = inla.spde.make.A(mesh, loc = pcoords)

# create spde
spde = inla.spde2.matern(mesh, alpha = 1.5)

# stack objects for inla function
stk.e = inla.stack(tag='est', A = list(A, 1),
    data = list(resp = y),
    effects = list(s = 1:spde$n.spde,
    data.frame(b0=1, x1 = coords1, x2 = coords2)))

stk.pred = inla.stack(tag='pred', A = list(Ap, 1),
    data = list(resp = NA),
    effects = list(s=1:spde$n.spde,
    data.frame(b0=1, x1=pcoords1, x2=pcoords2)))

stk.full = inla.stack(stk.e, stk.pred)

# formula for model
form = resp ~ 0 + b0 + x1 + x2 + f(s, model=spde)

p.res = inla(form, data = inla.stack.data(stk.full),
             control.predictor = list(A=inla.stack.A(stk.full),
             compute=TRUE, cdf=c(0)),
             control.compute = list(config = TRUE),
             control.family=list(hyper=list(prec=list(initial=1e2, fixed=TRUE)))
             )

p.res2 = inla(form, data = inla.stack.data(stk.full),
             control.predictor = list(A=inla.stack.A(stk.full),
             compute=TRUE, cdf=c(0)),
             control.compute = list(config = TRUE),
             control.family=list(hyper=list(theta=list(initial=20, fixed=TRUE)))
             )

--
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 http://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.

<myinla.R>