Re: [r-inla] INLA Linear Constraints

[INLA help]

On Fri, 2012-03-16 at 12:24 -0700, Dong Liang wrote:
> Hi Group,
>
> I was running a disease mapping model on simulated case-control data.
> The model includes covariates and a latent spatial field with the SPDE
> prior. This model was fitted using INLA under a minute.
>
> Then I put some linear constraints on the latent variables. But the
> program crashed. The marginal likelihood optimization routine diverged
> to infinity. I tried to change diagonal, strategy, int.strategy but
> have not identified the problem.
>
> Please advise the R code and data that generated the scenario. Any
> help will be appreciated.

the problem is a bit 'badly scaled'; using more O(1) coordinates make it
run more nicely. this is due to internal stuff in the code.

I attach a revised code, where the rescaling of the coordinates are just

borders = borders / 10^4
work[, "X"] = work[, "X"] /10^4
work[, "Y"] = work[, "Y"] /10^4

the constrained case now runs, although I strongly suspect there is very
little information in the data about the spatial component, and that is
why these data are 'difficult'.

--
INLA help <he...@r-inla.org>
R-INLA

[Finn Lindgren]

Somehow, the original message seems to have disappeared, but here are some further comments:

The mesh you are using is very ill-defined, as you have simply constructed a pure Delaunay triangulation of your input points, that are very clustered, making for poor FEM numerics.  As a minimum, you should use "refine=TRUE" (to get minimum triangle angle >=21 degrees) and "cutoff=1000" ( distances closer than this are likely too short for the large-scale spatial model to handle, so don't even try; data closer together than this will be treated as having the same true spatial field value, unless you introduce an "A matrix" in the observation equation).  Also with your data scale (metres? maybe not.), to reduce the edge effects where you actually have data, I suggest setting at least "offset=20000", as "100" is _really_ small...

I would also suggest trying the newer spde2 model, as it will set a prior for log(kappa) that is appropriate for the size of your domain; this may eliminate the need for rescaling the domain.  The default parametrisation of the spde2 model is theta1=log(tau) and theta2=log(kappa).

I couldn't test the resulting code, since you didn't include your linear constraint generator, which depends on the size and location data of the mesh, which will now be slightly different (3515 instead of 3757 vertices).  The changed bits:

mesh.dummy = inla.mesh.create(loc=borders,extend=list(n=16,offset=10000),refine=FALSE)

boundary = inla.mesh.boundary(mesh.dummy)

mesh = inla.mesh.create(loc=cbind(work[,'X'],work[,'Y'], 0),boundary=boundary,

                        refine=TRUE, cutoff=1000)

nmesh = mesh$n

mesh.idx = mesh$idx$loc

## Definition of the SPDE model

spde = inla.spde2.matern(mesh=mesh, alpha=2)

/Finn L

[Dong Liang]

Thanks for the comments. After applying the changes, the constraint case now runs, starting from the un-constraint fit. But the posterior mean of the spatial component was biased toward the null. I tried to start the constraint model from scratch, but the optimization still diverged to infinity. I think this data set is 'hard' for spatial modeling.

This data were simulated based on residential locations. Thus, clustering is expected since people tends to live close to each other. SPDE is ideal for health application since it can handle clustered point patterns on large spatial domain effectively.

[INLA help]

On Sun, 2012-03-18 at 17:48 -0700, Dong Liang wrote:
>
> Thanks for the comments. After applying the changes, the constraint
> case now runs, starting from the un-constraint fit. But the posterior
> mean of the spatial component was biased toward the null. I tried to
> start the constraint model from scratch, but the optimization still
> diverged to infinity. I think this data set is 'hard' for spatial
> modeling.

yes, this is a recurrent theme.

if you have

f(i, ..., constr=TRUE)

is means that

sum(x) = 0,

say, in the model, hence we would like that

sum( E(x | y) ) = 0

so that

sum(result$summary.random$i$mean)

should be (close to zero). This is not necessarily so in INLA. (Well, it
will be not that far off but not zero.) This beacuse, the marginal
expectation E( x_i | y) is computed from the marginal of x_i|y which
again is computed under the model where sum(x) = 0. However, the
estimate of x_i|y is done 'independently' for each i, so we do not 'see
the sum(x)=0 constraint' a posteriori, so even though sum(x)=0 apriori,
there is not 'joint distribution' a posteriori that will ensure that
sum(x)=0 also aposteriori. We only do the best guess we can for each i.

Only if strategy="gaussian", you will get that sum(E(x|y)) = 0 (or close
to zero), as the sum(x)=0 constraint survive after conditioning on
data.

Similar for the more general constraints, Ax=e.

> This data were simulated based on residential locations. Thus,
> clustering is expected since people tends to live close to each other.
> SPDE is ideal for health application since it can handle clustered
> point patterns on large spatial domain effectively.
>

binary data are known to be difficult as there is little information in
it, more informative priors on the SPDE model would help from INLA's
point of view.

[Dong Liang]

Thanks for your advices. I have changed my code accordingly. In addtion, I fixed an error in my  code so that linear constraints match the non-normal observations. However,

the constraint SPDE model still diverged: optimizer jumped back and forth on a ridge toward

infinity. I have tried setting diagonal option in control.inla , prec and prec.intercept in control.fixed. Neither fixed the problem. This is puzzling to me since these linear constraints were supposed to improve the parameterization of the model, and the computing. Please advise the attached code puzzle.R

 

Regards,

Dong

[Finn Lindgren]

Den 20 mar 2012 kl. 20:38 skrev Dong Liang <dong.m...@gmail.com>:

Thanks for your advices. I have changed my code accordingly. In addtion, I fixed an error in my  code so that linear constraints match the non-normal observations. However,

the constraint SPDE model still diverged: optimizer jumped back and forth on a ridge toward

infinity. I have tried setting diagonal option in control.inla , prec and prec.intercept in control.fixed. Neither fixed the problem. This is puzzling to me since these linear constraints were supposed to improve the parameterization of the model, and the computing. Please advise the attached code puzzle.R

Constraints in the SPDE models on meshes do not necessarily behave as intended.  I'll take a look at your code.

Finn

<puzzle.R>