Using the “spTimer” Package to Model Spatio-Temporal Data in R

The “spTimer” package uses three Bayesian models to fit Spatio-Temporal Data. The data may be given at sparse spatial stations, where observations at each station are considered time series. The package can model the residual spatio-temporal variation to measure uncertainty. It also gives flexibility to customize covariance function selection, the hyper-parameters of the prior distributions and the tuning parameters for the implemented MCMC algorithms.

One of the models is represented as an Independent Gaussian Process (GP):

Zlt = Olt + ϵlt, (2)

Olt = Xltβ + ηlt, (3)

for each l = 1,…,r and t = 1,…,Tl, where we assume that ϵlt and ηlt are independent and each are normally distributed with their respective parameters. Let O denote all the random effects, Olt, l = 1,…,r and t = 1,…,Tl. Let θ = (β, σ2ϵ, σ2η, φ, ν) denote all the parameters of this model and let π(θ) denote the prior distribution that we shall specify later.

Nugget Effect

Throughout, the notation ϵlt = (ϵl(s1,t),…, ϵl(sn,t))′ will be used to denote the so called nugget effect or the pure error term assumed to be independently normally distributed N(0,σ2ϵ In), where σ2ϵ is the unknown pure error variance and In is the identity matrix of order n. The spatio-temporal random effects will be denoted by ηlt = (ηl(s1, t),…,ηl(sn, t))′ and these will be assumed to follow N(0,Ση) independently in time, where Ση = σ2ηSη, σ2η is the site invariant spatial variance and Sη is the spatial correlation matrix obtained from the, often used, general Matern correlation function.

The hierarchical structure of the spatio-temporal models can be represented in three stages:

[data|process,parameter] ~ π(Z/O, θ)

[process|parameter] ~ π(O| θ)

[parameter] ~ π(θ)

The implementation can be done using the following code:


##libraries

library(spTimer)
library(sp)

## Reading the data and setting up the spatial coordinates

X <- read.csv(file="Production1.csv")

n<-10

days<-data.frame(X[,3])

X<-data.frame(X[,4])

N<-67

M<-2

colnames(X)<-c("BOPD")

L <- read.csv(file="Wells.csv")

Coordinates<-data.frame(L[1:n,8:9])

colnames(Coordinates)<-c("Latitude","Longitude")

Coordinates2<-data.frame(Coordinates,Coordinates)

Gk<-Fourier_basisfunctions(M,Coordinates2,n,N)

X<-data.frame(X,Gk)

time.inf<-spT.time(t.series=c(67), segments=1)

#Fitting the Spatio-Temporal model

fittedgibbs<-spT.Gibbs(formula= BOPD ~ ., data=X, model="GP", time.data=time.inf, coords=Coordinates,newdata=NULL, priors=NULL, initials=NULL, nItr=1000, nBurn=500, report=1, tol.dist=0.00005, distance.method="maximum",cov.fnc="matern", scale.transform="NONE", spatial.decay=spT.decay(distribution="FIXED"), annual.aggrn="NONE")

plot(fittedgibbs, residuals=TRUE)

#Extract parameters from the fitted model
fittedphi<-data.frame(fittedgibbs$phip)
fittednu<-data.frame(fittedgibbs$nup)
fittedsig2e<-data.frame(fittedgibbs$sig2ep)
fittedsig2eta<-data.frame(fittedgibbs$sig2etap)

#Plotting the Sigma parameter across iterations (Markov Chain)
plot(seq(1,500,1),fittedsig2e[,1],type="l")

#Plotting the Sigmaeta parameter across iterations (Markov Chain)
plot(seq(1,500,1),fittedsig2eta[,1],type="l")

# Create a grid for predicting
new.coords <- spT.grid.coords(Longitude=c(max(Coordinates[,2]), min(Coordinates[,2])),Latitude=c(max(Coordinates[,1]), min(Coordinates[,1])),by=c(60,60))

new.coords2<-data.frame(new.coords,new.coords)

# Basis functions for the new locations

Gk<-Fourier_basisfunctions(M,new.coords2,3600,N)

# Predicting on the new locations

pred.gp <-
predict(fittedgibbs, newdata=Gk, newcoords=new.coords)

# We can extract prediction samples and parameters from each iteration

str(pred.gp)

 

1 thought on “Using the “spTimer” Package to Model Spatio-Temporal Data in R

  1. Philipp Reply

    Hi! Great post and really interesting topic! Would you have those csv files at hand? Thanks!

Leave a Reply

Your email address will not be published. Required fields are marked *