-
Notifications
You must be signed in to change notification settings - Fork 8
Expand file tree
/
Copy pathspatialProcess.html
More file actions
435 lines (388 loc) · 16.1 KB
/
Copy pathspatialProcess.html
File metadata and controls
435 lines (388 loc) · 16.1 KB
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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"><html xmlns="http://www.w3.org/1999/xhtml"><head><title>R: Estimates a spatial process model.</title>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<link rel="stylesheet" type="text/css" href="R.css" />
</head><body>
<table width="100%" summary="page for spatialProcess {fields}"><tr><td>spatialProcess {fields}</td><td style="text-align: right;">R Documentation</td></tr></table>
<h2>
Estimates a spatial process model.
</h2>
<h3>Description</h3>
<p>For a given covariance function estimates the nugget (sigma^2) and process variance (rho) and the
range parameter (theta) by restricted maximum likelihood and then computes the
spatial model with these estimated parameters.
Other parameters of the covariance
are kept fixed and need to be specified.
</p>
<h3>Usage</h3>
<pre>
spatialProcess(x, y, weights = rep(1, nrow(x)), Z = NULL,
mKrig.args = list(m = 2),
cov.function = "stationary.cov", cov.args = list(Covariance = "Matern",
smoothness = 1), theta = NULL, theta.start = NULL, lambda.start = 0.5,
theta.range = NULL,
abstol = 1e-04, na.rm = TRUE, verbose = FALSE, REML = FALSE, ...)
## S3 method for class 'spatialProcess'
summary(object, ...)
## S3 method for class 'spatialProcess'
print(x, digits = 4, ...)
## S3 method for class 'spatialProcessSummary'
print(x, digits = 4, ...)
## S3 method for class 'spatialProcess'
plot(x, digits = 4, which = 1:4, ...)
</pre>
<h3>Arguments</h3>
<table summary="R argblock">
<tr valign="top"><td><code>x</code></td>
<td>
<p>Observation locations</p>
</td></tr>
<tr valign="top"><td><code>y</code></td>
<td>
<p>Observation values</p>
</td></tr>
<tr valign="top"><td><code>weights</code></td>
<td>
<p>Weights for the error term (nugget) in units of reciprocal variance.</p>
</td></tr>
<tr valign="top"><td><code>Z</code></td>
<td>
<p>A matrix of extra covariates for the fixed part of spatial model.
E.g. elevation for fitting climate data over space. </p>
</td></tr>
<tr valign="top"><td><code>mKrig.args</code></td>
<td>
<p>Arguments passed to the mKrig function.</p>
</td></tr>
<tr valign="top"><td><code>cov.function</code></td>
<td>
<p>A character string giving the name of the covariance function
for the spatial component.</p>
</td></tr>
<tr valign="top"><td><code>cov.args</code></td>
<td>
<p> A list specifying parameters and other components of the covariance function.</p>
</td></tr>
<tr valign="top"><td><code>theta</code></td>
<td>
<p>If not NULL the range parameter for the covariance is fixed at this value. </p>
</td></tr>
<tr valign="top"><td><code>theta.start</code></td>
<td>
<p>Starting value for MLE fitting of the scale (aka range) parameter. If omitted the starting value is taken from a grid search ove theta.</p>
</td></tr>
<tr valign="top"><td><code>lambda.start</code></td>
<td>
<p>Starting value for MLE fitting of the lambda parameter. Note lambda
is the ratio of the nugget variance to the process variance. In code variables this is
<code>sigma^2</code> divided by <code>rho</code>. </p>
</td></tr>
<tr valign="top"><td><code> theta.range </code></td>
<td>
<p>A range for the ML search to estimate theta. Default is based on quantiles of the location pairwise distances.</p>
</td></tr>
<tr valign="top"><td><code>na.rm</code></td>
<td>
<p>If TRUE NAs are removed from the data. </p>
</td></tr>
</table>
<table summary="R argblock">
<tr valign="top"><td><code>object</code></td>
<td>
<p>A spatialProcess object returned from the spatialProcess function.</p>
</td></tr>
</table>
<table summary="R argblock">
<tr valign="top"><td><code>REML</code></td>
<td>
<p> If TRUE the parameters are found by restricted maximum likelihood.</p>
</td></tr>
<tr valign="top"><td><code>verbose</code></td>
<td>
<p>If TRUE print out intermediate information for debugging.</p>
</td></tr>
<tr valign="top"><td><code>...</code></td>
<td>
<p>Any other arguments that will be passed to the <code>mKrig</code> function and interpreted
as additional arguments to the covariance function. E.g. <code>smoothness</code> for the
Matern covariance.
</p>
</td></tr>
<tr valign="top"><td><code>abstol</code></td>
<td>
<p>The absolute tolerance bound used to judge convergence. This is applied
to the difference in log likelihood values. </p>
</td></tr>
<tr valign="top"><td><code>digits</code></td>
<td>
<p>Number of significant digits in printed summary</p>
</td></tr>
<tr valign="top"><td><code>which</code></td>
<td>
<p>The vector 1:4 or any subset of 1:4, giving the plots to draw.
See the description ofthese plots below.</p>
</td></tr>
</table>
<h3>Details</h3>
<p>This function makes many choices for the user in terms of defaults and it is
important to be aware of these.
The spatial model is
</p>
<p>Y.k= P(x.k) + Z(x.k)%*%d2 + g(x.k) + e.k
</p>
<p>where ".k" means subscripted by k, Y.k is the dependent variable
observed at location x.k. P is a low degree polynomial (default is a
linear function in the spatial coordinates) and Z is a matrix of covariates (optional) that enter as a linear model the fixed part. g is a mean zero,
Gaussian stochastic process with a marginal variance of rho and a
scale (or range) parameter, theta. The measurement errors, e.k, are
assumed to be uncorrelated, normally distributed with mean zero and
standard deviation sigma. If weights are supplied then the variance of e is assumed to be <code>sigma^2/ weights</code>.
</p>
<p>Perhaps the most important aspect of this function is that
the range (theta), nugget (sigma**2) and process variance (rho) parameters
for the covariance are estimated by restricted maximum
likelihood and this is the model that is then used for spatial
prediction. Geostatistics usaually refers to sigma**2 + rho as the
"sill" and often these parameters are estimated by variogram fitting rather
than maximum likelihood. To be consistent with spline models and to focus
on the key part of model we reparametrize as lambda= sigma**2/
rho and rho. Thinking about h as the spatial signal and e as the noise lambda
can be interpreted as the noise to signal variance ratio in this spatial
context.(See the comparision with fitting the geoR model in the
examples section.)
</p>
<p>The likelihood and the cross valdiation function
can be concentrated to only depend on lambda and theta and so
in reported the optimiztation of these two criterion we focus
on this form of the parameters. Once lambda and theta are
found, the MLE for rho has a closed form and of course then
sigma is then determined from lambda and rho.
</p>
<p>Often the lambda
parameter is difficult to interpret when covariates and a
linear function of the coordinates is included and also when
the range becomes large relative to the size of the spatial
domain. For this reason it is convenient to report the
effective degrees of freedom (also referred to trA in R code and
the output summaries) associated with the predicted
surface or curve. This measure has a one to one relationship
with lamdba and is easier to interpret. For example an eff
degrees of freedom that is very small suggests that the
surface is rwell represented by a low ordoer
polynomial. Degrees of freedom close to the number of
locations indicates a surface that is close to interpolating
the observations and suggests a small or zero value for the
nugget variance.
</p>
<p>The default covariance model is assumed to follow a Matern
with smoothness set to 1.0. This is implementd using the
<code>stationary.cov</code> covariance that can take a argument for
the form of the covariance, a sill and range parameters and
possibily additional parameter might comtrol the shape.
</p>
<p>See the example below how to switch to another model. (Note
that the exponential is also part of the Matern family with
smoothness set to .5. )
</p>
<p>The parameter estimation is done by <code>MLESpatialProcess</code>
and the returned list from this function is added to the Krig
output object that is returned by this function. The estimate
is a version of maximum likelihood where the observations are
transfromed to remove the fixed linear part of the model. If
the user just wants to fix the range parameter theta then
<code>Krig</code> can be used.
</p>
<p>NOTE: The defaults for the <code>optim</code> function used in MLESpatialProcess are:
</p>
<pre>
list(method = "BFGS",
control=list(fnscale = -1,
ndeps = rep(log(1.1),length(cov.params.start)+1),
abstol = abstol,
maxit = 20))
</pre>
<p>There is always a hazard in providing a simple to use method that
makes many default choices for the spatial model. As in any analysis
be aware of these choices and try alternative models and parameter
values to assess the robustness of your conclusions. Also examine the
residuals to check the adequacy of the fit. See the examples below for
some help in how to do this easily in fields. Also see quilt.plot to
get an quick plot of a spatial field to discern obvious spatial paterns.
</p>
<p><strong>summary</strong> method forms a list with class <code>spatialProcessSummary</code> that hasa subset of information from the output object and also creates a table of the estimates of the linear parameters in the fixed part of the model.
With replacated fields there is an option to estimate different linear parameters for each field
( <code> collapseFixedEffect = FALSE </code> ) and in this case a table is not created because there is more than one estimate. See (<code>Omega</code> and <code>fixedEffectsCov</code>) in the <code>mKrig</code> object to build the standard errors.
</p>
<p><strong>plot</strong> method provides a panel of 4 diagnositic plots of the fit.
Use <code>set.panel(2,2)</code> to see all 4 at once. The third plot gives the likelihood and
GCV functions as a function of lambda evaluated at the global MLE for theta.
This is based on the gird evaluations in the component MLEInfo$MLEProfileLambda.
The fourth
plot is a profile likelihood trace for theta having maximized over lambda and is based on the component MLEInfo$MLEGrid.
</p>
<p><strong>print</strong> method prints the <code>spatialProcessSummary</code> object of the fit, adding some details and explanations.
</p>
<h3>Value</h3>
<p>An object of classes <code>mKrig</code> and <code>SpatialProcess</code>. The main difference
from mKrig is an extra component, <code>MLEInfo</code> that has the results of the profile likelihood grid
evaluation over theta ( having maximizing lamdba), joint maximization over theta and lambda,
and a grid evaluation over lambda with theta fixed at its MLE.
</p>
<h3>Author(s)</h3>
<p>Doug Nychka
</p>
<h3>See Also</h3>
<p>Tps, MLESpatialProcess, mKrigMLEGrid, mKrigMLEJoint, plot.Krig, predict.mKrig, predictSE.mKrig
</p>
<h3>Examples</h3>
<pre>
data( ozone2)
# x is a two column matrix where each row is a location in lon/lat
# coordinates
x<- ozone2$lon.lat
# y is a vector of ozone measurements at day 16 a the locations.
y<- ozone2$y[16,]
obj<- spatialProcess( x, y)
# summary of model
summary( obj)
# diagnostic plots
set.panel(2,2)
plot(obj)
# plot 1 data vs. predicted values
# plot 2 residuals vs. predicted
# plot 3 criteria to select the smoothing
# parameter lambda = sigma^2 / rho
# the x axis has transformed lambda
# in terms of effective degrees of freedom
# to make it easier to interpret
# Note that here the GCV function is minimized
# while the REML is maximzed.
# plot 4 the log profile likelihood used to
# determine theta.
#
# predictions on a grid
surface( obj)
#(see also predictSurface for more control on evaluation grid
# and plotting)
#
## Not run:
# working with covariates and filling in missing station data
# using an ensemble method
# see the example under help(sim.spatialProcess) to see how to
# handle a conditional simulation on a grid of predictions with
# covariates.
data(COmonthlyMet)
fit1E<- spatialProcess(CO.loc,CO.tmin.MAM.climate, Z=CO.elev,
theta.range= c(.25, 2.0) )
set.panel( 2,2)
plot( fit1E)
# conditional simulation at missing data
notThere<- is.na(CO.tmin.MAM.climate )
xp <- CO.loc[notThere,]
Zp <- CO.elev[notThere]
infill<- sim.spatialProcess( fit1E, xp=xp,
Z= Zp, M= 10)
#
# interpretation is that these infilled values are all equally plausible
# given the observations and also given the estimated covariance model
#
# for extra credit one could now standardized the infilled values to have
# conditional mean and variance from the exact computations
# e.g. predict( fit1E, xp=CO.loc[!good,], Z= CO.elev[!good])
# and predictSE(fit1E, xp=CO.loc[!good,], Z= CO.elev[!good])
# with these standardization one would still preserve the correlations
# among the infilled values that is also important for considering them as a
# multivariate prediction.
# conditional simulation on a grid but not using the covariate of elevation
fit2<- spatialProcess(CO.loc,CO.tmin.MAM.climate,
theta.range= c(.25, 2.0) )
# note larger range parameter
# create 2500 grids using handy function
gridList <- fields.x.to.grid( fit2$x, nx=50,ny=50)
xGrid<- make.surface.grid( gridList)
ensemble<- sim.spatialProcess( fit2, xp=xGrid, M= 5)
# this is an "n^3" computation so increasing the grid size
# can slow things down for computation
image.plot( as.surface( xGrid, ensemble[1,]))
set.panel()
## End(Not run)
## Not run:
data( ozone2)
# x is a two column matrix where each row is a location in lon/lat
# coordinates
x<- ozone2$lon.lat
# y is a vector of ozone measurements at day 16 a the locations.
y<- ozone2$y[16,]
# a comparison to using an exponential and Wendland covariance function
# and great circle distance -- just to make range easier to interpret.
obj <- spatialProcess( x, y,
Distance = "rdist.earth")
obj2<- spatialProcess( x, y,
cov.args = list(Covariance = "Exponential"),
Distance = "rdist.earth" )
obj3<- spatialProcess( x, y,
cov.args = list(Covariance = "Wendland",
dimension = 2,
k = 2),
Distance = "rdist.earth")
# obj2 could be also be fit using the argument:
# cov.args = list(Covariance = "Matern", smoothness=.5)
#
# Note very different range parameters - BTW these are in miles
# but similar nugget variances.
obj$pars
obj2$pars
obj3$pars
# since the exponential is Matern with smoothness == .5 the first two
# fits can be compared in terms of their likelihoods
# the REML value is slightly higher for obj verses obj2 (598.4 > 596.7)
# these are the _negative_ log likelihoods so suggests a preference for the
# exponential model
#
# does it really matter in terms of spatial prediction?
set.panel( 3,1)
surface( obj)
US( add=TRUE)
title("Matern sm= 1.0")
surface( obj2)
US( add=TRUE)
title("Matern sm= .5")
surface( obj3)
US( add=TRUE)
title("Wendland k =2")
# prediction standard errors
# these take a while because prediction errors are based
# directly on the Kriging weight matrix
# see mKrig for an alternative.
set.panel( 2,1)
out.p<- predictSurfaceSE( obj, nx=40,ny=40)
surface( out.p)
US( add=TRUE)
title("Matern sm= 1.0")
points( x, col="magenta")
#
out.p<- predictSurfaceSE( obj, nx=40,ny=40)
surface( out.p)
US( add=TRUE)
points( x, col="magenta")
title("Matern sm= .5")
## End(Not run)
set.panel(1,1)
## Not run:
### comparison with GeoR
data(ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
good<-!is.na(y)
x1<- x[good,]
y1<- y[good]
obj<- spatialProcess( x, y, mKrig.args= list(m=1), smoothness = .5 )
library( geoR)
ml.n <- likfit(coords= x1, data=y1, ini = c(570, 3), nug = 50)
# compare to
stuffFields<- obj$MLEInfo$MLEJoint$summary[c(1,3,4,5)]
stuffGeoR<- c( ml.n$loglik, ml.n$phi, sqrt(ml.n$nugget),ml.n$sigmasq)
test.for.zero( max(stuffFields/stuffGeoR), 1, tol=.004)
## End(Not run)
</pre>
<hr /><div style="text-align: center;">[Package <em>fields</em> version 9.9 <a href="00Index.html">Index</a>]</div>
</body></html>