In aaarchmiller/uhcplots: Used-Habitat Calibration Plots (UHC Plots)

UHC plots for validating step-selection functions fit to moose location data

In this vignette, we illustrate how UHC plots can be used to validate step-selection functions fit to animal telemetry data. We will use the moose data set from Fieberg et al. (in review), which is also available as part of the uhcplots package. We can access the data and explore the first 6 records as follows:

library(uhcplots)
data("moose12687")
head(moose12687)

The data set consists of 689 GPS locations from 2013 and from 2014 from a single collared moose. The data were subsampled to achieve a consistent 4.25-hour fix rate. A total of 10 available locations were generated for each used location by randomly selecting 10 step lengths and 10 turn angles to project the animal forward in time from the previous location (see Street, G., J. Fieberg, A.R. Rodgers, M. Carstensen, R. Moen, S.A. Moore, S.K. Windels, and J.D. Forester. In press. Habitat functional response mitigates reduced foraging opportunity: implications for animal fitness and space use. Landscape Ecology for a more detailed description). The group of (1 used, 10 available) locaitons at each time step forms a stratum.

We defined resource availability at used and available locations as the proportional cover of 4 land cover types within a 50 m radius buffer: deciduous forest (decid50), mixedwood forest (mixed50), coniferous forest (conif50) and treed wetlands (treedwet50). Lastly, the variable step records the distance between each location from the previous location (after diving through by 1000).

The steps for generating a UHC plot are identical to those from our first example, except that we:

1. Fit a conditional logistic regression model using the clogit function in the survival library rather than a logistic regression model.
2. We take a stratified random sample of locations rather than a simple random sample of locations in step 4c.

We will use the uhcsimstrat function rather than the uhcsim function to generate our predicted distributions of habitat covariates at the used locations in the test data set.

# load the survival library for the clogit function
library(survival)
library(dplyr)

The data set includes observations from 2 years (2013, 2014). We will treat the 2013 data as training data & 2014 as test data.

table(moose12687$year)
mdat.train<-subset(moose12687, year==2013)
mdat.test<-subset(moose12687, year==2014)

Here we explore UHC plots associated with the model containing deciduous forest (decid50), mixedwood forest (mixed50), coniferous forest (conif50) and treed wetlands (treedwet50). We also include step (for a justification,k see Forester et al. (2009) and Avgar et al. (2016)).

# set up 
textplot1 <- (expression(y %~% decid50 + mixed50 + conif50  + treedwet50 + 
                          step + strata(stratum)))
form1a <- (presence ~ decid50 + mixed50 + conif50  + treedwet50 + 
           step + strata(stratum))
form2a <- ~decid50 + mixed50 + conif50 + treedwet50 + step -1

Fit SSF models

We begin by fitting a step-selection function to the data for this animal for the year 2013 (training data).

ssf.train.full <- clogit(form1a, data=mdat.train)
summary(ssf.train.full)

UHC plots

We then:

1. Create simulation envelopes for the environmental characteristics at the observed locations in the test data using the uhcsimstrat function

2. Get density estimates for the environmental characteristics associated with the observed locations in the test data and also associated with the randomly chosen locations generated by the uhcsimstrat function. To do this, we will again use the uhcdenscalc function

design.mat.test.full <- model.matrix(form2a, data=mdat.test)
z <- model.matrix(~decid50 + mixed50 + conif50 + treedwet50 + step -1, 
               data = mdat.test)[,-5]
xchoice.full <- uhcsimstrat(nsims = 1000,
                      xmat = design.mat.test.full, 
                      stratum = mdat.test$stratum, 
                      fit_ssf = ssf.train.full,
                      z = z)    
denshats.decid.full <- uhcdenscalc(rand_sims = xchoice.full[,,1], 
                         dat = z[mdat.test$presence==1,1], 
                         avail = z[mdat.test$presence==0,1]) 
denshats.mixed.full <- uhcdenscalc(rand_sims=xchoice.full[,,2], 
                          dat=z[mdat.test$presence==1,2], 
                          avail=z[mdat.test$presence==0,2])  
denshats.conif.full <- uhcdenscalc(rand_sims=xchoice.full[,,3], 
                          dat=z[mdat.test$presence==1,3], 
                          avail=z[mdat.test$presence==0,3]) 
denshats.treedwt.full <- uhcdenscalc(rand_sims=xchoice.full[,,4], 
                           dat=z[mdat.test$presence==1,4],
                           avail=z[mdat.test$presence==0,4])

1. Then, we create UHC plots again using the uhcdensplot function illustrating:

2. the distribution of covariates at the available points in the test data (red dashed lines)

3. the distribution of covariates at the used points in the test data (solid black lines)
4. the predicted distribution of covariates at the used points in the test data (gray bands)

par(mfrow=c(1,4), mar=c(4,2,2,2), oma=c(3, 4, 7, 0), bty="L")

# First, the density plots
uhcdensplot(densdat = denshats.decid.full$densdat, 
            densrand = denshats.decid.full$densrand, 
            includeAvail = TRUE, 
            densavail = denshats.decid.full$densavail) 
    mtext(outer=F, side=2, line=3, "Density")
    mtext(outer=F, side=3, line=1, "Deciduous", cex=1)
uhcdensplot(densdat = denshats.mixed.full$densdat,
            densrand = denshats.mixed.full$densrand,
            includeAvail = TRUE, 
            densavail = denshats.mixed.full$densavail) 
    mtext(outer=F, side=3, line=1, "Mixedwood", cex=1)
uhcdensplot(densdat = denshats.conif.full$densdat, 
            densrand = denshats.conif.full$densrand,
            includeAvail = TRUE, 
            densavail = denshats.conif.full$densavail) 
    mtext(outer=F, side=3, line=1, "Conifer", cex=1)
uhcdensplot(densdat = denshats.treedwt.full$densdat, 
            densrand = denshats.treedwt.full$densrand, 
            includeAvail = TRUE, 
            densavail = denshats.treedwt.full$densavail) 
    mtext(outer=F, side=3, line=1, "Treed Wetlands", cex=1)
    mtext(outer=T, side=3, line=3,  textplot1, cex=1)

aaarchmiller/uhcplots documentation built on May 10, 2019, 2:05 a.m.