#2-d example
fit<- Tps(ozone$x, ozone$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(ozone$x, ozone$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<-predict.surface( 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<-predict.surface( fit,df=10)
#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<- predict.se( 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<- predict.se( out, xgrid)
lines( xgrid,fhat + 1.96* SE, col="red", lty=2)
lines(xgrid, fhat - 1.96*SE, col="red", lty=2)
# Adding a covariate to the fixed part of model
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 coeficient
# it is easy to get the smooth surface separate from the elevation.
out.p<-predict.surface( 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)
#
# however the correct way to do this is evaluate the estimate
# on a grid where you also have elevations
# To evaluate with elevations on a grid takes some more
# work because we need to use the same grid where elevations are available.
#
# An elevation DEM from the PRISM climate data product (4km resolution)
data(RMelevation)
grid.list<- list( x=RMelevation$x, y= RMelevation$y)
# now create single big 2 column matrix with all the lon.lat locations in this grid
bigX<- make.surface.grid( grid.list)
# and create a big vector of the elevations in the grid
bigZ<- c( RMelevation$z)
#NOTE: nrow( bigX) = 69938 = 289X242, length( bigZ)= 69938
fit<- predict(out, x= bigX, Z= bigZ)
# fit is just a big vector
# now coerce back into grid format (x,y,z components)
# this function works because bigX has some attributes set that indicate the parent grid.
fit<- as.surface( bigX, fit)
#
# decomposing into the different pieces:
# this is just the spatial model without elevation
fit0<- predict(out0, x=bigX)
fit0<- as.surface( bigX, fit0)
# this is the linear fixed part of the second spatial model: lon,lat and elevation
fit.fixed<- predict( out, x= bigX, just.fixed=TRUE, Z= bigZ)
fit.fixed<- as.surface( bigX, fit.fixed)
# This is the smooth part but also with the linear lon lat terms.
fit.smooth<- predict( out, x= bigX, drop.Z=TRUE)
fit.smooth<- as.surface( bigX, fit.smooth)
#
set.panel( 2,2)
image.plot( fit0)
title(" first spatial model (w/o elevation)")
image.plot( fit.fixed)
title(" fixed part of second model")
image.plot( fit.smooth)
title("smoothed component second model")
US( add=TRUE)
image.plot( fit)
title("full prediction second model")
set.panel()
###
### 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)
set.panel( 1,2)
surface( out2)
#
# now use great circle distance for this smooth
# note the different "theta" for the taper support ( there are
# about 70 miles in one degree of latitude).
#
fastTps( RMprecip$x,RMprecip$y,m=2,lambda= 1e-2,lon.lat=TRUE, theta=210) -> out3
print( out3) # note the effective degrees of freedom is different.
surface(out3)
set.panel()
#
# 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")
#
#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)Run the code above in your browser using DataLab