constr=TRUE and diagonal=1e-06 features in control.group option for rw1 and rw2 models

[Yimer Wasihun]

Dear All,

I am always thankful for the great help you have.

My work is on spatiotemporal modeling of geostatistical spatiotemporal data using SPDE approach. Let me show part of my code and I will proceed with the questions,

Formula <- y ~ -1 + Intercept + f(field, model=spde, group=field.group, diagonal=1e-06,

          control.group=list(model="rw1", scale.model=TRUE, adjust.for.con.comp=TRUE))   

In previous discussions, it has been mentioned that for intrinsic models (which are almost singular) we need to specify constr=TRUE and add diagonal=1e-06 by default (and very small values in general) to make them non-singular. 

Q1) Is the above formula correct for specifying the diagonal value to be added on the precision matrix of rw1? 

Q2) When we use rw1 or rw2 model under control.group,  how is it possible to include the constr=TRUE option in r-inla?

I was thinking diagonal=1e-06 and constr=TRUE to be included inside the control.group option, but there are no such features in r-inla. 

Using the above formula, I tried to use diagonal=1e-06, diagonal=1e-05, diagonal=1e-04 and diagonal=1e-03 , and the spatial predictions are not the same for all these diagonal values.  So, I am afraid if this is because the diagonal values have been misplaced, and constr=TRUE option wasn't included. 

I really appreciate if you have some inputs on this thing.

Kind regards,

Yimer

[Håvard Rue]

On Wed, 2016-03-23 at 19:40 -0700, Yimer Wasihun wrote:

My work is on spatiotemporal modeling of geostatistical spatiotemporal data using SPDE approach. Let me show part of my code and I will proceed with the questions,

Formula <- y ~ -1 + Intercept + f(field, model=spde, group=field.group, diagonal=1e-06,

          control.group=list(model="rw1", scale.model=TRUE, adjust.for.con.comp=TRUE))   

In previous discussions, it has been mentioned that for intrinsic models (which are almost singular) we need to specify constr=TRUE and add diagonal=1e-06 by default (and very small values in general) to make them non-singular. 

Q1) Is the above formula correct for specifying the diagonal value to be added on the precision matrix of rw1? 

Hi,

this is correct, although do you not need 'adjust.for.con.comp' as it is not used (its only for 

model = besag).

Formula <- y ~ -1 + Intercept + f(field, model=spde, group=field.group, diagonal=1e-06,

          control.group=list(model="rw1", scale.model=TRUE, adjust.for.con.comp=TRUE))   

and it is not clear if you really need 'diagonal=', or if you want the constr=F|T.

The 'diagonal' term is automatically added when you specify an intrinsic model, like 'rw2', for which also constr=T is default. 

This to avoid confounding with the intercept. Also, due to how the 'constr=T' is handled internally, ie the algorithm used to compute it,

then the model has to be proper BEFORE the constr=T is accounted for. So we deal with this by adding a small term on the diagonal. 

In your case, then the 'spde' model is likely already proper, hence you do not need diagonal=, unless you experience

numerical instabilities. 

be aware, that the group model rw1, is intrinsic, which normally, in this setting does not create trouble, as the level for each

location is determined by data, hence the 'time' defined by 'group' is proper. if you experience trouble here (read HUGE stderrors), 

then you may revisit this issue. 

I would guess that you can remove 'diagonal' and 'adjust.for.con.comp'.

Q2) When we use rw1 or rw2 model under control.group,  how is it possible to include the constr=TRUE option in r-inla?

yes, if you have 

f(space, group = time, constr=T)

this means that the (weighted) sum over space is zero, for every timepoint.  so this is replicated to |time| constraints. 

Using the above formula, I tried to use diagonal=1e-06, diagonal=1e-05, diagonal=1e-04 and diagonal=1e-03 , and the spatial predictions are not the same for all these diagonal values.  So, I am afraid if this is because the diagonal values have been misplaced, and constr=TRUE option wasn't included. 

try to set it to zero, which should be the correct option.   if you run into trouble, let me know

Best,

H

-- 
Håvard Rue
Department of Mathematical Sciences
Norwegian University of Science and Technology
N-7491 Trondheim, Norway
Voice:  +47-7359-3533    URL  : http://www.math.ntnu.no/~hrue  
Mobile: +47-9260-0021    Email: havar...@math.ntnu.no

R-INLA: www.r-inla.org

[Yimer Wasihun]

Dear Håvard,

Thank you for the detailed explanation and your quick reply.  I needed to know on both ‘constr=TRUE’ and ‘diagonal=‘  options, and sorry for the unclarity on the previous post. And your answers are helpful to understand how to use them.

 

• Based on your suggestion, I go through my code and re-write the formula by including diagonal=0 option.
  

Formula <- y ~ -1 + Intercept + f(field, model=spde, group=field.group, diagonal=0, control.group=list(model="rw1", scale.model=TRUE))

So, in this case, I found the following trouble. 

Function: GMRFLib_factorise_sparse_matrix_TAUCS(), Line: 829, Thread: 0

Fail to factorize Q. I will try to fix it

