Here is a “solution” for practical 4. As with any data analysis, there is no correct answer, but this shows how I would approach this analysis. The analysis here is conditional on selecting a detection function and DSM in the previous exercises; I’ve shown a variety of models selected in the previous solutions to show the differences between models.

Much of the text below is as in the exercise itself, so it should be relatively easy to navigate.

Additional text and code is highlighted using boxes like this.

Now we’ve fitted some models, let’s use the predict functions and the data from GIS to make predictions of abundance.

Loading the packages and data

## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-7. For overview type 'help("mgcv-package")'.
## Loading required package: mrds
## This is mrds 2.1.14
## Built: R 3.2.0; ; 2015-07-30 10:07:19 UTC; unix
## This is dsm 2.2.11
## Built: R 3.2.2; ; 2015-10-23 20:20:41 UTC; unix
# colourblind-friendly colourschemes
# to load and save raster data
## Loading required package: sp
## Warning: no function found corresponding to methods exports from 'raster'
## for: 'overlay'
## Attaching package: 'raster'
## The following object is masked from 'package:nlme':
##     getData
# models with only spatial terms
# models with all covariates

Loading prediction data

Before we can make predictions we first need to load the covariates into a “stack” from their files on disk using the stack() function from raster. We give stack() a vector of locations to load the rasters from. Note that in RStudio you can use tab-completion for these locations and avoid some typing.At this point we arbitrarily choose the time periods of the SST, NPP and EKE rasters (2 June 2004, or Julian date 153).

predictorStack <- stack(c("Covariates_for_Study_Area/Depth.img",
                          "Covariates_for_Study_Area/Global/DT\ all\ sat/MSLA_ke/2004/MSLA_ke_2004153.img"

We need to rename the layers in our stack to match those in the model we are going to use to predict. If you need a refresher on the names that were used there, call summary() on the DSM object.

names(predictorStack) <- c("Depth","SST","NPP", "DistToCAS", "EKE")

Now these are loaded, we can coerce the stack into something dsm can talk to using the function. Note we need the xy=TRUE to ensure that x and y are included in the prediction data. We also set the offset value – the area of each cell in our prediction grid.

predgrid <-, xy=TRUE)
predgrid$off.set <- (10*1000)^2

We can then predict for the model dsm_nb_xy_ms:

pp <- predict(dsm_nb_xy_ms, predgrid)

This is just a list of numbers – the predicted abundance per cell. We can sum these to get the estimated abundance for the study area:

sum(pp, na.rm=TRUE)
## [1] 1589.216

Because we predicted over the whole raster grid (including those cells without covariate values – e.g. land), some of the values in pp will be NA, so we can ignore them when we sum by setting na.rm=TRUE. We need to do this again when we plot the data too.

We can also plot this to get a spatial representation of the predictions:

# assign the predictions to the prediction grid data.frame
predgrid$Nhat_nb_xy <- pp
# remove the NA entries (because of the grid structure of the raster)
predgrid_plot <- predgrid[!$Depth),]
# plot!
p <- ggplot(predgrid_plot) +
      geom_tile(aes(x=x, y=y, fill=Nhat_nb_xy, width=10*1000, height=10*1000)) +
      coord_equal() + 

Copy the chunk above and make predictions for the other models you saved in the previous exercises. In particular, compare the models with only spatial terms to those with environmental covariates included.

We now want to plot the predictions for the 8 models that we have so far (4 negative binomial, 4 Tweedie; 4 spatial only, 4 with environmental covariats; 4 with bivariate smooths, 4 with additive spatial effects).

Duplicating the above code is a bit tiresome and can be prone to errors, so let’s nerd-out pretty heavily for this solution and show how we can use the ldply() function from plyr to do the same task many times.

# make a function that makes the predictions, adds them to a column named Nhat 
# and adds a column called "model" that stores the model name, then returns the
# data.frame.
make_pred_dat <- function(model_name, predgrid){
  # we use get() here to grab the object with the name of its argument
  predgrid[["Nhat"]] <- predict(get(model_name), predgrid)
  predgrid[["model"]] <- model_name

# load plyr and apply to a list of the names of the models, make_pred_dat returns
# a data.frame (hence this is an "ld" function: list->data.frame) that it then binds
# together
big_predgrid <- ldply(list("dsm_nb_xy", "dsm_nb_x_y", "dsm_nb_xy_ms", "dsm_nb_x_y_ms",
                           "dsm_tw_xy", "dsm_tw_x_y", "dsm_tw_xy_ms", "dsm_tw_x_y_ms"),
                       make_pred_dat, predgrid=predgrid_plot)

# make the plot, facetting using the model column
p <- ggplot(big_predgrid) +
      geom_tile(aes(x=x, y=y, fill=Nhat, width=10*1000, height=10*1000)) +
      coord_equal() + 
      facet_wrap(~model, nrow=2)+

Note here that the _ms models have the environmental covariates, the others are spatial-only.

We can also use plyr to help calculate overall abundance…

# ddply