Variable evaluates to NAN or INF

[Facundo Muñoz]

Sorry to keep bothering...

I'm getting this error:
file: smtp-taucs.c hgid: 336855f14e99 date: Mon Oct 03 07:47:35 2011 +0200
Function: GMRFLib_build_sparse_matrix_TAUCS(), Line: 654, Thread: 0
Variable evaluates to NAN or INF. idx=(2,2). I will try to fix it...
several times, and finally it crashes saying that the Newton-Raphson
optimizer did not converge.

Now, I have this binomial response (Ntrials=1), and a latent field with
a factor (4 levels), a linear covariate and a spatial (spde) model.

I tried to increase the prec.intercept to 1 (which fixed similar
numerical problems with other models) and even to 10, without any change
with this one.

I also tried to increase the cutoff of the mesh (was 500), and with
cutoff=2000 it stops crashing, but instead I get DIC=INF.

Looking into local.dic and local.p.eff it turns out that one specific
point is taking value INF.
Tried to take it away from the dataset and it crashes again.

�Is there any other parameter that I can tune?

Thanks
�acu.-

[Finn Lindgren]

On Oct 7, 3:43 am, Facundo Muñoz <famu...@gmail.com> wrote:
> Now, I have this binomial response (Ntrials=1), and a latent field with
> a factor (4 levels), a linear covariate and a spatial (spde) model.

The Binomial modems are notoriously difficult, since they easily
result in bimodal distributions, and this is a possible cause of this
problems. But there is one parameter that you need to check; I'm
working on the interface for the SPDE models, and have not written all
the needed documentation yet.

inla.spde.create tries to set a sensible prior for the spatial range,
via the kappa parameter, based on the size of the domain, but it
cannot do the same for the variance rescaling parameter tau. Look in
spde$f$hyper.default
And either change the theta.T specification, or override it with a
hyper=list(theta.T=...)
parameter to your f(model=spde, ...) call.
This may or may not help in your specific situation, but is something
that you likely need to do anyway.

This bit will be easier once I have time to add helper code for it, as
well as write full documentation and tutorials...

Finn

> I tried to increase the prec.intercept to 1 (which fixed similar
> numerical problems with other models) and even to 10, without any change
> with this one.
>
> I also tried to increase the cutoff of the mesh (was 500), and with
> cutoff=2000 it stops crashing, but instead I get DIC=INF.
>
> Looking into local.dic and local.p.eff it turns out that one specific
> point is taking value INF.
> Tried to take it away from the dataset and it crashes again.
>

> Is there any other parameter that I can tune?
>
> Thanks

> acu.-

[Facundo Muñoz]

Thaks Finn, I'll try.

But spde$f$hyper.default is a list comprising theta1, theta2 and theta3,
each of these having an initial value and two parameters.
theta1 have a positive initial value (7.83) while for the other two this
initial value is negative.

I'm looking forward for those tutorials and documentation.
I know it takes a lot of time. I would be happy to help, as far as I can.

Thank you
�acu.-

El 08/10/11 03:31, Finn Lindgren escribi�:

[Finn Lindgren]

Ah, yes, I forgot some of the details of the code...

The model (for R^2 and nu=1, the only case the interface currently covers fully) in the stationary case is

(kappa^2 - Delta ) ( tau x(s) ) = W(s)

where kappa controls the spatial scale (approximate range = sqrt(8)/kappa ), and tau is a rescaling parameter.  The resulting variance of the field is sigma^2 = 1/(4*pi*kappa^2*tau^2)

In the stationary case, only theta1 and theta2 are used, with
theta1 = log(tau)
theta2 = 2*log(kappa)

inla.spde.create attempts to put a sensible prior on theta2 based on the size of the domain.
For the rescaling parameter tau,
theta1 = -0.5*log(4*pi) - 0.5*theta2 - log(sigma)
so a sensible prior mean for theta1 depends on both the prior for theta2 and on the standard deviation of the field.  But inla.spde.create does not know anything about the range of the data values, and instead uses sigma=1.
However, using the above relation, you can modify the prior by subtracting your own log(sigma); something like

