tests/simLgcp.R

In geostatsp: Geostatistical Modelling with Likelihood and Bayes

library('geostatsp')
inlaSetThreads(2)

# exclude this line to use the RandomFields package
options(useRandomFields = FALSE)

mymodel = c(mean=-1.5, variance=1, 
				range=2, shape=2)

myraster = rast(nrows=15,ncols=15,xmin=0,xmax=10,ymin=0,ymax=10)

# some covariates, deliberately with a different resolution than myraster
covA = covB = myoffset = rast(ext(myraster), 10, 10)
values(covA) = as.vector(matrix(1:10, 10, 10))
values(covB) = as.vector(matrix(1:10, 10, 10, byrow=TRUE))
values(myoffset) = round(seq(-1, 1, len=ncell(myoffset)))

myCovariate = list(a=covA, b=covB, offsetFooBar = myoffset)

set.seed(0)
myLgcp=simLgcp(mymodel, myCovariate, 
    betas=c(a=-0.1, b=0.25), 
	offset='offsetFooBar',
	rasterTemplate=myraster)

if(inlaAvailable()) {
res = lgcp(data=myLgcp$events, 
		formula = ~ a + b + offset(offsetFooBar),
		grid=squareRaster(myoffset, 15), 
		covariates=myCovariate,
		prior=list(sd=0.2, range=0.4))

if(length(res$parameters)) {
res$parameters$summary[,c(1,3,5)]

lgcpRoc =  spatialRoc(
  fit=res, 
	rr=c(2,3), 
	truth=myLgcp, 
	random=FALSE)

head(lgcpRoc)

plot(lgcpRoc[,'onemspec'] , 
	lgcpRoc[,'2'], 
	type='o', 
	xlim=c(0,1), ylim=c(0,1),
	ylab='sensitivity', xlab='1-specificity'
)
}
}