Simulation1

Download

### Simulation study similar to that in section 3 of
### "Geostatistical inference under preferential sampling" Diggle et al. 2009

####Load the necessary libraries
library(INLA)
library(RandomFields)

### Define the model parameters the grid size and the
### number of data points to be simulated
model <- "exponential"
mean <- 4
variance <- 1.5
nugget <- 0
scale <- 0.15
beta=2
sd.data=0

nrow=50
ncol=50
x <- seq(0, 1, length.out=nrow)
y <- seq(0, 1, length.out=ncol)
n.data = 100
n = nrow*ncol

###simulate the latent field

S0 <- GaussRF(x=x, y=y, model=model, grid=TRUE,
param=c(mean, variance, nugget, scale))
S0.vec = inla.matrix2vector(S0)

### map to probability
p = exp(beta*S0.vec) / max(beta*exp(S0.vec))

### SAMPLE THE DATA USING DIFFERENT DESIGNS
###Completely random sampling design

### sample points
points.rand = rep(0,n)
points.rand[ sample(1:n, size=n.data, prob=rep(1/n,n),replace = FALSE)] = 1

### observe the marks
y.rand = rep(NA,n)
pp.rand = which(points.rand == 1)
y.rand[pp.rand] = rnorm(n.data, mean = S0.vec[pp.rand], sd=sd.data)

###Preferential sample design
### sample points
points.pref = rep(0,n)
points.pref[ sample(1:n, size=n.data, prob=p, replace = FALSE) ] = 1

### observe the marks
y.pref = rep(NA,n)
pp.pref = which(points.pref == 1)
y.pref[pp.pref] = rnorm(n.data, mean = S0.vec[pp.pref], sd=sd.data)

### Cluster design

###simulate an independent realization of S
S1 = GaussRF(x=x, y=y, model=model, grid=TRUE,
param=c(mean, variance, nugget, scale))
S1.vec = inla.matrix2vector(S1)
p1 = exp(beta*S1.vec) / max(beta*exp(S1.vec))

### sample points
points.clust = rep(0,n)
points.clust[ sample(1:n, size=n.data, prob=p1, replace = FALSE) ] = 1

### observe the marks
y.clust = rep(NA,n)
pp.clust = which(points.clust == 1)
y.clust[pp.clust] = rnorm(n.data, mean = S0.vec[pp.clust], sd=sd.data)


###Plot the data locations
loc.pref=rep(NA,n)
loc.pref[pp.pref]=1
loc.rand=rep(NA,n)
loc.rand[pp.rand]=1
loc.clust=rep(NA,n)
loc.clust[pp.clust]=1
dev.new()
par(mfrow=c(2,2))
image(S0,col=grey(seq(0,1,len=256)))
image(inla.vector2matrix(loc.pref,nrow,ncol),add=T)
title("Preferential sampling")

image(S0,col=grey(seq(0,1,len=256)))
image(inla.vector2matrix(loc.rand,nrow,ncol),add=T)
title("Random sampling")

image(S0,col=grey(seq(0,1,len=256)))
image(inla.vector2matrix(loc.clust,nrow,ncol),add=T)
title("Cluster sampling")

### Fit a simple geostatistical model as in equation (1) in Diggle et
### al as a model for the latent random field we use first a RW2 and
### then a matern field

### Prepare the data sets
data.simple.pref=list(y=y.pref,ii=seq(1,n))
data.simple.rand=list(y=y.rand,ii=seq(1,n))
data.simple.clust=list(y=y.clust,ii=seq(1,n))

###define the model formula (here we use a RW2 prior)
formula.simple = y~f(ii,model="rw2d",nrow=nrow,ncol=ncol)

###fit the three models
simple.model.pref = inla(formula.simple,family="gaussian",data=data.simple.pref,
control.predictor = list(compute = TRUE), verbose=T)

simple.model.rand = inla(formula.simple,family="gaussian",data=data.simple.rand,
control.predictor = list(compute = TRUE),verbose=T)

simple.model.clust = inla(formula.simple,family="gaussian",data=data.simple.clust,
control.predictor = list(compute = TRUE),verbose=T)


