DATASETS NEEDED:

coalash.dt (Coal Ash dataset)
www4.stat.ncsu.edu/~fuentes/coalash.dat

davis.dat  (topographic heights)
www4.stat.ncsu.edu/~fuentes/davis.dat

ozone.dat (ozone values)
www4.stat.ncsu.edu/~fuentes/ozone.dat



software R with library geoR, fields and akima



PART II


Spatial prediction:


##EXERCISE 1: Estimating covariance parameters
##EXERCISE 2: Moving neighborhood  kriging
##EXERCISE 3: ANISTROPY 
##EXERCISE 4: Maps and Images
########################################################################






# download your data using read.table
# to learn more about read.table  type
# > ? read.table

coal.ash<-read.table('coalash.dat')
coal.m<-as.matrix(coal.ash)
davis.dat<-read.table('davis.dat')
davis.m<-as.matrix(davis.dat)


###################################

##EXERCISE 1: Estimating model parameters
###A. Bayesian estimates of model parameters:


#Function:   krige.bayes
# You obtain posteriors for the trend, range, partial  sill and nugget
# Recommended (default) priors: 
# uniform for beta, 1/sigma^2 for the partial  sill
# and discrete uniform for the range  and nugget.
#

# aniso.pars are the parameters for anisotropy: stretching and rotation angle.
# 
# the default is with nugget 0, but you 
# can also get a posterior for the nugget when you give a prior to the nugget,
# by saying nugget.prior="uniform"
# 


#Additional information:
#----------------------
# This function computes (Bayesian) estimates of a spatial linear 
#model and performs Bayesian and/or kriging prediction in a set of 
#locations specified by the user.
# Priors options for mean and/or covariance parameters 
# coords       : data coordinates (vector for 1d or matrix for 2d data)
#  data         : vector with data values
# locations    : coordinates of points to be estimated (vector for 1d or 
# matrix for 2d data). If locations='no' only model parameters estimates are 
#returned if the full bayesian model is considered.
# trend.d      : trend data in data locations. Default is constant trend. 
#The options '1st' or '2nd' builts a first or second degree polinomial.
# trend.l      : trend data in locations to be estimated. Default 
#is constant trend. The options '1st' or '2nd' builts a first or second 
#degree polinomial.

#
#Example:
#nugget is fixed to zero or anyother value (here 0):
#


coal.by1<-krige.bayes( coords=coal.m[,2:3],data=coal.m[,4],
prior=prior.control(phi.discrete=seq(0, 3, l=21)))


# see the histograms of the posteriors
# for the mean, the partial sill and the range:

X11()
par(mfrow=c(1,3))

hist(coal.by1$posterior$sample$beta)

hist(coal.by1$posterior$sample$sigmasq)

hist(coal.by1$posterior$sample$phi)


#
#

# to see your function type:  coal.by1$call

#posterior for MEAN
coal.by1$posterior$beta$summary

# mean   median mode.cond 
# 9.740055 9.740496  9.758056

#posterior for partial sill

coal.by1$posterior$sigmasq

# mean   median mode.cond 
# 9.740055 9.740496  9.758056


#posterior for range
coal.by1$posterior$phi

#     mean    median      mode 
# 0.8645331 0.7894737 0.7894737



###B. Bayesian prediction:

#locations where we want Bayesian prediction:

loci<-matrix(c(2,3.5,4,5.5),
ncol=2,byrow=T)

 loci
     [,1] [,2] 
[1,]    2  3.5
[2,]    4  5.5

# this is the same as before for krige.bayes, but now we
# specify where you want to predict, by using the locations argument:

davis.bpred<-
krige.bayes( coords=davis.m[,1:2],data=davis.m[,3],
locations=loci,prior=prior.control(phi.discrete=seq(0, 3, l=21)))


davis.bpred$predictive$simulations

davis.bpred$predictive$mean


davis.bpred$predictive$variance

# 803.0342 735.7671
# 369.3989 461.2960


