Add category x space x survey interactions to model comparison

[athowes]

Can you use three options with R-INLA's group option somehow to do this? I guess worst case use the generic to define the space x survey interaction and then put an IID group for category on top of it

[athowes]

Think it's going to have to be like:

inla.rgeneric.besagxar1.model <- function(cmd = c("graph", "Q", "mu", "initial","log.norm.const", "log.prior", "quit"), theta = NULL)

for the Besag spatial and AR1 temporal interaction.

The rest of them have IID parts in so can be made without

[athowes]

Jeff's repo here could be helpful https://github.com/jeffeaton/inla-sandbox/blob/master/inla-internals.md

[athowes]

Rather than try to write a manual AR1 x Besag, I think it'll be easier to write a "manual" Besag x IID then use the group option for AR1. With proper scaling, I think that Besag x IID specified via the group option should be equivalent to a Besag which has copies of the adjacency graph like this:

This could be tested.

[athowes]

e.g. the adjacency matrix looks like

[athowes]

Need to fix this bug:

 Error in inla.interpret.formula(formula, data.same.len = data.same.len,  : 
  
The covariate `area_idx' are used in the following f()-terms: 4, 6
Only one is allowed since the naming convention of the f()-terms depends on unique names of the covariates

You can resolve this as follows:
Not allowed: formula = y ~ f(x, <etc>) + f(x, <etc>) + <etc>
Fix: data$xx = data$x; formula = y ~ f(x, <etc>) + f(xx, <etc>) + <etc> 

Can be accomplished using area_idx_copy and sur_idx_copy

[athowes]

Crashed at MOZ

Begin fitting Model 8x.


	GMRFLib version 3.0-0-snapshot, has recived error no [12]
	Reason    : The Newton-Raphson optimizer did not converge
	Message   : Condition `lambda < 1.0 / lambda_lim' is not TRUE
	Function  : GMRFLib_init_GMRF_approximation_store__intern
	File      : approx-inference.c
	Line      : 2953
	GitID     : file: approx-inference.c  1f6a39183ef43d8ef33f10ff3f04fd13f8432758 - Mon Feb 22 21:27:50 2021 +0300

Segmentation fault (core dumped)
[ failed     ]  2021-10-02 19:54:29
[ elapsed    ]  Failed report in 25.26072 mins
 Error in inla.inlaprogram.has.crashed() : 
  The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
  If this does not help, please contact the developers at <help@r-inla.org>. Begin fitting Model 8x.


	GMRFLib version 3.0-0-snapshot, has recived error no [12]
	Reason    : The Newton-Raphson optimizer did not converge
	Message   : Condition `lambda < 1.0 / lambda_lim' is not TRUE
	Function  : GMRFLib_init_GMRF_approximation_store__intern
	File      : approx-inference.c
	Line      : 2953
	GitID     : file: approx-inference.c  1f6a39183ef43d8ef33f10ff3f04fd13f8432758 - Mon Feb 22 21:27:50 2021 +0300

Segmentation fault (core dumped)
[ failed     ]  2021-10-02 19:54:29
[ elapsed    ]  Failed report in 25.26072 mins
 Error in inla.inlaprogram.has.crashed() : 
  The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
  If this does not help, please contact the developers at <help@r-inla.org>. 

[athowes]

Based on early model comparison plots, looks to me like the interactions are improving on things slightly: probably worth keeping in model comparison

[athowes]

• Note that if these are included there will need to be some work on getting them to have the right constraints (see issue Prevent interactions stealing the category-specific intercept variance #42)
• One option: write a custom R-INLA model; another option: custom constraints
• After: rerun model comparison

[athowes]

See https://github.com/jeffeaton/inla-sandbox/blob/c77d854f43a160841ae7172698e7453e96f08f78/inla-internals.R#L164 for additional information about the interactions between R-INLA's group and extraconstr / constr options.

Background

Currently, space x survey x category random effects are being implemented via (here is Model 6x which has Besag x IID x IID):

f(area_idx_copy, model = "besag", graph = adjM, scale.model = TRUE, group = sur_cat_idx,
  control.group = list(model = "iid"), constr = TRUE, hyper = tau_pc(x = 0.001, u = 2.5, alpha = 0.01))

The challenge here is how to get the right sum-to-zero constraints.

How to do it with simpler random effects

Using the default option constr = TRUE places a sum-to-zero constraint for each value of group. So for simpler age x category random effects

f(area_idx, model = "iid", group = cat_idx, control.group = list(model = "iid"),
            constr = TRUE, hyper = tau_pc(x = 0.001, u = 2.5, alpha = 0.01))

then we get the correct sum-to-zero within each category:

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

What we need for the interactions

What is required for the interaction random effects \delta_{itk} is:

sum_{it} \delta_{itk} = 0 \forall k = 1, ..., K.

What I have currently is:

sum_{i} \delta_{itk} = 0 \forall t, k

If there are K categories and T time points then I am getting K x T sum-to-zero constraints: not only does the sum over it have to be zero for all k = 1, ..., K but so does the sum over each i for all t.

One way to implement this would be to get group = cat_idx rather than group = sur_cat_idx. This would be OK for this particular model but would need a custom model for Besag x AR1.

What's required for each of the interaction models

Model	Space x Survey x Category	Proposed strategy
5x	IID x IID x IID	Done
6x	Besag x IID x IID	Use a repeat_matrix to put Besag x IID in the f and let category be group
8x	IID x AR1 x IID	Not sure.
9x	Besag x AR1 x IID	Not sure.

[athowes]

Got this error for 6x:

Error in inla(formula, data = df_model, family = "xPoisson", control.predictor = list(link = 1),  : 
  In f(area_sur_idx): 'covariate' must match 'values',  and both must either be 'numeric', or 'factor'/'character'. 

EDIT: fixed below.

[athowes]

Jeff suggests to use separate f terms for each category and share precision parameters between terms using copy as a potential workaround.

[athowes]

There is a discussion on the R-INLA Google group about this here: https://groups.google.com/g/r-inla-discussion-group/c/ypFiyve-LLY

[athowes]

I wrote copy in the above commit but I mean replicate.

replicate allows you to define random effects u_1, ..., u_K such that u_k ~ f(\theta), where \theta is shared between the random effects. For more information, see this section of Virgilio Gómez-Rubio's Gitbook.

I believe what I have implemented is close to being correct, but not totally. Going to try running through the countries to look at model fit criteria when using this almost correct interaction term!

• LSO gives convergence issue for R-INLA on Model 8x:

Begin fitting Model 8x.
	GitID: file: error-handler.c  35027aac8a4d2d47029a01aa0d7ced1e20be6ce8 - Tue Nov 16 16:36:14 2021 +0300
	Error:12 Reason: The Newton-Raphson optimizer did not converge
	Message: Condition `lambda < 1.0 / lambda_lim' is not TRUE
	Line:2536 Function: GMRFLib_init_GMRF_approximation_store__intern

Resources for rgeneric

• https://inla.r-inla-download.org/r-inla.org/doc/vignettes/rgeneric.pdf
• https://github.com/jeffeaton/inla-sandbox/blob/master/arima110-inla-rgeneric.R
• https://github.com/athowes/information-smoothing/blob/90967cfce15db3d6b4ba70f84405bc0c853c1bdc/src/example_three-areas/three-areas.Rmd#L149

Downstream issues to fix

• Something with the information criteria? Probably because wrong number of models. (Fixed here.)

id <- orderly::orderly_run("process_information-criteria")

[ failed     ]  2021-11-29 18:02:48
[ elapsed    ]  Failed report in 3.214426 secs
Error in max_idx[[i]] : subscript out of bounds
In addition: There were 22 warnings (use warnings() to see them)

• process_variance-proportions runs but needs to be updated.

[athowes]

#' BWA errored straight away :'(
#' MOZ crashed on 8x
#' MWI crashed on 5x