bayesGeostatExact: Simple Bayesian spatial linear model with fixed semivariogram...

Description Usage Arguments Value Author(s) Examples

View source: R/bayesRegression.R

Description

Given a observation coordinates and fixed semivariogram parameters the bayesGeostatExact function fits a simple Bayesian spatial linear model.

Usage

1
2
3
4
5
  bayesGeostatExact(formula, data = parent.frame(), n.samples,
                     beta.prior.mean, beta.prior.precision,
                     coords, cov.model="exponential", phi, nu, alpha,
                     sigma.sq.prior.shape, sigma.sq.prior.rate,
                     sp.effects=TRUE, verbose=TRUE, ...)

Arguments

formula

for a univariate model, this is a symbolic description of the regression model to be fit. See example below.

data

an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which spLM is called.

n.samples

the number of posterior samples to collect.

beta.prior.mean

beta multivariate normal mean vector hyperprior.

beta.prior.precision

beta multivariate normal precision matrix hyperprior.

coords

an n x 2 matrix of the observation coordinates in R^2 (e.g., easting and northing).

cov.model

a quoted key word that specifies the covariance function used to model the spatial dependence structure among the observations. Supported covariance model key words are: "exponential", "matern", "spherical", and "gaussian". See below for details.

phi

the fixed value of the spatial decay.

nu

if cov.model is "matern" then the fixed value of the spatial process smoothness must be specified.

alpha

the fixed value of the ratio between the nugget tau.sq and partial-sill sigma.sq parameters from the specified cov.model.

sigma.sq.prior.shape

sigma.sq (i.e., partial-sill) inverse-Gamma shape hyperprior.

sigma.sq.prior.rate

sigma.sq (i.e., partial-sill) inverse-Gamma 1/scale hyperprior.

sp.effects

a logical value indicating if spatial random effects should be recovered.

verbose

if TRUE, model specification and progress of the sampler is printed to the screen. Otherwise, nothing is printed to the screen.

...

currently no additional arguments.

Value

An object of class bayesGeostatExact, which is a list with the following tags:

p.samples

a coda object of posterior samples for the defined parameters.

sp.effects

a matrix that holds samples from the posterior distribution of the spatial random effects. The rows of this matrix correspond to the n point observations and the columns are the posterior samples.

args

a list with the initial function arguments.

Author(s)

Sudipto Banerjee [email protected],
Andrew O. Finley [email protected]

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
## Not run: 

data(FBC07.dat)
Y <- FBC07.dat[1:150,"Y.2"]
coords <- as.matrix(FBC07.dat[1:150,c("coord.X", "coord.Y")])

n.samples <- 500
n = length(Y)
p = 1

phi <- 0.15
nu <- 0.5

beta.prior.mean <- as.matrix(rep(0, times=p))
beta.prior.precision <- matrix(0, nrow=p, ncol=p)

alpha <- 5/5

sigma.sq.prior.shape <- 2.0
sigma.sq.prior.rate <- 5.0

##############################
##Simple linear model with
##the default exponential
##spatial decay function
##############################
set.seed(1)
m.1 <- bayesGeostatExact(Y~1, n.samples=n.samples,
                          beta.prior.mean=beta.prior.mean,
                          beta.prior.precision=beta.prior.precision,
                          coords=coords, phi=phi, alpha=alpha,
                          sigma.sq.prior.shape=sigma.sq.prior.shape,
                          sigma.sq.prior.rate=sigma.sq.prior.rate)



print(summary(m.1$p.samples))

##Requires MBA package to
##make surfaces
library(MBA)
par(mfrow=c(1,2))
obs.surf <-
  mba.surf(cbind(coords, Y), no.X=100, no.Y=100, extend=T)$xyz.est
image(obs.surf, xaxs = "r", yaxs = "r", main="Observed response")
points(coords)
contour(obs.surf, add=T)

w.hat <- rowMeans(m.1$sp.effects)
w.surf <-
  mba.surf(cbind(coords, w.hat), no.X=100, no.Y=100, extend=T)$xyz.est
