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):

Z_{lt} = O_{lt} + ϵ_{lt}, (2)

O_{lt} = X_{lt}β + η_{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 eﬀects, 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 eﬀect 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 eﬀects 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)

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