---
title: "Inference for Totals and Weighted Sums from Finite Spatial Populations"
author:
  - Matt Higham, Jay M. Ver Hoef, Bryce M. Frank
date: "`r Sys.Date()`"
output:
  prettydoc::html_pretty:
    theme: architect
    toc: true
vignette: >
  %\VignetteIndexEntry{Inference for Totals and Weighted Sums from Finite Spatial Populations}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```
```{css, echo = FALSE}
caption, .caption{
  font-style:italic;
  margin-top:0.5em;
  margin-bottom:0.5em;
  width:99%;
  text-align: left;
}
body{
  font-size: 12pt;
}
p {line-height: 1.7;}
tr:nth-child(even) {background-color: #f2f2f2;}
th {
    background-color: #4CAF50;
    color: black;
}
p {
padding-top: 7px;
padding-bottom: 7px;
}

```

```{r, echo = FALSE, message = FALSE}
########################################################################
########################################################################
########################################################################
#        Introduction
########################################################################
########################################################################
########################################################################
```

# Introduction

The `sptotal` package was developed for predicting a weighted sum, most commonly a mean or total, from a finite number of sample units in a fixed geographic area. Estimating totals and means from a finite population is an important goal for both academic research and management of environmental data. One naturally turns to classical sampling methods, such as simple random sampling or stratified random sampling. Classical sampling methods depend on probability-based sample designs and are robust. Very few assumptions are required because the probability distribution for inference comes from the sample design, which is known and under our control.  For design-based methods, sample plots are chosen at random, they are measured or counted, and inference is obtained from the probability of sampling those units randomly based on the design (e.g., Horwitz-Thompson estimation). As an alternative, we will use model-based methods, specifically geostatistics, to accomplish the same goals. Geostatistics does not rely on a specific sampling design. Instead, when using geostatistics, we assume the data were produced by a stochastic process with parameters that can be estimated. The relevant theory is given by Ver Hoef (2008). The `sptotal` package puts much of the code and plots in Ver Hoef (2008) in easily accessible, convenient functions.

In the `sptotal` package, our goal is to estimate some linear function of all of the sample units, call it $\tau(\mathbf{z}) = \mathbf{b}^\prime \mathbf{z}$, where $\mathbf{z}$ is a vector of the realized values for all the sample units and $\mathbf{b}$ is a vector of weights. By "realized," we mean that whatever processes produced the data have already happened, and that, if we had enough resources, we could measure them all, obtaining a complete census. If $\tau(\mathbf{z})$ is a population total, then every element of $\mathbf{b}$ contains a $1$. Generally, $\mathbf{b}$ can contain any set of weights that we would like to multiply times each value in a population, and then these are summed, yielding a weighted sum. 

The vector $\mathbf{b}$ contains the weights that we would apply if we could measure or count every observation, but, because of cost consideration, we usually only have a sample.

```{r, echo = FALSE, message = FALSE}
########################################################################
########################################################################
########################################################################
#        Data
########################################################################
########################################################################
########################################################################
```

# Data 

Prior to using the `sptotal` package, the data needs to be in `R` in the proper format. For this package, we assume that your data set is a `data.frame()` object, described below.

## Data Frame Structure

Data input for the `sptotal` package is a `data.frame`. The basic information required to fit a spatial linear model, and make predictions, are the response variable, covariates, the x- and y-coordinates, and a column of weights. You can envision your whole population of possible samples as a `data.frame` organized as follows,

```{r, echo = FALSE, message = FALSE, fig.align="center", cache = FALSE}
old.par = par(mar = c(0,0,0,0))
plot(c(0,1), c(0,1), type = 'n', bty = 'n', xaxt = 'n', yaxt = 'n', 
  xlab = '', ylab = '')
rect(0, .7, .08, 1, col = rgb(228/255,26/255,28/255))
rect(0, 0, .08, .68, col = 'white')
rect(.1, .7, .48, 1, col = rgb(55/255,126/255,184/255))
rect(.1, 0, .48, .68, col = rgb(55/255,126/255,184/255, alpha = .3))
rect(.5, .7, .58, 1, col = rgb(77/255,175/255,74/255))
rect(.6, .7, .68, 1, col = rgb(77/255,175/255,74/255))
rect(.5, 0, .58, .68, col = rgb(77/255,175/255,74/255, alpha = .3))
rect(.6, 0, .68, .68, col = rgb(77/255,175/255,74/255, alpha = .3))
rect(.7, .7, .78, 1, col = rgb(152/255,78/255,163/255))
rect(.7, 0, .78, .68, col = rgb(152/255,78/255,163/255, alpha = .3))
rect(.8, .7, 1, 1, col = rgb(255/255,127/255,0/255))
rect(.8, 0, 1, .68, col = rgb(255/255,127/255,0/255, alpha = .3))
par(old.par)
```

where the red rectangle represents the column of the response variable, and the top part, colored in red, are observed locations, and the lower part, colored in white, are the unobserved values. To the right, colored in blue, are possibly several columns containing covariates thought to be _predictive_ for the response value at each location. Covariates must be known for both observed and unobserved locations, and the covariates for unobserved locations are shown as pale blue below the darker blue covariates for observed locations above. It is also possible that there are no available covariates. 

The `data.frame` must have x- and y-coordinates, and they are shown as two columns colored in green, with the coordinates for the unobserved locations shown as pale green below the darker green coordinates for the observed locations above. The `data.frame` can have a column of weights. If one is not provided, we assume a column of all ones so that the prediction is for the population total. The column of weights is purple, with weights for the observed locations a darker shade, above the lighter shade of purple representing weights for unsampled locations. Finally, the `data.frame` may contain columns that are not relevant to predicting the weighted sum. These columns are represented by the orange color, with the sampled locations a darker shade, above the unsampled locations with the ligher shade.

Of course, the data do not have to be in exactly this order, either in terms of rows or columns.  Sampled and unsampled rows can be intermingled, and columns of response variable, covariates, coordinates, and weights can be also be intermingled. The figure above is an idealized graphic of the data. However, this figure helps envision how the data are used and illustrate the goal. We desire a weighted sum, where the weights (in the purple column) are multiplied with the response variable (red/white) column, and then summed.  Because some of the response values are unknown (the white values in the response column), covariates and spatial information (obtained from the x- and y-coordinates) are used to _predict_ the unobserved (white) values. The weights (purple) are then applied to both the observed response values (red), and the predicted response values (white), to obtain a weighted sum. Because we use predictions for unobserved response values, it is important to assess our uncertainty, and the software provides both an estimate of the weighted sum, mean, or total for the response variable as well as its estimated prediction variance.  

## Simulated Data Creation

To demonstrate the package, we created some simulated data so they are perfectly behaved, and we know exactly how they were produced. Here, we give a brief description before using the main features of the package. To get started, install the package

```{r, eval = FALSE}
install.packages("sptotal")
```

and then type

```{r}
library(sptotal)
```

Type 

```{r}
data(simdata)
```

and then `simdata` will be available in your workspace. To see the structure of `simdata`, type

```{r}
head(simdata)
```

Here, we see that `simdata` is a data frame with 400 records. The spatial coordinates are `numeric` variables in columns named `x` and `y`.  We created 7 continuous covariates, `X1` through `X7`.  The variables `X1` through `X5` were all created using the `rnorm()` function, so they are all standard normal variates that are independent between and within variable.  Variables `X6` and `X7` were independent from each other, but spatially autocorrelated within, each with a variance parameter of 1, an autocorrelation range parameter of 0.2 from an exponential model, and a small nugget effect of 0.01.  The variables `F1` and `F2` are factor variables with 3 and 5 levels, respectively. The variable `Z` is the response. Data were simulated from the model

\begin{align*}
Z_i = 10 & + 0 \cdot X1_i + 0.1 \cdot X2_i + 0.2 \cdot X3_i + 0.3 \cdot X4_i + \\ & 0.4 \cdot X5_i + 
    0.4 \cdot X6_i + 0.1 \cdot X7_i + F1_i + F2_i + \delta_i + \varepsilon_i
\end{align*}

where factor levels for `F1` have effects $0, 0.4, 0.8$, and factor levels for `F2` have effects $0, 0.1, 0.2, 0.3, 0.4$. The random errors $\{\delta_i\}$ are spatially autocorrelated from an exponential model,

$$
\textrm{cov}(\delta_i,\delta_j) = 2*\exp(-d_{i,j})
$$  

where $d_{i,j}$ is Euclidean distance between locations $i$ and $j$. In geostatistics terminology, this model has a partial sill of 2 and a range of 1. The random errors $\{\varepsilon_i\}$ are independent with variance 0.02, and this variance is called the nugget effect. Two columns with weights are included, `wts1` contains 1/400 for each row, so the weighted sum will yield a prediction of the overall mean. The column `wts2` contains a 1 for 25 locations, and 0 elsewhere, so the weighted sum will be a prediction of a total in the subset of 25 locations.

The spatial locations of `simdata` are in a $20 \times 20$ grid uniformly spaced in a box with sides of length 1,

```{r, fig.width = 5, fig.height = 5, fig.align = "center", message = FALSE, cache = FALSE}
require(ggplot2)
ggplot(data = simdata, aes(x = x, y = y)) + geom_point(size = 3) +
  geom_point(data = subset(simdata, wts2 == 1), colour = "red",
    size = 3)
```

The locations of the 25 sites where `wts2` is equal to one are shown in red.

We have simulated the data for the whole population. This is convenient, because we know the true means and totals. In order to compare with the prediction from the `sptotal` package, let's find the true population total

```{r}
sum(simdata[ ,'Z'])
```

as well as the total in the subset of 25 sites

```{r}
sum(simdata[ ,'wts2'] * simdata[ ,'Z'])
```

However, we will now sample from this population to provide a more realistic setting where we can measure only a part of the whole population. In order to make results reproducible, we use the `set.seed` command, along with `sample`.  The code below will replace some of the response values with `NA` to represent the unsampled sites.

```{r}
set.seed(1)
# take a random sample of 100
obsID <- sample(1:nrow(simdata), 100)
simobs <- simdata
simobs$Z <- NA
simobs[obsID, 'Z'] <- simdata[obsID, 'Z']
```

We now have a data set where the whole population is known, `simdata`, and another one, `simobs`, where 75% of the response variable of the population has been replaced by `NA`.  Next we show the sampled sites as solid circles, while the missing values are shown as open circles, and we use red again to show the sites within the small area of 25 locations.

```{r, fig.width = 5, fig.height = 5, fig.align = "center", cache = FALSE}
ggplot(data = simobs, aes(x = x, y = y)) +
  geom_point(shape = 1, size = 2.5, stroke = 1.5) +
  geom_point(data = subset(simobs, !is.na(Z)), shape = 16, size = 3.5) +
  geom_point(data = subset(simobs, !is.na(Z) & wts2 == 1), shape = 16,
    colour = "red", size = 3.5) +
  geom_point(data = subset(simobs, is.na(Z) & wts2 == 1), shape = 1,
    colour = "red", size = 2.5, stroke = 1.5)
```

We will use the `simobs` data to illustrate use of the `sptotal` package.

```{r, echo = FALSE, message = FALSE}
########################################################################
########################################################################
########################################################################
#        Using the sptotal package
########################################################################
########################################################################
########################################################################
```

# Using the `sptotal` Package

After your data is in a similar format to `simobs`, using the `sptotal` package occurs in two primary stages. In the first, we fit a spatial linear model. This stage estimates spatial regression coefficients and spatial autocorrelation parameters. In the second stage, we predict the unsampled locations for the response value, and create a prediction for the weighted sum (e.g. the total) of all response variable values, both observed and predicted. To show how the package works, we demonstrate on ideal, simulated data. Then, we give a realistic example on moose data and a second example on lakes data to provide further insight and documentation. The moose example also has a section on data preparation steps.

## Fitting a Spatial Linear Model: `slmfit`

We continue with our use of the simulated data, `simobs`, to illustrate fitting the spatial linear model. The spatial model-fitting function is `slmfit` (spatial-linear-model-fit), which uses a formula like many other model-fitting functions in `R` (e.g., the `lm()` function). To fit a basic spatial linear model we use

```{r}
slmfit_out1 <- slmfit(formula = Z ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + F1 + F2, 
  data = simobs, xcoordcol = 'x', ycoordcol = 'y',
  CorModel = "Exponential")
```

The documentation describes the arguments in more detail, but as mentioned earlier, the linear model includes a formula argument, and the `data.frame` that is being used as a data set. We also need to include which columns contain the $x$- and $y$-coordinates, which are arguments to `xcoordcol` and `ycoordcol`, respectively. In the above example, we specify `'x'` and `'y'` as the column coordinates arguments since the names of the coordinate columns in our simulated data set are `'x'` and `'y'`. We also need to specify a spatial autocorrelation model, which is given by the `CorModel` argument. As with many other linear model fits, we can obtain a summary of the model fit,

```{r}
summary(slmfit_out1) 
```

The output looks similar to the `summary` of a standard `lm` object, but there is some extra output at the end that gives our fitted covariance parameters. 

## Prediction: `predict`

After we have obtained a fitted spatial linear model, we can use the `predict.slmfit` function to construct a data frame of predictions for the unsampled sites. By default, the `predict` function assumes that we are predicting the population total and outputs this predicted total, the prediction variance for the total, as well as the original data frame input into `slmfit` appended with site-by-site predictions and site-by-site prediction variances. We name this object `pred_obj` in the chunk below.

```{r, results = "hide"}
pred_obj <- predict(slmfit_out1)
```

The output of the `predict` function is too lengthy to print here and can be cumbersome to read. To obtain useful output, we have a couple of choices. If you are familiar with `R`, then you might want to construct your own summary output, maps, variograms, etc. One of the outputs of `predict` is a data frame that is appended with the site by site predictions that would facilitate construction of, for example, a map with predictions. This is particularly useful if you have a shapefile. The site-by-site predictions for density are given by the variable `name_of_response_pred_density` while the site-by-site predictions for counts are given by `name_of_response_pred_count`. These two columns will only differ if you have provided a column for areas of each site. 

```{r, results = "hide"}
prediction_df <- pred_obj$Pred_df
prediction_df[ ,c("x", "y", "Z", "Z_pred_density")]
```

## Examining results: `check.variogram`, `get.predinfo`, and `get.predplot`

Alternatively, you could take advantage of some pre-specified package output functions `check.variogram`, `get.predinfo`, and `get.predplot`. 

First, `check.variogram` constructs an empirical variogram of the model residuals and overlays the REML-fitted variogram model with a line:

```{r, warning = FALSE}
check.variogram(slmfit_out1)
```

Note that the fitted line may not appear to fit the empirical variogram perfectly for a couple of reasons. First, only pairs of points that have a distance between 0 and one-half the maximum distance are shown. Second, the fitted model is estimated using REML, which may give different results than using weighted least squares.

We can also examine a histogram of the residuals as well as a histogram of the cross-validation (leave-one-out) residuals:

```{r, fig.height = 3}
qplot(residuals(slmfit_out1), bins = 20) + xlab("Residuals")
qplot(residuals(slmfit_out1, cross.validation = TRUE), bins = 20) +
  xlab("CV Residuals")
```

There is still one somewhat large cross-validation residual for an observed count that is larger than what would be predicted from a model without that particular count. The cause of this somewhat large residual can be attributed to random chance because we know that the data was simulated to follow all assumptions. 

If we are satisfied with the fit of the variogram, we can examine our prediction and obtain a confidence interval for this prediction using `get.predinfo`. While the `check.variogram` function takes the fitted spatial linear model as input, the `get.predinfo` function uses the object from the `predict` function, which we previously named `pred_obj`.

```{r}
get.predinfo(pred_obj, conf_level = 0.90)
```

We predict a total of 4817 units in this simulated region with 90% confidence bounds of (4779, 4856). The confidence interval is fairly small because we simulated data that were highly correlated, increasing precision in prediction for unobserved sites. You can see that the prediction of the total is close to the true value of 4834.326, and the true value is within the confidence interval.

Finally, to get a basic plot of the predictions, we can use the `get.predplot` function. 

```{r}
get.predplot(pred_obj)
```

The map shows the distribution of the response across sampled and unsampled sites. Its purpose is simply to give the user a quick idea of the distribution of the response. For example, we see from the plot that the predicted response is low in the upper-right region of the graph, is high in the middle of the region and in the upper-left corner of the region, and is low again at the lower portion of the area of interest. However, using the prediction data frame generated from the `predict` function, you can use `ggplot2` or any other plotting package to construct your own map that may be more useful in your context.

###Prediction for a Small Area of Interest

Spatial prediction can be used to estimate means and totals over finite populations from geographic regions, but can also be used for the special case of estimating a mean or total in a small area of interest. The term small area estimation refers to making an inference on a smaller geographic area within the overall study area. There may be few or no samples within that small area, so that estimation by classical sampling methods may not be possible or variances become exceedingly large. 

If we want to predict a quantity other than the population total, then we need to specify the column in our data set that has the appropriate prediction weights in a `wtscol` argument. For example, we might want to predict the total for a small area of interest. if we want to predict the total for the 25 sites in coloured in red, then we can use

```{r}
pred_obj2 <- predict(slmfit_out1, wtscol = "wts2")
get.predinfo(pred_obj2)
```

Recall that the true total for this small area was 273.4. We see that this is close to our prediction of 282.2 and is also within the bounds of our confidence interval. 

# Real Data Examples

## Moose Abundance from Aerial Surveys

The simulated data example assumes that the coordinates are a Transverse Mercator projection (TM), that the vector of the response is numeric and has `NA` values for sites that were not sampled, and that the areas of each site sampled are all the same. If this isn't the case for the data set you are working with, the following moose abundance example can help prepare your data for the functions in `sptotal`.

For an example with real data, we consider a data set on moose abundance in Alaska obtained from [Alaska Department of Fish and Game, Division of Wildlife Conservation](https://www.adfg.alaska.gov/index.cfm?adfg=divisions.wcoverview). Each observation corresponds to a moose counted at a particular site, but operational constraints do not permit all sites to be counted. Begin by loading the data into `R`. Unlike the simulated data, `AKmoose` is an `sp` object. In order to use the functions in this package, we need to extract the coordinates and relevant data from the `sp` object and put this information into a rectangular data frame. The easy part is getting the data; the more difficult part is getting the coordinates into a usable form. 

### Spatial Coordinates

Our goal is to append the x and y-coordinates to the data frame with the survey data. 

```{r, message = FALSE}
data(AKmoose)
```

We want to get the centroids of each of the sites and combine the centroids with the survey data. To obtain centroids, we use the `rgeos` package.

```{r, message = FALSE}
require(rgeos)
centroids <- data.frame(ID = AKmoose@data,
  x = rgeos::gCentroid(AKmoose, byid=TRUE)@coords[ ,'x'],
  y = rgeos::gCentroid(AKmoose, byid=TRUE)@coords[ ,'y'])
```

Next, for most spatial prediction, we want to use a transverse Mercator projection instead of latitude-longitude so that physical distance between two sites is accurately represented. The `LLtoTM()` function in this package provides a convenient way to convert latitude/longitude coordinates into user-defined transverse Mercator coordinates. 

```{r}
xy <- LLtoTM(mean(centroids$x), centroids$y, centroids$x)$xy
```

Transverse Mercator is based on minimizing distortion from a central meridian: the first argument specifies that the central meridian should be the mean longitude value in the data set. The second and third arguments to this function give the centroid latitude values and centroid longitude values. Finally, we add the transformed coordinates to our data frame with the survey data. We first extract the survey data from our `sp` object using

```{r}
moose_df <- AKmoose@data ## name the data set moose_df
head(moose_df) ## look at the first 6 observations
```

We see that, in addition to the `total` column, which has counts of moose, the data set also has `strat`, a covariate that is either `L` for Low or `M` for medium, and `surveyed`, which is a `0` if the site wasn't sampled and a `1` if the site was sampled. 

And then the tranformed coordinates can be added to the survey data frame.

```{r}
moose_df$x = xy[ ,'x']
moose_df$y = xy[ ,'y']
```

It might be helpful to compare the latitude and longitude coordinates of the original data frame to the tranformed coordinates in the new data frame to make sure that the transformation seems reasonable:

```{r, results = "hide"}
cbind(moose_df$x, moose_df$y, centroids$x, centroids$y)
```

Now, the `moose_df` data frame is in a more workable form for the `sptotal` package. However, there are still a couple of issues involving how the count data is stored and which sites were sampled that may be somewhat common in real data sets, which we address next.

### Count Vector Specifications

Let's look specifically at the counts in this moose data set in the `total` column:

```{r}
head(moose_df)
str(moose_df$total)
```

The first issue is that our original `sp` object had `total` as a factor, which `R` treats as a categorical variable. `total` should be numeric, and, in fact, the variable `surveyed` has the same issue. If we were to keep `total` as a factor and try to run `slmfit`, we would get a convenient error message, reminding us to make sure that our response variable is numeric, not a factor or character:

```{r, error = TRUE, purl = FALSE}
slmfit_out_moose <- slmfit(formula = total ~ strat, 
  data = moose_df, xcoordcol = 'x', ycoordcol = 'y',
  CorModel = "Exponential")
```

We first want to convert these two columns into numeric variables instead of factors. There are packages that can help with this conversion, like `dplyr` and `forcats`, but we opt for base `R` functions here.

```{r}
moose_df$surveyed <- as.numeric(levels(moose_df$surveyed))[moose_df$surveyed]
moose_df$total <- as.numeric(levels(moose_df$total))[moose_df$total]
```

This may not be an issue with the data frame you are working with, in which case you can ignore the above code. The `str()` command will tell you whether your variables are coded as factors or numeric.

After conversion to numeric variables, note that the first 6 observations for the `total` variable are all 0, but, the first two sites and the fourth, fifth, and sixth sites weren't actually sampled. Without some modification to this variable, `sptotal` wouldn't be able to differentiate between zeroes that were zero due to a site really having 0 counts or 0 density at the site and zeroes that were zero due to the site not being sampled. The following code converts the `total` variable on sites that were __not__ surveyed (`surveyed` = `0`) to `NA`.

```{r}
moose_df$total[moose_df$surveyed == 0] <- NA
head(moose_df)
```

### Fitting the Model and Obtaining Predictions

Now that

* we have x and y coordinates in TM format,

* our response variable is numeric and not a factor, and

* the column with our counts has `NA` values for sites that were not surveyed,

we can proceed to use the functions in `sptotal` in a similar way to how the functions were used for the simulated data. To get a sense of the data, we first give a plot of the raw observed counts:

```{r}
ggplot(data = moose_df, aes(x = x, y = y)) +
  geom_point(aes(colour = total), size = 4) +
  scale_colour_viridis_c() +
  theme_bw()
```

where the grey circles are sites that have not been sampled. 

```{r, results = "hide", fig.keep = "none"}
slmfit_out_moose <- slmfit(formula = total ~ strat, 
  data = moose_df, xcoordcol = 'x', ycoordcol = 'y',
  CorModel = "Exponential")
summary(slmfit_out_moose)
check.variogram(slmfit_out_moose)
qplot(residuals(slmfit_out_moose), bins = 20) + xlab("Residuals")
qplot(residuals(slmfit_out_moose, cross.validation = TRUE), bins = 20) +
  xlab("CV Residuals")

pred_moose <- predict(slmfit_out_moose)
get.predinfo(pred_moose)
get.predplot(pred_moose)
```

We obtain a predicted total of 1596 animals with 90% lower and upper confidence bounds of 921 and 2271 animals, respectively. Unlike the simulation setting, there is no "true total" we can compare our prediction to, because, in reality, not all sites were sampled!

### Sites with Different Areas

Finally, throughout all of the above analyses, we have assumed that the areas of each site were equal. Though this assumption is not accurate for the moose data, due to slightly differing areas based on differing latitudes and longitudes, the assumption approximately holds so that any differences in the prediction that incorporates area is negligible. But, suppose we had sites with very different areas. To showcase how to incorporate site area into the functions in this package, let's first create a "fake" area variable that has the first 700 sites in the region have an area of 1 square kilometer and has the last 160 sites in the region have an area of 2 square kilometeres. 

```{r}
moose_df$fake_area <- c(rep(1, 700), rep(2, 160))
```

For a spatial model, it makes much more sense to use density as the response variable instead of raw counts if the areas of the sites in the model are drastically different. By supplying an `areacol` argument to `slmfit`, the function converts counts to densities, and then gives regression parameters and covariance parameters for the _density_.

```{r}
slmfit_out_moose_area <- slmfit(formula = total ~ strat, 
  data = moose_df, xcoordcol = 'x', ycoordcol = 'y',
  CorModel = "Exponential", areacol = 'fake_area')
summary(slmfit_out_moose_area)
```

The `predict` function then keeps track of the `areacol` argument and gives output in the data frame that pertains to both `counts` and `densities`:

```{r}
pred_obj_area <- predict(slmfit_out_moose_area)
head(pred_obj_area$Pred_df[ ,c("total_pred_density", "total_pred_count",
  "fake_area")])
tail(pred_obj_area$Pred_df[ ,c("total_pred_density", "total_pred_count",
  "fake_area")])
```

Note that, for the first 6 observations, which have an area of 1, the `total_pred_density` and `total_pred_count` columns are identical, while, for the last 6 observations, which have an area of 2, the `total_pred_density` column is half that of the `total_pred_count` column.

The `get.predinfo` function assumes that we are interested in predicting the total count since we did not specify a column of weights so our prediction for the number of moose in this fake region is 1556 moose with a standard error of 393.62.

```{r}
get.predinfo(pred_obj_area)
```

If sites have differing areas, the `get.predplot` function doesn't make much sense to use because each site is represented by the same-sized dot. Here, it would be helpful to import the data frame with the predicted counts and densities into a shapefile so that you are able to construct your own graphics that reflect the different-sized sites.

## Mean Dissolved Organic Carbon from National Lakes Data

As another example, we took data from the [National Aquatic Resource Surveys](https://www.epa.gov/national-aquatic-resource-surveys/data-national-aquatic-resource-surveys). With concerns about global warming, the earth's capacity to store carbon is of great interest, and dissolved organic carbon (DOC) is an estimate of a lake's ability to store carbon.  We will estimate the mean mg/L for DOC from a sample of lakes.  If the total lake volume could be calculated (we will not attempt that), then the total dissolved carbon in a population of lakes could be estimated. We will examine DOC in lakes from the 2012 surveys.  We combined [site data](https://www.epa.gov/sites/production/files/2016-12/nla2012_wide_siteinfo_08232016.csv), [DOC data](https://www.epa.gov/sites/production/files/2016-12/nla2012_waterchem_wide.csv), and [habitat metrics](https://www.epa.gov/sites/production/files/2016-12/nla2012_wide_phabmet_10202016.csv) to create a data set of 1206 lakes in the conterminous United States. 

To access the data, type

```{r}
data(USlakes)
```

and to get a feeling for the data, create a histogram

```{r}
ggplot(data = USlakes, aes(x = log(DOC_RESULT))) + geom_histogram(bins = 20)
```

Even on the log scale, there appears to be some outliers with very high values, and these may be the result of errors in collection or lab analysis.  We will eliminate lakes that have log(DOC) values $>$ 5.

```{r}
lakes <- USlakes[log(USlakes$DOC_RESULT) < 5, ]
```

Our new data set has dimensions

```{r}
nrow(lakes)
```

sites, so we have eliminated 2 sites. To visualize our data some more, we make a bubble plot,

```{r}
plot(USlakes$XCOORD, USlakes$YCOORD, pch = 19, 
  cex = 2 * log(lakes$DOC_RESULT) / max(log(lakes$DOC_RESULT)))
```

and it appears that there is spatial patterning. 

We also have covariates that may help in prediction:

* `ELEVATION`: Elevation at lake coordinates (LAT_DD_N83, LON_DD_N83) from NHD Digital Elevation Map layer
* `RVFPUNDWOODY_RIP`: riparian zone and vegetation: fraction of understory with nonwoody cover present in the riparian zone
* `FCIBIG_LIT`: Fish cover: index of fish cover due to large structures in the littoral zone
* `RVFCGNDBARE_RIP`: riparian zone and vegetation: fraction of ground lacking cover in the riparian zone
* `RVFCGNDWOODY_RIP`: riparian zone and vegetation: fraction of ground cover by woody vegetation in the riparian zone

It is always important to make scatterplots of your variables to look for obvious relationships. 

```{r}
ggplot(data = lakes, aes(x = RVFPUNDWOODY_RIP, y = log(DOC_RESULT))) +
  geom_jitter(width = 0.02)
```

It looks like there might be a negative relationship between riparian nonwoody-understory cover and DOC.

### Creating a Subsample Data Set

We have the whole population of lakes, but, with budget cuts, it is likely that this whole population will not always be surveyed in its entirety.  So, we will ask the question, "If we sample from this population, can we still get a fairly precise estimate of the mean DOC?"

We will do the same thing that we did with the simulated data, and take a random sample of 500 lakes.  Also, because we want the mean, and not a total, we will create a weights column for the `lakeobs` data set, with each element $1/N$, where, here, $N = 1204$.

```{r}
set.seed(2)
LakeObsID <- sample(1:nrow(lakes), 500)
lakeobs <- lakes
lakeobs$DOC_RESULT <- NA
lakeobs[LakeObsID, 'DOC_RESULT'] <- lakes[LakeObsID, 'DOC_RESULT']
lakeobs$wts <- 1 / nrow(lakeobs)
```

### Fitting the Model and Making Predictions

Even though data are skewed, let's try it without taking log of response variable. Note that the mean of log-tranformed  variables is not equal to the log of the mean of set of variables. So if we want a total on the untransformed scale, it would be a mistake to transform the data first, model it, make predictions, sum the predictions, and then exponentiate.  It is much simpler to leave the data untransformed and rely on robustness of the method.  Let's see how this works. 

```{r}
slmfitout_exp_lakes <- slmfit(formula = DOC_RESULT ~ ELEVATION + RVFPUNDWOODY_RIP + 
  FCIBIG_LIT + RVFCGNDBARE_RIP + RVFCGNDWOODY_RIP, data = lakeobs, 
  xcoordcol = 'XCOORD', ycoordcol = 'YCOORD', CorModel = "Exponential")
summary(slmfitout_exp_lakes)
```

We see that all covariates are highly significant. There is substantial autocorrelation because the range parameter is very large, and the partial sill is about six times that of the nugget effect. We fit the model again, but this time with the spherical autocorrelation model.

```{r}
slmfitout_sph_lakes <- slmfit(formula = DOC_RESULT ~ ELEVATION + RVFPUNDWOODY_RIP + 
  FCIBIG_LIT + RVFCGNDBARE_RIP + RVFCGNDWOODY_RIP, data = lakeobs, 
  xcoordcol = 'XCOORD', ycoordcol = 'YCOORD', CorModel = "Spherical")
summary(slmfitout_sph_lakes)
```

We can use AIC to compare the use of the two autocorrelation models.

```{r}
AIC(slmfitout_exp_lakes)
AIC(slmfitout_sph_lakes)
```

Based on AIC, there is not much difference in fit between the two structures. We will use the exponential covariance structure going forward.

```{r}
pred_exp_lakes <- predict(slmfitout_exp_lakes,  wtscol = "wts")
get.predinfo(pred_exp_lakes, conf_level = 0.95)
mean(lakes$DOC_RESULT)
```

We can see that the prediction, 7.975, is close to the true value, 7.65, and that the confidence interval is quite narrow, and it does contain the true value. If a standard error of 0.196, yielding a coefficient of variation of 0.196/7.975 = 0.0245, is acceptable, then sampling 500 lakes could save money and still provide a useful result on DOC.

# Statistical Background

An alternative to a sampling-based appraoch is to assume that the data were generated by a stochastic process and use model-based approaches. It is assumed that the response variable is a realization of a spatial stochastic process. Geostatistical models and methods are used (for a review, see Cressie, 1993). Geostatistics was developed for point samples. If the samples are very small relative to the population size, an infinite population is assumed. In classical geostatistics, the average value over an area can be predicted using methods such as block kriging. Thus it appears that this is closely related to small area estimation, but where samples come from point locations rather than a finite set of sample units. While there is a large literature on geostatistics and block kriging methods, they have been developed for infinite populations. This package is designed for the case where we have a finite collection of plots and we assume that the data were produced by a spatial stochastic process. Detailed developements are given in Ver Hoef (2001, 2008). Comparisons to classical sampling methods can be found in Ver Hoef (2002), and applications in forestry are contained in Ver Hoef and Temesgen (2013) and Temesgen and Ver Hoef (2015).

# Citation

To cite this package, type

```{r, results = "hide"}
citation("sptotal")
```

# References

Cressie, N. 1993. *Statistics for Spatial Data, Revised Edition* John Wiley and Sons, NY.

Temesgen, H. and Ver Hoef, J.M. 2015. Evaluation of the Spatial Linear Model,
Random Forest and Gradient Nearest-Neighbour Methods for Imputing Potential Pro-
ductivity and Biomass of the Pacific Northwest Forests. *Forestry* **88(1)**: 131–142.

Ver Hoef, J.M. 2001. Predicting Finite Populations from Spatially Correlated Data.
*2000 Proceedings of the Section on Statistics and the Environment of the American
Statistical Association*, pgs. 93 – 98.

Ver Hoef, J.M. 2002. Sampling and Geostatistics for Spatial Data. *Ecoscience* **9**:
152–161.

Ver Hoef, J. M. 2008. Spatial Methods for Plot-Based Sampling of Wildlife Populations. *Environmental and Ecological Statistics* **15**: 3-13.

Ver Hoef, J.M. and Temesgen, H. 2013. A Comparison of the Spatial Linear Model
to Nearest Neighbor (k-NN) Methods for Forestry Applications. *PloS ONE* **8(3)**:
e59129.