image(w.surf, xaxs = "r", yaxs = "r", main="Estimated random effects")
points(coords)
contour(w.surf, add=T)


##############################
##Simple linear model with
##the matern spatial decay
##function. Note, nu=0.5 so
##should produce the same
##estimates as m.1
##############################
set.seed(1)
m.2 <- bayesGeostatExact(Y~1, n.samples=n.samples,
                          beta.prior.mean=beta.prior.mean,
                          beta.prior.precision=beta.prior.precision,
                          coords=coords, cov.model="matern",
                          phi=phi, nu=nu, alpha=alpha,
                          sigma.sq.prior.shape=sigma.sq.prior.shape,
                          sigma.sq.prior.rate=sigma.sq.prior.rate)

print(summary(m.2$p.samples))

##############################
##This time with the
##spherical just for fun
##############################
m.3 <- bayesGeostatExact(Y~1, n.samples=n.samples,
                          beta.prior.mean=beta.prior.mean,
                          beta.prior.precision=beta.prior.precision,
                          coords=coords, cov.model="spherical",
                          phi=phi, alpha=alpha,
                          sigma.sq.prior.shape=sigma.sq.prior.shape,
                          sigma.sq.prior.rate=sigma.sq.prior.rate)

print(summary(m.3$p.samples))

##############################
##Another example but this
##time with covariates
##############################
data(FORMGMT.dat)

n = nrow(FORMGMT.dat)
p = 5 ##an intercept an four covariates

n.samples <- 50

phi <- 0.0012

coords <- cbind(FORMGMT.dat$Longi, FORMGMT.dat$Lat)
coords <- coords*(pi/180)*6378

beta.prior.mean <- rep(0, times=p)
beta.prior.precision <- matrix(0, nrow=p, ncol=p)

alpha <- 1/1.5

sigma.sq.prior.shape <- 2.0
sigma.sq.prior.rate <- 10.0

m.4 <-
  bayesGeostatExact(Y~X1+X2+X3+X4, data=FORMGMT.dat, n.samples=n.samples,
                     beta.prior.mean=beta.prior.mean,
                     beta.prior.precision=beta.prior.precision,
                     coords=coords, phi=phi, alpha=alpha,
                     sigma.sq.prior.shape=sigma.sq.prior.shape,
                     sigma.sq.prior.rate=sigma.sq.prior.rate)

print(summary(m.4$p.samples))