hyper=spde$f$hyper.default
hyper$theta1$initial = hyper$theta1$initial - log(sigma)
hyper$theta1$param[1] = hyper$theta1$param[1] - log(sigma)
formula = y ~ f(model=spde, hyper=hyper)

where sigma is you prior guess for the field standard deviation.

(In the future interface, there will be an option to use a different parameterisation such that log(sigma) is a parameter instead of log(tau).)

/Finn

[Finn Lindgren]

...and you can also try changing the prior precisions for theta1 and theta2 (theta3 is only used for non-stationary models), i.e. hyper$theta1$param[2] and hyper$theta2$param[2] . 

/Finn

[Facundo Muñoz]

That's valuable info!!

It works. Although I had to specify both a lower mean and variance for
the prior of theta2, up to a point where the spatial effect is
negligible a priori, haha.

It is worth noting here in the discussion group another approach
suggested by Håvard in a private mail.
As an alternative to tuning the prior of the spatial effect, you can use
the weight parameter in the f() specification. For example,
weight = rep(0.001, n)
f(idx, weight, model="spde"...)

Finally, I also could fix the numerical problems by specifying a more
informative prior for the fixed effects (leaving the spatial effect as
by default).
Namely, using
control.fixed = list(prec.intercept = .1, prec=.1)

So, thanks a lot to Håvard and Finn for their invaluable work and help,
and hope this will help someone else running through the same issues.
ƒacu.-

El 09/10/11 12:47, Finn Lindgren escribió:

[Facundo Muñoz]

... as a curiosity: the priors on theta_i are Normal distributions?

Thanks
ƒacu.-

El 09/10/11 12:35, Finn Lindgren escribió:

[Finn Lindgren]

 Yes, in the current version they are independent Normal distributions.
In the version under development they will be jointly multivariate
Normal (allowing a bit more flexibility in the parameterisation).
/Finn

[Victor]

I have been using the new build spde2 and despite trying the above approach, initializing theta1, the problem failed to resolve. here is extract (long format),

