Reproducibility of Knorr-Held spatiotemporal models

[vit...@osu.edu]

I was experimenting with the Knorr-Held 

spacetime models, as shown in Ch 7.1 of

Blangiardo+Cameletti. I was a bit surprised to

see that while the Type-I and Type-IV results

reproduced exactly (at least as they are reported

by INLA's summary function), the Type-II and

Type-III did not; and this persisted even after

setting the random seed before running the models.

The differences between runs can be quite substantial.

Should I be concerned about this? Is there something about

the Type-II and Type-III models that would lead one

to expect this behavior, when the Type-I and Type-IV

models don't show it?

I'm running INLA version 20.03.17 on MS Windows-10 (64-bit).

Here's the script, mostly truncated from the

version supplied by Blangiardo + Cameletti

### R code for Chapter 7.1 - Georgia Data

# This does the 4 types of space-time interactions a la Knorr-Held

### Last update: 18/08/2014

# modified by pav for working directory and other folders

# This version just does the Type I-IV models, without any

# intervening stuff or graphics

###################################################

### Set working directory and load packages

###################################################

#remove(list=ls())

my.dir <- "c:/work/R/blangiardo/"

require(INLA)

inla.setOption(scale.model.default=FALSE)

set.seed(12345)

require(maptools)

###################################################

### Code for Section 7.1.1

###################################################

data <- read.csv(paste(my.dir,"LBWGeorgiaData/data.final.csv",sep=""))

georgia <- readShapePoly(paste(my.dir,"LBWGeorgiaData/co13_d00.shp",sep=""))

# Remove duplicates (area 88 and 94)

names <- georgia$NAME

georgia <- georgia[-c(99,100,105,137),]

data.georgia <- attr(georgia, "data")

low.vector <- as.vector(as.matrix(data[,2:12])) #by column

E.vector <- as.vector(as.matrix(data[,13:23]))  #by column

year <- numeric(0)

for(i in 1:11){ 

  year<- append(year,rep(i,dim(data)[1]))

}

# county replicates the county names 11 times

county <- as.factor(rep(data[,1],11))

data <- data.frame(y= low.vector, E= E.vector, ID.area=as.numeric(county), ID.area1=as.numeric(county), year=year,

                   ID.year = year, ID.year1=year, ID.area.year = seq(1,length(county)))

# Spatial graph

Georgia.adj <- paste(my.dir,"LBWGeorgiaData/Georgia.graph",sep="")

# Temporal graph

Temp.adj <- paste(my.dir,"LBWGeorgiaData/Temp.graph",sep="")

data$ID.area.int <- data$ID.area # PAV I've added these to data

data$ID.year.int <- data$ID.year

###################################################

#--- Type I interaction ---#

formula.intI<- y ~ + f(ID.area,model="bym", graph=Georgia.adj) +

  f(ID.year,model="rw2") + f(ID.year1,model="iid") + f(ID.area.year,model="iid")

mod.intI <- inla(formula.intI,family="poisson",data=data,E=E, 

                 control.predictor=list(compute=TRUE), 

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

summary(mod.intI)

#--- Type II interaction ---#

formula.intII<- y ~ f(ID.area,model="bym",graph=Georgia.adj) +

  f(ID.year,model="rw2") + f(ID.year1,model="iid") +

  f(ID.area.int,model="iid", group=ID.year.int,control.group=list(model="rw2")) 

mod.intII <- inla(formula.intII,family="poisson",data=data,E=E, 

                  control.predictor=list(compute=TRUE), 

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

summary(mod.intII)

#--- Type III interaction ---#

formula.intIII<- y ~  f(ID.area,model="bym",graph=Georgia.adj) +

  f(ID.year,model="rw2") + 

  f(ID.year1,model="iid") +

  f(ID.year.int,model="iid", group=ID.area.int,control.group=list(model="besag", graph=Georgia.adj))

mod.intIII <- inla(formula.intIII,family="poisson",data=data,E=E, 

                   control.predictor=list(compute=TRUE), 

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

summary(mod.intIII)

#--- Type IV interaction ---#

formula.intIV<- y ~ f(ID.area,model="bym",graph=Georgia.adj) +

  f(ID.year,model="rw2") +

  f(ID.year1,model="iid") + 

  f(ID.area.int,model="besag", graph=Georgia.adj,group=ID.year.int,control.group=list(model="rw2"))

mod.intIV <- inla(formula.intIV,family="poisson",data=data,E=E, 

                  control.predictor=list(compute=TRUE), 

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

summary(mod.intIV)

[Helpdesk]

first, can you upgrade to a recent version ? (and I guess you use R-4.0)

in what way are they not reproduced ?

> --
> 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/25f76426-3aa4-497c-a61e-7c1fbcde8da8n%40googlegroups.com
> .

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