sum-to-zero constraint

[Reza Hajipour]

Hello

I want to add a sum-to-zero constraint on the random effects of a spatio-temporal model. Is it enough to put  constr=TRUE   as follows:

formula <- O ~ X +
  f(ID.area, model="besag", graph=.... , constr=TRUE, rankdef=1) +
  f(ID.year, model="rw1", constr=TRUE, rankdef=1)

[Esmail Abdul Fattah]

Hi Reza,

Yes.

Best Regards,

Esmail

--
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 view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/37d02205-d9b8-4d4e-8ea0-8fd729e78bcbn%40googlegroups.com.

[Helpdesk (Haavard Rue)]

yes, this gives a sum-to-zero constraint. for the besag and rw1 I would
also add scale.model=TRUE

for 'besag' and 'rw1' which are intrinsice, then constr=TRUE is default.

you add more constraints with the argument 'extraconstr' if needed.

> --
> 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 view this discussion on the web, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/37d02205-d9b8-4d4e-8ea0-8fd729e78bcbn%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org

[Reza Hajipour]

Thank you.

What about the "BYM" model where we have two random effects, one with ICAR prior and one with gaussian prior, does it give sum-to-zero as follows:

f(ID.area, model="BYM", graph=.... , constr=TRUE, rankdef=1) 

[Helpdesk (Haavard Rue)]

in this case the constr=T will only impose a constraint on the 'besag'
part of the model, not the 'iid'

[Reza Hajipour]

Thank you.

The constraint on 'besag' (ICAR) part is enough.

I have another random effect part in the model, interaction between space and time, specifically the interaction between 'besag' (ICAR) and ' rw1' (or rw2). This is specified by the Kronecker product of the besag and  rw. In this case how can we make sum-to-zero constraint?

Based on the  Blangiardo, M. & Cameletti, M. (2015), Spatial and spatio-temporal Bayesian models with R-INLA, page 244, the coding is as follows:

> ID.area.int <- data$ID.area
> ID.year.int <- data$ID.year

> f(ID.area.int, model="besag", graph=... , group=ID.year.int, control.group=list(model="rw2"))

Where should we add the  constr=T?

[Elias T. Krainski]

Please refer to this thread: https://groups.google.com/g/r-inla-discussion-group/c/iEhzLk9rLX8/m/RfS4SgbdAwAJ

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/9b0e3cf7-cb03-4667-bfd6-141f8b496d0an%40googlegroups.com.

[Reza Hajipour]

sum-to-zero constraints ensure the identifiability, if yes, is it necessary to have this constraint on 'iid' parts, as follows?

formula <- O ~ X +

  f(ID.area, model="besag", graph=.... , constr=TRUE) +

  f(ID.area1, model="iid", constr=TRUE) +
  f(ID.year, model="rw2", constr=TRUE) +
  f(ID.year1, model="iid", constr=TRUE)

As you mentioned the constr=TRUE is default, if I don't add it to f(...) parts of the above formula, the conster is True and the sum-to-zero constraint is done for all random effects?

[Elias T. Krainski]

Hi, 

You do not always need constr = TRUE. It is not needed for 'iid' model, or any other model with proper joint distribution. But sometimes you want it, for example for the case of having a long range correlation that interacts with some other model component. You only need constraint in intrinsic models if another model component includes its null space. 

Please take a look at the example below, where I illustrate the fact that you do not always need to add a constraint while working with intrinsic models.

library(INLA)

n <- 30
x <- scale(cumsum(rnorm(n)))

plot(x)

y <- rpois(n, exp(2 + x))

fit1 <- inla(
    y ~ f(i, model = 'rw1'),
    data = list(i = 1:n, y = y),
    family = 'poisson'
)

## The intercept captures the overall mean
fit1$summary.fix

## posterior mean for the constrained random effect is around zero
summary(fit1$summary.random$i$mean)