name=[INLA.Model]
                strategy=[default]
        store results in directory=[C:Temp/Rtmpys8p9F/file2d0732a0/results.files]
                output:
                        cpo=[1]
                        dic=[1]
                        kld=[1]
                        mlik=[1]
                        q=[0]
                        graph=[0]
                        hyperparameters=[1]
                        summary=[1]
                        return.marginals=[1]
                        nquantiles=[4]  [ 0.025 0.5 0.95 0.975 ]
                        ncdf=[0]  [ ]
        parse section=[2] name=[Predictor] type=[PREDICTOR]
        inla_parse_predictor ...
                section=[Predictor]
                dir=[predictor]
                PRIOR->name=[loggamma]
                PRIOR->from_theta=[function (x) <<NEWLINE>>exp(x)]
                PRIOR->to_theta = [function (x) <<NEWLINE>>log(x)]
                PRIOR->PARAMETERS=[1, 1e-005]
                initialise log_precision[11]
                fixed=[1]
                user.scale=[1]
                n=[3621]
                m=[21267]
                ndata=[21267]
                compute=[1]
                Aext=[C:/Temp/Rtmpys8p9F/file2d0732a0/data.files/file481d7b54]
                AextPrecision=[3.269e+006]
                output:
                        summary=[1]
                        return.marginals=[1]
                        nquantiles=[4]  [ 0.025 0.5 0.95 0.975 ]
                        ncdf=[0]  [ ]
        parse section=[1] name=[INLA.Data1] type=[DATA]
        inla_parse_data [section 1]...
                tag=[INLA.Data1]
                family=[GAUSSIAN]
                link=[IDENTITY]
                likelihood=[GAUSSIAN]
                file->name=[C:Temp/Rtmpys8p9F/file2d0732a0/data.files/file5dc458b]
                file->name=[C:Temp/Rtmpys8p9F/file2d0732a0/data.files/file400d44ff]
                read n=[10284] entries from file=[C:d0732a0/data.files/file5dc458b]
                        0/3428  (idx,a,y,d) = (1, 1, 2.20727, 1)
                        1/3428  (idx,a,y,d) = (20, 1, -1.#INF, 1)
                        2/3428  (idx,a,y,d) = (34, 1, -1.#INF, 1)
                        3/3428  (idx,a,y,d) = (39, 1, 3.44202, 1)
                        4/3428  (idx,a,y,d) = (80, 1, -1.#INF, 1)
                use variant [0]
                        bit 0 is off
                        bit 1 is off
                        bit 2 is off
                        bit 3 is off
                initialise log_precision[4]
                fixed=[0]
                PRIOR->name=[loggamma]
                PRIOR->from_theta=[function (x) <<NEWLINE>>exp(x)]
                PRIOR->to_theta = [function (x) <<NEWLINE>>log(x)]
                PRIOR->PARAMETERS=[1, 5e-005]
        parse section=[3] name=[field] type=[FFIELD]
        inla_parse_ffield...
                section=[field]
                dir=[random.effect00000001]
                model=[spde2]
                PRIOR0->name=[mvnorm]
                PRIOR0->from_theta=[function (x) <<NEWLINE>>x]
                PRIOR0->to_theta = [function (x) <<NEWLINE>>x]
                PRIOR0->PARAMETERS0[0]=[-3.56078]
                PRIOR0->PARAMETERS0[1]=[-0.189643]
                PRIOR0->PARAMETERS0[2]=[0.1]
                PRIOR0->PARAMETERS0[3]=[0]
                PRIOR0->PARAMETERS0[4]=[0]
                PRIOR0->PARAMETERS0[5]=[0.1]
                constr=[0]
                diagonal=[0]
                si=[0] (if possible)
                id.names=<not present>
                compute=[1]
                nrep=[1]
                ngroup=[7]
                read covariates from file = tmpys8p9F/file2d0732a0/data.files/file3d4e7858] 4/3621  (idx,y) = (4, 0.001)
                spde2.prefix = [C:/Users/Temp/Rtmpys8p9F/file2d0732a0/data.files/fileb26704b/file373d6050.]
                spde2.transform = [identity]
                ntheta = [2]
                initialise theta[0]=[-3.56078]
                fixed[0]=[0]
                initialise theta[1]=[-0.189643]
                fixed[1]=[0]
                computed/guessed rank-deficiency = [0]
                output:
                        summary=[1]
                        return.marginals=[1]
                        nquantiles=[4]  [ 0.025 0.5 0.95 0.975 ]
                        ncdf=[0]  [ ]
                group.type = ar1
                initialise group_rho_intern[2]
                group.fixed=[0]
                GROUP.PRIOR->name=[normal]
                GROUP.PRIOR->from_theta=[function (x) <<NEWLINE>>log((1 + x)/(1 - x))]
                GROUP.PRIOR->to_theta = [function (x) <<NEWLINE>>2 * exp(x)/(1 + exp(x)) - 1]
                GROUP.PRIOR->GROUP.PARAMETERS=[0, 0.15]
        Index table: number of entries[2], total length[40785]
                tag                            start-index     length
                Predictor                               0      24888
                field                               24888      15897
        parse section=[4] name=[INLA.Parameters] type=[INLA]
        inla_parse_INLA...
                section[INLA.Parameters]
                        lincomb.derived.only = [Yes]
                        lincomb.derived.correlation.matrix = [No]
                global_node.factor = 2.000
                global_node.degree = 2147483647
                reordering = -1
        diagonal (expert emergency) = 0.001
Contents of ai_param 00000000057BC910
        Optimiser: DEFAULT METHOD
                Option for domin-BFGS: epsx = 0.01
                Option for domin-BFGS: epsf = 1e-005 (rounding error)
                Option for domin-BFGS: epsg = 0.01
                Option for GSL-BFGS2: tol  = 0.1
                Option for GSL-BFGS2: step_size = 1
                Option for GSL-BFGS2: epsx = 0.01
                Option for GSL-BFGS2: epsf = 0.001
                Option for GSL-BFGS2: epsg = 0.01
                Restart: 0
                Mode known: No
        Gaussian approximation:
                abserr_func = 0.0005
                abserr_step = 0.0005
                optpar_fp = 0
                optpar_nr_step_factor = -0.1
        Gaussian data: Yes
        Strategy:       Use the Gaussian approximation
        Fast mode:      On
        Use linear approximation to log(|Q +c|)? Yes
                Method:  Compute the derivative exact
        SI directory:   <NONE>
        Parameters for improved approximations
                Number of points evaluate:       9
                Step length to compute derivatives numerically:  0.000122
                Cutoff value to construct local neigborhood:     0.000100
        Log calculations:        On
        Log calculated marginal for the hyperparameters:         On
        Integration strategy:    Use points from Central Composite Design (CCD)
                f0 (CCD only):   1.100000
                dz (GRID only):  1.000000
                Adjust weights (GRID only):      On
                Difference in log-density limit (GRID only):     2.500000
                Skip configurations with (presumed) small density (GRID only):   On
        Gradient is computed using Forward difference with step-length 0.010000
        Hessian is computed using Central difference with step-length 0.100000
        Hessian matrix is forced to be a diagonal matrix? [No]
        Compute effective number of parameters? [Yes]
        Perform a Monte Carlo error-test? [No]
        Interpolator [Auto]
        CPO required diff in log-density [3]
        Stupid search mode:
                Status     [On]
                Max iter   [1000]
                Factor     [1.01]
        Numerical integration of hyperparameters:
                Maximum number of function evaluations [100000]
                Relative error ....................... [1e-005]
                Absolute error ....................... [1e-006]
        To stabalise the numerical optimisation:
                Minimum value of the -Hesssian [0]
        CPO manual calculation[No]

inla_build: check for unused entries in[C:/Users/Temp/Rtmpys8p9F/file2d0732a0/Model.ini]
inla_INLA...
        Strategy = [DEFAULT]
        Size is [40785] and strategy [LARGE] is chosen
        Size of full graph=[40785]
        Found optimal reordering=[amdc] nnz(L)=[5047128] and use_global_nodes(user)=[no]
        List of hyperparameters:
                theta[0] = [Log precision for the Gaussian observations]
                theta[1] = [Theta1 for field]
                theta[2] = [Theta2 for field]

        file: smtp-taucs.c  hgid: 6bab309d6ee2  date: Fri Jun 08 15:13:40 2012 +0200
        Function: GMRFLib_build_sparse_matrix_TAUCS(), Line: 669, Thread: 0

[Finn Lindgren]

Hi,
I suspect that your mesh is ill-conditioned. Please make sure that
you have a "good" triangulation; this often requires the use of the
"cutoff" parameter to inla.mesh.create or inla.mesh.create.helper,
that determines how close measurements are allowed to be and still
have their own triangulation node.
/Finn

--
Finn Lindgren
email: finn.l...@gmail.com
skype: finn.lindgren / +46-46-2880644

[Victor]

Many thanks,

i will try changing the parameters for mesh and specifying the cutoff - will report back if i have any success

"good" conditioned here i presume notably the interior angle with respect to distance? and in other words i have measurements that are very close in space which should not from a triangulation and these rather makes my matrix ill conditioned

best
Victor

[Finn Lindgren]

On 27 June 2012 09:41, Victor <vale...@gmail.com> wrote:
> "good" conditioned here i presume notably the interior angle with respect to
> distance?

That is one part of it, yes; setting min.angle=26 is usually ok (the
triangulation algorithm is theoretically guaranteed to terminate for
min.angle<=21, but values as high as 30 often work in practice;
however, since there is a tradeoff between "nice" triangles and "few"
triangles, it's likely not worthwhile to use very large
min.angle-values).

> and in other words i have measurements that are very close in
> space which should not from a triangulation and these rather makes my matrix
> ill conditioned

Precisely; the size ratio between large and small triangles shouldn't
be too large. Closely clustered sites with empty space in between are
among the worst offenders.

/Finn