This is actually the main reason why I previously added diagonal=1e-06. This problem will go away when I add  diagonal=1e-06. But I guess, in previous formula, the diagonal=1e-06 is not added on the precision matrix of either "rw1" or "spde" model separately, rather on the joint precision matrix. Keeping diagonal=0 and replacing rw1 by ar1 model also get ride of the problem. I therefore suspect this can be related to the use of rw1, but not to the spde model like as you explained. 

 And then I tried the following 

 Formula <- y ~ -1 + Intercept + f(field, model=spde, group=field.group, diagonal=0, control.group=list(model="rw1", scale.model=TRUE, diagonal=1e-06))

but I found the following problem 

 Error in inla.check.control(control.group) : 

  Name `diagonal' in control-argument `control.group', is void.

Q1) So would you mind please to help me how to get rid of this? 

• About your suggestion on the constr=TRUE option, I tried it like as follows. 
  

Formula <- y ~ -1 + Intercept + f(field, model=spde, group=field.group, diagonal=0, constr=TRUE, control.group=list(model="rw1", scale.model=TRUE))

Here also, I found the following problem:

Error in INLA::f(field, model = spde, group = field.group, constr = TRUE,  : 

  Option 'constr=TRUE' is disabled for model='spde2' and 'spde' and 'spde3'

So, I removed constr=TRUE option here, and this okay because I don’t need to put this constraint for the spde model. What I previously needed was to constrain the rw1 and rw2 models, but now you already told me that constr=TRUE is the default for these models under the ’control.group’ feature. This clarify my confusion. 

Finally, the only worry is about the “diagonal=“ option.  So, I will really appreciate your inputs on this. Thank you so much for everything. 

Kind regards,

Yimer

[Håvard Rue]

On Thu, 2016-03-24 at 09:27 -0700, Yimer Wasihun wrote:
> 
> Function: GMRFLib_factorise_sparse_matrix_TAUCS(), Line: 829, Thread:
> 0
> Fail to factorize Q. I will try to fix it
> This is actually the main reason why I previously added diagonal=1e-
> 06. This problem will go away when I add  diagonal=1e-06. But I
> guess, in previous formula, the diagonal=1e-06 is not added on the
> precision matrix of either "rw1" or "spde" model separately, rather
> on the joint precision matrix. Keeping diagonal=0 and replacing rw1
> by ar1 model also get ride of the problem. I therefore suspect this
> can be related to the use of rw1, but not to the spde model like as
> you explained. 

hmm, themn you should add it. 

>  And then I tried the following 
>  Formula <- y ~ -1 + Intercept + f(field, model=spde,
> group=field.group, diagonal=0, control.group=list(model="rw1",
> scale.model=TRUE, diagonal=1e-06))
> but I found the following problem 
>  Error in inla.check.control(control.group) : 
>   Name `diagonal' in control-argument `control.group', is void.

there is no option to add it there. its eqv, since the diagonal term is
 SPDE_ii * RW1_ii. 

> About your suggestion on the constr=TRUE option, I tried

>   Option 'constr=TRUE' is disabled for model='spde2' and 'spde' and
> 'spde3'

you have to add it when you create the spde-model, this because the
correct way is to say that constr=T is 

\sum_i C_i x_i = 0

which is what is done when you create the spde, while constr=T in the
f() does \sum_i x_i = 0.

> Finally, the only worry is about the “diagonal=“ option.  So, I will
> really appreciate your inputs on this. Thank you so much for
> everything. 

you have to think if what you estimate is 'identifiable', in any case
you'll see that if you have large marginal stdev's. 

[Håvard Rue]

On Thu, 2016-03-24 at 09:27 -0700, Yimer Wasihun wrote:

> Dear Håvard,
>


the JSS article discuss in fact your questions, 

http://www.r-inla.org/examples/tutorials/spde-tutorial-from-jss

[Elias T Krainski]

On 24/03/16 17:27, Yimer Wasihun wrote:

Formula <- y ~ -1 + Intercept + f(field, model=spde, group=field.group, diagonal=0, control.group=list(model="rw1", scale.model=TRUE))

The precision for this space-time model is the Kronecker product between the precision matrix from the SPDE model and the structure matrix defined by the temporal neighbourhood, Q (x) R, and has rank deficiency equals the number of nodes in the mesh.

The diagonal added (in the way Håvard Rue described) is enough to make the model proper, however, as he warned, you may see HUGE standard errors.

The sum-to-zero constraints, see [1] and [2], is the usual way to solve this problem. In this case the number of such constraints is equal to the number of mesh nodes, as it needs sum-to-zero constraint for each time series. Since the computational time grows squared on the number of such constraints, it is not a cheap solution. The use of 'ar1' model is computationally convenient.

Elias

[Elias T Krainski]

On 24/03/16 20:35, Elias T Krainski wrote:
> The sum-to-zero constraints, see [1] and [2], is the usual way to
> solve this problem.

the references
[1] http://www.ncbi.nlm.nih.gov/pubmed/10960871
[2] http://onlinelibrary.wiley.com/doi/10.1002/env.1065/abstract

Elias

[Yimer Wasihun]

Dear Havard and Elias,

Thank you so much for your help both of you. I went through the papers you suggested me to read, and they were very helpful to understand the importance of linear constraints when we use rank deficiency latent models. 

Kind regards,

Yimer