## Now, take the intercept ou of the linear predictor
## Then, the constraint in the random effect is no longer needed
fit2 <- inla(
    y ~ 0 + f(i, model = 'rw1', constr = FALSE),
    data = list(i = 1:n, y = y),
    family = 'poisson'
)

## unconstrained random effect is around the intercept (not fitted separately)
summary(fit2$summary.random$i$mean)

## However, the fitted values are similar (not exactly the same as the priors are different and there are some small numerical differences while doing the computations)
summary(fit1$summary.fitted.values$mean)
summary(fit2$summary.fitted.values$mean)

regards,

Elias

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/aebac503-c764-430e-8a96-308782ab0b72n%40googlegroups.com.

[Reza Hajipour]

Thank you so much. It was so helpful for me.

As you explained and based on the articles I have studied, the models with improper joint distribution like ICAR, leading to an identifiability problem with the model intercept.

My model also includes the spatial varying coefficient covariates (along with fixed effect intercept) as follow:

formula <- O ~ X +   f(ID.area1, X, model="besag", graph=TR.adj, constr=T) +

  f(ID.area2, model="besag", graph=TR.adj, constr=T) +

  f(ID.area3, model="iid", constr=F) 

In this case, does the ICAR specification of covariate 'X' again lead to an identifiability problem with the model intercept (does it interact with the model intercept) and is the constraint necessary or not?

[Elias T. Krainski]

 For this case, you have to consider either having "+ X" or not. If you have it, the constraint is needed, otherwise not. This will lead to having either (beta + beta_i) X_i (where beta_i should be constrained) or just beta_i X_i.

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/6c55e4cf-be57-4d01-99f4-39c80055f79an%40googlegroups.com.

[Reza Hajipour]

It's great. Thank you.

Yes, I have both fixed slop and random slop covariates. In this case, as you mentioned, 'beta' and 'beta_i' interact, so we need the constraint.

I am fitting a spatio-temporal model with varying coefficients (random slop) covariates. If we have evidence that covariates vary across regions over time, is the model as follows?

formula <- O ~ X +   

f(ID.area1, X, model="besag", graph=TR.adj, constr=T) +

f(ID.year1, X, model="rw2", constr=T) +

f(ID.area2, model="besag", graph=TR.adj, constr=T) +

f(ID.area3, model="iid", constr=F) +

f(ID.year2, model="rw2", constr=T) +

f(ID.year3, model="iid", constr=F)

Is (beta + beta_i) X_it enough or we need (beta + beta_i + beta_t) X_it? I mean is f(ID.year1, X, model="rw2", constr=T) necessary or not?

[Elias T. Krainski]

This fits (beta + beta_i + beta_t) X_it (although the indexing here is not very precise) and having constraints in both 'rw2' and 'besag' is fine. But this is not yet a spatio-temporal model for the regression coefficient, where you would have beta_it.

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/cc5c2d20-8d6b-419e-87a2-a6550ca37b9cn%40googlegroups.com.

[FU peter]

Hello

When using ar1 as a time effect, do the following models need to add additional constraints?

Do ID.year and ID.year1 need to add constr=TRUE?

Elias T. Krainski <eliask...@gmail.com> 于2024年3月7日周四 19:30写道：

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/CAA0dk0EnU67DjBUo567gGDEZMzsr7bF8eF2Ex6bt9xBiQBw6Fg%40mail.gmail.com.

[FU peter]

Hello

When using ar1 as a time effect, do the following models need to add additional constraints?

Do ID.year and ID.year1 need to add constr=TRUE?

st2<- O ~ 1+f(ID.city,model="bym",graph=g,scale.model = T,constr=T,hyper=pc1)+

  f(ID.year,model="ar1",hyper=pc2,constr=T)+

  f(ID.year1, model = "iid",hyper=pc2,constr=T)+

  f(ID.city1,model="iid",hyper=pc2,group=ID.year2,control.group = list(model="ar1"))

