Prevent interactions stealing the category-specific intercept variance

[athowes]

I used to not include category random effects in the model, but then I thought that it would be strange not to have first order effects if I'm including interactions. Anyway, it turns out that the category effects are basically not doing anything. This is a bit of an issue because I think there is sufficient flexibility in any of the interaction random effects to subsume the category random effects. And since we care about the proportion variance this is bad, since it seems somewhat arbitrary which of the interaction effects actually takes the category part. I think the solution is some kind of constraint on all of the interaction random effects to stop them absorbing the category random effects. I'm already constraining them to sum-to-zero as in here so I think it'll have to be something more complicated. Possibly constrain to sum to zero within each i, t or a, not only over all ik, tk or ak.

In other words, say for the age-category random effects, what is required is:

sum_a \alpha_{ak} = 0 \forall k = 1, ..., K

whereas what is there at the moment is

sum_{k = 1}^K \sum_a \alpha_{ak} = 0

[athowes]

From here, impose additional linear constraints on the latent field u via Au = e using extraconstr = list(A = A, e = e).

A should be a matrix which has ncol = length(u) and nrow equal to the number of constraints desired. e should length equal to the number of constraints required.

So in this case, say for the age groups 15-19, 20-24 and 25-29, we would have something like:

A = (1, 1, 1, 0, 0, 0, 0, 0, 0) (0, 0, 0, 1, 1, 1, 0, 0, 0) (0, 0, 0, 0, 0, 0, 1, 1, 1)

e = (0) (0) (0)

[athowes]

Looks good for this being done / doable:

> tmp1$fit$summary.random$age_cat_idx$mean[1:3] %>% sum()
[1] -0.04291764
> tmp2$fit$summary.random$age_cat_idx$mean[1:3] %>% sum()
[1] 3.814697e-05

where tmp1 has the usual sum-to-zero constraint and tmp2 has age-category specific sum-to-zero constraints.

[athowes]

Looks plausible to me that with the group option being used that the constraint is applied to each group already. This is with the extraconstr on age x category and regular on the area x category (using group option):

> tmp$fit$summary.random$area_idx$mean[1:33] %>% sum()
[1] -0.0006370672

[athowes]

The only place that the extraconstr are required (or the least bad option) are for the interaction random effects. Here is an example of the Besag x AR1 interactions:

#' ALERT! Constraints not properly being applied here!
#' Going to need a custom extraconstr to make this one correct
#' Or the other option is to put area_sur_idx, but then a different model is needed

#' Model 9x: category random effects (IID), age x category random effects (IID),
#' space x category random effects (Besag), survey x category random effects (AR1)
#' space x survey x category random effects (Besag x AR1)
formula9x <- update(
formula9,
. ~ . + f(area_cat_idx, model = "besag", graph = interaction_adjM, scale.model = TRUE, group = sur_idx,
control.group = list(model = "ar1", hyper = list(rho = list(prior = "pc.cor1", param = c(0, 0.75)))),
constr = TRUE, hyper = tau_pc(x = 0.001, u = 2.5, alpha = 0.01))
)

The issue is that here group = sur_idx rather than cat_idx. Putting group = cat_idx would require to have f(area_sur_idx, ...) but then there would have to be a model = "besag-ar1" which there isn't (unless you defined it yourself using generic0)

It might be that for grouped models the only constraints you can specify are within-group, in which case this workaround for "besag-ar1" not existing isn't going to work and in order to get the interactions with proper constraints it would require a custom implementation.

[athowes]

This is working fine looking at what comes out of the variance comparison script. Haven't done this for the interaction random effects but I'll close this and leave a note on that issue to fix if I end up including them.