Detecting unstable precision estimates

[Coralie Williams]

Hi everyone,

 

I'm running simulations using INLA to fit phylogenetic models with one covariate and two random effects (species IID and phylogenetic using a correlation matrix). In approximately 2 per cent of runs, the model returns large random effect variance estimates >10 (e.g., either for species or phylo or for both), but without any convergence warnings. I’ve attached the model output and the code of one model with simulated data that shows a large variance estimate for the phylo random effect even though I have used a penalised complexity (PC) prior on the phylo effect.

 

To improve the estimates, I am thinking of:

·        Changing the prior to param = c(1, 0.1), given I simulated the variances of the random components to be lower than 1.

·        Adding a PC prior to the species random effect.

·        Centering the simulated covariate (mean = 0)

 

My question:
Are there built-in INLA diagnostics or object components (beyond `model_inla$ok` and `model_inla$mode$mode.status`) that I can query to automatically flag runs with unstable or extreme precision estimates? I’d like a single logical flag or numeric indicator in the fitted object for use in my simulation, without plotting or inspecting marginals manually.

 

Thank you in advance for your help and insights.

Coralie

  

Session info:

R 4.4.2

R-INLA 24.06.27

[Helpdesk (Haavard Rue)]

its plausible to think that gaussian data is the most easiest case, but its not.
the unknown scale of the data makes this (and similar) likelihood sensitive to
non-centering issues, and secondly, also initial values. its almost impossible
to write code that cannot be broken by scaling and offset data and covariates.
also the posterior *is* multimodal, as several interpreations of the data
exists.

in your case, then centering of covariate 'x' will help. this just shift the
burden onto the intercept, which is 'centered' using a fixed offset. for the
'random effect' I would suggest adding constr=T to center them, as your priors
is pretty vague, which implicite suggest that you allow them to fly away from
zero.

I would guess this formula will be pretty stable. if you run into issues, then
can fix a higher initial values for 'iid' and 'generic0' adding

hyper = list(prec=list(initial=8)))

so that its higher than the initial value for the likelihood (which 4, which can
be lowered if you know roughtly the scale)

let me know how this goes

best
H


dat$x <- dat$x - mean(dat$x)
dat$off <- mean(dat$y)
model_inla <- inla(yi ~ offset(off) + x + f(species,
#hyper = pcprior,
model = "iid",
constr = TRUE) +
f(phylo, ## this needs to be a numeric to work
model = "generic0",
constr = TRUE,
hyper=pcprior,
Cmatrix = phylo.prec.mat),
family = "gaussian",
data = dat)

> --
> 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, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/5196ecb9-d758-4adc-8417-d1cf3c219e9fn%40googlegroups.com
> .

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

[Coralie Williams]

Thank you very much for your response, it was helpful. I tested these approaches both in combination and on their own using the same set of simulated datasets.

Using the default priors reduced the number of very large random effect variance estimates. Scaling the covariate and setting different starting values also helped, but to a smaller extent. Adding constr=T did not seem to make a difference in these datasets. 

For the purpose of my simulation study, I’ve decided to use default priors and scale the covariate, which reduced extreme cases to fewer than 0.05% across 10,000 runs (and I will exclude these cases from the results).

Thanks again for your insights, they have been helpful beyond this project!
Coralie