%\VignetteIndexEntry{An introduction to the package geoR}
%\VignetteDepends{}
%\VignetteKeywords{spatial}
%\VignettePackage{geoR}
\documentclass[12pt,a4paper]{article}
\usepackage[latin1]{inputenc}
%\usepackage[brazil,brazilian]{babel}
\usepackage[dvips]{graphicx}
%\usepackage{psfig,epsfig,psfrag}
\usepackage{color,enumerate,longtable}
\usepackage{texnames,lastpage,verbatim,fancyvrb}
%% possible conflicts with Rd.sty
%\usepackage{url,html}
%\usepackage{heqn}
%\usepackage{amssymb,amsmath,amsfonts}
%
% Page size
%
%\usepackage{setspace,calc}
%\setlength{\voffset}{-1.5in}
%\setlength{\hoffset}{-1in}
%\setlength{\topmargin}{30mm}
%\setlength{\oddsidemargin}{40mm}
%\setlength{\textwidth}{210mm-10mm-\oddsidemargin}
%\setlength{\textheight}{297mm-10mm-\topmargin-\headheight-\headsep}
%
% Sweave settings
%
% Setting default size for \includegraphics 
% setting global defaults for figures: directory, names and sizes  
%%%%% \ SweaveOpts{prefix.string=figures/fig,width=4,height=4}
\SweaveOpts{width=4,height=4}
% emulating Sweave: Sweave style for R input code and output
\DefineVerbatimEnvironment{Rin}{Verbatim}{fontshape=sl}
\DefineVerbatimEnvironment{Rout}{Verbatim}{}
\newcommand{\Rcmd}[1]{\textsl{\texttt{#1}}}

% loading Defs from Rd.sty
\usepackage{/usr/local/lib/R/share/texmf/Rd}
% \newcommand{\strong}[1]{{\normalfont\fontseries{b}\selectfont #1}}
% \let\pkg=\strong
% \newcommand{\code}[1]{\texttt{\small #1}}
% \newcommand{\R}{{\normalfont\textsf{R}}{}}
% \newcommand{\sQuote}[1]{`#1'}
% \newcommand{\dQuote}[1]{``#1''}

\title{\pkg{geoR} : Package for Geostatistical Data Analysis\\An illustrative session}
\author{Paulo J. Ribeiro Jr. \& Peter J. Diggle}
\date{Last update: \today}

\begin{document}

\maketitle

%% Sweave settings for includegraphics default plot size (Sweave dafault is 0.8)  
%% notice this must be after begin{document}
%%% \setkeys{Gin}{width=0.9\textwidth}

<<label="R settings", echo=FALSE>>= 
options("width"=70)
options(SweaveHooks=list(fig=function() par(mar=c(3,3,1,0.5), mgp=c(2,1,0))))
options(geoR.messages = FALSE)
@

\section{Introduction} 
The package \pkg{geoR}  provides functions for geostatistical data analysis using the  software \R.
This document illustrates some (but not all!)
of the capabilities of the package.

The objective is to familiarise the reader with the 
\pkg{geoR}'s commands for data analysis 
and show some of the graphical outputs which can be produced.

The commands used here are just illustrative, providing basic examples of the package handling.
We did not attempt to perform a definitive analysis of the data-set used throughout the exemples neither to cover all the details of the package capability.
      
In what follows:
\begin{itemize}
\item the \R{} commands are shown in \Rcmd{slanted typewriter fonts like this}, 
\item the corresponding output, if any, is shown in \texttt{typewriter fonts like this}.
\end{itemize}

Typically, default arguments are used for the function calls and the user 
is encouraged to inspect other arguments of the functions using the \command{args} and \command{help} functions.
For instance, to see all the arguments for the function \code{variog} type
\code{args(variog)} and/or \code{help(variog)}.

% this must be changed:
The commands shown in this page are also 
available in the file \url{http://www.est.ufpr.br/geoR/geoRdoc/geoRintro.R}.

We refer to the \pkg{geoR} documentation (\url{http://www.est.ufpr.br/geoR/geoR.html#help})
for more details on the functions included in the package \pkg{geoR}.

        
\section{Starting a Session and Loading Data}

\subsection{Loading the \pkg{package}}        
After starting an \R{} session, load \pkg{geoR} with the command \code{library} (or \code{require}).
If the package is loaded correctly a message will be displayed.
<<label="Loading package">>=   
library(geoR)
@

If the installation directory for the package is the default location for \R{} packages, type:
\begin{Rin}
> library(geoR, lib.loc="PATH_TO_geoR")
\end{Rin}
%<<eval=FALSE>>=
%> library(geoR, lib.loc="PATH_TO_geoR")
%@ 
where \option{"PATH\_TO\_geoR"} is the path to the directory where \pkg{geoR} was installed. 

\subsection{Using data}        
Typically, data are stored as an object (a list) of the class \code{geodata}. 
An object of this class contains two obligatory elements: the coordinates of data locations 
as first element (\code{\$coords}) and the data values as second element (\code{\$data}), 
which can be a vector or a matrix.  
Objets of the class \code{geodata} may have other elements such as covariates and coordinates of the borders of the study area.

We refer to the documentation for the functions \command{as.geodata} and \command{read.geodata} for more information on how to import/convert data and on the definitions for the class  \code{"geodata"}.
% add this as an appendix?
Check \url{http://www.est.ufpr.br/geoR/importingASCII.html} for information on  how to read data from an ASCII (text) file.

There a a few data-sets included in the package distribution.
For the examples included in this document we use the data set \code{s100} included in the \pkg{geoR} distribution.
To load this data type:
<<label="sourcing data">>=
data(s100)
@ 
The list of all data-sets included in the package is given by   
\code{data(package='geoR')}.


\section{Exploratory Tools}
        
A quick summary for the \code{geodata} object
can be obtained using a method for \command{summary}
which will return a information on the coordinates and data values.
<<label="summary of geodata">>=
summary(s100)
@ 

The elements of the list \code{\$covariate}, \code{\$borders} and/or \code{units.m} will be 
also summarized if present in the \code{geodata} object.

\subsection{Plotting data locations and values}

The function \command{plot.geodata} shows a $2\times2$ display with data locations (top plots) and data {\it versus} coordinates (bottom plots).
For an object of the class \code{geodata} the command \command{plot(s100)} produce the output shown in Figure~\ref{fig:plot}.
\begin{figure}[h!]
\centering
<<echo=F,fig=T,width=6,height=6>>=
plot(s100)
@
\label{fig:plot} 
\caption{Plot produced by \texttt{plot.geodata}.}
\end{figure}

%Notice that the top-right plot is produced using the package \scatterplot3d}.
%If this package is not installed a histogram of the data will replace this plot.

The function \code{points.geodata} produces a plot showing the data locations. Alternatively, points indicating the data locations can be added to a current plot.
There are options to specify point sizes, patterns and colors, which can be set to be proportional to the data values or specified quantiles.
Some examples of graphical outputs are illustrated by the commands below and corresponding plots as shown in Figure~\ref{fig:points}.
%\footnote{Notice that we start by saving the current graphical parameters in the object \code{par.ori} to restore them at the end}.
<<label="points.geodata",eval=F>>=
par(mfrow=c(2,2))
points(s100,xlab="Coord X",ylab="Coord Y")        
points(s100,xlab="Coord X",ylab="Coord Y", pt.divide="rank.prop")
points(s100,xlab="Coord X",ylab="Coord Y", cex.max=1.7, col=gray(seq(1, 0.1, l=100)), pt.divide="equal")
points(s100,pt.divide="quintile",xlab="Coord X",ylab="Coord Y")
@ 
\begin{figure}[h!]
\centering
<<echo=F,fig=T,width=6,height=6>>=
par(mfrow=c(2,2))
points(s100,xlab="Coord X",ylab="Coord Y")        
points(s100,xlab="Coord X",ylab="Coord Y", pt.divide="rank.prop")
points(s100,xlab="Coord X",ylab="Coord Y", cex.max=1.7, col=gray(seq(1, 0.1, l=100)), pt.divide="equal")
points(s100,pt.divide="quintile",xlab="Coord X",ylab="Coord Y")
@ 
\label{fig:points}
\caption{Plot produced by \texttt{points.geodata}.}
\end{figure}

\subsection{Empirical variograms}

Empirical variograms are calculated using the function \code{variog}. There are options for the {\it classical} or {\it modulus} estimator. 
Results can be returned as variogram clouds, binned or smoothed variograms.
There are methods for \command{plot} to facilitate the display of the results as shown in Figure~\ref{fig:variog}. 
<<label="variograms">>=
cloud1 <- variog(s100, option = "cloud", max.dist=1)
cloud2 <- variog(s100, option = "cloud", estimator.type = "modulus", max.dist=1)
bin1 <- variog(s100, uvec=seq(0,1,l=11))
bin2  <- variog(s100, uvec=seq(0,1,l=11), estimator.type= "modulus")
par(mfrow=c(2,2))
plot(cloud1, main = "classical estimator")
plot(cloud2, main = "modulus estimator")
plot(bin1, main = "classical estimator")
plot(bin2, main = "modulus estimator")
@ 
\begin{figure}[h!]
\centering
<<echo=F,fig=T,width=6,height=6>>=
par(mfrow=c(2,2))
plot(cloud1, main = "classical estimator")
plot(cloud2, main = "modulus estimator")
plot(bin1, main = "classical estimator")
plot(bin2, main = "modulus estimator")
@ 
\caption{Ploting of \texttt{variog} results.}
\label{fig:variog}
\end{figure}

Several results are returned by the function \command{variog}. 
The first three are the more important ones
and contains the distances, the estimated semivariance and the number of pairs for each bin.
<<>>=
names(bin1)
@ 

Furthermore, the points of the variogram clouds can be grouped into classes of distances ("bins") and displayed with a box-plot for each bin as shown in Figure~\ref{fig:boxplot}. This can be used as an exploratory tool to access variogram results.
<<>>=
bin1 <- variog(s100,uvec = seq(0,1,l=11), bin.cloud = T)
bin2 <- variog(s100,uvec = seq(0,1,l=11), estimator.type = "modulus", bin.cloud = T)
par(mfrow=c(1,2))
plot(bin1, bin.cloud = T, main = "classical estimator")
plot(bin2, bin.cloud = T, main = "modulus estimator")
@ 
\setkeys{Gin}{width=0.9\textwidth}
\begin{figure}[h!]
\centering
<<echo=F,fig=T,width=8,height=4>>=
par(mfrow=c(1,2))
plot(bin1, bin.cloud = T, main = "classical estimator")
plot(bin2, bin.cloud = T, main = "modulus estimator")
@ 
\label{fig:boxplot}
\caption{Box-plots for each of the variogram bins.}
\end{figure}

Directional variograms can also be computed by the function
\code{variog} using the arguments
\option{direction} and \option{tolerance}. For example, to compute a
variogram for the direction 60 degrees with the default tolerance
angle (22.5 degrees) the command would be:
<<>>=
vario60 <- variog(s100, max.dist = 1, direction=pi/3) 
@ 

For a quick computation in four directions we can use the function \code{variog4} which by default computes variogram for the direction angles $0^o$, $45^o$, $90^o$ and $135^o$.
<<>>=
vario.4 <- variog4(s100, max.dist = 1)
@ 

The Figure~\ref{fig:directional} show the directional variograms obtained with the commands above.
\setkeys{Gin}{width=0.9\textwidth}
\begin{figure}[h!]
\centering
<<fig=T,width=8,height=4>>=
par(mfrow=c(1,2), mar=c(3,3,1.5,.5))
plot(vario60)
title(main = expression(paste("directional, angle = ", 60 * degree)))       
plot(vario.4, lwd=2)
@ 
\caption{Directional variograms}
\label{fig:directional}
\end{figure}


\section{Parameter Estimation}

Theoretical and empirical variograms can be plotted and visually compared. For example, the left panel in Figure~\ref{fig:variomodel} shows the theoretical variogram model used to simulate the data \code{s100} and two estimated variograms.
<<eval=F>>=
bin1 <- variog(s100, uvec = seq(0,1,l=11))
plot(bin1)
lines.variomodel(cov.model = "exp", cov.pars = c(1,0.3), nugget = 0, max.dist = 1,  lwd = 3)
smooth <- variog(s100, option = "smooth", max.dist = 1, n.points = 100, kernel = "normal", band = 0.2)
lines(smooth, type ="l", lty = 2)
legend(0.4, 0.3, c("empirical", "exponential model", "smoothed"), lty = c(1,1,2), lwd = c(1,3,1))
@ 

In practice we usually don't know the true parameters which have top be estimated by some method. 
In the package \pkg{geoR} the model parameters can be estimated:
\begin{enumerate}
\item{\textit{\dQuote{by eye}}: trying different models over empirical variograms (using the function \code{lines.variomodel}),}
\item{\textit{by least squares fit of empirical variograms}: with options for ordinary (OLS) and weighted (WLS) least squares (using the function \code{variofit}),}
\item{\textit{by likelihood based methods}: with options for maximum likelihood (ML) and restricted maximum likelihood (REML) (using the function \code{likfit}),}
\item{\textit{Bayesian methods}: are also implemented and will be presented in Section 5 (using the function \code{krige.bayes}).}
\end{enumerate}

Fitting \dQuote{by eye} consists of drawing curves of theoretical variogram functions over an empirical variogram, changing the variogram model and/or its parameters and, at last, choosing one of them.
The following commands show how to add a line with a variogram model to a variogram plot.
Three different variogram models are used. 
<<eval=F>>=
plot(variog(s100, max.dist=1))
lines.variomodel(cov.model="exp", cov.pars=c(1,.3), nug=0, max.dist=1)
lines.variomodel(cov.model="mat", cov.pars=c(.85,.2), nug=0.1, kappa=1,max.dist=1, lty=2)
lines.variomodel(cov.model="sph", cov.pars=c(.8,.8), nug=0.1,max.dist=1, lwd=2)
@ 

\setkeys{Gin}{width=0.9\textwidth}
\begin{figure}[h!]
\centering
<<echo=F,fig=T,width=8,height=4>>=
par(mfrow=c(1,2))
plot(bin1)
lines.variomodel(cov.model = "exp", cov.pars = c(1,0.3), nugget = 0, max.dist = 1,  lwd = 3)
smooth <- variog(s100, option = "smooth", max.dist = 1, n.points = 100, kernel = "normal", band = 0.2)
lines(smooth, type ="l", lty = 2)
legend(0.4, 0.3, c("empirical", "exponential model", "smoothed"), lty = c(1,1,2), lwd = c(1,3,1))

plot(variog(s100, max.dist=1))
lines.variomodel(cov.model="exp", cov.pars=c(1,.3), nug=0, max.dist=1)
lines.variomodel(cov.model="mat", cov.pars=c(.85,.2), nug=0.1, kappa=1,max.dist=1, lty=2)
lines.variomodel(cov.model="sph", cov.pars=c(.8,.8), nug=0.1,max.dist=1, lwd=2)
@ 
\caption{Theoretical variogram curves added to empirical variograms.}
\label{fig:variomodel}
\end{figure}

When using the parameter estimation functions \code{variofit} and \code{likfit} 
the nugget effect parameter can either be estimated or set to a fixed value. 
The same applies for smoothness, anisotropy and transformation parameters. 
Options for taking trends into account are also included. 
Trends can be specified as polynomial functions of the coordinates and/or 
linear functions of given covariates.

An example call to \code{likfit} is given below.
Methods for \code{print()} and \code{summary()} have been written to summarize the resulting objects.
<<echo=F>>=
options(geoR.messages = TRUE)
@ 
<<label="Maximum Likelihood">>=
ml <- likfit(s100, ini = c(1,0.5))
ml
summary(ml)
@ 

The commands below shows how to fit models by using different methods, with options for fixed or estimated nugget parameter.  Notice there are other features not illustrated here such as estimation of trends, anisotropy, smoothness and Box-Cox transformation parameter.
Notice in the call above that the functions show some messages while they are running --- and we don't want to see them in the following calls.
To prevent this we can set the argument \code{messages = FALSE} at each function call or,
to set it globally for all functions use \code{options()} as follows.
<<>>=
options(geoR.messages = FALSE)
@ 

\begin{itemize}
\item Fitting models with nugget fixed to zero
<<>>=
ml <- likfit(s100, ini = c(1,0.5), fix.nugget = T)
reml <- likfit(s100, ini = c(1,0.5), fix.nugget = T, method = "RML")
ols <- variofit(bin1, ini = c(1,0.5), fix.nugget = T, weights="equal")
wls <- variofit(bin1, ini = c(1,0.5), fix.nugget = T)
@ 

\item Fitting models with a fixed value for the nugget
<<>>=
ml.fn <- likfit(s100, ini = c(1,0.5), fix.nugget = T, nugget = 0.15)
reml.fn <- likfit(s100, ini = c(1,0.5), fix.nugget = T, nugget = 0.15, method = "RML")
ols.fn <- variofit(bin1,ini = c(1,0.5), fix.nugget = T, nugget = 0.15, weights="equal")
wls.fn <- variofit(bin1, ini = c(1,0.5), fix.nugget = T, nugget = 0.15)
@ 

\item Fitting models estimated nugget
<<>>=
ml.n <- likfit(s100, ini = c(1,0.5), nug = 0.5)
reml.n <- likfit(s100, ini = c(1,0.5), nug = 0.5, method = "RML")
ols.n <- variofit(bin1, ini = c(1,0.5), nugget=0.5, weights="equal")
wls.n <- variofit(bin1, ini = c(1,0.5), nugget=0.5)
@ 
\end{itemize}

\setkeys{Gin}{width=\textwidth}

Now, the comands for plotting fitted models against empirical variogram as show in Figure~\ref{fig:fit} are:
<<eval=F>>=
par(mfrow = c(1,3))
plot(bin1, main = expression(paste("fixed ", tau^2 == 0)))
lines(ml, max.dist = 1)
lines(reml, lwd = 2, max.dist = 1)
lines(ols, lty = 2, max.dist = 1)
lines(wls, lty = 2, lwd = 2, max.dist = 1)
legend(0.5, 0.3, legend=c("ML","REML","OLS","WLS"),lty=c(1,1,2,2),lwd=c(1,2,1,2), cex=0.7)

plot(bin1, main = expression(paste("fixed ", tau^2 == 0.15)))
lines(ml.fn, max.dist = 1)
lines(reml.fn, lwd = 2, max.dist = 1)
lines(ols.fn, lty = 2, max.dist = 1)
lines(wls.fn, lty = 2, lwd = 2, max.dist = 1)
legend(0.5, 0.3, legend=c("ML","REML","OLS","WLS"), lty=c(1,1,2,2), lwd=c(1,2,1,2), cex=0.7)

plot(bin1, main = expression(paste("estimated  ", tau^2)))
lines(ml.n, max.dist = 1)
lines(reml.n, lwd = 2, max.dist = 1)
lines(ols.n, lty = 2, max.dist = 1)
lines(wls.n, lty =2, lwd = 2, max.dist = 1)
legend(0.5, 0.3, legend=c("ML","REML","OLS","WLS"), lty=c(1,1,2,2), lwd=c(1,2,1,2), cex=0.7)
@ 
\setkeys{Gin}{width=\textwidth}
\begin{figure}
\centering
<<echo=F,fig=T,width=6,height=2>>=
par(mfrow = c(1,3))
plot(bin1, main = expression(paste("fixed ", tau^2 == 0)))
lines(ml, max.dist = 1)
lines(reml, lwd = 2, max.dist = 1)
lines(ols, lty = 2, max.dist = 1)
lines(wls, lty = 2, lwd = 2, max.dist = 1)
legend(0.5, 0.3, legend=c("ML","REML","OLS","WLS"),lty=c(1,1,2,2),lwd=c(1,2,1,2), cex=0.7)

plot(bin1, main = expression(paste("fixed ", tau^2 == 0.15)))
lines(ml.fn, max.dist = 1)
lines(reml.fn, lwd = 2, max.dist = 1)
lines(ols.fn, lty = 2, max.dist = 1)
lines(wls.fn, lty = 2, lwd = 2, max.dist = 1)
legend(0.5, 0.3, legend=c("ML","REML","OLS","WLS"), lty=c(1,1,2,2), lwd=c(1,2,1,2), cex=0.7)

plot(bin1, main = expression(paste("estimated  ", tau^2)))
lines(ml.n, max.dist = 1)
lines(reml.n, lwd = 2, max.dist = 1)
lines(ols.n, lty = 2, max.dist = 1)
lines(wls.n, lty =2, lwd = 2, max.dist = 1)
legend(0.5, 0.3, legend=c("ML","REML","OLS","WLS"), lty=c(1,1,2,2), lwd=c(1,2,1,2), cex=0.7)
@ 
\label{fig:fit}
\caption{Empirical variogram and models fitted by different methods.}
\end{figure}

Two kinds of variogram {\it envelopes}
computed by simulation are 
illustrated in the figure below. 

The plot on the left-hand side shows an envelope based on 
permutations  of the data values across the locations, 
i.e. envelopes built under 
the assumption of no spatial correlation.           
The envelopes shown on the right-hand side are based on 
simulations from a given set of model parameters, in this example the parameter 
estimates from the WLS variogram fit. 
This envelope shows the variability of the empirical variogram. 
<<label="Variogram Envelops">>=
env.mc <- variog.mc.env(s100, obj.var=bin1)
env.model <- variog.model.env(s100, obj.var=bin1, model=wls)
@

\setkeys{Gin}{width=0.9\textwidth}
\begin{figure}[h!]
\centering
<<fig=T,width=8,height=4>>=
par(mfrow=c(1,2))
plot(bin1, envelope=env.mc)
plot(bin1, envelope=env.model)
@ 
\label{fig:envelop}
\caption{Variogram envelopes.}
\end{figure}

Profile likelihoods (1-D and 2-D) are computed by the function \code{proflik}.
Here we show the profile likelihoods for the covariance parameters of the model without nugget effect previously fitted by \code{likfit}. \\
\strong{WARNING: RUNNING THE NEXT COMMAND CAN BE TIME-CONSUMING}
%\begin{Rin}
\setkeys{Gin}{width=\textwidth}
\begin{figure}
\centering
<<fig=T,width=6,height=2>>=
prof <- proflik(ml, geodata = s100, sill.val = seq(0.48, 2, l=11),
range.val = seq(0.1, 0.52, l=11), uni.only = FALSE)
par(mfrow=c(1,3))
plot(prof, nlevels=16)
@ 
\label{fig:profile}
\caption{Profile likelihoods for the model parameters}
\end{figure}
%\end{Rin}

\section{Cross-Validation}

The function \code{xvalid} performs 
cross-validation either using the {\it leaving-one-out} strategy 
or using a different set of locations provided
by the user through the argument \option{location.xvalid}.

For the first strategy, data points are removed one by one and
predicted by kriging using the remaining data.
The commands below illustrates cross-validation for the models
fitted by maximum likelihood and weighted least squares.
In the following two calls the model parameters remains the
same for the prediction at each location. 
<<label="Cross-validation">>=
xv.ml <- xvalid(s100, model=ml)
xv.wls <- xvalid(s100, model=wls)
@ 

Graphical results are shown for the cross-validation results where the
leaving-one-out strategy combined with the wls estimates for the parameters was used.
Cross-validation residuals are obtained subtracting the observed data minus the predicted value.
Standardised residuals are obtained dividing by the square root of the prediction variance (\sQuote{kriging variance}).
By default the 10 plots shown in the Figure~\ref{fig:xvalid} are produced but the user can restrict the choice using the function arguments.
<<eval=F>>=
par(mfcol=c(5,2),mar=c(3,3,1,.5), mgp=c(1.5,.7,0))
plot(xv.wls)
@ 
\setkeys{Gin}{width=0.7\textwidth}
\begin{figure}[t]
\centering
<<echo=F,fig=T,width=5,height=6>>=
par(mfcol=c(5,2),mar=c(3,3,1,.5), mgp=c(1.5,.7,0))
plot(xv.wls)
@ 
\caption{Cross-validations results}
\label{fig:xvalid}
\end{figure}

A variation of this method is illustrated by
the next two calls where the model parameters are re-estimated each
time a point is removed from the data-set.\\
\strong{WARNING: RUNNING THE NEXT COMMAND CAN BE TIME-CONSUMING}
<<eval=F>>=
xvR.ml <- xvalid(s100, model=ml, reest=TRUE)
xvR.wls <- xvalid(s100, model=wls, reest=TRUE, variog.obj=bin1)
@ 

\section{Spatial Interpolation}
        
Conventional geostatistical spatial interpolation ({\it kriging}) can be performed with options for:
\begin{enumerate}
\item{\textit{Simple kriging}}
\item{\textit{Ordinary kriging}}
\item{\textit{Trend (universal) kriging}}
\item{\textit{External trend kriging}}
\end{enumerate}

There are additional options for Box-Cox transformation (and back transformation of the results) and anisotropic models.
Simulations can be drawn from the resulting predictive 
distributions if requested.  

As a first example consider the prediction at four locations labeled {\it 1, 2, 3, 4} and indicated in the figure below.
\setkeys{Gin}{width=0.5\textwidth}
\begin{figure}[h!]
\centering
<<fig=T>>=
plot(s100$coords, xlim=c(0,1.2), ylim=c(0,1.2), xlab="Coord X", ylab="Coord Y")
loci <- matrix(c(0.2, 0.6, 0.2, 1.1, 0.2, 0.3, 1.0, 1.1), ncol=2)
text(loci, as.character(1:4), col="red")
polygon(x=c(0,1,1,0), y=c(0,0,1,1), lty=2)
@ 
\label{fig:loc}
\caption{Data locations and points to be predicted}
\end{figure}

The command to perform {\it ordinary kriging} using the parameters estimated by weighted least squares with nugget fixed to zero would be:
<<label="Kriging">>=
kc4 <- krige.conv(s100, locations = loci, krige = krige.control(obj.m = wls))
@ 

The output is a list including the predicted values  (\code{kc4\$predict}) and the kriging variances (\code{kc4\$krige.var}).

Consider now a second example. The goal is to perform prediction
on a grid covering the area and to display the results. 
Again, we use ordinary kriging. 
The commands commands below defines a grid of locations and performs the prediction at those locations.
<<>>=
pred.grid <-  expand.grid(seq(0,1, l=51), seq(0,1, l=51))
kc <- krige.conv(s100, loc = pred.grid, krige = krige.control(obj.m = ml))
@ 
A method for the function \command{image} can be used for displaying predicted values as shown in the next Figure, as well as  other prediction results returned by \command{krige.conv}.
\begin{figure}[h!]
\centering
<<fig=T>>=
image(kc, loc = pred.grid, col=gray(seq(1,0.1,l=30)), xlab="Coord X", ylab="Coord Y")
@ 
\label{fig:krige}
\caption{Map of the kriging estimates}
\end{figure}

\section{BAYESIAN ANALYSIS}

Bayesian analysis for Gaussian models is implemented by the function \code{krige.bayes}. 
It can be performed for different "degrees of uncertainty", meaning that model parameters can be treated as fixed or random.

As an example consider a model without nugget and including uncertainty in the mean, sill and range parameters. 
Prediction at the four locations indicated above is performed by typing a command like: 

\strong{WARNING: RUNNING THE NEXT COMMAND CAN BE TIME-CONSUMING}
<<label="Bayesian">>=
bsp4 <- krige.bayes(s100, loc = loci, prior = prior.control(phi.discrete = seq(0,5,l=101), phi.prior="rec"), output=output.control(n.post=5000))
@ 

Histograms showing posterior distribution for the model parameters can be plotted by typing:
\setkeys{Gin}{width=\textwidth}
\begin{figure}[h!]
\centering
<<fig=T,width=6,height=2>>=
par(mfrow=c(1,3), mar=c(3,3,1,.5), mgp=c(2,1,0))
hist(bsp4$posterior$sample$beta, main="", xlab=expression(beta), prob=T)
hist(bsp4$posterior$sample$sigmasq, main="", xlab=expression(sigma^2), prob=T)
hist(bsp4$posterior$sample$phi, main="", xlab=expression(phi), prob=T)
@ 
\label{fig:histpost}
\caption{Histograms of samples from posterior distribution}
\end{figure}
Using summaries of these posterior distributions (means, medians or modes)
we can check the "estimated Bayesian variograms"  
against the empirical variogram, as shown in the next figure. 
Notice that it is also possible to compare these estimates
with other fitted variograms such as the ones computed in Section 3.
\setkeys{Gin}{width=0.45\textwidth}
\begin{figure}[h!]
\centering
<<fig=T>>=
plot(bin1, ylim = c(0,1.5))
lines(bsp4, max.dist = 1.2, summ = mean)
lines(bsp4, max.dist = 1.2, summ = median, lty = 2)
lines(bsp4, max.dist = 1.2, summ = "mode", post="par",lwd = 2, lty = 2)
legend(0.25, 0.4, legend = c("variogram posterior mean", "variogram posterior median", "parameters posterior mode"), lty = c(1,2,2), lwd = c(1,1,2), cex = 0.8)
@ 
\label{fig:bayesvario}
\caption{Variogram models based on the posterior distributions}
\end{figure}

% must say something about options for the posterior for the entire variogram

The next figure shows predictive distributions at the four selected locations. 
Dashed lines show Gaussian distributions with mean and variance given by results of ordinary kriging obtained 
in Section 4. 
The full lines correspond to the Bayesian prediction. The plot shows results of density estimation 
using samples from the predictive distributions. 
<<eval=F>>=
par(mfrow=c(2,2))
for(i in 1:4){
  kpx <- seq(kc4$pred[i] - 3*sqrt(kc4$krige.var[i]), kc4$pred[i] +3*sqrt(kc4$krige.var[i]), l=100)
  kpy <- dnorm(kpx, mean=kc4$pred[i], sd=sqrt(kc4$krige.var[i]))
  bp <- density(bsp4$predic$sim[i,])
  rx <- range(c(kpx, bp$x))
  ry <- range(c(kpy, bp$y))
  plot(cbind(rx, ry), type="n", xlab=paste("Location", i), ylab="density", xlim=c(-4, 4), ylim=c(0,1.1))
  lines(kpx, kpy, lty=2)
  lines(bp)
}
@ 
\setkeys{Gin}{width=0.8\textwidth}
\begin{figure}[h!]
\centering
<<echo=F,fig=T,width=6,height=6>>=
par(mfrow=c(2,2))
for(i in 1:4){
  kpx <- seq(kc4$pred[i] - 3*sqrt(kc4$krige.var[i]), kc4$pred[i] +3*sqrt(kc4$krige.var[i]), l=100)
  kpy <- dnorm(kpx, mean=kc4$pred[i], sd=sqrt(kc4$krige.var[i]))
  bp <- density(bsp4$predic$sim[i,])
  rx <- range(c(kpx, bp$x))
  ry <- range(c(kpy, bp$y))
  plot(cbind(rx, ry), type="n", xlab=paste("Location", i), ylab="density", xlim=c(-4, 4), ylim=c(0,1.1))
  lines(kpx, kpy, lty=2)
  lines(bp)
}
@ 
\label{fig:locbayes}
\caption{Predictive distributions at the four selected locations}
\end{figure}

Consider now, under the same model assumptions, obtaining simulations from the predictive distributions on a grid of points covering the area. The commands to define the grid and perform Bayesian prediction are:
<<>>=
pred.grid <-  expand.grid(seq(0,1, l=31), seq(0,1, l=31))
#WARNING: RUNNING THE NEXT COMMAND CAN BE TIME-CONSUMING
bsp <- krige.bayes(s100, loc = pred.grid, prior = prior.control(phi.discrete = seq(0,5,l=51)), 
output=output.control(n.predictive=2))
@ 

Maps with the summaries and simulations of the predictive distribution can be plotted as follows.
<<eval=F>>=
par(mfrow=c(2,2), mar=c(3,3,1,.5), mgp=c(1.5,.7,0))
image(bsp, loc = pred.grid, main = "predicted", col=gray(seq(1,0.1,l=30)))
image(bsp, val ="variance", loc = pred.grid, 
main = "prediction variance", col=gray(seq(1,0.1,l=30)))
image(bsp, val = "simulation", number.col = 1, loc = pred.grid, 
main = "a simulation from\nthe predictive distribution", col=gray(seq(1,0.1,l=30)))
image(bsp, val = "simulation", number.col = 2,loc = pred.grid, 
main = "another simulation from \n the predictive distribution", col=gray(seq(1,0.1,l=30)))
@ 
\setkeys{Gin}{width=0.8\textwidth}
\begin{figure}[h!]
\centering
<<echo=F,fig=T,width=6,height=6>>=
par(mfrow=c(2,2), mar=c(3,3,1,.5), mgp=c(1.5,.7,0))
image(bsp, loc = pred.grid, main = "predicted", col=gray(seq(1,0.1,l=30)))
image(bsp, val ="variance", loc = pred.grid, 
main = "prediction variance", col=gray(seq(1,0.1,l=30)))
image(bsp, val = "simulation", number.col = 1, loc = pred.grid, 
main = "a simulation from\nthe predictive distribution", col=gray(seq(1,0.1,l=30)))
image(bsp, val = "simulation", number.col = 2,loc = pred.grid, 
main = "another simulation from \n the predictive distribution", col=gray(seq(1,0.1,l=30)))
@ 
\label{fig:mapbayes}
\caption{Maps obtained from the predictive distribution}
\end{figure}

\strong{Note:} Further examples for the function \code{krige.bayes} are given in the file
\url{http://www.est.ufpr.br/geoR/geoRdoc/examples.krige.bayes.R"} .

\section{Simulation of Gaussian Random Fields}

The function \code{grf} generates simulations of Gaussian random fields on regular or irregular sets of locations.
Some of its functionality is illustrated by the next commands.  
\setkeys{Gin}{width=\textwidth}
\begin{figure}[h!]
\centering
<<fig=T,width=8,height=4>>=
par(mfrow = c(1,2))
sim1 <- grf(100, cov.pars=c(1, .25))
points.geodata(sim1, main="simulated locations and values")
plot(sim1, max.dist=1, main="true and empirical variograms") 
@ 
\label{fig:grf1}
\caption{Simulation results produced with the function \texttt{grf}.}
\end{figure}

\begin{figure}[h!]
\centering
<<fig=T,width=8,height=4>>=
par(mfrow=c(1,2))
sim2 <- grf(441, grid="reg", cov.pars=c(1, .25))
image(sim2, main="a small-ish simulation", col=gray(seq(1, .1, l=30)))
sim3 <- grf(40401, grid="reg", cov.pars=c(10, .2), met="circ")
image(sim3, main="simulation on a fine grid", col=gray(seq(1, .1, l=30)))
@ 
\label{fig:grf2}
\caption{Simulations with different resolutions}
\end{figure}

\strong{NOTE:} we recommend the package 
\pkg{RandomFields} (\url{http://cran.r-project.org/src/contrib/PACKAGES.html#RandomFields}) for a more comprehensive implementation for simulation of Gaussian Random Fields. 


\section{Citing \pkg{geoR}}

The function \code{cite.geoR} shows information on how to cite geoR in publications, and we remember the function \code{citation} shows how to cite \R. 

<<>>=
cite.geoR()
citation()
@ 

\end{document}