Elias T. Krainski <eliask...@gmail.com> 于2024年3月7日周四 19:30写道：

This fits (beta + beta_i + beta_t) X_it (although the indexing here is not very precise) and having constraints in both 'rw2' and 'besag' is fine. But this is not yet a spatio-temporal model for the regression coefficient, where you would have beta_it.

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/CAA0dk0EnU67DjBUo567gGDEZMzsr7bF8eF2Ex6bt9xBiQBw6Fg%40mail.gmail.com.

[Elias T. Krainski]

The 'iid' models do not need constr = TRUE.

[FU peter]

Thank you.

When fitting Knorr-Held models, using "rw1" and "besag" in the group did not take into account the required constraints. Is it the same for using "ar1" in the group?

Elias T. Krainski <eliask...@gmail.com> 于2024年5月12日周日 14:33写道：

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/CAA0dk0HBG8FN4Kkk9zGUSxnDtnG-JXsVaOPrOrfRxX4dRWiX-Q%40mail.gmail.com.

[Elias T. Krainski]

Yes.

[Helpdesk (Haavard Rue)]

the recommentation is to use 'bym2' instead of 'bym', and for this, constr=T is
default.

you do not need, formally, the constr=T, for ID.year and ID.year1, but due to
confounding, sometimes this is done. especially for shorter `in time' data.
constr=T also make less informative priors 'less wrong'.

I guess you see the potential confounding for ar1 and iid for year, if the years
are few...

> > > > > https://groups.google.com/d/msgid/r-inla-discussion-group/6c55e4cf-be57-4d01-99f4-39c80055f79an%40googlegroups.com
> > > > > .

> > > --
> > > 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 view this discussion on the web, visit

> > > https://groups.google.com/d/msgid/r-inla-discussion-group/cc5c2d20-8d6b-419e-87a2-a6550ca37b9cn%40googlegroups.com

[FU peter]

Thank you.

I'm using this code on the official book, and the example data has 18 years of data, is this considered few?

If AR1 doesn't need additional constraints, is it appropriate to use ‘groups’ to model spatio-temporal interactions?

Helpdesk (Haavard Rue) <he...@r-inla.org> 于2024年5月15日周三 15:48写道：

[Helpdesk (Haavard Rue)]

for stability, I would add constr=T as it separate effects and avoid confounding

[FU peter]

Thank you.

I fitted the spatio-temporal interaction model using this code

st2<- O ~ 1+f(ID.city,model=“bym”,graph=g,scale.model = T,constr=T,hyper=pc1)+

  f(ID.year,model=“ar1”,hyper=pc2,constr=T)+

  f(ID.year1, model = “iid”,hyper=pc2,constr=T)+

  f(ID.city1, model=“iid”,hyper=pc2,constr=T,group=ID.year2,control.group = list(model=“ar1”))

 The spatio-temporal interaction effects extracted from this are consistent with “time-varying”, but I can't extract the correct temporal and spatial effects

I use generic0 to impose constraints on the spatio-temporal interaction terms

st3<- O ~ 1+f(ID.city,model=“bym”,graph=g,scale.model = T,constr=T,hyper=pc1)+

  f(ID.year,model=“ar1”,hyper=pc2,constr=T)+

  f(ID.year1, model = “iid”,hyper=pc2,constr=T)+

  f(cityyear,model=“generic0”, Cmatrix=R2, rankdef=rank.def2,

    constr=TRUE, hyper=list(prec=list(prior=pc2)),

    extraconstr=list(A=A.constr2, e=e2) 

  )

I can extract the correct time and space effects from this, but the spatio-temporal interaction effects do not match the “time change”

How can I do this?

Helpdesk (Haavard Rue) <he...@r-inla.org> 于2024年5月17日周五 20:36写道：

[Ismael Ahmad Rab]

How can you also add the constraint on the iid part