Problems with SPDE model

[Stephan Thies]

Hi everyone,

I already posted into this group last week to ask for help with a dynamic spatio-temporal AR model. I currently fit an SPDE model following the instructions that I found in several documents and the paper by Cameletti et al (2012). However INLA always crahses when I enlarge my dataset beyond a (to me) unknown threshold and Windows gives the error: "inla.exe has stopped working". After I select "close program" I get the additional warning message from the verbose output:

running command '"\\icnas1.cc.ic.ac.uk/sft16/R/win-library/3.3/INLA/bin/windows/64bit/inla.exe"  -b -s -v  "C:/Users/SFT16/AppData/Local/Temp/RtmpQB7vaS/fileb857be763/Model.ini"' had status 253

I attached the full verbose output of a working model (currently 10 stations, mesh with 25 knots and 478 days) and two verbose outputs of models that do not work (also 10 stations, mesh with 25 knots and 479 days). Do you have any idea why INLA crashes?

Furthermore, I am slightly confused about the results of the model that works. It is a hierarchical model of the form:

Z(s,t)= Y(s,t) + e(s,t); e(s,t) ~ iid Gaussian
Y(s,t)= rho * Y(s,t-1) + eps(s,t); eps(s,t) ~ temporal independent, spatially coloured Gaussian

The empirical autocorrelation in the observed stationwise data z usually lies around 0.7. However, the posterior of the AR coefficient rho is estimated by the model to be very close to 1. Can such a posterior be the result of a wrong SPDE approximation of the GF and is this an unusual estimate?

Thank you for your help!
Stephan

[Finn Lindgren]

Hi,

Many of us were at the same conference last week, so didn't have a chance to look at this fully.

You should probably check if using  inla.spde2.pcmatern instead inla.spde2.matern works better for the small model, as it looks like your parameter estimates walk away in a weird direction. The pcmatern model has a better prior distribution for the parameters.  Can you please state formula code for the model as well, as I suspect there might be some issue there as well. (Also include any inla.stack code etc, and the inla() call itself)

Finn

--
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 post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.

<Crash2.txt>

<Working Model.txt>

<Crash1.txt>

[Stephan Thies]

Hi,

thank you - I will try out the inla.spde2.pcmatern function tomorrow. My code for the model construction is the following:

# Contruct mesh:
ana.mesh<-inla.mesh.2d(loc = coordinates(MS.loc[c(new.control.stations, predict.station)]),
                       loc.domain = coordinates(LDN.union@polygons[[1]]@Polygons[[1]]),
                       max.edge = c(50*1000))

# Mesh to spde
matern.alpha<-1.5
ana.spde <- inla.spde2.matern(mesh = ana.mesh, alpha = matern.alpha)

# A
ana.A <- inla.spde.make.A(mesh=ana.mesh,
                          loc=coordinates(MS.loc[data.ana$code]),
                          group=data.ana$time2,
                          n.group = max(data.ana$time2))

# Produce indices
ana.indices<-inla.spde.make.index("field", n.spde=ana.spde$n.spde, n.group = max(data.ana$time2))

# Prepare data stack for inla:
ana.stack <- inla.stack(tag='ana',
                        data=list(logNOX=data.ana$nox.na),
                        A=list(ana.A,1),
                        effects=list(c(list(intercept = rep(1, ana.mesh$n)),ana.indices),
                                     list(ws=data.ana$ws,
                                          solar = data.ana$solar,
                                          rain = data.ana$rain,
                                          temp=data.ana$temp,
                                          bp = data.ana$bp,
                                          rhum = data.ana$rhum,
                                          code = data.ana$code,
                                          trend = data.ana$trend,
                                          i = data.ana$ID)))

inla.data.ana<-inla.stack.data(ana.stack, SPDE = ana.spde)

# Generate summary list
model.simp<-list()

#Simple Model without covariates (except station specific intercept):
model.simp$formula <- logNOX ~ 0 + code + f(field, model= SPDE, group=field.group, control.group=list(model="ar1")) + f(i, model = "iid")

model.simp$inla <- inla(model.simp$formula,
                        data = inla.data.ana ,
                        family = "gaussian",
                        control.fixed=list(expand.factor.strategy='inla'),
                        control.predictor=list(A=inla.stack.A(ana.stack), compute = T),
                        control.family = list(hyper = list(theta = list(param=c(1,20), fixed = T))),
                        debug = F,
                        verbose = T)

Am Dienstag, 1. August 2017 16:48:42 UTC+1 schrieb Finn Lindgren:

Hi,

Many of us were at the same conference last week, so didn't have a chance to look at this fully.

You should probably check if using  inla.spde2.pcmatern instead inla.spde2.matern works better for the small model, as it looks like your parameter estimates walk away in a weird direction. The pcmatern model has a better prior distribution for the parameters.  Can you please state formula code for the model as well, as I suspect there might be some issue there as well. (Also include any inla.stack code etc, and the inla() call itself)

Finn

On 1 Aug 2017, at 16:12, Stephan Thies <stephan...@gmail.com> wrote:

Hi everyone,

I already posted into this group last week to ask for help with a dynamic spatio-temporal AR model. I currently fit an SPDE model following the instructions that I found in several documents and the paper by Cameletti et al (2012). However INLA always crahses when I enlarge my dataset beyond a (to me) unknown threshold and Windows gives the error: "inla.exe has stopped working". After I select "close program" I get the additional warning message from the verbose output:

running command '"\\icnas1.cc.ic.ac.uk/sft16/R/win-library/3.3/INLA/bin/windows/64bit/inla.exe"  -b -s -v  "C:/Users/SFT16/AppData/Local/Temp/RtmpQB7vaS/fileb857be763/Model.ini"' had status 253

I attached the full verbose output of a working model (currently 10 stations, mesh with 25 knots and 478 days) and two verbose outputs of models that do not work (also 10 stations, mesh with 25 knots and 479 days). Do you have any idea why INLA crashes?

Furthermore, I am slightly confused about the results of the model that works. It is a hierarchical model of the form:

Z(s,t)= Y(s,t) + e(s,t); e(s,t) ~ iid Gaussian
Y(s,t)= rho * Y(s,t-1) + eps(s,t); eps(s,t) ~ temporal independent, spatially coloured Gaussian

The empirical autocorrelation in the observed stationwise data z usually lies around 0.7. However, the posterior of the AR coefficient rho is estimated by the model to be very close to 1. Can such a posterior be the result of a wrong SPDE approximation of the GF and is this an unusual estimate?

Thank you for your help!
Stephan

--
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-group+unsub...@googlegroups.com.

To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.

<Crash2.txt>

<Working Model.txt>

<Crash1.txt>

[Finn Lindgren]

Is the model practically identifiable? Specifically, how does "code" and "i"/"ID" vary across the observations?

Finn

To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.

[Stephan Thies]

I think it should be. "code" is a station specific factor which allows for a station specific intercept (in the smaller model, which works, these coefficients are estimated in a sensible manner).
"i/ID" varies for each observation and the corresponding iid effect is intended to mimic the random noise eps(s,t) since the noise in the control.family argument is fixed to a relatively small value. (Somewhere in the INLA group I read that this is a way to do predictions that include observational noise) However inla also crashes when I remove the iid effect on "i/ID" and allow the noise to be estimated.