#Bayesian predictive distributions:
par(mfrow=c(2,1)

hist(davis.bpred$predictive$simulations[1,])

hist(davis.bpred$predictive$simulations[2,])


#
#
# If you want a nugget prior, instead of nugget=0,
#say  tausq.rel.discrete= c(1,2,5,6)
# or any other vector with the values of the discrete uniform prior
# of nugget/partial sill. Example:

coal.b2<-krige.bayes(coords=coal.m[,2:3],data=coal.m[,4],
prior=prior.control(tausq.rel.prior="uniform",
             tausq.rel.discrete=seq(0,0.5,l=6),
             phi.discrete=seq(0,3,l=25)))


X11()

par(mfrow=c(1,4))

hist(coal.b2$posterior$sample$beta)

hist(coal.b2$posterior$sample$sigmasq)

hist(coal.b2$posterior$sample$phi)
hist(coal.b2$posterior$sample$tausq.rel)



############################################################
##EXERCISE 2: Moving neighbourhood  kriging

###"Kriging performed in global neighbourhood":
#      cov.pars : covariance parameters vector (partial sill,range)
#      m0 : defines the type of kriging:
#                                 'sk': simple kriging (no trend)
#                                'ok': ordinary kriging (constant trend)
#                                'kt': kriging with a trend model(universal)
#      kappa : kappa (smoothing) is the smoothing
# parameter for Matern or powered exponential covariance function

coal.k<-ksline(coords=coal.m[,2:3],data=coal.m[,4],locations=loci,
cov.pars=c(3,1),nugget=0,
cov.model="matern",
kappa=.5,m0="ok")

##coal.k
  coal.k$predict

  [1] 10.07305 11.22320

  coal.k$krige.var

  [1] 2.684456 1.349229

  coal.k$beta:
         [,1] 
  [1,] 9.752618

   coal.k$message:
  [1] "Kriging performed in global neighborhood"


##Universal Kriging:
#
# the value of trend here is 2 which means we have a polynomial
# of degree 2 for a two dimensional problem, therefore
# we need to estimate 5 parameters

coal.uk<-ksline(coords=coal.m[,2:3],data=coal.m[,4],locations=loci,
cov.pars=c(3,1),nugget=0,
cov.model="matern",
kappa=.5,m0="kt",trend=2)

coal.uk$predict
[1] 10.77988 11.24153

coal.uk$krige.var
[1] 3.023352 1.349461



coal.uk$beta
      coefficients:
[1,]  1.124853e+01   (1)
[2,] -2.090557e-01   (x)
[3,] -1.381999e-02   (y)
[4,] -2.817476e-03   (x^2)
[5,] -9.021021e-04   (y^2)
[6,]  6.312104e-03   (x*y)


### kriging performed in moving neighborhood:

coal.wind <-ksline(coords=coal.m[,2:3],data=coal.m[,4],locations=loci,
cov.pars=c(3,1),nugget=0,
cov.model="matern",
kappa=.5,m0="ok",nwin=5)
#nwin number of closest neighbors

coal.wind 
  $predict:
  [1] 10.43828 10.84857

   $krige.var:
   [1] 3.136432 1.357889
  

###########################################3
# linear trend not subregions:
# when trend is one implies a linear trend in two dimensions
# therefore we have 3 parameters to estimate

coal.wind1 <-ksline(coords=coal.m[,2:3],data=coal.m[,4],locations=loci,
cov.pars=c(3,1),nugget=0,
cov.model="matern",
kappa=.5,m0="kt",trend=1,nwin="full")

> coal.wind1
coal.wind1$predict
[1] 10.62524 11.23732

coal.wind1$krige.var
[1] 2.774784 1.349341

coal.wind1$beta
             [,1] 
[1,]  10.86216801
[2,]  -0.16286289
[3,]   0.01077067


### kriging performed in moving neighborhood with linear trend

coal.wind2 <-ksline(coords=coal.m[,2:3],data=coal.m[,4],locations=loci,
cov.pars=c(3,1),nugget=0,
cov.model="matern",
kappa=.5,m0="kt",trend=1,nwin=5)

coal.wind2
coal.wind2$predict
[1]  9.669875 10.876154

coal.wind2$krige.var
[1] 8.350462 1.362241



t(coal.wind2$beta[1:2,])
         [,1]      [,2]      [,3] 
[1,] 7.467317 0.4058152 0.3430783
[2,] 7.030745 0.4716316 0.3130314





######################
## EXERCISE 3: ANISOTROPY.

#The function coords.aniso transform your coordinates
# into the coordinates in the new space where we have isotropy
# first one is  the rotation angle  (in radians) and the second is lambda
# where lambda is
# the stretching parameter  
# and it is greater than 1, if lambda is 1 there is not
# stretching
#
# 


new.coord<-coords.aniso(coal.m[,2:3],aniso.pars=c(.2,3))

# do kriging for anisotropic case
 coal.wind2 <-ksline(coords=coal.m[,2:3],data=coal.m[,4],locations=loci,
cov.pars=c(3,1),nugget=0,
cov.model="exponential",
m0="kt",trend=1,aniso.pars=c(.2,3))


#################################
##EXERCISE 4: Maps and Images
#use library akima and fields

X11()
 ozone.dat<-read.table('ozone.dat') 
       rx <- range(ozone.dat[,1])
       ry <- range(ozone.dat[,2])

      
#to obtain map of us within the limits rx (long) and ry (lat):
      
 US(xlim = rx, ylim = ry, lwd = 2, col = 1,add=F)

  

    
       surf<- interp(ozone.dat[,1], ozone.dat[,2], ozone.dat[,3])
#to create an image:
       image(surf$x,surf$y,surf$z, add = F,graphics.reset=T,col=topo.colors(12))
#to add points in the image:
       text(ozone.dat[,1:2], labels=ozone.dat[,3])
       
US(xlim = rx, ylim = ry, lwd = 2, col = 1,add=T)

#
#

#to create a 3-d plot:
        persp(surf$x,surf$y,surf$z)

# to get a contour:
        contour(surf$x,surf$y,surf$z)


#to add a legend to the image:
#in the fields library

 image.plot(surf$x,surf$y,surf$z, add = F,graphics.reset=T,col=topo.colors(12))

#
#DESCRIPTION of image:
#       Creates an image, under some graphics devices,  of  shades
#       of gray or colors that represent a third dimension.
#
#USAGE:
#       image(x, y, z, zlim = range(z), add = F)
#
#REQUIRED ARGUMENTS:
#x:     vector containing x coordinates of grid over which  z'  is
#       evaluated.   The  values  should  be  in increasing order;
#       missing values are not accepted.
#y:     vector of grid y coordinates.  The  values  should  be  in
#       increasing order; missing values are not accepted.
#z:    is a matrix, z[i,j]'  is
#       evaluated at x[i]', y[j]').  The rows of z' are indexed by
#       x',  and  the  columns  by  y'.  Missing values (NA's) are
#       allowed.
#####
#
# interp: Interpolates the value  of  the  third  variable  onto  an
# evenly spaced grid of the first two variables. default is 40x40
# grid


#to save the plot as a postscript file
postscript(file="example.ps")
 image.plot(surf$x,surf$y,surf$z, add = F,graphics.reset=T,col=topo.colors(12))
dev.off()