Primary programmer: John Fieberg | Secondary programmer: Althea A. ArchMiller
Date: 20160805
Complete code and analysis to create UHC plots for moose SSFs prepared in support of: Fieberg et al. Species distribution models: Predictive snipers or shots in the dark?
Minnesota moose telemetry data was collected from a single moose in 2013 (training data) and 2014 (test data). See Street et al. (2016) for full description of data development.
Three SSF were fit using these 2013 NLCD cover classes:
SSF Models. (All models included step length divided by 1000 (“step”) to scale the magnitude of the regression coefficient to that of the NLCD classes.)
Set seed
set.seed(0193842)
Load libraries
library(survival) # for running the clogit (standard conditional logit)
library(MASS) # for mvrnorm
library(dplyr) # for data manipulations
library(ezknitr) # for creating html file
library(ggplot2)
Error: package or namespace load failed for 'ggplot2'
library(KernSmooth)
uhcplots functions from GitHub
library(devtools)
install_github("aaarchmiller/uhcplots")
Skipping install of 'uhcplots' from a github remote, the SHA1 (a10e7099) has not changed since last install.
Use `force = TRUE` to force installation
library(uhcplots)
Load and process data
Read in the dataset for moose #12687
mdat <- read.csv("examples/data/moose12687.csv",
stringsAsFactors = F,
header = T)
# or use package's data
mdat <- moose12687
mdat$step <- mdat$step/1000 # divide step lengths by 1000
The data set includes observations from 2 years (2013, 2014). We will treat the 2013 data as training data & 2014 as test data
table(mdat$year)
2013 2014
7579 7579
mdat.train<-subset(mdat, year==2013)
mdat.test<-subset(mdat, year==2014)
#################################################################
Full model
outfig1 <- "Paper/graphics/ssfplota.pdf"
outreg1 <- "examples/Regmods/SSF1_full.R"
panlabs1 <- c("A)","B)", "C)", "D)")
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
Reduced model 1
outfig2 <- "Paper/graphics/ssfplotb.pdf"
outreg2 <- "examples/Regmods/SSF2_reduc.R"
panlabs2 <- c("E)","F)", "G)", "H)")
textplot2 <- (expression(y %~% decid50 + conif50 + treedwet50 +
step + strata(stratum)))
form1b <- (presence ~ decid50 + conif50 + treedwet50 + step + strata(stratum))
form2b <- ~decid50 + conif50 + treedwet50 + step -1
Reduced model 2
outfig3 <- "Paper/graphics/ssfplotc.pdf" # Reduced model 2
outreg3 <- "examples/Regmods/SSF3_reduc.R"
panlabs3 <- c("I)","J)", "K)", "L)")
textplot3 <- (expression(y %~% mixed50 + step + strata(stratum)))
form1c <- (presence ~ mixed50 + step + strata(stratum))
form2c <- ~mixed50 + step -1
#################################################################
Full Model
ssf.train.full <- clogit(form1a, data=mdat.train)
# Save regression models (to make a table later)
save(ssf.train.full, file=outreg1)
ssf.train.full
Call:
clogit(form1a, data = mdat.train)
coef exp(coef) se(coef) z p
decid50 0.489 1.631 0.329 1.49 0.14
mixed50 1.382 3.982 0.328 4.21 2.5e-05
conif50 -0.296 0.744 0.377 -0.78 0.43
treedwet50 0.396 1.486 0.313 1.26 0.21
step -6325.695 0.000 247.919 -25.52 < 2e-16
Likelihood ratio test=1427 on 5 df, p=0
n= 7579, number of events= 689
Reduced model 1
ssf.train.reduc1 <- clogit(form1b, data=mdat.train)
# Save regression models (to make a table later)
save(ssf.train.reduc1, file=outreg2)
ssf.train.reduc1
Call:
clogit(form1b, data = mdat.train)
coef exp(coef) se(coef) z p
decid50 -0.598 0.550 0.189 -3.16 0.0016
conif50 -1.365 0.255 0.266 -5.13 2.9e-07
treedwet50 -0.696 0.499 0.160 -4.35 1.4e-05
step -6435.575 0.000 248.037 -25.95 < 2e-16
Likelihood ratio test=1407 on 4 df, p=0
n= 7579, number of events= 689
Reduced model 2
ssf.train.reduc2 <- clogit(form1c, data=mdat.train)
# Save regression models (to make a table later)
save(ssf.train.reduc2, file=outreg3)
ssf.train.reduc2
Call:
clogit(form1c, data = mdat.train)
coef exp(coef) se(coef) z p
mixed50 1.026 2.791 0.158 6.5 7.9e-11
step -6384.514 0.000 248.190 -25.7 < 2e-16
Likelihood ratio test=1414 on 2 df, p=0
n= 7579, number of events= 689
#################################################################
Full Model
[1]
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)
[2]
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])
Reduced Model 1
[1]
design.mat.test.reduc1 <- model.matrix(form2b, data=mdat.test)
z <- model.matrix(~decid50 + mixed50 + conif50 + treedwet50 + step -1,
data = mdat.test)[,-5]
xchoice.reduc1 <- uhcsimstrat(nsims = 1000,
xmat = design.mat.test.reduc1,
stratum = mdat.test$stratum,
fit_ssf = ssf.train.reduc1,
z = z)
[2]
denshats.decid.reduc1 <- uhcdenscalc(rand_sims = xchoice.reduc1[,,1],
dat = z[mdat.test$presence==1,1],
avail = z[mdat.test$presence==0,1])
denshats.mixed.reduc1 <- uhcdenscalc(rand_sims=xchoice.reduc1[,,2],
dat=z[mdat.test$presence==1,2],
avail=z[mdat.test$presence==0,2])
denshats.conif.reduc1 <- uhcdenscalc(rand_sims=xchoice.reduc1[,,3],
dat=z[mdat.test$presence==1,3],
avail=z[mdat.test$presence==0,3])
denshats.treedwt.reduc1 <- uhcdenscalc(rand_sims=xchoice.reduc1[,,4],
dat=z[mdat.test$presence==1,4],
avail=z[mdat.test$presence==0,4])
Reduced Model 2
[1]
design.mat.test.reduc2 <- model.matrix(form2c, data=mdat.test)
z <- model.matrix(~decid50 + mixed50 + conif50 + treedwet50 + step -1,
data = mdat.test)[,-5]
xchoice.reduc2 <- uhcsimstrat(nsims = 1000,
xmat = design.mat.test.reduc2,
stratum = mdat.test$stratum,
fit_ssf = ssf.train.reduc2,
z = z)
[2]
denshats.decid.reduc2 <- uhcdenscalc(rand_sims = xchoice.reduc2[,,1],
dat = z[mdat.test$presence==1,1],
avail = z[mdat.test$presence==0,1])
denshats.mixed.reduc2 <- uhcdenscalc(rand_sims=xchoice.reduc2[,,2],
dat=z[mdat.test$presence==1,2],
avail=z[mdat.test$presence==0,2])
denshats.conif.reduc2 <- uhcdenscalc(rand_sims=xchoice.reduc2[,,3],
dat=z[mdat.test$presence==1,3],
avail=z[mdat.test$presence==0,3])
denshats.treedwt.reduc2 <- uhcdenscalc(rand_sims=xchoice.reduc2[,,4],
dat=z[mdat.test$presence==1,4],
avail=z[mdat.test$presence==0,4])
Uses uhcdensplot function
Full model
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(side=3, line=1, panlabs1[1], cex=1.2, ad=0)
mtext(outer=F, side=2, line=3, "Density")
mtext(outer=F, side=3, line=1, "Deciduous", cex=1.4)
legend(0.3, 10.5, c("Available", "Used", "Predicted"),
lty=c(1, 2, 1), lwd=c(2,3, 10),
col=c("Black", "red", "gray"), bty="n", cex=1.8)
uhcdensplot(densdat = denshats.mixed.full$densdat,
densrand = denshats.mixed.full$densrand,
includeAvail = TRUE,
densavail = denshats.mixed.full$densavail)
mtext(side=3, line=1, panlabs1[2], cex=1.2, ad=0)
mtext(outer=F, side=3, line=1, "Mixedwood", cex=1.4)
uhcdensplot(densdat = denshats.conif.full$densdat,
densrand = denshats.conif.full$densrand,
includeAvail = TRUE,
densavail = denshats.conif.full$densavail)
mtext(side=3, line=1, panlabs1[3], cex=1.2, ad=0)
mtext(outer=F, side=3, line=1, "Conifer", cex=1.4)
uhcdensplot(densdat = denshats.treedwt.full$densdat,
densrand = denshats.treedwt.full$densrand,
includeAvail = TRUE,
densavail = denshats.treedwt.full$densavail)
mtext(side=3, line=1, panlabs1[4], cex=1.2, ad=0)
mtext(outer=F, side=3, line=1, "Treed Wetlands", cex=1.4)
mtext(outer=T, side=3, line=3, textplot1, cex=1.8)
plot.uhc.full <- recordPlot()
quartz_off_screen
2
Reduced model 1
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.reduc1$densdat,
densrand = denshats.decid.reduc1$densrand,
includeAvail = TRUE,
densavail = denshats.decid.reduc1$densavail)
mtext(side=3, line=1, panlabs2[1], cex=1.2, ad=0)
mtext(outer=F, side=2, line=3, "Density")
mtext(outer=F, side=3, line=1, "Deciduous", cex=1.4)
legend(0.3, 10.5, c("Available", "Used", "Predicted"),
lty=c(1, 2, 1), lwd=c(2,3, 10),
col=c("Black", "red", "gray"), bty="n", cex=1.8)
uhcdensplot(densdat = denshats.mixed.reduc1$densdat,
densrand = denshats.mixed.reduc1$densrand,
includeAvail = TRUE,
densavail = denshats.mixed.reduc1$densavail)
mtext(side=3, line=1, panlabs2[2], cex=1.2, ad=0)
mtext(outer=F, side=3, line=1, "Mixedwood", cex=1.4)
uhcdensplot(densdat = denshats.conif.reduc1$densdat,
densrand = denshats.conif.reduc1$densrand,
includeAvail = TRUE,
densavail = denshats.conif.reduc1$densavail)
mtext(side=3, line=1, panlabs2[3], cex=1.2, ad=0)
mtext(outer=F, side=3, line=1, "Conifer", cex=1.4)
uhcdensplot(densdat = denshats.treedwt.reduc1$densdat,
densrand = denshats.treedwt.reduc1$densrand,
includeAvail = TRUE,
densavail = denshats.treedwt.reduc1$densavail)
mtext(side=3, line=1, panlabs2[4], cex=1.2, ad=0)
mtext(outer=F, side=3, line=1, "Treed Wetlands", cex=1.4)
mtext(outer=T, side=3, line=3, textplot2, cex=1.8)
plot.uhc.reduc1 <- recordPlot()
quartz_off_screen
2
Reduced model 2
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.reduc2$densdat,
densrand = denshats.decid.reduc2$densrand,
includeAvail = TRUE,
densavail = denshats.decid.reduc2$densavail)
mtext(side=3, line=1, panlabs3[1], cex=1.2, ad=0)
mtext(outer=F, side=2, line=3, "Density")
mtext(outer=F, side=3, line=1, "Deciduous", cex=1.4)
legend(0.3, 10.5, c("Available", "Used", "Predicted"),
lty=c(1, 2, 1), lwd=c(2,3, 10),
col=c("Black", "red", "gray"), bty="n", cex=1.8)
uhcdensplot(densdat = denshats.mixed.reduc2$densdat,
densrand = denshats.mixed.reduc2$densrand,
includeAvail = TRUE,
densavail = denshats.mixed.reduc2$densavail)
mtext(side=3, line=1, panlabs3[2], cex=1.2, ad=0)
mtext(outer=F, side=3, line=1, "Mixedwood", cex=1.4)
uhcdensplot(densdat = denshats.conif.reduc2$densdat,
densrand = denshats.conif.reduc2$densrand,
includeAvail = TRUE,
densavail = denshats.conif.reduc2$densavail)
mtext(side=3, line=1, panlabs3[3], cex=1.2, ad=0)
mtext(outer=F, side=3, line=1, "Conifer", cex=1.4)
uhcdensplot(densdat = denshats.treedwt.reduc2$densdat,
densrand = denshats.treedwt.reduc2$densrand,
includeAvail = TRUE,
densavail = denshats.treedwt.reduc2$densavail)
mtext(side=3, line=1, panlabs3[4], cex=1.2, ad=0)
mtext(outer=F, side=3, line=1, "Treed Wetlands", cex=1.4)
mtext(outer=T, side=3, line=3, textplot3, cex=1.8)
plot.uhc.reduc2 <- recordPlot()
quartz_off_screen
2
#################################################################
\(f^u(z)-\hat{f}^u(z)_i\)
Uses uhcdiffdensplot function
Full Model
par(mfrow=c(3,4), mar=c(2,2,2,2), oma=c(3, 0, 5, 0), bty="L")
uhcdiffdensplot(densdat = denshats.decid.full$densdat,
densrand = denshats.decid.full$densrand)
mtext(side=3, line=1, panlabs1[1], cex=1, ad=0)
mtext(outer=F, side=2, line=3, "Observed-Predicted")
uhcdiffdensplot(densdat = denshats.mixed.full$densdat,
densrand = denshats.mixed.full$densrand)
mtext(side=3, line=1, panlabs1[2], cex=1, ad=0)
uhcdiffdensplot(densdat = denshats.conif.full$densdat,
densrand = denshats.conif.full$densrand)
mtext(side=3, line=1, panlabs1[3], cex=1, ad=0)
uhcdiffdensplot(densdat = denshats.treedwt.full$densdat,
densrand = denshats.treedwt.full$densrand)
mtext(side=3, line=1, panlabs1[4], cex=1, ad=0)
# **Reduced Model 1**
uhcdiffdensplot(densdat = denshats.decid.reduc1$densdat,
densrand = denshats.decid.reduc1$densrand)
mtext(side=3, line=1, panlabs2[1], cex=1, ad=0)
mtext(outer=F, side=2, line=3, "Observed-Predicted")
uhcdiffdensplot(densdat = denshats.mixed.reduc1$densdat,
densrand = denshats.mixed.reduc1$densrand)
mtext(side=3, line=1, panlabs2[2], cex=1, ad=0)
uhcdiffdensplot(densdat = denshats.conif.reduc1$densdat,
densrand = denshats.conif.reduc1$densrand)
mtext(side=3, line=1, panlabs2[3], cex=1, ad=0)
uhcdiffdensplot(densdat = denshats.treedwt.reduc1$densdat,
densrand = denshats.treedwt.reduc1$densrand)
mtext(side=3, line=1, panlabs2[4], cex=1, ad=0)
# **Reduced Model 2**
uhcdiffdensplot(densdat = denshats.decid.reduc2$densdat,
densrand = denshats.decid.reduc2$densrand)
mtext(side=3, line=1, panlabs3[1], cex=1, ad=0)
mtext(outer=F, side=2, line=3, "Observed-Predicted")
uhcdiffdensplot(densdat = denshats.mixed.reduc2$densdat,
densrand = denshats.mixed.reduc2$densrand)
mtext(side=3, line=1, panlabs3[2], cex=1, ad=0)
uhcdiffdensplot(densdat = denshats.conif.reduc2$densdat,
densrand = denshats.conif.reduc2$densrand)
mtext(side=3, line=1, panlabs3[3], cex=1, ad=0)
uhcdiffdensplot(densdat = denshats.treedwt.reduc2$densdat,
densrand = denshats.treedwt.reduc2$densrand)
mtext(side=3, line=1, panlabs3[4], cex=1, ad=0)
References:
Street, G.M., Fieberg, J., Rodgers, A.R., Carstensen, M., Moen, R., Moore, S.A., et al. (2016). Habitat functional response mitigates reduced foraging opportunity: Implications for animal fitness and space use. Landscape Ecology, 1-15.
Session Information
sessionInfo()
R version 3.3.0 (2016-05-03)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.11.6 (El Capitan)
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] dplyr_0.5.0 MASS_7.3-45 survival_2.39-5
[4] uhcplots_0.1.0 devtools_1.12.0 SDMTools_1.1-221
[7] ENMeval_0.2.1 rJava_0.9-8 dismo_1.1-1
[10] raster_2.5-8 sp_1.2-3 knitr_1.14
[13] KernSmooth_2.23-15 ezknitr_0.6
loaded via a namespace (and not attached):
[1] Rcpp_0.12.7 RColorBrewer_1.1-2 formatR_1.4
[4] git2r_0.15.0 plyr_1.8.4 R.methodsS3_1.7.1
[7] R.utils_2.4.0 iterators_1.0.8 tools_3.3.0
[10] digest_0.6.10 tibble_1.2 evaluate_0.10
[13] memoise_1.0.0 gtable_0.2.0 lattice_0.20-34
[16] Matrix_1.2-7.1 foreach_1.4.3 DBI_0.5-1
[19] curl_2.2 withr_1.0.2 stringr_1.1.0
[22] httr_1.2.1 grid_3.3.0 R6_2.2.0
[25] magrittr_1.5 splines_3.3.0 scales_0.4.1
[28] codetools_0.2-15 assertthat_0.1 mime_0.5
[31] colorspace_1.2-7 stringi_1.1.2 lazyeval_0.2.0
[34] doParallel_1.0.10 munsell_0.4.3 markdown_0.7.7
[37] R.oo_1.20.0
Spun with: ezspin(file = “examples/example_3_moose.R”, out_dir =“examples/output”, fig_dir = “figures”, keep_md = FALSE)