multiple likelihoods including bgev

[TheCorinna1994]

Dear all,

I try to run a 2 year coregionalization model with the following likelihoods bgev, gaussian and bgev, which looks like this

inla(formula4_max, family=c("bgev","gaussian","bgev"),
                              data=inla.stack.data(stack_all_early,spde=spde.nonstat),
                              control.family = list(hyper=hyper.bgev, control.bgev=control.bgev),
                              control.predictor = list(A=inla.stack.A(stack_all_early),compute=TRUE,link=1),
                             control.compute = list(dic=TRUE,cpo=TRUE, waic=TRUE, config=TRUE,openmp.strategy="pardiso.parallel"),
                             control.inla=list(int.strategy="eb"),verbose=TRUE)

Unfortunately I get the following error: 

Error in inla.core(formula = formula, family = family, contrasts = contrasts,  :

  Argument 'control.family' does not match length(family)= 3

If I choose a PC prior for precision of errors for the gaussian distribution (eprec<-list(hyper=list(prec=list(prior="pc.prec", param=c(1,0.01)))) ) it still yields the same error statement.

Do you have any idea how to circumvent this problem?

All the best, Corinna

[TheCorinna1994]

Maybe it is due to the inla.mdata in the formula since this is not applicable to the gaussian distribution.

 formula4_max<-inla.mdata(predictor, spread.x_early, tail.x_early) ~-1+  Intercept1+Intercept2+Intercept3+
  f(field_rain, model=spde.nonstat, group = field_rain.group, control.group = list(model="ar1", hyper=rho))+
  f(field_temp, model=spde.nonstat, group = field_temp.group, control.group = list(model="ar1", hyper=rho))+
  f(field_max, model=spde.nonstat, group = field_max.group, control.group = list(model="ar1", hyper=rho))+
  f(field_rain_copy,copy="field_rain", group = field_rain_copy.group, fixed=FALSE, hyper=hyper)+
  f(field_rain_copy2, copy="field_rain", group = field_rain_copy2.group, fixed=FALSE, hyper=hyper)+
  f(field_temp_copy, copy="field_temp", group = field_temp_copy.group, fixed=FALSE, hyper=hyper)

Still I dont know how to deal with this problem.

Thank you!

[Finn Lindgren]

Since you have three "family" values, you should make control.family
be a list of three elements, one for each of the three;
I think
control.family = list(list(hyper=hyper.bgev, control.bgev =
control.bgev), list(), list(hyper=hyper.bgev, control.bgev =
control.bgev))
is what you need (if you want the same settings for both bgev observations).

> --
> 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/803f3a6d-3417-41b0-ae5c-9b257a1850b2n%40googlegroups.com.



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

[TheCorinna1994]

Thank. But now I get a new error statement, saying: Error in inla.create.data.frame(y.orig=yy, mf=mf, E=E, scale=scale, : ncol(Y)==1 is not TRUE

Do you maybe have another advice for me please?

Thank you so much!

[Finn Lindgren]

That does look more related to mdata… I have no advice at this time… I’ll soon work on mdata for inlabru and then I’ll have more insights, but I’m busy with exam marking first.

Finn

On 3 Jan 2024, at 15:08, TheCorinna1994 <corinn...@hotmail.com> wrote:



To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/a3406d55-b7aa-4cb5-9e9f-92f61a820c1cn%40googlegroups.com.

[Håvard Rue]

On Wed, 2024-01-03 at 06:08 -0800, TheCorinna1994 wrote:
> inla.mdata(predictor, spread.x_early, tail.x_early) ~

this says one family only, so it crash with

family=c("bgev","gaussian","bgev"),

with that family, then one could use


Y <- list(inla.mdata(predictor, spread.x_early, tail.x_early),
y.gaussian,
inla.mdata(predictor, spread.x_early, tail.x_early))

with obvious changes in the variable names and you also need to pad with
NA's, similar to

> Y
[,1] [,2] [,3]
[1,] 1 NA NA
[2,] NA 2 NA
[3,] NA NA 3

using family=rep("gaussian", 3)

however, the use of 'inla.mdata' does not play with inla.stack so you
have to do deconstruct that one or wait for Finn (its has been on the
TODO list for quite some time, but is actually not that straight
forward...)

Best
H
--
Håvard Rue
hr...@r-inla.org

[Helpdesk (Haavard Rue)]

You have to do this without stack, essentially. start easier, make the
mesh nodes at the observational location so you can bypass the A-matrix.
which this works, you can add that in if you like.

On Fri, 2024-01-05 at 04:29 -0800, TheCorinna1994 wrote:
> Thank you for your answer.
> Maybe one last question: I defined a stack for each dataset according
> to https://becarioprecario.bitbucket.io/spde-gitbook/ch-manipula.html
>
> stack_rain_early<-
> inla.stack(data=list(predictor=cbind(data.rain$Rain_max,NA,NA)),
>                                A=list(A_rain_early),
>                                effects=list(c(index_rain,
> list(Intercept1=1))),
>                                    tag="rain_stack_early" )
>
> stack_temp_early<-inla.stack(data=list(predictor2=cbind(NA,
> data.temp$Temp_mean, NA)),
>                        A=list(A_temp_early),
>                        effects=list(c(index_temp, index_rain_copy
> ,list(Intercept2=1))),
>                        tag="temp_stack_early")
>
> stack_stations_early<-inla.stack(data=list(predictor=cbind(NA,NA,
> data.mon$Rain_max)),
>                            A=list(A_stations_early,1),
>                            effects=list(c(index_max, index_temp_copy,
> index_rain_copy2,list(Intercept3=1)),
>                                         list(covars.data_early)),
>                            tag="stations_stack_early")
>
> stack_all_early<-inla.stack(stack_rain_early,stack_temp_early,
> stack_stations_early)
>
> Then I define different spread and tail parameters (for the rain
> dataset I do not have elevation dependencies):
>
> n.rain.early<-length(data.rain$Rain_max)
> null.matrix_early<-matrix(nrow=n.rain_early, ncol=0)
> spread.x.rain_early<-null.matrix_early
> tail.x.rain_early<-null.matrix_early
>
> spread.x_early<-data.mon$Elevation
> tail.x_early<-data.mon$Elevation
>
> and I rewrite the formula
>
> formula4_max<-list(inla.mdata(predictor, spread.x.rain_early,
> tail.x.rain_early),predictor,inla.mdata(predictor, spread.x_early,
> tail.x_early))  ~-1+  Intercept1+Intercept2+Intercept3+

>   f(field_rain, model=spde.nonstat, group = field_rain.group,
> control.group = list(model="ar1", hyper=rho))+
>   f(field_temp, model=spde.nonstat, group = field_temp.group,
> control.group = list(model="ar1", hyper=rho))+
>   f(field_max, model=spde.nonstat, group = field_max.group,
> control.group = list(model="ar1", hyper=rho))+
>   f(field_rain_copy,copy="field_rain", group = field_rain_copy.group,
> fixed=FALSE, hyper=hyper)+
>   f(field_rain_copy2, copy="field_rain", group =
> field_rain_copy2.group, fixed=FALSE, hyper=hyper)+
>   f(field_temp_copy, copy="field_temp", group = field_temp_copy.group,
> fixed=FALSE, hyper=hyper)
>

> but the error is now: Error in inla.core(formula=formula,
> family=family, contrasts=contrasts, : MPredictor==ny is not TRUE
>
> I thought I deconstructed the stacks already.... maybe you have a
> quick answer. Otherwise I appreciate your help and would wait for more
> information at some later point.
>
> All the best,
> Corinna

he...@r-inla.org

[TheCorinna1994]

Thank you so much. I will definitely try.-

All the best,

Corinna