# Multi species single season with a list of models
# Co-occurence of Grizzly and black bears as measured by camera
#
# Data from
# Ladle A, Steenweg R, Shepherd B, Boyce MS (2018)
# The role of human outdoor recreation in shaping patterns of grizzly bear-black
# bear co-occurrence.
# PLOS ONE 13(2): e0191730. https://doi.org/10.1371/journal.pone.0191730
# Ladle A, Steenweg R, Shepherd B, Boyce MS (2018)
# Data from: The role of human outdoor recreation in shaping patterns of
# grizzly bear-black bear co-occurrence.
# Dryad Digital Repository. https://doi.org/10.5061/dryad.s81k3
# Briefly,
# Cameras placed on trails in summer 2014
# Classified images as type of bear and type of human activity.
# 2020-06-24 CJS real$rho and real$nu now in derived$
# 2018-12-26 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
library(ggplot2)
library(readxl)
library(reshape2)
library(RPresence)
# get the data. Sites x # of detections grouped into 4 day intervals
gb.data <- read.csv(file.path("..","Raw data","four_day_grizzly.csv"), header=FALSE, strip.white=TRUE, as.is=TRUE)
bb.data <- read.csv(file.path("..","Raw data","four_day_black.csv"), header=FALSE, strip.white=TRUE, as.is=TRUE)
# Notice that in some days, more than 1 detection recorded.
range(gb.data, na.rm=TRUE)
## [1] 0 6
apply(gb.data, 1, max, na.rm=TRUE)
## [1] 2 0 0 6 0 0 0 0 0 1 1 1 0 3 1 1 1 0 0 1 0 1 0 3 1 0 0 2 0 0 0 0 0 0 0 0 0
## [38] 1 0 1 0 1 0 1 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 2
## [75] 0 3 0 3 0 1 2 1 0 2 1 4 1 0 1 0 0 0 0 0 0 0 1 1 1 0 0 1 0 1 2 0 2 2 2 1 0
## [112] 0 1 0 1 0 0 1 4 2 1 1 1 0 0 0 1 0 0 0 1 1 1 0 2 2 2 1 2 0 2 1 1 0 1 1 0 1
## [149] 0 0 0 0 0 1 1 1 2 0 0 0 1 1 1 2 1 0 1 0 1 0 2 2 2 0 0 0 0 0 0 1 2 0
# Convert history to compressed format
# 0=neither species detected,
# 1=only species A detected,
# 2=only species B detected,
# 3=both species detected
# we need to ignore multple detections of each species in a camera
input.history <- 1*(bb.data>=1) + 2*(gb.data>=1)
input.history
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18
## [1,] 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [4,] 2 2 2 2 2 2 2 2 NA NA NA NA NA NA NA NA NA NA
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [8,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [9,] 0 0 0 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA
## [10,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0
## [11,] 2 0 0 0 0 2 0 0 2 0 0 0 0 0 0 0 0 0
## [12,] 0 2 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0
## [13,] 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [14,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0
## [15,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0
## [16,] 0 2 0 0 2 0 0 NA NA NA NA NA NA NA NA NA NA NA
## [17,] 0 0 0 2 0 2 0 0 0 0 0 0 0 0 0 0 0 0
## [18,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [19,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [20,] 1 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0
## [21,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [22,] 2 0 0 0 2 0 0 0 0 0 0 0 2 0 0 0 0 0
## [23,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [24,] 2 0 0 0 0 2 0 0 0 0 2 2 0 0 0 0 0 0
## [25,] 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [26,] 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## [27,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [28,] 0 0 2 0 1 0 0 0 0 0 0 0 0 2 2 0 NA NA
## [29,] 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA
## [30,] 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA
## [31,] 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA
## [32,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [33,] 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [34,] 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [35,] 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA
## [36,] 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA
## [37,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [38,] 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0
## [39,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [40,] 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 1 0
## [41,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [42,] 0 0 2 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0
## [43,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [44,] 1 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 1 0
## [45,] 0 1 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0
## [46,] 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [47,] 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## [48,] 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0
## [49,] 0 1 1 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA
## [50,] 0 0 0 0 2 1 0 0 0 0 0 0 0 0 0 0 0 0
## [51,] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [52,] 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [53,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [54,] 0 0 0 0 0 0 1 0 0 2 0 0 0 0 0 0 0 0
## [55,] 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [56,] 0 0 0 1 0 0 0 0 0 0 0 2 0 0 0 0 0 0
## [57,] 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## [58,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [59,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [60,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [61,] 0 3 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0
## [62,] 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1
## [63,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [64,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [65,] 0 2 0 0 2 2 0 0 0 2 0 2 0 0 NA NA NA NA
## [66,] 0 0 0 1 0 0 0 0 0 0 NA NA NA NA NA NA NA NA
## [67,] 0 1 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA
## [68,] 0 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA
## [69,] 1 1 1 0 1 0 1 0 0 0 0 1 1 1 1 0 0 1
## [70,] 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
## [71,] 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
## [72,] 0 0 0 0 0 0 0 0 0 0 0 0 NA NA NA NA NA NA
## [73,] 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0
## [74,] 2 0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [75,] 0 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 1 0
## [76,] 3 2 2 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0
## [77,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
## [78,] 2 1 2 2 2 0 0 0 0 0 1 0 0 0 NA NA NA NA
## [79,] 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## [80,] 0 2 0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [81,] 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 NA NA
## [82,] 1 0 1 0 1 1 0 0 0 0 0 0 0 2 0 NA NA NA
## [83,] 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [84,] 3 3 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## [85,] 2 0 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA
## [86,] 2 0 2 2 0 0 0 0 2 0 0 2 0 0 2 2 2 3
## [87,] 1 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0
## [88,] 0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0 NA NA
## [89,] 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [90,] 0 0 0 0 1 0 0 0 0 1 NA NA NA NA NA NA NA NA
## [91,] 0 1 1 0 0 0 1 NA NA NA NA NA NA NA NA NA NA NA
## [92,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [93,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 NA
## [94,] 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [95,] 0 0 0 0 0 0 0 0 0 NA NA NA NA NA NA NA NA NA
## [96,] 1 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0
## [97,] 0 0 0 0 0 0 0 0 0 2 NA NA NA NA NA NA NA NA
## [98,] 2 2 0 0 2 0 0 NA NA NA NA NA NA NA NA NA NA NA
## [99,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 NA NA
## [100,] 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [101,] 0 0 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA
## [102,] 2 0 0 0 0 2 2 2 NA NA NA NA NA NA NA NA NA NA
## [103,] 1 1 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 1
## [104,] 2 0 0 1 0 0 NA NA NA NA NA NA NA NA NA NA NA NA
## [105,] 2 1 0 0 0 0 0 0 2 0 2 0 2 0 NA NA NA NA
## [106,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [107,] 2 0 2 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [108,] 2 0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [109,] 2 2 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [110,] 0 0 0 0 2 0 0 0 0 0 0 0 0 0 NA NA NA NA
## [111,] 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 1 0
## [112,] 0 0 0 0 0 0 0 0 0 0 0 0 NA NA NA NA NA NA
## [113,] 0 0 2 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA
## [114,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [115,] 0 1 0 1 0 0 0 0 0 2 0 0 NA NA NA NA NA NA
## [116,] 1 1 1 1 1 1 1 0 0 0 1 1 0 NA NA NA NA NA
## [117,] 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [118,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0
## [119,] 3 1 2 2 2 0 0 2 2 0 2 0 0 2 0 0 0 0
## [120,] 0 0 2 0 2 0 0 2 0 NA NA NA NA NA NA NA NA NA
## [121,] 1 0 0 0 2 3 0 NA NA NA NA NA NA NA NA NA NA NA
## [122,] 2 2 0 0 0 0 0 0 0 0 0 0 0 0 NA NA NA NA
## [123,] 2 0 0 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA
## [124,] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 NA NA
## [125,] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 NA NA
## [126,] 0 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA
## [127,] 0 1 2 2 1 1 1 0 0 0 1 0 0 0 0 0 NA NA
## [128,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [129,] 0 0 0 0 0 0 0 0 0 NA NA NA NA NA NA NA NA NA
## [130,] 0 0 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA
## [131,] 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 NA NA NA
## [132,] 0 0 0 0 2 0 0 0 0 1 0 0 0 0 2 0 2 0
## [133,] 0 0 0 0 0 0 0 0 2 0 0 0 NA NA NA NA NA NA
## [134,] 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [135,] 0 0 2 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0
## [136,] 2 2 0 0 0 2 0 0 0 0 2 0 0 0 0 0 0 0
## [137,] 2 2 2 0 0 0 2 0 0 0 2 2 0 0 0 0 0 0
## [138,] 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0
## [139,] 2 2 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0
## [140,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [141,] 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [142,] 0 2 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## [143,] 1 2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## [144,] 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0
## [145,] 0 2 2 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## [146,] 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0
## [147,] 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## [148,] 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0
## [149,] 0 0 0 0 0 0 0 0 0 0 0 0 0 NA NA NA NA NA
## [150,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [151,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [152,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [153,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [154,] 0 0 2 0 0 0 1 1 0 0 0 0 NA NA NA NA NA NA
## [155,] 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [156,] 1 1 2 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0
## [157,] 2 2 1 0 0 0 1 0 0 2 2 2 0 0 0 0 0 0
## [158,] 0 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA
## [159,] 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## [160,] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [161,] 1 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0
## [162,] 0 1 0 0 1 0 0 0 0 0 0 2 0 0 1 0 0 0
## [163,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0
## [164,] 0 2 0 2 0 0 0 0 0 0 0 0 0 0 2 0 0 0
## [165,] 1 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0
## [166,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [167,] 1 1 2 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
## [168,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [169,] 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2
## [170,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
## [171,] 0 2 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [172,] 2 0 0 0 0 0 NA NA NA NA NA NA NA NA NA NA NA NA
## [173,] 2 2 2 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0
## [174,] 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## [175,] 0 1 1 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0
## [176,] 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## [177,] 0 0 0 1 0 0 1 1 1 1 1 1 0 0 0 0 0 0
## [178,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [179,] 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 0
## [180,] 0 0 2 0 0 NA NA NA NA NA NA NA NA NA NA NA NA NA
## [181,] 2 0 2 0 0 0 0 0 0 0 0 0 NA NA NA NA NA NA
## [182,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Get the site covariates
site.covar <- read.csv(file.path("..","Raw data","covariates.csv"), header=TRUE, strip.white=TRUE, as.is=TRUE)
names(site.covar)
## [1] "Cam_name" "X" "Y" "Land_type" "droad" "dstream"
## [7] "elev" "tpi" "NDVI"
names(site.covar) <- make.names(names(site.covar))
names(site.covar)
## [1] "Cam_name" "X" "Y" "Land_type" "droad" "dstream"
## [7] "elev" "tpi" "NDVI"
head(site.covar)
## Cam_name X Y Land_type droad
## 1 amethyst lake 415655 5838411 415655 5838411 Jasper 12759.7998
## 2 ancient wall 376992 5919664 376992 5919664 Jasper 47406.6016
## 3 astoria 427308 5839183 427308 5839183 Jasper 1134.0400
## 4 bess pass 343359 5911520 343359 5911520 Crown 76867.2031
## 5 brewster skywalk c2 481109 5789822 481109 5789822 Jasper 44.7214
## 6 brewster skywalk c2b 481109 5789822 481107 5789822 Jasper 49.2443
## dstream elev tpi NDVI
## 1 55.0000 1994 -212.058 7102.968
## 2 53.1507 1760 -373.933 6576.896
## 3 640.3120 1722 -259.626 6440.972
## 4 530.0000 1695 0.000 6789.304
## 5 101.2420 1825 -608.371 5605.156
## 6 97.0824 1825 -608.371 5605.156
# Standardize Elevation covariates
elev.mean <- mean(site.covar$elev)
elev.s <- sd (site.covar$elev)
site.covar$elev.std <- (site.covar$elev - elev.mean)/elev.s
Nsites <- nrow(input.history)
Nvisits<- ncol(input.history)
# create the pao file
bear.pao <- createPao(input.history,
unitcov=site.covar,
title="bear multi species - co-occurance")
summary(bear.pao)
## paoname=pres.pao
## title=bear multi species - co-occurance
## Naive occ=0.6813187
## naiveR =0.4450549
## nunits nsurveys nseasons nsurveyseason nmethods
## "182" "18" "1" "18" "1"
## nunitcov nsurvcov
## "10" "1"
## unit covariates : Cam_name X Y Land_type droad dstream elev tpi NDVI elev.std
## survey covariates: SURVEY
# Use the psiAB parameterization
# Parameters of the model are
# psiA - occupancy probability of species A
# psiBA - occupancy probability of species B if species A is present
# psiBa - occupancy probability of species B if species A is absent
#
# If species are independent thatn psiBA = psiBa.
# Alternatively, nu = odds(B|A)/odds(B|a) = 1.
#
# Detection parameters
# pA - probability of detection if A is alone in the site
# pB - probability of detection if B is alone in the site
# rA - probability of detecting A only given both on the site
# rBA - probability of detecting B given that A was detected and both are on the site
# rBa - probability of detecting B given that A not detected and both are on the site
# Ifspecies do not interact, then
# rBA = rBa
# and
# rho = odds(rBA)/odds(rBa) = 1
#
# The following special variables are available for modelling PSI
# SP species effect
# INT interaction of effect on presence of species B when species A was, or was not present
#
# Model for PSI.... impact on parameters
# psi~1 psiA=psiBA=psiBa (1 parameter)
# psi~SP psiA psiBA=psiBa (2 parameters)
# psi~SP+INT psiA psiBA psiBa (3 parameters)
#
# The following special variables are available for p
# SP species effect
# INT_o is a detection-level interaction where the OCCUPANCY of one species
# changes the detection probability of the other species
# INT_d is a detection-level interaction where DETECTION of one species changes the
# detection probability of the other species in the same survey.
#
# Model for p.... impact on parameters
# p~1 pA = pB = rA = rBA = rBa (1 parameter)
# p~SP pA=rA, pB=rBA=rBa (2 parameters)
# p~SP+INT_o pA=rA, pB, rBA=rBa (3 parameters)
# p~SP+INT_o+SP:INT_o pA, rA, pB, rBA=rBa (4 parameters)
# p~SP+INT_o+SP:INT_o+INT_d pA, rA, pB, rBA, rBa (5 parameters)
# define the list of models to fit
# Notice the commas between the column and the placement of the quotes
model.list.csv <- textConnection("
p, psi
~SP+INT_o+SP:INT_o + INT_d, ~SP+INT
~SP+INT_o+SP:INT_o + INT_d, ~SP
~SP+INT_o+SP:INT_o , ~SP+INT
~SP+INT_o+SP:INT_o, ~SP + elev.std+SP:elev.std+INT+INT:elev.std")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list
## p psi
## 1 ~SP+INT_o+SP:INT_o + INT_d ~SP+INT
## 2 ~SP+INT_o+SP:INT_o + INT_d ~SP
## 3 ~SP+INT_o+SP:INT_o ~SP+INT
## 4 ~SP+INT_o+SP:INT_o ~SP + elev.std+SP:elev.std+INT+INT:elev.std
# fit the models
model.fits <- plyr::alply(model.list, 1, function(x,detect.pao){
cat("\n\n***** Starting ", unlist(x), "\n")
fit <- RPresence::occMod(model=list(as.formula(paste("psi",x$psi)),
as.formula(paste("p" ,x$p ))),
data=detect.pao,type="so.2sp.1", randint=10)
fit
},detect.pao=bear.pao)
##
##
## ***** Starting ~SP+INT_o+SP:INT_o + INT_d ~SP+INT
##
##
## ***** Starting ~SP+INT_o+SP:INT_o + INT_d ~SP
##
##
## ***** Starting ~SP+INT_o+SP:INT_o ~SP+INT
##
##
## ***** Starting ~SP+INT_o+SP:INT_o ~SP + elev.std+SP:elev.std+INT+INT:elev.std
# Look the output from a specific model
check.model <- 4
names(model.fits[[check.model]])
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "npar" "aic" "beta" "real"
## [11] "derived" "warnings"
model.fits[[check.model]]$beta
## $psi
## est se
## A1_psiA 0.748239 0.281195
## A2_psiBA 0.145727 0.746871
## A3_psiA.elev.std_psiA -1.000606 0.279150
## A4_psiBa -0.597586 1.747908
## A5_psiBA.SP2:elev.std_psiBA -0.970599 0.817819
## A6_psiBa.elev.std:INT2_psiBa 10.246286 5.942839
##
## $psi.VC
## A1 A2 A3 A4 A5 A6
## A1 0.079070 -0.145933 -0.024289 0.178608 0.082469 0.329030
## A2 -0.145933 0.557816 0.080219 -0.922256 -0.342742 -0.342424
## A3 -0.024289 0.080219 0.077925 -0.127854 -0.142763 0.033377
## A4 0.178608 -0.922256 -0.127854 3.055182 0.564698 1.474631
## A5 0.082469 -0.342742 -0.142763 0.564698 0.668827 -0.696643
## A6 0.329030 -0.342424 0.033377 1.474631 -0.696643 35.317330
##
## $p
## est se
## B1_pA[1] -4.772841 0.856238
## B2_pB[1] 2.975097 0.896673
## B3_rA[1] 2.847030 0.836951
## B4_rBA[1] -3.786278 0.913364
##
## $p.VC
## B1 B2 B3 B4
## B1 0.733144 -0.757301 -0.712042 0.762779
## B2 -0.757301 0.804022 0.733891 -0.808261
## B3 -0.712042 0.733891 0.700487 -0.747203
## B4 0.762779 -0.808261 -0.747203 0.834234
##
## $VC
## A1 A2 A3 A4 A5 A6 B1
## A1 0.079070 -0.145933 -0.024289 0.178608 0.082469 0.329030 -0.008140
## A2 -0.145933 0.557816 0.080219 -0.922256 -0.342742 -0.342424 -0.183313
## A3 -0.024289 0.080219 0.077925 -0.127854 -0.142763 0.033377 -0.042126
## A4 0.178608 -0.922256 -0.127854 3.055182 0.564698 1.474631 0.295266
## A5 0.082469 -0.342742 -0.142763 0.564698 0.668827 -0.696643 0.210166
## A6 0.329030 -0.342424 0.033377 1.474631 -0.696643 35.317330 -0.561425
## B1 -0.008140 -0.183313 -0.042126 0.295266 0.210166 -0.561425 0.733144
## B2 0.015459 0.184249 0.052723 -0.340038 -0.234071 0.635433 -0.757301
## B3 0.006617 0.165834 0.039935 -0.245849 -0.196825 0.493199 -0.712042
## B4 -0.012218 -0.185518 -0.052563 0.317151 0.237333 -0.567216 0.762779
## B2 B3 B4
## A1 0.015459 0.006617 -0.012218
## A2 0.184249 0.165834 -0.185518
## A3 0.052723 0.039935 -0.052563
## A4 -0.340038 -0.245849 0.317151
## A5 -0.234071 -0.196825 0.237333
## A6 0.635433 0.493199 -0.567216
## B1 -0.757301 -0.712042 0.762779
## B2 0.804022 0.733891 -0.808261
## B3 0.733891 0.700487 -0.747203
## B4 -0.808261 -0.747203 0.834234
model.fits[[check.model]]$dmat
## $psi
## a1 a2 a3 a4 a5 a6
## psiA "1" "0" "elev.std_psiA" "0" "0" "0"
## psiBA "1" "1" "elev.std_psiBA" "0" "SP2:elev.std_psiBA" "0"
## psiBa "1" "1" "elev.std_psiBa" "1" "SP2:elev.std_psiBa" "elev.std:INT2_psiBa"
##
## $p
## b1 b2 b3 b4
## pA[1] "1" "0" "0" "0"
## pA[2] "1" "0" "0" "0"
## pA[3] "1" "0" "0" "0"
## pA[4] "1" "0" "0" "0"
## pA[5] "1" "0" "0" "0"
## pA[6] "1" "0" "0" "0"
## pA[7] "1" "0" "0" "0"
## pA[8] "1" "0" "0" "0"
## pA[9] "1" "0" "0" "0"
## pA[10] "1" "0" "0" "0"
## pA[11] "1" "0" "0" "0"
## pA[12] "1" "0" "0" "0"
## pA[13] "1" "0" "0" "0"
## pA[14] "1" "0" "0" "0"
## pA[15] "1" "0" "0" "0"
## pA[16] "1" "0" "0" "0"
## pA[17] "1" "0" "0" "0"
## pA[18] "1" "0" "0" "0"
## pB[1] "1" "1" "0" "0"
## pB[2] "1" "1" "0" "0"
## pB[3] "1" "1" "0" "0"
## pB[4] "1" "1" "0" "0"
## pB[5] "1" "1" "0" "0"
## pB[6] "1" "1" "0" "0"
## pB[7] "1" "1" "0" "0"
## pB[8] "1" "1" "0" "0"
## pB[9] "1" "1" "0" "0"
## pB[10] "1" "1" "0" "0"
## pB[11] "1" "1" "0" "0"
## pB[12] "1" "1" "0" "0"
## pB[13] "1" "1" "0" "0"
## pB[14] "1" "1" "0" "0"
## pB[15] "1" "1" "0" "0"
## pB[16] "1" "1" "0" "0"
## pB[17] "1" "1" "0" "0"
## pB[18] "1" "1" "0" "0"
## rA[1] "1" "0" "1" "0"
## rA[2] "1" "0" "1" "0"
## rA[3] "1" "0" "1" "0"
## rA[4] "1" "0" "1" "0"
## rA[5] "1" "0" "1" "0"
## rA[6] "1" "0" "1" "0"
## rA[7] "1" "0" "1" "0"
## rA[8] "1" "0" "1" "0"
## rA[9] "1" "0" "1" "0"
## rA[10] "1" "0" "1" "0"
## rA[11] "1" "0" "1" "0"
## rA[12] "1" "0" "1" "0"
## rA[13] "1" "0" "1" "0"
## rA[14] "1" "0" "1" "0"
## rA[15] "1" "0" "1" "0"
## rA[16] "1" "0" "1" "0"
## rA[17] "1" "0" "1" "0"
## rA[18] "1" "0" "1" "0"
## rBA[1] "1" "1" "1" "1"
## rBA[2] "1" "1" "1" "1"
## rBA[3] "1" "1" "1" "1"
## rBA[4] "1" "1" "1" "1"
## rBA[5] "1" "1" "1" "1"
## rBA[6] "1" "1" "1" "1"
## rBA[7] "1" "1" "1" "1"
## rBA[8] "1" "1" "1" "1"
## rBA[9] "1" "1" "1" "1"
## rBA[10] "1" "1" "1" "1"
## rBA[11] "1" "1" "1" "1"
## rBA[12] "1" "1" "1" "1"
## rBA[13] "1" "1" "1" "1"
## rBA[14] "1" "1" "1" "1"
## rBA[15] "1" "1" "1" "1"
## rBA[16] "1" "1" "1" "1"
## rBA[17] "1" "1" "1" "1"
## rBA[18] "1" "1" "1" "1"
## rBa[1] "1" "1" "1" "1"
## rBa[2] "1" "1" "1" "1"
## rBa[3] "1" "1" "1" "1"
## rBa[4] "1" "1" "1" "1"
## rBa[5] "1" "1" "1" "1"
## rBa[6] "1" "1" "1" "1"
## rBa[7] "1" "1" "1" "1"
## rBa[8] "1" "1" "1" "1"
## rBa[9] "1" "1" "1" "1"
## rBa[10] "1" "1" "1" "1"
## rBa[11] "1" "1" "1" "1"
## rBa[12] "1" "1" "1" "1"
## rBa[13] "1" "1" "1" "1"
## rBa[14] "1" "1" "1" "1"
## rBa[15] "1" "1" "1" "1"
## rBa[16] "1" "1" "1" "1"
## rBa[17] "1" "1" "1" "1"
## rBa[18] "1" "1" "1" "1"
names(model.fits[[check.model]]$real)
## [1] "psiA" "psiBA" "psiBa" "pA" "pB" "rA" "rBA" "rBa"
model.fits[[check.model]]$real$psiA[1:5,]
## est se lower upper
## unit1 0.2425 0.0945 0.1046 0.4674
## unit2 0.4520 0.0792 0.3058 0.6069
## unit3 0.4902 0.0746 0.3488 0.6332
## unit4 0.5175 0.0715 0.3796 0.6527
## unit5 0.3880 0.0867 0.2366 0.5647
model.fits[[check.model]]$real$psiBA[1:5,]
## est se lower upper
## unit1 0.0561 0.0583 0.0068 0.3400
## unit2 0.2769 0.1216 0.1043 0.5573
## unit3 0.3414 0.1250 0.1485 0.6064
## unit4 0.3912 0.1260 0.1856 0.6445
## unit5 0.1858 0.1092 0.0525 0.4843
model.fits[[check.model]]$real$psiBa[1:5,]
## est se lower upper
## unit1 1.0000 0.0000 0.0015 1
## unit2 0.9997 0.0018 0.0263 1
## unit3 0.9989 0.0056 0.0403 1
## unit4 0.9973 0.0121 0.0537 1
## unit5 1.0000 0.0003 0.0122 1
model.fits[[check.model]]$derived$nu[1:5,]
## est se lower upper
## unit1 7.371492e-09 8.540258e-08 1.013353e-18 53.62287
## unit2 1.187424e-04 7.229182e-04 7.803536e-10 18.06842
## unit3 5.725349e-04 2.988507e-03 2.063643e-08 15.88435
## unit4 1.750777e-03 8.076598e-03 2.072230e-07 14.79189
## unit5 8.053426e-06 6.118529e-05 2.747405e-12 23.60688
model.fits[[check.model]]$real$pA[1:5,]
## est se lower upper
## unit1_1-1 0.0084 0.0071 0.0016 0.0433
## unit2_1-1 0.0084 0.0071 0.0016 0.0433
## unit3_1-1 0.0084 0.0071 0.0016 0.0433
## unit4_1-1 0.0084 0.0071 0.0016 0.0433
## unit5_1-1 0.0084 0.0071 0.0016 0.0433
model.fits[[check.model]]$real$pB[1:5,]
## est se lower upper
## unit1_1-1 0.1421 0.0183 0.1099 0.1819
## unit2_1-1 0.1421 0.0183 0.1099 0.1819
## unit3_1-1 0.1421 0.0183 0.1099 0.1819
## unit4_1-1 0.1421 0.0183 0.1099 0.1819
## unit5_1-1 0.1421 0.0183 0.1099 0.1819
model.fits[[check.model]]$real$rA[1:5,]
## est se lower upper
## unit1_1-1 0.1272 0.0108 0.1074 0.15
## unit2_1-1 0.1272 0.0108 0.1074 0.15
## unit3_1-1 0.1272 0.0108 0.1074 0.15
## unit4_1-1 0.1272 0.0108 0.1074 0.15
## unit5_1-1 0.1272 0.0108 0.1074 0.15
model.fits[[check.model]]$real$rBA[1:5,]
## est se lower upper
## unit1_1-1 0.0608 0.0071 0.0483 0.0764
## unit2_1-1 0.0608 0.0071 0.0483 0.0764
## unit3_1-1 0.0608 0.0071 0.0483 0.0764
## unit4_1-1 0.0608 0.0071 0.0483 0.0764
## unit5_1-1 0.0608 0.0071 0.0483 0.0764
model.fits[[check.model]]$real$rBa[1:5,]
## est se lower upper
## unit1_1-1 0.0608 0.0071 0.0483 0.0764
## unit2_1-1 0.0608 0.0071 0.0483 0.0764
## unit3_1-1 0.0608 0.0071 0.0483 0.0764
## unit4_1-1 0.0608 0.0071 0.0483 0.0764
## unit5_1-1 0.0608 0.0071 0.0483 0.0764
model.fits[[check.model]]$derived$rho[1:5,]
## est se lower upper
## unit1_1-1 1 0 1 1
## unit2_1-1 1 0 1 1
## unit3_1-1 1 0 1 1
## unit4_1-1 1 0 1 1
## unit5_1-1 1 0 1 1
names(model.fits[[check.model]]$derived)
## [1] "nu" "rho"
# Model averaging
aic.table <- RPresence::createAicTable(model.fits)
aic.table$table
## Model
## 4 psi(SP P elev.std P SP T elev.std P INT P INT T elev.std)p(SP P INT_o P SP T INT_o)
## 3 psi(SP P INT)p(SP P INT_o P SP T INT_o)
## 1 psi(SP P INT)p(SP P INT_o P SP T INT_o P INT_d)
## 2 psi(SP)p(SP P INT_o P SP T INT_o P INT_d)
## AIC neg2ll npar warn.conv warn.VC DAIC modlike wgt
## 4 2392.828 2372.828 10 5.76 0 0.0000 1e+00 0.9997
## 3 2409.173 2395.173 7 7.00 0 16.3442 3e-04 0.0003
## 1 2424.142 2408.142 8 7.00 0 31.3132 0e+00 0.0000
## 2 2425.193 2411.193 7 7.00 0 32.3644 0e+00 0.0000
names(aic.table)
## [1] "table" "models" "ess"
ma.psiA<- RPresence::modAvg(aic.table, param="psiA")
ma.psiA$parameter <- "psiA"
ma.psiA <- cbind(ma.psiA, site.covar)
ma.psiBA<- RPresence::modAvg(aic.table, param="psiBA")
ma.psiBA$parameter <- "psiBA"
ma.psiBA <- cbind(ma.psiBA, site.covar)
ma.psiBa<- RPresence::modAvg(aic.table, param="psiBa")
ma.psiBa$parameter <- "psiBa"
ma.psiBa <- cbind(ma.psiBa, site.covar)
ma.psiB <- site.covar
ma.psiB$parameter <- 'psiB'
ma.psiB$est <- ma.psiBA$est * ma.psiA$est +
ma.psiBa$est * (1-ma.psiA$est)
# too hard to compute the se
psi.all <- plyr::rbind.fill(ma.psiA, ma.psiBA, ma.psiBa, ma.psiB)
head(psi.all)
## est se lower_0.95 upper_0.95 parameter
## 1 0.2426393 0.09485377 0.1043374 0.4683937 psiA
## 2 0.4520801 0.07933756 0.3057571 0.6071842 psiA
## 3 0.4902693 0.07470921 0.3486771 0.6334401 psiA
## 4 0.5175615 0.07158990 0.3794996 0.6529933 psiA
## 5 0.3880981 0.08688928 0.2364156 0.5650799 psiA
## 6 0.3880981 0.08688928 0.2364156 0.5650799 psiA
## Cam_name X Y Land_type droad
## 1 amethyst lake 415655 5838411 415655 5838411 Jasper 12759.7998
## 2 ancient wall 376992 5919664 376992 5919664 Jasper 47406.6016
## 3 astoria 427308 5839183 427308 5839183 Jasper 1134.0400
## 4 bess pass 343359 5911520 343359 5911520 Crown 76867.2031
## 5 brewster skywalk c2 481109 5789822 481109 5789822 Jasper 44.7214
## 6 brewster skywalk c2b 481109 5789822 481107 5789822 Jasper 49.2443
## dstream elev tpi NDVI elev.std
## 1 55.0000 1994 -212.058 7102.968 1.8858767
## 2 53.1507 1760 -373.933 6576.896 0.9404520
## 3 640.3120 1722 -259.626 6440.972 0.7869215
## 4 530.0000 1695 0.000 6789.304 0.6778340
## 5 101.2420 1825 -608.371 5605.156 1.2030700
## 6 97.0824 1825 -608.371 5605.156 1.2030700
ggplot(data=psi.all, aes(x=elev, y=est))+
ggtitle("Estimated occupancy parameters vs. elevation")+
geom_point()+
geom_ribbon( aes(ymin=lower_0.95, ymax=upper_0.95),alpha=0.2)+
facet_wrap(~parameter, ncol=2)
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
# Why is psiBa so bad? This may be a case of complete separation.
# The detection probabilities are all fairly large so that over 5 visits, there is almost certain detection
# if the species is present. Consequently, the observed occupancy is likely a pretty good indication of the
# actual occupancy.
RPresence::modAvg(aic.table, param="pA")[1,]
## est se lower_0.95 upper_0.95
## unit1_1-1 0.008479604 0.008619334 0.001145378 0.0599582
RPresence::modAvg(aic.table, param="pB")[1,]
## est se lower_0.95 upper_0.95
## unit1_1-1 0.1421618 0.01884126 0.1090716 0.1832263
RPresence::modAvg(aic.table, param="rA")[1,]
## est se lower_0.95 upper_0.95
## unit1_1-1 0.1271818 0.01085396 0.1073835 0.1500168
RPresence::modAvg(aic.table, param="rBA")[1,]
## est se lower_0.95 upper_0.95
## unit1_1-1 0.06080363 0.00710522 0.04828086 0.07631401
RPresence::modAvg(aic.table, param="rBa")[1,]
## est se lower_0.95 upper_0.95
## unit1_1-1 0.06080363 0.007105256 0.04828081 0.0763141
# Compute the observed occupance of each species at each camera
# Don't forget that cameras record number of times bear are seen
obs.occ <- data.frame(elev = site.covar$elev,
oogb = pmin(1, apply(gb.data,1,max, na.rm=TRUE)),
oobb = pmin(1, apply(bb.data,1,max, na.rm=TRUE)))
head(obs.occ)
## elev oogb oobb
## 1 1994 1 0
## 2 1760 0 0
## 3 1722 0 0
## 4 1695 1 0
## 5 1825 0 0
## 6 1825 0 0
ggplot(data=obs.occ, aes(x=elev, y=oobb))+
ggtitle("Observed occupancy of bb given gb")+
geom_point()+
geom_smooth(method="glm", method.args = list(family = quasibinomial(link = 'logit')))+
facet_wrap(~oogb, ncol=1)
## `geom_smooth()` using formula 'y ~ x'
# and look at when the species are reversed
ggplot(data=obs.occ, aes(x=elev, y=oogb))+
ggtitle("Observed occupancy of gb given bb")+
geom_point()+
geom_smooth(method="glm", method.args = list(family = quasibinomial(link = 'logit')))+
facet_wrap(~oobb, ncol=1)
## `geom_smooth()` using formula 'y ~ x'
