MCMC algorithm for fitting the model.
spate.mcmc(y,coord=NULL,lengthx=NULL,lengthy=NULL,Sind=NULL,n=NULL,
IncidenceMat=FALSE,x=NULL,SV=c(rho0=0.2,sigma2=0.1,
zeta=0.25,rho1=0.2,gamma=1,alpha=0.3,muX=0,muY=0,tau2=0.005),
betaSV=rep(0,dim(x)[1]),RWCov=NULL,parh=NULL,tPred=NULL,
sPred=NULL,P.rho0=Prho0,P.sigma2=Psigma2,P.zeta=Pzeta,P.rho1=Prho1,
P.gamma=Pgamma,P.alpha=Palpha,P.mux=Pmux,P.muy=Pmuy,P.tau2=Ptau2,
lambdaSV=1,sdlambda=0.01,P.lambda=Plambda,DataModel="Normal",
DimRed=FALSE,NFour=NULL,indEst=1:9,Nmc=10000,BurnIn =1000,
path=NULL,file=NULL,SaveToFile=FALSE,PlotToFile=FALSE,
FixEffMetrop=TRUE,saveProcess=FALSE,Nsave=200,seed=NULL,
Padding=FALSE,adaptive=TRUE,NCovEst=500,BurnInCovEst=500,
MultCov=0.5,printRWCov=FALSE,MultStdDevLambda=0.75,
Separable=FALSE,Drift=!Separable,Diffusion=!Separable,
logInd=c(1,2,3,4,5,9),nu=1,plotTrace=TRUE,
plotHist=FALSE,plotPairs=FALSE,trueVal=NULL,
plotObsLocations=FALSE,trace=TRUE,monitorProcess=FALSE,
tProcess=NULL,sProcess=NULL)The function returns a 'spateMCMC' object with, amongst others, the following entries
Matrix containing samples from the posterior of the hyper-parameters and the regression coefficient
Array with samples from the posterior of the spatio-temporal process
(Estimated) proposal covariance matrix
Observed data in an T x N matrix with columns and rows corresponding to time and space (observations on a grid stacked into a vector), respectively. By default, at each time point, the observations are assumed to lie on a square grid with each axis scaled so that it has unit length.
If specified, this needs to be a matrix of dimension N x 2 with coordinates of the N observation points. Observations in 'y' can either be on a square grid or not. If not, the coordinates of each observation point need to be specified in 'coord'. According to these coordinates, each observation location is then mapped to a grid cell. If 'coord' is not specified, the observations in 'y' are assumed to lie on a square grid with each axis scaled so that it has unit length.
Use together with 'coord' to specify the length of the x-axis. This is usefull if the observations lie in a rectangular area instead of a square. The length needs to be at least as large as the largest x-distance in 'coord.
Use together with 'coord' to specify the length of the y-axis. This is usefull if the observations lie in a rectangular area instead of a square. The length needs to be at least as large as the largest y-distance in 'coord.
Vector of indices of grid cells where observations are made, in case, the observation are not made at every grid cell. Alternatively, the coordinates of the observation locations can be specfied in 'coord'.
Number of point per axis of the square into which the points are mapped. In total, the process is modeled on a grid of size n*n.
Logical; if 'TRUE' an incidence matrix relating the latent process to observation locations is used. This is only recommended to use when the observations are relatively low-dimensional and when the latent process is modeled in a reduced dimensional space as well.
Covariates in an array of dimensions p x T X N, where p denotes the number of covariates, T the number of time points, and N the number of spatial points.
Starting values for parameters. Parameters for the SPDE in the following order: rho_0, sigma^2, zeta, rho_1, gamma, alpha, mu_x, mu_y, tau^2. rho_0 and sigma^2 are the range and marginal variance of the Matern covariance funtion for the innovation term epsilon. zeta is the damping parameter. rho_1, gamma, and alpha parametrize the diffusion matrix with rho_1 being a range parameter, gamma and alpha determining the amount and the direction, respectively, of anisotropy. mu_x and mu_y are the two components of the drift vector. tau^2 denotes the nugget effect or measurment error.
Starting values for regression coefficients.
Covariance matrix of the proposal distribution used in the random walk Metropolis-Hastings step for the hyper-parameters.
Only used in prediction mode. If 'parh' is not 'NULL', this indicates that 'spate.mcmc' is used for making predictions at locations (tPred,sPred) instead of applying the traditional MCMC algorithm. In case 'parh' is not 'NULL', it is a Npar x Nsim matrix containing Nsim samples from the posterior of the Npar parameters. This argument is used by the wrapper function 'spate.predict'.
Time points where predictions are made.This needs to be a vector if predictions are made at multiple times. For instance, if T is the number of time points in the data 'y', then tPred=c(T+1, T+2) means that predictions are made at time 'T+1' and 'T+2'. This argument is used by the wrapper function 'spate.predict'.
Vector of indices of grid cells (positions of locations in the stacked spatial vector) where predictions are made. This argument is used by the wrapper function 'spate.predict'.
Function specifying the prior for rho0.
Function specifying the prior for sigma2.
Function specifying the prior for zeta.
Function specifying the prior for rho1.
Function specifying the prior for gamma.
Function specifying the prior for alpha.
Function specifying the prior for mux.
Function specifying the prior for muy.
Function specifying the prior for tau2.
Starting value for transformation parameter lambda in the Tobit model.
Standard deviation of the proposal distribution used in the random walk Metropolis-Hastings step for lambda.
Function specifying the prior for lambda.
Specifies the data model. "Normal" or "SkewTobit" are available options.
Logical; if 'TRUE' dimension reduction is applied. This means that not the full number (n*n) of Fourier functions is used but rather only a reduced dimensional basis of dimension 'NFour'.
If 'DimRed' is 'TRUE', this specifies the number of Fourier functions.
A vector of numbers specifying which for which parameters the posterior should be computed and which should be held fix (at their starting value). If the corresponding to the index of rho_0, sigma^2, zeta, rho_1, gamma, alpha, mu_x, mu_y, tau^2 is present in the vector, the parameter will be estimated otherwise not. Default is indEst=1:9 which means that one samples from the posterior for all parameters.
Number of MCMC samples.
Length of the burn-in period.
Path, in case plots and / or the spateMCMC object should be save in a file.
File name, in case plots and / or the spateMCMC object should be save in a file.
Indicates whether the spateMCMC object should be save in a file.
Indicates whether the MCMC output analysis plots should be save in a file.
The fixed effects, i.e., the regression coefficients, can either be sampled in a Gibbs step or updated together with the hyperparameters in the Metropolis-Hastings step. The latter is the default and recommended option since correlations between fixed effects and the random process can result in slow mixing.
Logical; if 'TRUE' samples from the posterior of the latent spatio-temporal process xi are saved.
Number of samples from the posterior of the latent spatio-temporal process xi that should be save.
Seed for random generator.
Indicates whether padding is applied or not. If the range parameters are large relative to the domain, this is recommended since otherwise spurious periodicity can occur.
Indicates whether an adaptive Metropolis-Hastings algorithm is used or not. If yes, the proposal covariance matrix 'RWCov' is adaptively estimated during the algorithm and tuning does not need to be done by hand.
Minimal number of samples to be used for estimating the proposal matrix.
Burn-in period for estimating the proposal matrix.
Numeric used as multiplier for the adaptively estimated proposal cocariance matrix 'RWCov' of the hyper-parameters. I.e., the estimated covariance matrix is multiplied by 'MultCov'.
Logical, if 'TRUE' the estimated proposal cocariance matrix is printed each time.
Numeric used as multiplier for the adaptively estimated proposal standard deviation of the Tobit transformation parameter lambda. I.e., the estimated standard deviation is multiplied by 'MultStdDevLambda'.
Indicates whether a separable model, i.e., no transport / drift and no diffusion, should be estimated.
Indicates whether a drift term should be included.
Indicates whether a diffusion term should be included.
Indicates which parameters are sampled on the log-scale. Default is logInd=c(1, 2, 3, 4, 5, 9) corresponding to rho_0, sigma2, zeta, rho_1, gamma, and tau^2.
Smoothness parameter of the Matern covariance function for the innovations. By default this equals 1 corresponding to the Whittle covariance function.
Indicates whether trace plots are made.
Indicates whether histograms of the posterior distributions are made.
Indicates whether scatter plots of the hyper-parameters and the regression coefficients are made.
In simulations, true values can be supplied for comparison with the MCMC output.
Logical; if 'TRUE' the observations locations are ploted together with the grid cells.
Logical; if 'TRUE' tracing information on the progress of the MCMC algorithm is produced.
Logical; if 'TRUE' in addition to the trace plots of the hyper-parameters, the mixing properties of the latent process xi=Phi*alpha is monitored. This is done by plotting the current sample of the process. More specifically, the time series at locations 'sProcess' and the spatial fieldd at time points 'tProcess'.
To be secified if 'monitorProcess=TRUE'. Time points at which spatial fields of the sampled process should be plotted.
To be secified if 'monitorProcess=TRUE'. Locations at which time series of the sampled process should be plotted.
Fabio Sigrist
##Specify hyper-parameters
par <- c(rho0=0.1,sigma2=0.2,zeta=0.5,rho1=0.1,gamma=2,alpha=pi/4,muX=0.2,muY=-0.2,tau2=0.01)
##Simulate data
spateSim <- spate.sim(par=par,n=20,T=20,seed=4)
w <- spateSim$w
##Below is an example to illustrate the use of the MCMC algorithm.
##In practice, more samples are needed for a sufficiently large effective sample size.
##The following takes a couple of minutes.
##Load the precomputed object some lines below to save time.
##spateMCMC <- spate.mcmc(y=w,x=NULL,SV=c(rho0=0.2,sigma2=0.1,
## zeta=0.25,rho1=0.2,gamma=1,alpha=0.3,muX=0,muY=0,tau2=0.005),
## RWCov=diag(c(0.005,0.005,0.05,0.005,0.005,0.001,0.0002,0.0002,0.0002)),
## Nmc=10000,BurnIn=2000,seed=4,Padding=FALSE,plotTrace=TRUE,NCovEst=500,
## BurnInCovEst=500,trueVal=par,saveProcess=TRUE)
##spateMCMC
##plot(spateMCMC.fit,true=par,postProcess=TRUE)
##Instead of waiting, you can also use this precomputed object
data("spateMCMC")
spateMCMC
plot(spateMCMC,true=par,medianHist=FALSE)
Run the code above in your browser using DataLab