maternGmrfPrec: Precision matrix for a Matern spatial correlation

In geostatsp: Geostatistical Modelling with Likelihood and Bayes

Precision matrix for a Matern spatial correlation

Description

Produces the precision matrix for a Gaussian random field on a regular square lattice, using a Markov random field approximation.

Usage

 
maternGmrfPrec(N, ...)
## S3 method for class 'dgCMatrix'
maternGmrfPrec(N, 
	param=c(variance=1, range=1, shape=1, cellSize=1),
  adjustEdges=FALSE,...) 
## Default S3 method:
maternGmrfPrec(N, Ny=N, 	
  param=c(variance=1, range=1, shape=1, cellSize=1),
  adjustEdges=FALSE, ...)
NNmat(N, Ny=N, nearest=3, adjustEdges=FALSE)
## S3 method for class 'SpatRaster'
NNmat(N, Ny=N, nearest=3, adjustEdges=FALSE)
## Default S3 method:
NNmat(N, Ny=N, nearest=3, adjustEdges=FALSE)

Arguments

N	Number of grid cells in the x direction, or a matrix denoting nearest neighbours.
Ny	Grid cells in the y direction, defaults to N for a square grid
param	Vector of model parameters, with named elements: scale, scale parameter for the correlation function; prec, precision parameter; shape, Matern differentiability parameter (0, 1, or 2); and cellSize, the size of the grid cells. Optionally, variance and range can be given in place of prec and scale, when the former are present and the latter are missing the reciprocal of the former are taken.
adjustEdges	If TRUE, adjust the precision matrix so it does not implicitly assume the field takes values of zero outside the specified region. Defaults to FALSE. Can be a character string specifying the parameters to use for the correction, such as 'optimal' or 'optimalShape', with TRUE equivalent to 'theo'
nearest	Number of nearest neighbours to compute
...	Additional arguments passed to maternGmrfPrec.dsCMatrix

Details

The numbering of cells is consistent with the terra package. Cell 1 is the top left cell, with cell 2 being the cell to the right and numbering continuing row-wise.

The nearest neighbour matrix N has: N[i,j]=1 if i=j; takes a value 2 if i and j are first ‘rook’ neighbours; 3 if they are first ‘bishop’ neighbours; 4 if they are second ‘rook’ neighbours; 5 if ‘knight’ neighbours; and 6 if third ‘rook’ neighbours.

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    0    0    6    0    0    0
[2,]    0    0    5    4    5    0    0
[3,]    0    5    3    2    3    5    0
[4,]    6    4    2    1    2    4    6
[5,]    0    5    3    2    3    5    0
[6,]    0    0    5    4    5    0    0
[7,]    0    0    0    6    0    0    0

Value

A sparse matrix dsCMatrix-class object, containing a precision matrix for a Gaussian random field or (from the NNmat function) a matrix denoting neighbours.

Examples

# produces the matrix above
	matrix(NNmat(11, 11, nearest=5)[,11*5+6],11, 11)

	params=c(range = 3,	shape=2, variance=5^2)
	
	myGrid = squareRaster(ext(0,20,0,10), 40)
		
	# precision matrix without adjusting for edge effects
	precMat =maternGmrfPrec(N=myGrid, param=params) 
	
	attributes(precMat)$info$precisionEntries
	
	midcell = cellFromRowCol(myGrid, 
		round(nrow(myGrid)/2), round(ncol(myGrid)/2)) # the middle cell
	edgeCell = cellFromRowCol(myGrid, 5,5)# cell near corner

# show precision of middle cell 
	precMid=matrix(precMat[,midcell], 
		nrow(myGrid), ncol(myGrid), byrow=TRUE)

	precMid[round(nrow(precMid)/2) + seq(-5, 5), 
		round(ncol(precMid)/2) + seq(-3, 3)]

	# and with the adjustment
	precMatCorr =maternGmrfPrec(
		N = myGrid, param=params, 
		adjustEdges=TRUE) 

	

# variance matrices
	varMat = Matrix::solve(precMat)
	varMatCorr = Matrix::solve(precMatCorr)

# compare covariance matrix to the matern
	xseq = seq(-ymax(myGrid), ymax(myGrid), len=1000)/1.5
	plot(xseq, matern(xseq, param=params),
	 type = 'l',ylab='cov', xlab='dist',
	 ylim=c(0, params["variance"]*1.1),
	 main="matern v gmrf")

	# middle cell
	varMid=matrix(varMat[,midcell], 
		nrow(myGrid), ncol(myGrid), byrow=TRUE)
	varMidCorr=matrix(varMatCorr[,midcell], 
		nrow(myGrid), ncol(myGrid), byrow=TRUE)
	xseqMid = yFromRow(myGrid) - yFromCell(myGrid, midcell)	
	points(xseqMid, varMid[,colFromCell(myGrid, midcell)], 
		col='red')
	points(xseqMid, varMidCorr[,colFromCell(myGrid, midcell)],
		 col='blue', cex=0.5)

	# edge cells
	varEdge=matrix(varMat[,edgeCell], 
	  nrow(myGrid), ncol(myGrid), byrow=TRUE)
	varEdgeCorr = matrix(varMatCorr[,edgeCell], 
	  nrow(myGrid), ncol(myGrid), byrow=TRUE)
	xseqEdge = yFromRow(myGrid) - yFromCell(myGrid, edgeCell)
	points(xseqEdge, 
		varEdge[,colFromCell(myGrid, edgeCell)], 
		pch=3,col='red')
	points(xseqEdge, 
	  varEdgeCorr[,colFromCell(myGrid, edgeCell)], 
	  pch=3, col='blue')
	
	legend("topright", lty=c(1, NA, NA, NA, NA), 
	  pch=c(NA, 1, 3, 16, 16),
		col=c('black','black','black','red','blue'),
		legend=c('matern', 'middle','edge','unadj', 'adj')
		)


	# construct matern variance matrix

	myraster = attributes(precMat)$raster
	covMatMatern = matern(myraster, param=params)
 
 	prodUncor = crossprod(covMatMatern, precMat)
 	prodCor = crossprod(covMatMatern, precMatCorr)

 	quantile(Matrix::diag(prodUncor),na.rm=TRUE)
 	quantile(Matrix::diag(prodCor),na.rm=TRUE)
 	
 	quantile(prodUncor[lower.tri(prodUncor,diag=FALSE)],na.rm=TRUE)	
 	quantile(prodCor[lower.tri(prodCor,diag=FALSE)],na.rm=TRUE)	

 	

geostatsp documentation built on Sept. 23, 2023, 1:06 a.m.