-
Notifications
You must be signed in to change notification settings - Fork 8
/
Tps.html
540 lines (462 loc) · 19.8 KB
/
Tps.html
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
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
<!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: Thin plate spline regression</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 Tps {fields}"><tr><td>Tps {fields}</td><td style="text-align: right;">R Documentation</td></tr></table>
<h2>
Thin plate spline regression
</h2>
<h3>Description</h3>
<p>Fits a thin plate spline surface to irregularly spaced data. The
smoothing parameter is chosen by generalized cross-validation. The assumed
model is additive Y = f(X) +e where f(X) is a d dimensional surface.
This function also works for just a single dimension and is a special case of a spatial process estimate
(Kriging). A "fast" version of this function uses a compactly supported Wendland covariance and computes the estimate for a fixed smoothing parameter.
</p>
<h3>Usage</h3>
<pre>
Tps(x, Y, m = NULL, p = NULL, scale.type = "range", lon.lat = FALSE,
miles = TRUE, method = "GCV", GCV = TRUE, ...)
fastTps(x, Y, m = NULL, p = NULL, theta, lon.lat=FALSE,
find.trA = TRUE, lambda=0, ...)
</pre>
<h3>Arguments</h3>
<table summary="R argblock">
<tr valign="top"><td><code>x</code></td>
<td>
<p>Matrix of independent variables. Each row is a location or a set of
independent covariates.
</p>
</td></tr>
<tr valign="top"><td><code>Y</code></td>
<td>
<p>Vector of dependent variables.
</p>
</td></tr>
<tr valign="top"><td><code>m</code></td>
<td>
<p>A polynomial function of degree (m-1) will be
included in the model as the drift (or spatial trend) component.
Default is the value such that 2m-d is greater than zero where d is the
dimension of x.
</p>
</td></tr>
<tr valign="top"><td><code>p</code></td>
<td>
<p>Polynomial power for Wendland radial basis functions. Default is 2m-d
where d is the dimension of x.
</p>
</td></tr>
<tr valign="top"><td><code>scale.type</code></td>
<td>
<p>The independent variables and knots are scaled to the specified
scale.type.
By default the scale type is "range", whereby
the locations are transformed
to the interval (0,1) by forming (x-min(x))/range(x) for each x.
Scale type of "user" allows specification of an x.center and x.scale by
the
user. The default for "user" is mean 0 and standard deviation 1. Scale
type of "unscaled" does not scale the data.
</p>
</td></tr>
<tr valign="top"><td><code>theta</code></td>
<td>
<p>The tapering range that is passed to the Wendland compactly
supported covariance. The covariance (i.e. the radial basis function) is
zero beyond range theta. The larger theta the closer this model will approximate the
standard thin plate spline.</p>
</td></tr>
<tr valign="top"><td><code>lon.lat</code></td>
<td>
<p>If TRUE locations are interpreted as lognitude and
latitude and great circle distance is used to find distances among
locations. The theta scale parameter for <code>fast.Tps</code> (setting the
compact support of the Wendland function) in this case is in units of
miles (see example and caution below). </p>
</td></tr>
<tr valign="top"><td><code>method</code></td>
<td>
<p> Determines what "smoothing" parameter should be
used. The default is to estimate standard GCV Other choices are:
GCV.model, GCV.one, RMSE, pure error and REML. The differences
are explained in the Krig help file.</p>
</td></tr>
<tr valign="top"><td><code>GCV</code></td>
<td>
<p>If TRUE the decompositions are done to efficiently evaluate the estimate, GCV function and likelihood at multiple values of lambda. </p>
</td></tr>
<tr valign="top"><td><code>miles</code></td>
<td>
<p>If TRUE great circle distances are in miles if FALSE
distances are in kilometers</p>
</td></tr>
<tr valign="top"><td><code>lambda</code></td>
<td>
<p>Smoothing parameter the ratio of error variance to
process variance, default is zero which corresponds to
interpolation. See fastTpsMLE to estimate this paramter
from the data.</p>
</td></tr>
<tr valign="top"><td><code>find.trA</code></td>
<td>
<p>If TRUE will estimate the effective degrees of freedom
using a simple Monte Carlo method. This will add to the computational
burden by approximately <code>NtrA</code> solutions of the linear system but
the cholesky decomposition is reused.</p>
</td></tr>
<tr valign="top"><td><code>...</code></td>
<td>
<p>For <code>Tps</code> any argument that is valid for the
<code>Krig</code> function. Some of the main ones are listed below.
</p>
<p>For <code>fastTps</code> any argument that is suitable for the <code>mKrig</code>
function see help on mKrig for these choices.
</p>
<p>Arguments for Tps:
</p>
<dl>
<dt>lambda</dt><dd><p> Smoothing parameter that is the ratio of the error
variance (sigma**2) to the scale parameter of the covariance
function. If omitted this is estimated by GCV. </p>
</dd>
<dt>Z</dt><dd><p>Linear covariates to be included in fixed part of the model
that are distinct from the default low order polynomial in
<code>x</code></p>
</dd>
<dt>df</dt><dd><p> The effective number of parameters for the fitted
surface. Conversely, N- df, where N is the total number of
observations is the degrees of freedom associated with the
residuals. This is an alternative to specifying lambda and much
more interpretable.</p>
</dd>
<dt>cost</dt><dd><p> Cost value used in GCV criterion. Corresponds to a
penalty for increased number of parameters. The default is 1.0 and
corresponds to the usual GCV.</p>
</dd>
<dt>weights</dt><dd><p> Weights are proportional to the reciprocal variance
of the measurement error. The default is no weighting i.e. vector
of unit weights. </p>
</dd>
<dt>nstep.cv</dt><dd><p> Number of grid points for minimum GCV search. </p>
</dd>
<dt>x.center</dt><dd><p> Centering values are subtracted from each column of
the x matrix. Must have scale.type="user".</p>
</dd>
<dt>x.scale</dt><dd><p> Scale values that divided into each column after
centering. Must have scale.type="user".</p>
</dd>
<dt>rho</dt><dd><p>Scale factor for covariance. </p>
</dd>
<dt>sigma2</dt><dd><p> Variance of errors or if weights are not equal to 1
the variance is sigma**2/weight.</p>
</dd>
<dt>verbose</dt><dd><p> If true will print out all kinds of intermediate
stuff. </p>
</dd>
<dt>mean.obj</dt><dd><p> Object to predict the mean of the spatial
process. </p>
</dd>
<dt>sd.obj</dt><dd><p> Object to predict the marginal standard deviation of
the spatial process. </p>
</dd>
<dt>null.function</dt><dd><p>An R function that creates the matrices for
the null space model. The default is fields.mkpoly, an R
function that creates a polynomial regression matrix with all
terms up to degree m-1. (See Details) </p>
</dd>
<dt>offset</dt><dd><p> The offset to be used in the GCV
criterion. Default is 0. This would be used when Krig/Tps is
part of a backfitting algorithm and the offset has to be
included to reflect other model degrees of freedom. </p>
</dd>
</dl>
</td></tr>
</table>
<h3>Details</h3>
<p>Both of these functions are special cases of using the
<code>Krig</code> and <code>mKrig</code> functions. See the help on each of these
for more information on the calling arguments and what is returned.
<code>Tps</code> makes use of the stable computations via eigen decompositions in <code>Krig</code>. <code>fastTps</code> follows the more standard computations for spatial statistics centered around the Cholesky decomposition in <code>mKrig</code>.
</p>
<p>A thin plate spline is the result of minimizing the residual sum of
squares subject to a constraint that the function have a certain
level of smoothness (or roughness penalty). Roughness is
quantified by the integral of squared m-th order derivatives. For one
dimension and m=2 the roughness penalty is the integrated square of
the second derivative of the function. For two dimensions the
roughness penalty is the integral of
</p>
<p>(Dxx(f))**22 + 2(Dxy(f))**2 + (Dyy(f))**22
</p>
<p>(where Duv denotes the second partial derivative with respect to u
and v.) Besides controlling the order of the derivatives, the value of
m also determines the base polynomial that is fit to the data.
The degree of this polynomial is (m-1).
</p>
<p>The smoothing parameter controls the amount that the data is
smoothed. In the usual form this is denoted by lambda, the Lagrange
multiplier of the minimization problem. Although this is an awkward
scale, lambda =0 corresponds to no smoothness constraints and the data
is interpolated. lambda=infinity corresponds to just fitting the
polynomial base model by ordinary least squares.
</p>
<p>This estimator is implemented by passing the right generalized covariance
function based on radial basis functions to the more general function
Krig. One advantage of this implementation is that once a Tps/Krig object
is created the estimator can be found rapidly for other data and smoothing
parameters provided the locations remain unchanged. This makes simulation
within R efficient (see example below). Tps does not currently support the
knots argument where one can use a reduced set of basis functions. This is
mainly to simplify the code and a good alternative using knots would be to use a
valid covariance from the Matern family and a large range parameter.
</p>
<p>CAUTION about <code>lon.lat=TRUE</code>: The option to use great circle distance
to define the radial basis functions (<code>lon.lat=TRUE</code>) is very useful
for small geographic domains where the spherical geometry is well approximated by a plane. However, for large domains the spherical distortion be large enough that the basis function no longer define a positive definite system and Tps will report a numerical error. An alternative is to switch to a three
dimensional thin plate spline the locations being the direction cosines. This will
give approximate great circle distances for locations that are close and also the numerical methods will always have a positive definite matrices.
</p>
<p>Here is an example using this idea for <code>RMprecip</code> and also some
examples of building grids and evaluating the Tps results on them:
</p>
<pre>
# a useful function:
dircos<- function(x1){
coslat1 <- cos((x1[, 2] * pi)/180)
sinlat1 <- sin((x1[, 2] * pi)/180)
coslon1 <- cos((x1[, 1] * pi)/180)
sinlon1 <- sin((x1[, 1] * pi)/180)
cbind(coslon1*coslat1, sinlon1*coslat1, sinlat1)}
# fit in 3-d to direction cosines
out<- Tps(dircos(RMprecip$x),RMprecip$y)
xg<-make.surface.grid(fields.x.to.grid(RMprecip$x))
fhat<- predict( out, dircos(xg))
# coerce to image format from prediction vector and grid points.
out.p<- as.surface( xg, fhat)
surface( out.p)
# compare to the automatic
out0<- Tps(RMprecip$x,RMprecip$y, lon.lat=TRUE)
surface(out0)
</pre>
<p>The function <code>fastTps</code> is really a convenient wrapper function that
calls <code>mKrig</code> with the Wendland covariance function. This is
experimental and some care needs to exercised in specifying the taper
range and power ( <code>p</code>) which describes the polynomial behavior of
the Wendland at the origin. Note that unlike Tps the locations are not
scaled to unit range and this can cause havoc in smoothing problems with
variables in very different units. So rescaling the locations <code> x<- scale(x)</code>
is a good idea for putting the variables on a common scale for smoothing.
This function does have the potential to approximate estimates of Tps
for very large spatial data sets. See <code>wendland.cov</code> and help on
the SPAM package for more background.
Also, the function <code>predictSurface.fastTps</code> has been made more efficient for the
case of k=2 and m=2. Also see the handy function <code><a href="mKrigMLE.html">fastTpsMLE</a></code> to estimate lambda by maximum likelihood.
</p>
<p>See also the mKrig function for handling larger data sets and also for an example
of combining Tps and mKrig for evaluation on a huge grid.
</p>
<h3>Value</h3>
<p>A list of class Krig. This includes the
fitted values, the predicted surface evaluated at the
observation locations, and the residuals. The results of the grid
search minimizing the generalized cross validation function are
returned in gcv.grid. Note that the GCV/REML optimization is
done even if lambda or df is given.
Please see the documentation on Krig for details of the returned
arguments.
</p>
<h3>References</h3>
<p>See "Nonparametric Regression and Generalized Linear Models"
by Green and Silverman.
See "Additive Models" by Hastie and Tibshirani.
</p>
<h3>See Also</h3>
<p><code><a href="Krig.html">Krig</a></code>,
<code><a href="mKrig.html">mKrig</a></code>,
<code><a href="spatialProcess.html">spatialProcess</a></code>,
<code><a href="summary.Krig.html">summary.Krig</a></code>,
<code><a href="predict.Krig.html">predict.Krig</a></code>,
<code><a href="predictSE.Krig.html">predictSE.Krig</a></code>,
<code><a href="predictSurface.html">predictSurface</a></code>,
<code><a href="predictSurface.html">predictSurface.fastTps</a></code>,
<code><a href="plot.Krig.html">plot.Krig</a></code>,
<code><a href="surface.Krig.html">surface.Krig</a></code>,
<code><a href="sreg.html">sreg</a></code>,
<code><a href="mKrigMLE.html">fastTpsMLE</a></code>
</p>
<h3>Examples</h3>
<pre>
#2-d example
fit<- Tps(ChicagoO3$x, ChicagoO3$y) # fits a surface to ozone measurements.
set.panel(2,2)
plot(fit) # four diagnostic plots of fit and residuals.
set.panel()
# summary of fit and estiamtes of lambda the smoothing parameter
summary(fit)
surface( fit) # Quick image/contour plot of GCV surface.
# NOTE: the predict function is quite flexible:
look<- predict( fit, lambda=2.0)
# evaluates the estimate at lambda =2.0 _not_ the GCV estimate
# it does so very efficiently from the Krig fit object.
look<- predict( fit, df=7.5)
# evaluates the estimate at the lambda values such that
# the effective degrees of freedom is 7.5
# compare this to fitting a thin plate spline with
# lambda chosen so that there are 7.5 effective
# degrees of freedom in estimate
# Note that the GCV function is still computed and minimized
# but the lambda values used correpsonds to 7.5 df.
fit1<- Tps(ChicagoO3$x, ChicagoO3$y,df=7.5)
set.panel(2,2)
plot(fit1) # four diagnostic plots of fit and residuals.
# GCV function (lower left) has vertical line at 7.5 df.
set.panel()
# The basic matrix decompositions are the same for
# both fit and fit1 objects.
# predict( fit1) is the same as predict( fit, df=7.5)
# predict( fit1, lambda= fit$lambda) is the same as predict(fit)
# predict onto a grid that matches the ranges of the data.
out.p<-predictSurface( fit)
image( out.p)
# the surface function (e.g. surface( fit)) essentially combines
# the two steps above
# predict at different effective
# number of parameters
out.p<-predictSurface( fit,df=10)
## Not run:
# predicting on a grid along with a covariate
data( COmonthlyMet)
# predicting average daily minimum temps for spring in Colorado
# NOTE to create an 4km elevation grid:
# data(PRISMelevation); CO.elev1 <- crop.image(PRISMelevation, CO.loc )
# then use same grid for the predictions: CO.Grid1<- CO.elev1[c("x","y")]
obj<- Tps( CO.loc, CO.tmin.MAM.climate, Z= CO.elev)
out.p<-predictSurface( obj,
grid.list=CO.Grid, ZGrid= CO.elevGrid)
image.plot( out.p)
US(add=TRUE, col="grey")
contour( CO.elevGrid, add=TRUE, levels=c(2000), col="black")
## End(Not run)
## Not run:
#A 1-d example with confidence intervals
out<-Tps( rat.diet$t, rat.diet$trt) # lambda found by GCV
out
plot( out$x, out$y)
xgrid<- seq( min( out$x), max( out$x),,100)
fhat<- predict( out,xgrid)
lines( xgrid, fhat,)
SE<- predictSE( out, xgrid)
lines( xgrid,fhat + 1.96* SE, col="red", lty=2)
lines(xgrid, fhat - 1.96*SE, col="red", lty=2)
#
# compare to the ( much faster) B spline algorithm
# sreg(rat.diet$t, rat.diet$trt)
# Here is a 1-d example with 95 percent CIs where sreg would not
# work:
# sreg would give the right estimate here but not the right CI's
x<- seq( 0,1,,8)
y<- sin(3*x)
out<-Tps( x, y) # lambda found by GCV
plot( out$x, out$y)
xgrid<- seq( min( out$x), max( out$x),,100)
fhat<- predict( out,xgrid)
lines( xgrid, fhat, lwd=2)
SE<- predictSE( out, xgrid)
lines( xgrid,fhat + 1.96* SE, col="red", lty=2)
lines(xgrid, fhat - 1.96*SE, col="red", lty=2)
## End(Not run)
# More involved example adding a covariate to the fixed part of model
## Not run:
set.panel( 1,3)
# without elevation covariate
out0<-Tps( RMprecip$x,RMprecip$y)
surface( out0)
US( add=TRUE, col="grey")
# with elevation covariate
out<- Tps( RMprecip$x,RMprecip$y, Z=RMprecip$elev)
# NOTE: out$d[4] is the estimated elevation coefficient
# it is easy to get the smooth surface separate from the elevation.
out.p<-predictSurface( out, drop.Z=TRUE)
surface( out.p)
US( add=TRUE, col="grey")
# and if the estimate is of high resolution and you get by with
# a simple discretizing -- does not work in this case!
quilt.plot( out$x, out$fitted.values)
#
# the exact way to do this is evaluate the estimate
# on a grid where you also have elevations
# An elevation DEM from the PRISM climate data product (4km resolution)
data(RMelevation)
grid.list<- list( x=RMelevation$x, y= RMelevation$y)
fit.full<- predictSurface( out, grid.list, ZGrid= RMelevation)
# this is the linear fixed part of the second spatial model:
# lon,lat and elevation
fit.fixed<- predictSurface( out, grid.list, just.fixed=TRUE, ZGrid= RMelevation)
# This is the smooth part but also with the linear lon lat terms.
fit.smooth<-predictSurface( out, grid.list, drop.Z=TRUE)
#
set.panel( 3,1)
fit0<- predictSurface( out0, grid.list)
image.plot( fit0)
title(" first spatial model (w/o elevation)")
image.plot( fit.fixed)
title(" fixed part of second model (lon,lat,elev linear model)")
US( add=TRUE)
image.plot( fit.full)
title("full prediction second model")
set.panel()
## End(Not run)
###
### fast Tps
# m=2 p= 2m-d= 2
#
# Note: theta =3 degrees is a very generous taper range.
# Use some trial theta value with rdist.nearest to determine a
# a useful taper. Some empirical studies suggest that in the
# interpolation case in 2 d the taper should be large enough to
# about 20 non zero nearest neighbors for every location.
fastTps( RMprecip$x,RMprecip$y,m=2,lambda= 1e-2, theta=3.0) -> out2
# note that fastTps produces an mKrig object so one can use all the
# the overloaded functions that are defined for the mKrig class.
# summary of what happened note estimate of effective degrees of
# freedom
print( out2)
## Not run:
set.panel( 1,2)
surface( out2)
#
# now use great circle distance for this smooth
# Here "theta" for the taper support is the great circle distance in degrees latitude.
# Typically for data analysis it more convenient to think in degrees. A degree of
# latitude is about 68 miles (111 km).
#
fastTps( RMprecip$x,RMprecip$y,m=2,lambda= 1e-2,lon.lat=TRUE, theta= 3.0) -> out3
print( out3) # note the effective degrees of freedom is different.
surface(out3)
set.panel()
## End(Not run)
## Not run:
#
# simulation reusing Tps/Krig object
#
fit<- Tps( rat.diet$t, rat.diet$trt)
true<- fit$fitted.values
N<- length( fit$y)
temp<- matrix( NA, ncol=50, nrow=N)
sigma<- fit$shat.GCV
for ( k in 1:50){
ysim<- true + sigma* rnorm(N)
temp[,k]<- predict(fit, y= ysim)
}
matplot( fit$x, temp, type="l")
## End(Not run)
#
#4-d example
fit<- Tps(BD[,1:4],BD$lnya,scale.type="range")
# plots fitted surface and contours
# default is to hold 3rd and 4th fixed at median values
surface(fit)
</pre>
<hr /><div style="text-align: center;">[Package <em>fields</em> version 9.9 <a href="00Index.html">Index</a>]</div>
</body></html>