Re: [r-inla] adding hierarchy level of Hyper-Hyperparameter

[Helpdesk]

it might be possible using a double-hack to by-pass the system and
creating a 'fake' hyperparmeter.

I need some more time to work out an example

On Thu, 2022-07-28 at 12:07 -0700, Amelie Liese wrote:
> Hi, 
>
> I am fitting an spde model and a generic car model on cosmological
> matter density field data and have two questions:
>
> 1. I would like to add an additional hierarchy level where the
> Hyperparameters of the latent field are conditioned on the
> cosmological parameters (the relationship between hyperparameters and
> cosmological parameters is investigated by a regression beforehand).
> Very simplified:
> Bildschirmfoto 2022-07-28 um 21.02.10.png
> In order to estimate the posterior marginals of the cosmological
> parameters
>
> Bildschirmfoto 2022-07-28 um 21.04.38.png
>
> Is there a way to implement this in the spde approach or for a generic
> car model?
>
>
> 2. How can I integrate different simulated fields in the generic car
> model inference? Is there something similar to the stacking function
> and group indexing from the spde approach?
>
>
>
>
>
>
>
> --
> 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/0da2dc57-c66b-4a7d-a802-c07f6f50ece9n%40googlegroups.com
> .

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

[Amelie Liese]

Thank you for your answer! I am wondering wether I can solve this problem by two steps:

1. To get posteriori distribution of cosmological parameters:

 sampling from marginal posterior of hyperparameters with a transformation  function. Where the transformation function is the regression function 

 cosmic.Para = beta_0 + beta_1*HP1 + … + eps

cosm.Para.samp <- inla.posterior.sample.eval(function(...) {Regressionformula},

                                              m1.samp)

2. Include the cosmological parameters in the prior:

For example, very simplified:

prec.prior <- list(prec = list(prior = "loggamma", param = c(cosm.para1, cosm.para2)),

                                                                                                            initial = 4, fixed = FALSE)

   f <- y ~ 1 + f(idx, model = "iid", hyper = prec.prior)

cosm.para1 and cosm.para2 are sampled from their prior distribution beforehand.

Is the implementation in inlabru the same?

Even if this approach would make sense I suppose that I couldn´t use all of the R-inla features for analysis? 

Regards,

Amelie Liese

[Finn Lindgren]

In inlabru you need to modify the component syntax slightly.
Instead of +f(...) you would use some label, so the component
definition would be +thelabel(idx,model="iid",hyper=prec.prior)

From this I take it that your "\theta" in the earlier model definition
is this precision parameter, so that it's just a single scalar
parameter?
Then I think the joint theta,cp1,cp2 model should be doable as an
generic model, with \theta parameterised as a quantile transformation
of a variable v.
For example, let the prior distributions be
v ~ N(0,1)
log_cp1 ~ LogGamma(...)
log_cp2 ~ LogGamma(...)
Then, when constructing Q, you'd just transform v into it's intended
Ga(cp1,cp2) distribution for the precision (theta):
Q <- Diagonal(n, x = qgamma(pnorm(v), shape = exp(log_cp1), rate =
exp(log_cp2))

For a more numerically stable construction, use the
inlabru::bru_forward_transformation() method:
inlabru::bru_forward_transformation(qgamma, v, shape = exp(log_cp1),
rate = exp(log_cp2))
This uses al the internal log=TRUE and log.p=TRUE options in the pnorm
and qgamma functions, so that the transformation is well-defined in
the distribution tails.

I'm not entirely convinced that you can actually get useful estimates
of cp1 and cp2 out of this model though; not due to inla limitations,
but due to lack of information in the data to actually identify them
in the model!

Finn

> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/8fcd3ca0-9414-4861-89d6-e575209fa555n%40googlegroups.com.



--
Finn Lindgren
email: finn.l...@gmail.com

[Virgilio Gómez-Rubio]

Hi,

For this, you could try to implement your new model using the rgeneric latent effect. Alternatively, you could combine INLA with a number of methods… You can check, for example, these papers (code available):

https://www.tandfonline.com/doi/pdf/10.1080/10618600.2022.2067551

https://arxiv.org/abs/2201.07493

In both cases, Importance sampling is combined with INLA to fit models that cannot otherwise be fit with INLA. These include modeling the hyperparameters using a linear predictor.

Best,

Virgilio

[Amelie Liese]

Hi Finn,

thank you. \theta is not a scalar but contains all the hyperparameters of the latent field (precision of lognormal observations, range of field, stdev of field). 

I guess that in that case I have to specify not only prec*Q but  prec*Q = prec*(I - rho*W)  where rho is the autocorrelation coefficient which I can compute from stdev and range of the latent field?

Yes, I understand that there is a problem with parameter identifiability. That´s why I will probably drop this part in my term paper. Can you give me any hint what could solve the problem of identifiability  in my case?

Kind regards,

Amelie Liese

[Amelie Liese]

Thank you for the link! Is the code for the paper of the second link also available? I could just find code for the first link.

Kind regards,

Amelie Liese

[Amelie Liese]

Hello Finn,

I would like to get back to my question from August. My question wasn’t clear enough.

I am fitting an SPDE model with r-inla (not inlabru) where the range hyperparameter  of the latent field should be conditioned on a cosmological parameter (Omega_M) wich has a flat prior.

range |  Omega_M ~ gamma(a*Omega_M, b*Omega_M)
With a, b being some factor.

I guess the code examples of your answer are not intended for a generic SPDE model ?
I think I have to use  inla.spde2.generic() function. However I could only find some little information about it and don’t know how to use it correctly. 

Thanks and regards,

Amelie Liese

Beyond code chunk may help to understand my problem:

range.prior <- list(range = list(prior = „loggamma“, param = c(a*Omega_M, b*Omega_M)), initial = 4, fixed = FALSE)

# Where I use the ML-Estimate for the range and sigma parameter as priors to generate the spde object

spde <- inla.spde2.pcmatern(mesh, alpha = 2, prior.sigma = c(sigma0_ML,0.5),prior.range = c(range0_ML,0.5))

 s.index <- inla.spde.make.index(name = "spatial.field", n.spde = spde$n.spde)

 s.stack <- inla.stack(data = list(y=obs), A = list(A),

                        effects = list(c(s.index, list(Intercept = 1))), tag = "s.data")

  s.stack.pred <- inla.stack(data = list(y=NA), A = list(1), 

                             effects = list(c(s.index, list(Intercept = 1))), tag = "s.pred")

  join.stack <- inla.stack(s.stack, s.stack.pred)

  formula = y ~ -1 + Intercept + f(spatial.field, model = spde, hyper = range.prior)

  res_inla <- inla(formula, data = inla.stack.data(join.stack, spde = spde), family = "lognormal",

                   control.predictor=list(A=inla.stack.A(join.stack), compute =TRUE),

                   control.compute = list(cpo = TRUE, dic = TRUE, waic = TRUE, config = TRUE),

                   control.inla = list(int.strategy = "grid"))