##Requires MBA package to
##make surfaces
library(MBA)
par(mfrow=c(1,2))
obs.surf <-
  mba.surf(cbind(coords, resid(lm(Y~X1+X2+X3+X4, data=FORMGMT.dat))),
                 no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(obs.surf, xaxs = "r", yaxs = "r", main="Observed response")
points(coords)
contour(obs.surf, add=T)

w.hat <- rowMeans(m.4$sp.effects)
w.surf <-
  mba.surf(cbind(coords, w.hat), no.X=100, no.Y=100, extend=TRUE)$xyz.est
image(w.surf, xaxs = "r", yaxs = "r", main="Estimated random effects")
contour(w.surf, add=T)
points(coords, pch=1, cex=1)



## End(Not run)

Example output

Loading required package: coda
Loading required package: magic
Loading required package: abind
Loading required package: Formula
-------------------------------------------------
	General model description
-------------------------------------------------
Model fit with 150 observations.
Number of covariates 1 (including intercept if specified).
Using the exponential spatial correlation model.

-------------------------------------------------
		Sampling
-------------------------------------------------
Sampled: 500 of 500, 100%
-------------------------------------------------
	Recovering spatial effects
-------------------------------------------------

Iterations = 1:500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 500 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

               Mean     SD Naive SE Time-series SE
(Intercept) -0.3157 0.6722  0.03006        0.03006
sigma.sq     5.0903 0.5919  0.02647        0.02647
tau.sq       5.0903 0.5919  0.02647        0.02647

2. Quantiles for each variable:

              2.5%     25%     50%    75%  97.5%
(Intercept) -1.723 -0.8045 -0.2886 0.1606 0.9405
sigma.sq     4.131  4.6694  5.0527 5.4517 6.2817
tau.sq       4.131  4.6694  5.0527 5.4517 6.2817

-------------------------------------------------
	General model description
-------------------------------------------------
Model fit with 150 observations.
Number of covariates 1 (including intercept if specified).
Using the matern spatial correlation model.

-------------------------------------------------
		Sampling
-------------------------------------------------
Sampled: 500 of 500, 100%
-------------------------------------------------
	Recovering spatial effects
-------------------------------------------------

Iterations = 1:500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 500 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

               Mean     SD Naive SE Time-series SE
(Intercept) -0.3157 0.6722  0.03006        0.03006
sigma.sq     5.0903 0.5919  0.02647        0.02647
tau.sq       5.0903 0.5919  0.02647        0.02647

2. Quantiles for each variable:

              2.5%     25%     50%    75%  97.5%
(Intercept) -1.723 -0.8045 -0.2886 0.1606 0.9405
sigma.sq     4.131  4.6694  5.0527 5.4517 6.2817
tau.sq       4.131  4.6694  5.0527 5.4517 6.2817

-------------------------------------------------
	General model description
-------------------------------------------------
Model fit with 150 observations.
Number of covariates 1 (including intercept if specified).
Using the spherical spatial correlation model.

-------------------------------------------------
		Sampling
-------------------------------------------------
Sampled: 500 of 500, 100%
-------------------------------------------------
	Recovering spatial effects
-------------------------------------------------

Iterations = 1:500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 500 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

               Mean     SD Naive SE Time-series SE
(Intercept) -0.3971 0.3274  0.01464        0.01464
sigma.sq     4.6919 0.5533  0.02475        0.02475
tau.sq       4.6919 0.5533  0.02475        0.02475

2. Quantiles for each variable:

              2.5%     25%     50%     75%  97.5%
(Intercept) -1.029 -0.6051 -0.4133 -0.1895 0.2561
sigma.sq     3.705  4.2761  4.6593  5.0797 5.9158
tau.sq       3.705  4.2761  4.6593  5.0797 5.9158

-------------------------------------------------
	General model description
-------------------------------------------------
Model fit with 1342 observations.
Number of covariates 5 (including intercept if specified).
Using the exponential spatial correlation model.

-------------------------------------------------
		Sampling
-------------------------------------------------
Sampled: 50 of 50, 100%
-------------------------------------------------
	Recovering spatial effects
-------------------------------------------------

Iterations = 1:50
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 50 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

                 Mean        SD  Naive SE Time-series SE
(Intercept) 10.918589 0.1945425 2.751e-02      3.895e-02
X1           0.061436 0.0011248 1.591e-04      1.591e-04
X2          -0.002175 0.0003358 4.749e-05      4.935e-05
X3          -0.051269 0.0104243 1.474e-03      1.474e-03
X4          -0.070742 0.0096346 1.363e-03      1.363e-03
sigma.sq     0.049152 0.0018755 2.652e-04      2.652e-04
tau.sq       0.032768 0.0012503 1.768e-04      1.768e-04

2. Quantiles for each variable:

                 2.5%       25%       50%       75%     97.5%
(Intercept) 10.545826 10.799551 10.900222 11.019453 11.281954
X1           0.059535  0.060639  0.061282  0.062114  0.063950
X2          -0.002733 -0.002371 -0.002206 -0.001965 -0.001486
X3          -0.071007 -0.057915 -0.051231 -0.045231 -0.033330
X4          -0.088345 -0.076629 -0.070915 -0.065610 -0.054462
sigma.sq     0.046031  0.047729  0.048875  0.050444  0.052985
tau.sq       0.030688  0.031819  0.032583  0.033629  0.035323

spBayes documentation built on July 20, 2017, 1:02 a.m.