Baseline (random effects only) areal model crash

19 views
Skip to first unread message

valentina ndolo

unread,
May 18, 2022, 5:24:43 AM5/18/22
to R-inla discussion group
Hi,

I am running a spatial-temporal model using areal data (sub-county level) and BYM correlation. When I run the full model with the covariates + spatial random effects + temporal trend, the model runs successfully. But when I remove the covariates to run a baseline (random effects only model), then the program crashes. 

Note: my dataset has 52,200 observations spanning 15 years. When I subset the data to include only 12 years instead of 15, the problem goes away.

Please see the code below for the full mode and baseline model:

# Add Covariates Plus SRF + Annual Trend
I5.b <- inla(formula = cases01 ~ Total.Population.std + Population.Density.std +  
               Dairy.Exotic.Cattle.std + Agricultural.Land.Area.std +
               Livestock.Prod.Households.std +
               f(County,
                 model = "bym",
                 graph = ken.adj,
                 hyper = HyperBYM,
                 adjust.for.con.comp = FALSE,
                 constr = TRUE,
                 scale.model = TRUE) +
               f(Year,model = "rw2") +
               f(Year1,model = "iid"),
             family = "binomial",
             data = train,
             control.compute = list(dic = TRUE, waic = TRUE, cpo=TRUE))

# Baseline model (without covariates)
I4.b <- inla(formula = cases01 ~ 1 +
               f(County,
                 model = "bym",
                 graph = ken.adj,
                 hyper = HyperBYM,
                 adjust.for.con.comp = FALSE,
                 constr = TRUE,
                 scale.model = TRUE) +
               f(Year,model = "rw2") +
               f(Year1,model = "iid"),
             family = "binomial",
             data = train,
             control.compute = list(dic = TRUE, waic = TRUE, cpo=TRUE))

Many thanks for your help.

Best,
Val.

Thierry Onkelinx

unread,
May 18, 2022, 5:53:05 AM5/18/22
to valentina ndolo, R-inla discussion group
Dear Val,

Why do you have Year both as rw2 and as iid?
I recommend you set sensible PC priors for the random effects. See INLA::inla.doc("pc.prec"). And example is given in INLA::inla.doc("rw1").

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry....@inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////




Op wo 18 mei 2022 om 11:24 schreef valentina ndolo <valenti...@gmail.com>:
--
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/cba22812-9cbf-452a-914b-c0c1b558650an%40googlegroups.com.

valentina ndolo

unread,
May 18, 2022, 6:10:11 AM5/18/22
to Thierry Onkelinx, R-inla discussion group
Hi Thierry,

I have year as both rw2 and iid to account for: the variation due to the temporal trend (rw2) and the unstructured variation between the years that is not explained by the temporal trend (iid).
I also set the following PC priors for the random effects:

# Set priors:
HyperBYM <- list(
  prec.unstruct=list(prior = "pc.prec",param = c(0.01, 0.5)),
  prec.spatial=list(prior = "pc.prec",param = c(0.01, 0.5)))

Best regards,
Val.

Thierry Onkelinx

unread,
May 18, 2022, 6:16:45 AM5/18/22
to valentina ndolo, R-inla discussion group
Dear Val.

You should also set prior for the Year random effect. I'd rather use only rw1 for Year than combining rw2 and iid. If you insist doing the latter, then you really need to pay attention to their priors. Both random effects can model similar trends.

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry....@inbo.be
Havenlaan 88 bus 73, 1000 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data. ~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////




Op wo 18 mei 2022 om 12:10 schreef valentina ndolo <valenti...@gmail.com>:
Reply all
Reply to author
Forward
0 new messages