<!--

-->


<title> R - PASI </title>


```{r, setup, include=FALSE}
##-----------------------------------------------------------------------------
## Knitr settings. Do not run.

opts_chunk$set(
    tidy=FALSE,
    cache=TRUE,
    results="hold",
    fig.path="figures/pasi",
    fig.width=9,
    fig.height=7,
    dpi=100,
    fig.align="center",
    fig.show="hold")
options(width=90)
options(geoR.messages=FALSE)
```

<center>

<h1><b>PASI</b></h1>

<h2>Pan-American Advanced Study Institute
on Spatio-Temporal Statistics</h2>

June 16-26, 2014

Búzios, RJ, Brazil


<h2><b><i>Practicals - geoR</b></i></h2>


[Paulo Justiniano Ribeiro Jr (I)](http://leg.ufpr.br/~paulojus/ "email: paulojus@ufpr.br"),   

(I) [LEG: Statistics and Geoinformation Lab]((http://leg.ufpr.br/)),
[DEST: Department of Statistics](http://www.est.ufpr.br/),
[UFPR: Federal University of Paraná](http://www.ufpr.br/)


</center>

<a id="toc"></a>**Table of contents**

  * [A (very) basic geoR session](#basic)
  * [Loading and exploring the data](#data)
  * [Model estimation](#estim)
  * [Model prediction](#predi)
  * [Bayesian inferance](#bayes)

*****************************************************************************

This is a brief overview of some functionalities in geoR.<br> 
We refer to the 
<a href="http://www.leg.ufpr.br/geoR/tutorials">geoR "Tutorials" page</a>
for further examples.

<a id="basic"></a>A (very) basic geoR session
-----------------------------------------------------------------------------

```{r, message=FALSE}
##-----------------------------------------------------------------------------
##
##        Pratical session: geostatistical analysis with geoR package
## 
##                       Paulo Justiniano Ribeiro Jr
## 
##-----------------------------------------------------------------------------
## Required package.
require(geoR)

## Loading a data-set included in the package
data(s100)
```

A quick data display with exploratory checks on spatial patterns, trends on the coordinates
and distribution of the data.
```{r, s100plot, fig.width=6, fig.height=6}
plot(s100, lowess=TRUE)
```

Another function allows for mode options for displaying the data.
```{r, s100points, fig.width=8, fig.height=4}
points(s100)
args(points.geodata)
```

The brief summary include information on the distances between data points.
```{r}
summary(s100)
```

There are several options for computing the empirical variogram.
Here we use default options except for the argument setting the maximum distance
to be considered between pairs of points. 
```{r, s100vario, fig.width=4, fig.height=4}
s100.v <- variog(s100, max.dist=1)
plot(s100.v)
```

Turning to the parameter estimation  we fit the model with a constant mean,
isotropic exponential correlation function and allowing for estimating the
variance of the conditionally independent observations (nugget variance)

```{r, s100-ML}
(s100.ml <- likfit(s100, ini=c(1, 0.15)))
summary(s100.ml)
```

Finally, we start the spatial prediction defining a grid of points.
The kriging function by default performs ordinary kriging 
It minimaly requires the data, prediction locations and estimated model parameters.
```{r, s100-kriging}
s100.gr <- expand.grid((0:100)/100, (0:100)/100)
s100.kc <- krige.conv(s100, locations=s100.gr, krige=krige.control(obj.model=s100.ml))
names(s100.kc)
```

If the locations forms a grid the predictions can be visualised
as image, contour or persp plots.
Image plots shows the predicted values (left) and the respective std errors.
```{r, s100kc, fig.width=8, fig.height=4}
par(mfrow=c(1,2), mar=c(3.5, 3.5, 0.5, 0.5))
image(s100.kc, col=gray(seq(1, 0.2, l=21)))
contour(s100.kc, nlevels=11, add=TRUE)
image(s100.kc, val = sqrt(krige.var), col=gray(seq(1, 0.2, l=21)))
```

It is worth notice that under the Gaussian model the prediction errors
depends on the coordinates of the data, and not to their values (only through 
model parameter estimates).
```{r, s100kcerr, fig.width=4, fig.height=4}
image(s100.kc, val = sqrt(krige.var), col=gray(seq(1, 0.2, l=21)), coords.data=TRUE)
```

[Back to table of contents](#toc)

<a id="data"></a>Loading and exploring the data
-----------------------------------------------------------------------------

There are some datasets included in geoR which can be listed as follows.
We use some of them to illustrate further options in the functions. 
```{r, eval=F}
##-----------------------------------------------------------------------------
## Datasets included in geoR.
data(package="geoR")
```

Some of them are of the class "geodata" which is convenient (but not compulsory!)
to use with geoR' functions.
Others are simply data.frames from which "geodata" can be created. 
"geodata" (S3) objects are just lists with a class "geodata" and 
must have elements "coords" and "data". Optional elements includes
"covariates", "units.m" among others (see help(as.geodata)).
```{r, some-datasets}
class(s100)
class(parana)
class(Ksat)
class(camg)
head(camg)
ctc020 <- as.geodata(camg, coords.col=c(1,2), data.col=7, covar.col=c(3,4))
```

The function calls shown above mostly use default options and it is
worth examine the arguments (and the help files). 
```{r, function-arguments}
args(plot.geodata)
args(points.geodata)
args(variog)
```

We point in particular to three arguments:
<ul>
<li><b>borders:</b> a matrix with coordinates of a polygon limiting the area</li>
<li><b>lambda:</b> parameter of the BoxCox transformation</li>
<li><b>trend:</b> a matrix, formula or special words ("1st", "2nd") specifying a mean model</li>
</ul>

Borders are conveniently dealt with if included in the <tt>geodata</tt> object.
As an example lets add borders to the s100 data
```{r, s100borders, fig.width=5, fig.height=5}
s100$borders <- matrix(c(0,0,1,0,1,1,0,1,0,0), ncol=2, byrow=TRUE)
points(s100)
```

If you don't have the coordinates of the border consider setting them manually by using "locator()"
over the "points()" plot.


The dataset Ksat with data on the saturated soil conductivity
has a clear asymetric behaviour and the BoxCox transformation suggests to
log transform the data.
```{r, ksatboxcox, fig.width=4, fig.height=4}
data(Ksat)
require(MASS)
summary(Ksat)
boxcox(Ksat)
```

Differences on exploratory plots are clear.
```{r, ksat-plot, fig.width=6, fig.height=6}
plot(Ksat, lowess=TRUE)
```

```{r, ksat0-plot, fig.width=6, fig.height=6}
plot(Ksat, lowess=TRUE, lambda=0)
```

We now use the Paraná rainfall data to illustrate
the usage and effect of trend argument. Notice also the data also include
the borders of the state
```{r, parana-data}
data(parana)
summary(parana)
```

It is clear from the exploratory plot that there is a gradient of the 
response along the L-W and N-S coordinates.  
```{r, parana-plot, fig.width=6, fig.height=6}
plot(parana, lowess=TRUE)
```

If the argument trend is provided a (ordinary) regression is fitted 
and the residuals are plotted. The options <tt>trend="1st"</tt>
defines a linear trend on the coordinates and
is equivalent to <tt>trend=~coords</tt>.
More generally a formula on a set of covariates can be used.
```{r, parana-trend-plot, fig.width=6, fig.height=6}
plot(parana, trend = "1st", lowess=TRUE)
```

The effect on the trend in the variograms are very much noticeable in this example.
```{r, parana-trend-vario, fig.width=8, fig.height=4}
par(mfrow=c(1,2))
parana.vario <- variog(parana, max.dist=400)
plot(parana.vario)
parana.variot <- variog(parana, trend="1st",  max.dist=400)
plot(parana.variot)
```

Next we show a simple Monte Carlo test using the empirical variogram.
The envelope reflects the assumption of no spatial correlation and
is given by permutations of the data across the sampling locations.
```{r, s100.vario.env, fig.width=4, fig.height=4}
s100.v <- variog(s100, max.dist=1)
s100.v.env <- variog.mc.env(s100, obj.variog=s100.v)
plot(s100.v, env=v.env)
```

To close this section we show a helper function which facilitate
building empirical variograms on different directions which may be used 
to investigate anisotropy.
```{r, s100.vario4, fig.width=4, fig.height=4}
s100.v4 <- variog4(s100, max.dist=1)
plot(s100.v4, env=s100.v.env, omni=TRUE)
```

[Back to table of contents](#toc)

<a id="estim"></a>Model estimation
-----------------------------------------------------------------------------
The are four methods for estimating the (covariance) model parameters implemented in geoR.
<ul>
<li><b>visual variogram fitting:</b> using the function eyefit().</li>
<li><b>variogram fitting:</b> using the function variofit() basically fitting a nnon-linear model to the empirical variogram. </li>
<li><b>likelihood (ML or REML):</b> using the function likfit().</li>
<li><b>Bayesian inference:</b> using the function krige.bayes().</li>
</ul>

We now illustrate the first three, leaving Bayesian inference for later.

eyefit() provides an interactive visual model fitting tool which can be used 
on its own to provide estimates, or to set reasonable initial values for 
the other methods. Try experimenting with the following.

```{r, s100-eyefit, eval=FALSE}
data(s100)
s100.v <- variog(s100, max.dist=1)
plot(s100)
eyefit(s100)
```

variofit() allows for several options.
In what follows we show variogram fits for the parana data
under constant mean and a trend  on the coordinates. By default a
exponential correlation function (which corresponds to the Martèrn model with kappa=0.5) is used. 
We also show fits for e spherical and Matèrn (with kappa=1.5) variogram models. 
```{r, parana-variofit, fig.width=8, fig.height=4}
parana.vfit.exp <- variofit(parana.vario)
parana.vfit.mat1.5  <- variofit(parana.vario, kappa=1.5)
parana.vfit.sph  <- variofit(parana.vario, cov.model="sph")
#
parana.vtfit.exp <- variofit(parana.variot)
parana.vtfit.mat1.5  <- variofit(parana.variot, kappa=1.5)
parana.vtfit.sph  <- variofit(parana.variot, cov.model="sph")
par(mfrow=c(1,2))
plot(parana.vario)
lines(parana.vfit.exp);lines(parana.vfit.mat1.5, col=2);lines(parana.vfit.sph, col=4)
plot(parana.variot)
lines(parana.vtfit.exp);lines(parana.vtfit.mat1.5, col=2);lines(parana.vtfit.sph, col=4)
```

Maximum likelihood estimates are obtained fitting the model to the 
data (and not to the empirical variogram). Models can be compared by
fitting measures.
```{r, parana-fits}
(parana.ml0 <- likfit(parana, ini=c(4500,50), nug=500))
(parana.ml1 <- likfit(parana, trend="1st", ini=c(1000,50), nug=100))
(parana.ml2 <- likfit(parana, trend="2nd", ini=c(1000,50), nug=100))
logLik(parana.ml0)
logLik(parana.ml1)
logLik(parana.ml2)
```

[Back to table of contents](#toc)


<a id="predi"></a>Model prediction
-----------------------------------------------------------------------------

There are two functions for spatial prediction in geoR.
<ul>
<li><b>krige.conv()</b> performing "plug-in" predictions, i.e. treating estimated parameters as if they were the true values</li>
<li><b>krige.bayes()</b> for Bayesian inference and prediction</li>
</ul>

For the former, the function takes: the data, coordinates of the locations,
kriging options (through krige.control()) and output options (through output.control()).
We use the Paraná rainfall data to illustrate more options for this function.

The first step is to define a grid of prediction locations.
```{r, include=FALSE}
options(warn=-1)
```

```{r, parana-grid, fig.width=6, fig.height=6}
parana.gr <- pred_grid(parana$borders, by=15)
points(parana)
points(parana.gr, pch=19, col=2, cex=0.25)
parana.gr0 <- locations.inside(parana.gr, parana$borders)
points(parana.gr0, pch=19, col=4, cex=0.25)
```

```{r, loadkrige, include=FALSE}
load("parana.kc")
load("parana.bayes")
```


The following kriging call corresponds to the "universal kriging".
The arguments "trend.d" and "trend.l" require the specification
(or terms in a formula) on both, data and prediction locations, respectively.

```{r, parana-kriging}
args(krige.control)
args(output.control)
KC <- krige.control(obj.m = parana.ml1, trend.d="1st", trend.l="1st")
OC <- output.control(simulations=TRUE, n.pred=1000,
		     quantile=c(0.10, 0.25, 0.5, 0.75, 0.90), threshold = 350) 
parana.kc <- krige.conv(parana, loc=parana.gr, krige = KC, output = OC) 
names(parana.kc)
```

From the output provided we picked the following in the next plot:
prediction means, medians, a simulation from the (joint) predictive distribution,
and the probability of exceeding the value of 350.

```{r, parana-krige, fig.width=6, fig.height=6}
par(mfrow=c(2,2), mar=c(3,3,0.5, 0.5))
image(parana.kc, col = terrain.colors(21), x.leg=c(500, 750), y.leg=c(0, 50))
image(parana.kc, val = parana.kc$quantile[,3], col=terrain.colors(21), 
	x.leg=c(500, 750), y.leg=c(0, 50))
image(parana.kc, val = parana.kc$simulation[,1], col=terrain.colors(21), 
	x.leg=c(500, 750), y.leg=c(0, 50))
image(parana.kc, val = 1-parana.kc$prob, col=terrain.colors(21), 
	x.leg=c(500, 750), y.leg=c(0, 50))
```

More generally, simulations can be used to derive
maps of quantities of interest. The top panels in next example calls 
extract minimum and maximum simulated values at each location. 
The bottom panels shows probabilities of exceeding 300 at each location,
and predictive distribution of the proportion of the area exceeding 300. 

```{r, parana-krige2, fig.width=6, fig.height=6}
par(mfrow=c(2,2), mar=c(3,3,0.5, 0.5))
image(parana.kc, val = apply(parana.kc$simulation, 1, min), col=terrain.colors(21), 
	x.leg=c(500, 750), y.leg=c(0, 50))
image(parana.kc, val = apply(parana.kc$simulation, 1, max), col=terrain.colors(21), 
	x.leg=c(500, 750), y.leg=c(0, 50))
image(parana.kc, val=apply(parana.kc$simulations, 1, function(x) mean(x > 300, na.rm=T)), 
                 x.leg=c(500, 750), y.leg=c(0,50), col=gray(seq(1, 0.2, length=21)))
hist(apply(parana.kc$simulations, 2, function(x) mean(x > 300, na.rm=T)), main="")
```

We now pick four particular locations for prediction
```{r, parana-kc4, fig.width=4, fig.height=4}
loc4 <- cbind(c(300, 480, 680, 244), c(460, 260, 170, 270))
parana.kc4 <- krige.conv(parana, loc=loc4, krige = KC, output = OC) 
points(parana)
points(loc4, col=2, pch=19)
```

```{r, parana-kc4-hist, fig.width=8, fig.height=2, echo=FALSE}
par(mfrow=c(1,4), mar=c(3.5,3.5, 0.5,0.5))
apply(parana.kc4$simulation, 1, 
      function(x){ hist(x, prob=T, xlim=range(parana$data), main=""); lines(density(x))})
```

```{r, include=FALSE}
options(warn=0)
```

[Back to table of contents](#toc)


<a id="bayes"></a>Bayesian Inference

For the Bayesian analysis both, inference on model parameters
and spatial predictions  (if prediction locations are provided) are performed with a call to the function krige.bayes(). In the next call we run an analysis 
setting a discrete prior for the parameter in the (exponential) correlation function
and fixing the "nugget" parameter to zero. 
More general call can set a discrete prior also for the relative nugget through the 
"tausq.rel" parameter.
```{r}
args(krige.bayes)
parana.bayes <- krige.bayes(parana, loc=parana.gr, 
                           model = model.control(trend.d="1st", trend.l="1st"),
		 	   prior= prior.control(phi.prior="rec", phi.disc=seq(0,150, by=15)),
			   output = OC) 
names(parana.bayes)
```

```{r, parana-bayes-prior-post, fig.width=4, fig.height=4}
names(parana.bayes$posterior)
plot(parana.bayes)
```

For predictions, the output is stored in the  obj$predictive are similar to the returned by
krige.conv() and can be explored and ploted using similar function calls.

```{r, parana-bayes-pred, fig.width=9, fig.height=3}
names(parana.bayes$predictive)
par(mfrow=c(1,3))
image(parana.bayes, col=terrain.colors(21))
image(parana.bayes, val = apply(parana.bayes$pred$simul, 1, quantile, prob=0.90),
      col=terrain.colors(21))
hist(apply(parana.bayes$pred$simul, 2, median), main="")
```
-----------------------------------------------------------------------------

<a id="cite"></a>How to cite geoR
```{r}
citation(package="geoR")
```

This page was produced with the <a href="http://yihui.name/knitr/">knitr</a>
using this <a href="http://www.leg.ufpr.br/geoR/PASI/geoR-PASI.Rmd">Rmarkdown</a> source file.
The HTML was generated with the following command
```{r eval=F}
require(knitr)
knit2html("geoR-PASI.Rmd")
```
and  <a href="http://www.leg.ufpr.br/geoR/PASI/geoR-PASI.R">a file with only the <b>R</b> code</a> 
can be extracted (tangle) with:
```{r, include=FALSE, eval=FALSE}
purl("geoR-PASI.Rmd")
```

[Back to table of contents](#toc)