###define the model formula (here we use a matern prior and fix the range to be 10)
formula.simple1 = y~f(ii,model="matern2d",nrow=nrow,ncol=ncol,initial=c(1,log(10)),
fixed=c(FALSE,TRUE))

###fit the three models
simple.model.pref1 = inla(formula.simple1,family="gaussian",data=data.simple.pref,
control.predictor = list(compute = TRUE), verbose=T)

simple.model.rand1 = inla(formula.simple1,family="gaussian",data=data.simple.rand,
control.predictor = list(compute = TRUE),verbose=T)

simple.model.clust1 = inla(formula.simple1,family="gaussian",data=data.simple.clust,
control.predictor = list(compute = TRUE),verbose=T)

### Model with preferential sampling
### as in eq (5) of Diggle et al.
### Prepare the data sets
ii=c(1:n,rep(NA,n))
jj = c(rep(NA,n), 1:n)
alpha = c(rep(0,n), rep(1,n))
mu = c(rep(1,n), rep(0,n))

### preferential sampling
yy.pref = matrix(NA,2*n,2)
yy.pref[1:n,1] = y.pref
yy.pref[n+1:n,2] = points.pref

data.pref.pref = list(yy=yy.pref,mu=mu,ii=ii,jj=jj,alpha=alpha)

### random sampling
yy.rand = matrix(NA,2*n,2)
yy.rand[1:n,1] = y.rand
yy.rand[n+1:n,2] = points.rand

data.pref.rand = list(yy=yy.rand,mu=mu,ii=ii,jj=jj,alpha=alpha)

### cluster sampling
yy.clust = matrix(NA,2*n,2)
yy.clust[1:n,1] = y.clust
yy.clust[n+1:n,2] = points.clust

data.pref.clust = list(yy=yy.clust,mu=mu,ii=ii,jj=jj,alpha=alpha)

### define the model formula (here with RW prior)
formula.pref = yy ~ alpha + mu + f(ii, model = "rw2d", nrow=nrow, ncol=ncol,
initial = 5, fixed=c(FALSE), param = c(1,0.001), constr=TRUE, bvalue=1) +
f(jj, copy="ii", fixed=FALSE, param=c(0,0.1)) -1

pref.model.pref = inla(formula.pref, family = c("gaussian", "poisson"),
control.family = list(list(initial = 4, fixed=FALSE), list()),
control.inla= list(strategy = "gaussian"),
data = data.pref.pref, verbose = TRUE)

pref.model.rand = inla(formula.pref, family = c("gaussian", "poisson"),
control.family = list(list(initial = 4, fixed=FALSE), list()),
control.inla= list(strategy = "gaussian"),
data = data.pref.rand, verbose = TRUE)

pref.model.clust = inla(formula.pref, family = c("gaussian", "poisson"),
control.family = list(list(initial = 4, fixed=FALSE), list()),
control.inla= list(strategy = "gaussian"),
data = data.pref.clust, verbose = TRUE)


### define the model formula (here with Matern prior, range is fixed to be 10)
formula.pref1 = yy ~ alpha + mu + f(ii, model = "matern2d", nrow=nrow, ncol=ncol,
initial = c(3, log(10)), fixed=c(FALSE,TRUE),
param = c(1,0.001, NA,NA), constr=TRUE) +
f(jj, copy="ii", fixed=FALSE, param=c(0,0.1)) -1

pref.model.pref1 = inla(formula.pref1, family = c("gaussian", "poisson"),
control.family = list(list(initial = 4, fixed=FALSE), list()),
control.inla= list(strategy = "gaussian"),
data = data.pref.pref, verbose = TRUE)

pref.model.rand1 = inla(formula.pref1, family = c("gaussian", "poisson"),
control.family = list(list(initial = 4, fixed=FALSE), list()),
control.inla= list(strategy = "gaussian"),
data = data.pref.rand, verbose = TRUE)

pref.model.clust1 = inla(formula.pref1, family = c("gaussian", "poisson"),
control.family = list(list(initial = 4, fixed=FALSE), list()),
control.inla= list(strategy = "gaussian"),
data = data.pref.clust, verbose = TRUE)

Comments