# 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.
# 2018-12-26 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
library(car)
## Loading required package: carData
library(ggplot2)
library(readxl)
library(reshape2)
library(RMark)
## This is RMark 2.2.5
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
# 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
## [36] 0 0 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
## [71] 0 0 1 2 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
## [106] 0 2 2 2 1 0 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
## [141] 2 1 1 0 1 1 0 1 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
## [176] 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.data <- (bb.data>=1) + 2*(gb.data>=1)
input.data
## 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
input.chistory <- data.frame(lapply(as.data.frame(input.data), as.character), stringsAsFactors=FALSE)
input.chistory <- data.frame(lapply(input.chistory,
car::recode, " '0'='00'; '2'='10';'1'='01';'3'='11';", as.numeric=FALSE, as.factor=FALSE),
stringsAsFactors=FALSE)
input.history <- data.frame(ch=apply(input.chistory, 1, paste, sep="",collapse=""), freq=1)
# Change any NA to . in the chapter history
select <- grepl("NA", input.history$ch)
input.history[ select,][1:5,]
## ch freq
## 2 0000NANANANANANANANANANANANANANANANA 1
## 4 1010101010101010NANANANANANANANANANA 1
## 9 0000000000000000NANANANANANANANANANA 1
## 16 00100000100000NANANANANANANANANANANA 1
## 28 00001000010000000000000000101000NANA 1
input.history$ch <- gsub("NA","..", input.history$ch, fixed=TRUE)
input.history[ select,][1:5,]
## ch freq
## 2 0000................................ 1
## 4 1010101010101010.................... 1
## 9 0000000000000000.................... 1
## 16 00100000100000...................... 1
## 28 00001000010000000000000000101000.... 1
# 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
# add the site covariate to the chapture hisotry
input.history <- cbind(input.history, site.covar)
head(input.history)
## ch freq
## 1 000000000000101000000000000000000000 1
## 2 0000................................ 1
## 3 000000000000000000000000000000000000 1
## 4 1010101010101010.................... 1
## 5 000000000000000000000000000000000000 1
## 6 000000000000000000000000000000000000 1
## 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
Nsites <- nrow(input.history)
Nvisits<- nchar(input.history$ch[1])
bear.data <- process.data(data=input.history,
model="2SpecConOccup") # this parameterization is more stable
summary(bear.data)
## Length Class Mode
## data 12 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 182 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 18 -none- numeric
## begin.time 1 -none- numeric
## age.unit 1 -none- numeric
## initial.ages 1 -none- numeric
## group.covariates 0 -none- NULL
## nstrata 1 -none- numeric
## strata.labels 0 -none- NULL
## counts 0 -none- NULL
## reverse 1 -none- logical
## areas 0 -none- NULL
## events 0 -none- NULL
# are there any time-specific covariates? If so add them to the ddl's
bear.ddl <- RMark::make.design.data(bear.data) # need for covariate predictions
# What are the parameter names for Single Season Single Species models
setup.parameters("2SpecConOccup", check=TRUE)
## [1] "PsiA" "PsiBA" "PsiBa" "pA" "pB" "rA" "rBA" "rBa"
# 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
setup.parameters("2SpecConOccup", check=TRUE)
## [1] "PsiA" "PsiBA" "PsiBa" "pA" "pB" "rA" "rBA" "rBa"
# Get the list of models. NOtice NO equal signs here
model.list.csv <- textConnection("
PsiA, PsiBA, PsiBa, pA, pB, rA, rBA, rBa
~elev_std, ~elev_std, ~elev_std, ~1, ~1, ~1, ~1, ~SHARE
~1, ~1, ~1, ~1, ~1, ~1, ~1, ~1
~1, ~1, ~1, ~1, ~1, ~1, ~1, ~SHARE
~1, ~1, ~SHARE, ~1, ~1, ~1, ~1, ~1
~1, ~1, ~SHARE, ~1, ~1, ~1, ~1, ~SHARE")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list$model.number <- 1:nrow(model.list)
model.list
## PsiA PsiBA PsiBa pA pB rA rBA rBa model.number
## 1 ~elev_std ~elev_std ~elev_std ~1 ~1 ~1 ~1 ~SHARE 1
## 2 ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1 2
## 3 ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~SHARE 3
## 4 ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~1 4
## 5 ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~SHARE 5
# fit the models
myobs <- ls()
myobs <- myobs[ grepl("m...",myobs,fixed=TRUE)]
cat("Removing ", myobs, "\n")
## Removing
rm(list=myobs)
model.fits <- plyr::dlply(model.list, "model.number", function(x,input.data, input.ddl){
cat("\n\n***** Starting ", unlist(x), "\n")
#browser()
model.parameters=list(
PsiA =list(formula=as.formula(eval(x$PsiA ))),
PsiBA =list(formula=as.formula(eval(x$PsiBA))),
PsiBa =list(formula=as.formula(eval(x$PsiBa))),
pA =list(formula=as.formula(eval(x$pA ))),
pB =list(formula=as.formula(eval(x$pB ))),
rA =list(formula=as.formula(eval(x$rA ))),
rBA =list(formula=as.formula(eval(x$rBA))),
rBa =list(formula=as.formula(eval(x$rBa)))
)
if(x$rBa == '~SHARE'){
model.parameters$rBA =list(formula=as.formula(eval(x$rBA)), share=TRUE)
model.parameters$rBa = NULL
}
if(x$PsiBa == '~SHARE'){
model.parameters$PsiBA =list(formula=as.formula(eval(x$PsiBA)), share=TRUE)
model.parameters$PsiBa = NULL
}
fit <- RMark::mark(input.data, ddl=input.ddl,
model="2SpecConOccup",
model.parameters=model.parameters
#,brief=TRUE,output=FALSE, delete=TRUE
#,invisible=TRUE,output=TRUE # set for debugging
)
mnumber <- paste("m...",formatC(x$model.number, width = 3, format = "d", flag = "0"),sep="")
assign( mnumber, fit, envir=.GlobalEnv)
#browser()
fit
},input.data=bear.data, input.ddl=bear.ddl)
##
##
## ***** Starting ~elev_std ~elev_std ~elev_std ~1 ~1 ~1 ~1 ~SHARE 1
##
## Note: only 9 parameters counted of 10 specified parameters
## AICc and parameter count have been adjusted upward
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~elev_std)PsiBA(~elev_std)PsiBa(~elev_std)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 10 (unadjusted=9)
## -2lnL: 2360.583
## AICc : 2381.87 (unadjusted=2379.6297)
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.4927662 0.2342841 0.0335693 0.9519630
## PsiA:elev_std 0.3729369 0.2104800 -0.0396040 0.7854778
## PsiBA:(Intercept) 4.4911960 1.0479925 2.4371306 6.5452614
## PsiBA:elev_std -0.7055059 0.8999841 -2.4694747 1.0584629
## PsiBa:(Intercept) -1.6751620 0.8887774 -3.4171658 0.0668418
## PsiBa:elev_std -3.1110476 1.4590367 -5.9707596 -0.2513357
## pA:(Intercept) 16.2869740 1186.2913000 -2308.8441000 2341.4180000
## pB:(Intercept) -1.4079082 0.1902395 -1.7807777 -1.0350387
## rA:(Intercept) -2.1926439 0.1152314 -2.4184976 -1.9667903
## rBA:(Intercept) -2.8737722 0.1206942 -3.1103327 -2.6372116
##
##
## Real Parameter PsiA
## 1
## 0.6207579
##
##
## Real Parameter PsiBA
## 1
## 0.988917
##
##
## Real Parameter PsiBa
## 1
## 0.1577372
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999
## 8 9 10 11 12 13 14
## 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999
## 15 16 17 18
## 0.9999999 0.9999999 0.9999999 0.9999999
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642
## 8 9 10 11 12 13 14
## 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642
## 15 16 17 18
## 0.1965642 0.1965642 0.1965642 0.1965642
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7 8
## 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413
## 9 10 11 12 13 14 15 16
## 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413
## 17 18
## 0.100413 0.100413
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654
## 8 9 10 11 12 13 14
## 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654
## 15 16 17 18
## 0.0534654 0.0534654 0.0534654 0.0534654
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654
## 8 9 10 11 12 13 14
## 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654
## 15 16 17 18
## 0.0534654 0.0534654 0.0534654 0.0534654
##
##
## ***** Starting ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1 2
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
##
## Npar : 8
## -2lnL: 2395.169
## AICc : 2412.002
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.9244542 0.2836842 0.3684331 1.4804753
## PsiBA:(Intercept) 2.4147570 0.7766575 0.8925083 3.9370057
## PsiBa:(Intercept) -0.9648749 0.5621798 -2.0667473 0.1369976
## pA:(Intercept) -0.5662061 0.6739858 -1.8872182 0.7548060
## pB:(Intercept) -0.8931236 0.3503614 -1.5798318 -0.2064153
## rA:(Intercept) -2.5318376 0.1777977 -2.8803210 -2.1833541
## rBA:(Intercept) -2.7247165 0.3960677 -3.5010092 -1.9484239
## rBa:(Intercept) -2.7027629 0.1673516 -3.0307720 -2.3747537
##
##
## Real Parameter PsiA
## 1
## 0.7159488
##
##
## Real Parameter PsiBA
## 1
## 0.9179457
##
##
## Real Parameter PsiBa
## 1
## 0.2759032
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.3621127 0.3621127 0.3621127 0.3621127 0.3621127 0.3621127 0.3621127
## 8 9 10 11 12 13 14
## 0.3621127 0.3621127 0.3621127 0.3621127 0.3621127 0.3621127 0.3621127
## 15 16 17 18
## 0.3621127 0.3621127 0.3621127 0.3621127
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.2904657 0.2904657 0.2904657 0.2904657 0.2904657 0.2904657 0.2904657
## 8 9 10 11 12 13 14
## 0.2904657 0.2904657 0.2904657 0.2904657 0.2904657 0.2904657 0.2904657
## 15 16 17 18
## 0.2904657 0.2904657 0.2904657 0.2904657
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.0736562 0.0736562 0.0736562 0.0736562 0.0736562 0.0736562 0.0736562
## 8 9 10 11 12 13 14
## 0.0736562 0.0736562 0.0736562 0.0736562 0.0736562 0.0736562 0.0736562
## 15 16 17 18
## 0.0736562 0.0736562 0.0736562 0.0736562
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.0615305 0.0615305 0.0615305 0.0615305 0.0615305 0.0615305 0.0615305
## 8 9 10 11 12 13 14
## 0.0615305 0.0615305 0.0615305 0.0615305 0.0615305 0.0615305 0.0615305
## 15 16 17 18
## 0.0615305 0.0615305 0.0615305 0.0615305
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.0628105 0.0628105 0.0628105 0.0628105 0.0628105 0.0628105 0.0628105
## 8 9 10 11 12 13 14
## 0.0628105 0.0628105 0.0628105 0.0628105 0.0628105 0.0628105 0.0628105
## 15 16 17 18
## 0.0628105 0.0628105 0.0628105 0.0628105
##
##
## ***** Starting ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~SHARE 3
##
## Note: only 6 parameters counted of 7 specified parameters
##
## AICc and parameter count have been adjusted upward
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 7 (unadjusted=6)
## -2lnL: 2394.819
## AICc : 2409.463 (unadjusted=2407.299)
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.6590358 0.2731688 0.1236249 1.1944467
## PsiBA:(Intercept) 4.2797053 0.8524949 2.6088153 5.9505952
## PsiBa:(Intercept) -0.5651372 0.4154503 -1.3794199 0.2491455
## pA:(Intercept) 20.9254000 0.0000000 20.9254000 20.9254000
## pB:(Intercept) -1.2123090 0.2684189 -1.7384102 -0.6862079
## rA:(Intercept) -2.2678667 0.1222930 -2.5075610 -2.0281725
## rBA:(Intercept) -2.8326722 0.1242855 -3.0762717 -2.5890727
##
##
## Real Parameter PsiA
## 1
## 0.6590438
##
##
## Real Parameter PsiBA
## 1
## 0.9863424
##
##
## Real Parameter PsiBa
## 1
## 0.3623596
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.2292927 0.2292927 0.2292927 0.2292927 0.2292927 0.2292927 0.2292927
## 8 9 10 11 12 13 14
## 0.2292927 0.2292927 0.2292927 0.2292927 0.2292927 0.2292927 0.2292927
## 15 16 17 18
## 0.2292927 0.2292927 0.2292927 0.2292927
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.0938194 0.0938194 0.0938194 0.0938194 0.0938194 0.0938194 0.0938194
## 8 9 10 11 12 13 14
## 0.0938194 0.0938194 0.0938194 0.0938194 0.0938194 0.0938194 0.0938194
## 15 16 17 18
## 0.0938194 0.0938194 0.0938194 0.0938194
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7 8
## 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584
## 9 10 11 12 13 14 15 16
## 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584
## 17 18
## 0.055584 0.055584
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7 8
## 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584
## 9 10 11 12 13 14 15 16
## 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584 0.055584
## 17 18
## 0.055584 0.055584
##
##
## ***** Starting ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~1 4
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
##
## Npar : 7
## -2lnL: 2422.278
## AICc : 2436.922
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 2.2730503 0.4424141 1.405919 3.1401820
## PsiBA:(Intercept) 2.2927617 0.4769166 1.358005 3.2275182
## pA:(Intercept) -0.7625181 0.3378840 -1.424771 -0.1002655
## pB:(Intercept) -0.8772258 0.2504671 -1.368141 -0.3863103
## rA:(Intercept) -2.8620846 0.1175235 -3.092431 -2.6317386
## rBA:(Intercept) -2.6529489 0.3767958 -3.391469 -1.9144290
## rBa:(Intercept) -2.9866548 0.1274735 -3.236503 -2.7368067
##
##
## Real Parameter PsiA
## 1
## 0.9066203
##
##
## Real Parameter PsiBA
## 1
## 0.9082758
##
##
## Real Parameter PsiBa
## 1
## 0.9082758
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.3180998 0.3180998 0.3180998 0.3180998 0.3180998 0.3180998 0.3180998
## 8 9 10 11 12 13 14
## 0.3180998 0.3180998 0.3180998 0.3180998 0.3180998 0.3180998 0.3180998
## 15 16 17 18
## 0.3180998 0.3180998 0.3180998 0.3180998
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7 8
## 0.293753 0.293753 0.293753 0.293753 0.293753 0.293753 0.293753 0.293753
## 9 10 11 12 13 14 15 16
## 0.293753 0.293753 0.293753 0.293753 0.293753 0.293753 0.293753 0.293753
## 17 18
## 0.293753 0.293753
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7 8 9
## 0.05406 0.05406 0.05406 0.05406 0.05406 0.05406 0.05406 0.05406 0.05406
## 10 11 12 13 14 15 16 17 18
## 0.05406 0.05406 0.05406 0.05406 0.05406 0.05406 0.05406 0.05406 0.05406
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.0658075 0.0658075 0.0658075 0.0658075 0.0658075 0.0658075 0.0658075
## 8 9 10 11 12 13 14
## 0.0658075 0.0658075 0.0658075 0.0658075 0.0658075 0.0658075 0.0658075
## 15 16 17 18
## 0.0658075 0.0658075 0.0658075 0.0658075
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.0480324 0.0480324 0.0480324 0.0480324 0.0480324 0.0480324 0.0480324
## 8 9 10 11 12 13 14
## 0.0480324 0.0480324 0.0480324 0.0480324 0.0480324 0.0480324 0.0480324
## 15 16 17 18
## 0.0480324 0.0480324 0.0480324 0.0480324
##
##
## ***** Starting ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~SHARE 5
##
## Note: only 5 parameters counted of 6 specified parameters
##
## AICc and parameter count have been adjusted upward
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 6 (unadjusted=5)
## -2lnL: 2431.749
## AICc : 2444.229 (unadjusted=2442.0898)
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 2.294224 0.3912768 1.527322 3.0611270
## PsiBA:(Intercept) 4.441123 0.8129395 2.847761 6.0344842
## pA:(Intercept) 20.102175 0.0000000 20.102175 20.1021750
## pB:(Intercept) -0.966923 0.2329031 -1.423413 -0.5104328
## rA:(Intercept) -2.639570 0.0875883 -2.811243 -2.4678967
## rBA:(Intercept) -3.077571 0.1169311 -3.306756 -2.8483865
##
##
## Real Parameter PsiA
## 1
## 0.9083976
##
##
## Real Parameter PsiBA
## 1
## 0.9883545
##
##
## Real Parameter PsiBa
## 1
## 0.9883545
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.2754942 0.2754942 0.2754942 0.2754942 0.2754942 0.2754942 0.2754942
## 8 9 10 11 12 13 14
## 0.2754942 0.2754942 0.2754942 0.2754942 0.2754942 0.2754942 0.2754942
## 15 16 17 18
## 0.2754942 0.2754942 0.2754942 0.2754942
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.0666348 0.0666348 0.0666348 0.0666348 0.0666348 0.0666348 0.0666348
## 8 9 10 11 12 13 14
## 0.0666348 0.0666348 0.0666348 0.0666348 0.0666348 0.0666348 0.0666348
## 15 16 17 18
## 0.0666348 0.0666348 0.0666348 0.0666348
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7 8
## 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042
## 9 10 11 12 13 14 15 16
## 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042
## 17 18
## 0.044042 0.044042
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7 8
## 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042
## 9 10 11 12 13 14 15 16
## 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042 0.044042
## 17 18
## 0.044042 0.044042
# examine individual model results
model.number <-1
summary(model.fits[[model.number]])
## Output summary for 2SpecConOccup model
## Name : PsiA(~elev_std)PsiBA(~elev_std)PsiBa(~elev_std)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 10 (unadjusted=9)
## -2lnL: 2360.583
## AICc : 2381.87 (unadjusted=2379.6297)
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.4927662 0.2342841 0.0335693 0.9519630
## PsiA:elev_std 0.3729369 0.2104800 -0.0396040 0.7854778
## PsiBA:(Intercept) 4.4911960 1.0479925 2.4371306 6.5452614
## PsiBA:elev_std -0.7055059 0.8999841 -2.4694747 1.0584629
## PsiBa:(Intercept) -1.6751620 0.8887774 -3.4171658 0.0668418
## PsiBa:elev_std -3.1110476 1.4590367 -5.9707596 -0.2513357
## pA:(Intercept) 16.2869740 1186.2913000 -2308.8441000 2341.4180000
## pB:(Intercept) -1.4079082 0.1902395 -1.7807777 -1.0350387
## rA:(Intercept) -2.1926439 0.1152314 -2.4184976 -1.9667903
## rBA:(Intercept) -2.8737722 0.1206942 -3.1103327 -2.6372116
##
##
## Real Parameter PsiA
## 1
## 0.6207579
##
##
## Real Parameter PsiBA
## 1
## 0.988917
##
##
## Real Parameter PsiBa
## 1
## 0.1577372
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999
## 8 9 10 11 12 13 14
## 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999 0.9999999
## 15 16 17 18
## 0.9999999 0.9999999 0.9999999 0.9999999
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642
## 8 9 10 11 12 13 14
## 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642 0.1965642
## 15 16 17 18
## 0.1965642 0.1965642 0.1965642 0.1965642
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7 8
## 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413
## 9 10 11 12 13 14 15 16
## 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413 0.100413
## 17 18
## 0.100413 0.100413
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654
## 8 9 10 11 12 13 14
## 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654
## 15 16 17 18
## 0.0534654 0.0534654 0.0534654 0.0534654
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654
## 8 9 10 11 12 13 14
## 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654 0.0534654
## 15 16 17 18
## 0.0534654 0.0534654 0.0534654 0.0534654
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed
## PsiA g1 a0 t1 0.6207579 0.0551546000 5.083915e-01 0.7215098
## PsiBA g1 a0 t1 0.9889170 0.0114862000 9.196152e-01 0.9985652
## PsiBa g1 a0 t1 0.1577372 0.1180795000 3.176330e-02 0.5167042
## pA g1 a0 t1 0.9999999 0.0001001955 6.586088e-302 1.0000000
## pB g1 a0 t1 0.1965642 0.0300439000 1.442071e-01 0.2621084
## rA g1 a0 t1 0.1004130 0.0104089000 8.177300e-02 0.1227341
## rBA g1 a0 t1 0.0534654 0.0061080000 4.268300e-02 0.0667816
## note
## PsiA g1 a0 t1
## PsiBA g1 a0 t1
## PsiBa g1 a0 t1
## pA g1 a0 t1
## pB g1 a0 t1
## rA g1 a0 t1
## rBA g1 a0 t1
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.4927662 0.2342841 0.0335693 0.9519630
## PsiA:elev_std 0.3729369 0.2104800 -0.0396040 0.7854778
## PsiBA:(Intercept) 4.4911960 1.0479925 2.4371306 6.5452614
## PsiBA:elev_std -0.7055059 0.8999841 -2.4694747 1.0584629
## PsiBa:(Intercept) -1.6751620 0.8887774 -3.4171658 0.0668418
## PsiBa:elev_std -3.1110476 1.4590367 -5.9707596 -0.2513357
## pA:(Intercept) 16.2869740 1186.2913000 -2308.8441000 2341.4180000
## pB:(Intercept) -1.4079082 0.1902395 -1.7807777 -1.0350387
## rA:(Intercept) -2.1926439 0.1152314 -2.4184976 -1.9667903
## rBA:(Intercept) -2.8737722 0.1206942 -3.1103327 -2.6372116
model.fits[[model.number]]$results$derived
## $`Species Interaction Factor`
## estimate se lcl ucl
## 1 1.467892 0.1130743 1.246267 1.689518
##
## $`Occupancy of Species B`
## estimate se lcl ucl
## 1 0.6736986 0.05248427 0.5639011 0.7672615
##
## $`Occupancy of Both Species`
## estimate se lcl ucl
## 1 0.613878 0.05527848 0.5016372 0.7151922
get.real(model.fits[[model.number]], "PsiA", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiA g1 a0 t1 1 1 0.6207579 0.0551546 0.5083915
## ucl fixed note group age time Age Time
## PsiA g1 a0 t1 0.7215098 1 0 1 0 0
get.real(model.fits[[model.number]], "PsiBA", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiBA g1 a0 t1 2 2 0.988917 0.0114862 0.9196152
## ucl fixed note group age time Age Time
## PsiBA g1 a0 t1 0.9985652 1 0 1 0 0
get.real(model.fits[[model.number]], "PsiBa", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiBa g1 a0 t1 3 3 0.1577372 0.1180795 0.0317633
## ucl fixed note group age time Age Time
## PsiBa g1 a0 t1 0.5167042 1 0 1 0 0
get.real(model.fits[[model.number]], "pA", se=TRUE)
## all.diff.index par.index estimate se
## pA g1 a0 t1 4 4 0.9999999 0.0001001955
## pA g1 a1 t2 5 4 0.9999999 0.0001001955
## pA g1 a2 t3 6 4 0.9999999 0.0001001955
## pA g1 a3 t4 7 4 0.9999999 0.0001001955
## pA g1 a4 t5 8 4 0.9999999 0.0001001955
## pA g1 a5 t6 9 4 0.9999999 0.0001001955
## pA g1 a6 t7 10 4 0.9999999 0.0001001955
## pA g1 a7 t8 11 4 0.9999999 0.0001001955
## pA g1 a8 t9 12 4 0.9999999 0.0001001955
## pA g1 a9 t10 13 4 0.9999999 0.0001001955
## pA g1 a10 t11 14 4 0.9999999 0.0001001955
## pA g1 a11 t12 15 4 0.9999999 0.0001001955
## pA g1 a12 t13 16 4 0.9999999 0.0001001955
## pA g1 a13 t14 17 4 0.9999999 0.0001001955
## pA g1 a14 t15 18 4 0.9999999 0.0001001955
## pA g1 a15 t16 19 4 0.9999999 0.0001001955
## pA g1 a16 t17 20 4 0.9999999 0.0001001955
## pA g1 a17 t18 21 4 0.9999999 0.0001001955
## lcl ucl fixed note group age time Age Time
## pA g1 a0 t1 6.586088e-302 1 1 0 1 0 0
## pA g1 a1 t2 6.586088e-302 1 1 1 2 1 1
## pA g1 a2 t3 6.586088e-302 1 1 2 3 2 2
## pA g1 a3 t4 6.586088e-302 1 1 3 4 3 3
## pA g1 a4 t5 6.586088e-302 1 1 4 5 4 4
## pA g1 a5 t6 6.586088e-302 1 1 5 6 5 5
## pA g1 a6 t7 6.586088e-302 1 1 6 7 6 6
## pA g1 a7 t8 6.586088e-302 1 1 7 8 7 7
## pA g1 a8 t9 6.586088e-302 1 1 8 9 8 8
## pA g1 a9 t10 6.586088e-302 1 1 9 10 9 9
## pA g1 a10 t11 6.586088e-302 1 1 10 11 10 10
## pA g1 a11 t12 6.586088e-302 1 1 11 12 11 11
## pA g1 a12 t13 6.586088e-302 1 1 12 13 12 12
## pA g1 a13 t14 6.586088e-302 1 1 13 14 13 13
## pA g1 a14 t15 6.586088e-302 1 1 14 15 14 14
## pA g1 a15 t16 6.586088e-302 1 1 15 16 15 15
## pA g1 a16 t17 6.586088e-302 1 1 16 17 16 16
## pA g1 a17 t18 6.586088e-302 1 1 17 18 17 17
get.real(model.fits[[model.number]], "pB", se=TRUE)
## all.diff.index par.index estimate se lcl
## pB g1 a0 t1 22 5 0.1965642 0.0300439 0.1442071
## pB g1 a1 t2 23 5 0.1965642 0.0300439 0.1442071
## pB g1 a2 t3 24 5 0.1965642 0.0300439 0.1442071
## pB g1 a3 t4 25 5 0.1965642 0.0300439 0.1442071
## pB g1 a4 t5 26 5 0.1965642 0.0300439 0.1442071
## pB g1 a5 t6 27 5 0.1965642 0.0300439 0.1442071
## pB g1 a6 t7 28 5 0.1965642 0.0300439 0.1442071
## pB g1 a7 t8 29 5 0.1965642 0.0300439 0.1442071
## pB g1 a8 t9 30 5 0.1965642 0.0300439 0.1442071
## pB g1 a9 t10 31 5 0.1965642 0.0300439 0.1442071
## pB g1 a10 t11 32 5 0.1965642 0.0300439 0.1442071
## pB g1 a11 t12 33 5 0.1965642 0.0300439 0.1442071
## pB g1 a12 t13 34 5 0.1965642 0.0300439 0.1442071
## pB g1 a13 t14 35 5 0.1965642 0.0300439 0.1442071
## pB g1 a14 t15 36 5 0.1965642 0.0300439 0.1442071
## pB g1 a15 t16 37 5 0.1965642 0.0300439 0.1442071
## pB g1 a16 t17 38 5 0.1965642 0.0300439 0.1442071
## pB g1 a17 t18 39 5 0.1965642 0.0300439 0.1442071
## ucl fixed note group age time Age Time
## pB g1 a0 t1 0.2621084 1 0 1 0 0
## pB g1 a1 t2 0.2621084 1 1 2 1 1
## pB g1 a2 t3 0.2621084 1 2 3 2 2
## pB g1 a3 t4 0.2621084 1 3 4 3 3
## pB g1 a4 t5 0.2621084 1 4 5 4 4
## pB g1 a5 t6 0.2621084 1 5 6 5 5
## pB g1 a6 t7 0.2621084 1 6 7 6 6
## pB g1 a7 t8 0.2621084 1 7 8 7 7
## pB g1 a8 t9 0.2621084 1 8 9 8 8
## pB g1 a9 t10 0.2621084 1 9 10 9 9
## pB g1 a10 t11 0.2621084 1 10 11 10 10
## pB g1 a11 t12 0.2621084 1 11 12 11 11
## pB g1 a12 t13 0.2621084 1 12 13 12 12
## pB g1 a13 t14 0.2621084 1 13 14 13 13
## pB g1 a14 t15 0.2621084 1 14 15 14 14
## pB g1 a15 t16 0.2621084 1 15 16 15 15
## pB g1 a16 t17 0.2621084 1 16 17 16 16
## pB g1 a17 t18 0.2621084 1 17 18 17 17
get.real(model.fits[[model.number]], "rA", se=TRUE)
## all.diff.index par.index estimate se lcl
## rA g1 a0 t1 40 6 0.100413 0.0104089 0.081773
## rA g1 a1 t2 41 6 0.100413 0.0104089 0.081773
## rA g1 a2 t3 42 6 0.100413 0.0104089 0.081773
## rA g1 a3 t4 43 6 0.100413 0.0104089 0.081773
## rA g1 a4 t5 44 6 0.100413 0.0104089 0.081773
## rA g1 a5 t6 45 6 0.100413 0.0104089 0.081773
## rA g1 a6 t7 46 6 0.100413 0.0104089 0.081773
## rA g1 a7 t8 47 6 0.100413 0.0104089 0.081773
## rA g1 a8 t9 48 6 0.100413 0.0104089 0.081773
## rA g1 a9 t10 49 6 0.100413 0.0104089 0.081773
## rA g1 a10 t11 50 6 0.100413 0.0104089 0.081773
## rA g1 a11 t12 51 6 0.100413 0.0104089 0.081773
## rA g1 a12 t13 52 6 0.100413 0.0104089 0.081773
## rA g1 a13 t14 53 6 0.100413 0.0104089 0.081773
## rA g1 a14 t15 54 6 0.100413 0.0104089 0.081773
## rA g1 a15 t16 55 6 0.100413 0.0104089 0.081773
## rA g1 a16 t17 56 6 0.100413 0.0104089 0.081773
## rA g1 a17 t18 57 6 0.100413 0.0104089 0.081773
## ucl fixed note group age time Age Time
## rA g1 a0 t1 0.1227341 1 0 1 0 0
## rA g1 a1 t2 0.1227341 1 1 2 1 1
## rA g1 a2 t3 0.1227341 1 2 3 2 2
## rA g1 a3 t4 0.1227341 1 3 4 3 3
## rA g1 a4 t5 0.1227341 1 4 5 4 4
## rA g1 a5 t6 0.1227341 1 5 6 5 5
## rA g1 a6 t7 0.1227341 1 6 7 6 6
## rA g1 a7 t8 0.1227341 1 7 8 7 7
## rA g1 a8 t9 0.1227341 1 8 9 8 8
## rA g1 a9 t10 0.1227341 1 9 10 9 9
## rA g1 a10 t11 0.1227341 1 10 11 10 10
## rA g1 a11 t12 0.1227341 1 11 12 11 11
## rA g1 a12 t13 0.1227341 1 12 13 12 12
## rA g1 a13 t14 0.1227341 1 13 14 13 13
## rA g1 a14 t15 0.1227341 1 14 15 14 14
## rA g1 a15 t16 0.1227341 1 15 16 15 15
## rA g1 a16 t17 0.1227341 1 16 17 16 16
## rA g1 a17 t18 0.1227341 1 17 18 17 17
get.real(model.fits[[model.number]], "rBA", se=TRUE)
## all.diff.index par.index estimate se lcl
## rBA g1 a0 t1 58 7 0.0534654 0.006108 0.042683
## rBA g1 a1 t2 59 7 0.0534654 0.006108 0.042683
## rBA g1 a2 t3 60 7 0.0534654 0.006108 0.042683
## rBA g1 a3 t4 61 7 0.0534654 0.006108 0.042683
## rBA g1 a4 t5 62 7 0.0534654 0.006108 0.042683
## rBA g1 a5 t6 63 7 0.0534654 0.006108 0.042683
## rBA g1 a6 t7 64 7 0.0534654 0.006108 0.042683
## rBA g1 a7 t8 65 7 0.0534654 0.006108 0.042683
## rBA g1 a8 t9 66 7 0.0534654 0.006108 0.042683
## rBA g1 a9 t10 67 7 0.0534654 0.006108 0.042683
## rBA g1 a10 t11 68 7 0.0534654 0.006108 0.042683
## rBA g1 a11 t12 69 7 0.0534654 0.006108 0.042683
## rBA g1 a12 t13 70 7 0.0534654 0.006108 0.042683
## rBA g1 a13 t14 71 7 0.0534654 0.006108 0.042683
## rBA g1 a14 t15 72 7 0.0534654 0.006108 0.042683
## rBA g1 a15 t16 73 7 0.0534654 0.006108 0.042683
## rBA g1 a16 t17 74 7 0.0534654 0.006108 0.042683
## rBA g1 a17 t18 75 7 0.0534654 0.006108 0.042683
## ucl fixed note group age time Age Time
## rBA g1 a0 t1 0.0667816 1 0 1 0 0
## rBA g1 a1 t2 0.0667816 1 1 2 1 1
## rBA g1 a2 t3 0.0667816 1 2 3 2 2
## rBA g1 a3 t4 0.0667816 1 3 4 3 3
## rBA g1 a4 t5 0.0667816 1 4 5 4 4
## rBA g1 a5 t6 0.0667816 1 5 6 5 5
## rBA g1 a6 t7 0.0667816 1 6 7 6 6
## rBA g1 a7 t8 0.0667816 1 7 8 7 7
## rBA g1 a8 t9 0.0667816 1 8 9 8 8
## rBA g1 a9 t10 0.0667816 1 9 10 9 9
## rBA g1 a10 t11 0.0667816 1 10 11 10 10
## rBA g1 a11 t12 0.0667816 1 11 12 11 11
## rBA g1 a12 t13 0.0667816 1 12 13 12 12
## rBA g1 a13 t14 0.0667816 1 13 14 13 13
## rBA g1 a14 t15 0.0667816 1 14 15 14 14
## rBA g1 a15 t16 0.0667816 1 15 16 15 15
## rBA g1 a16 t17 0.0667816 1 16 17 16 16
## rBA g1 a17 t18 0.0667816 1 17 18 17 17
get.real(model.fits[[model.number]], "rBa", se=TRUE)
## all.diff.index par.index estimate se lcl
## rBa g1 a0 t1 76 7 0.0534654 0.006108 0.042683
## rBa g1 a1 t2 77 7 0.0534654 0.006108 0.042683
## rBa g1 a2 t3 78 7 0.0534654 0.006108 0.042683
## rBa g1 a3 t4 79 7 0.0534654 0.006108 0.042683
## rBa g1 a4 t5 80 7 0.0534654 0.006108 0.042683
## rBa g1 a5 t6 81 7 0.0534654 0.006108 0.042683
## rBa g1 a6 t7 82 7 0.0534654 0.006108 0.042683
## rBa g1 a7 t8 83 7 0.0534654 0.006108 0.042683
## rBa g1 a8 t9 84 7 0.0534654 0.006108 0.042683
## rBa g1 a9 t10 85 7 0.0534654 0.006108 0.042683
## rBa g1 a10 t11 86 7 0.0534654 0.006108 0.042683
## rBa g1 a11 t12 87 7 0.0534654 0.006108 0.042683
## rBa g1 a12 t13 88 7 0.0534654 0.006108 0.042683
## rBa g1 a13 t14 89 7 0.0534654 0.006108 0.042683
## rBa g1 a14 t15 90 7 0.0534654 0.006108 0.042683
## rBa g1 a15 t16 91 7 0.0534654 0.006108 0.042683
## rBa g1 a16 t17 92 7 0.0534654 0.006108 0.042683
## rBa g1 a17 t18 93 7 0.0534654 0.006108 0.042683
## ucl fixed note group age time Age Time
## rBa g1 a0 t1 0.0667816 1 0 1 0 0
## rBa g1 a1 t2 0.0667816 1 1 2 1 1
## rBa g1 a2 t3 0.0667816 1 2 3 2 2
## rBa g1 a3 t4 0.0667816 1 3 4 3 3
## rBa g1 a4 t5 0.0667816 1 4 5 4 4
## rBa g1 a5 t6 0.0667816 1 5 6 5 5
## rBa g1 a6 t7 0.0667816 1 6 7 6 6
## rBa g1 a7 t8 0.0667816 1 7 8 7 7
## rBa g1 a8 t9 0.0667816 1 8 9 8 8
## rBa g1 a9 t10 0.0667816 1 9 10 9 9
## rBa g1 a10 t11 0.0667816 1 10 11 10 10
## rBa g1 a11 t12 0.0667816 1 11 12 11 11
## rBa g1 a12 t13 0.0667816 1 12 13 12 12
## rBa g1 a13 t14 0.0667816 1 13 14 13 13
## rBa g1 a14 t15 0.0667816 1 14 15 14 14
## rBa g1 a15 t16 0.0667816 1 15 16 15 15
## rBa g1 a16 t17 0.0667816 1 16 17 16 16
## rBa g1 a17 t18 0.0667816 1 17 18 17 17
# collect models and make AICc table
model.set <- RMark::collect.models( type="2SpecConOccup")
model.set
## model
## 1 PsiA(~elev_std)PsiBA(~elev_std)PsiBa(~elev_std)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
## 3 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
## 2 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
## 4 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
## 5 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
## npar AICc DeltaAICc weight Deviance
## 1 10 2381.870 0.00000 9.999987e-01 2360.583
## 3 7 2409.463 27.59293 1.019228e-06 2394.819
## 2 8 2412.002 30.13202 2.863613e-07 2395.169
## 4 7 2436.922 55.05233 1.110550e-12 2422.278
## 5 6 2444.229 62.35915 2.876617e-14 2431.749
names(model.set)
## [1] "m...001" "m...002" "m...003" "m...004" "m...005"
## [6] "model.table"
model.set$model.table
## PsiA PsiBA PsiBa pA pB rA rBA rBa
## 1 ~elev_std ~elev_std ~elev_std ~1 ~1 ~1 ~1
## 3 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## 2 ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## 4 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## 5 ~1 ~1 ~1 ~1 ~1 ~1
## model
## 1 PsiA(~elev_std)PsiBA(~elev_std)PsiBa(~elev_std)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
## 3 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
## 2 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
## 4 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
## 5 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
## npar AICc DeltaAICc weight Deviance
## 1 10 2381.870 0.00000 9.999987e-01 2360.583
## 3 7 2409.463 27.59293 1.019228e-06 2394.819
## 2 8 2412.002 30.13202 2.863613e-07 2395.169
## 4 7 2436.922 55.05233 1.110550e-12 2422.278
## 5 6 2444.229 62.35915 2.876617e-14 2431.749
# model averaged values at average covariate value
PsiA.ma <- RMark::model.average(model.set, param="PsiA")
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
PsiA.ma
## par.index estimate se fixed note group age time
## PsiA g1 a0 t1 1 0.6207579 0.05515463 1 0 1
## Age Time
## PsiA g1 a0 t1 0 0
# make predictions as a function of distance
# We need to know the indices to the parameters
# Look at the all.diff.index values to get thei correct index number
get.real(model.fits[[model.number]], "PsiA", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiA g1 a0 t1 1 1 0.6207579 0.0551546 0.5083915
## ucl fixed note group age time Age Time
## PsiA g1 a0 t1 0.7215098 1 0 1 0 0
PsiA.ma.p <- RMark::covariate.predictions(model.set, indices=1, data=input.history)$estimates
PsiA.ma.p$parameter <- 'PsiA'
get.real(model.fits[[model.number]], "PsiBa", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiBa g1 a0 t1 3 3 0.1577372 0.1180795 0.0317633
## ucl fixed note group age time Age Time
## PsiBa g1 a0 t1 0.5167042 1 0 1 0 0
PsiBa.ma.p <- RMark::covariate.predictions(model.set, indices=3, data=input.history)$estimates
PsiBa.ma.p$parameter <- 'PsiBa'
get.real(model.fits[[model.number]], "PsiBA", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiBA g1 a0 t1 2 2 0.988917 0.0114862 0.9196152
## ucl fixed note group age time Age Time
## PsiBA g1 a0 t1 0.9985652 1 0 1 0 0
PsiBA.ma.p <- RMark::covariate.predictions(model.set, indices=2, data=input.history)$estimates
PsiBA.ma.p$parameter <- 'PsiBA'
# Not easy to get the se of dervied parameters at the covariate values
PsiB.ma.p <- data.frame(estimate=PsiBA.ma.p$estimate*PsiA.ma.p$estimate +
PsiBa.ma.p$estimate*(1-PsiA.ma.p$estimate),
elev=PsiA.ma.p$elev, stringsAsFactors=FALSE)
PsiB.ma.p$parameter <- "PsiB"
plotdata <- plyr::rbind.fill(PsiA.ma.p,
PsiBA.ma.p,
PsiBa.ma.p,
PsiB.ma.p)
# It is clear that the fit is not satisfactory with a huge se for PsiBa and very steep decline.
# Try reversing the species and you get a much better fit
ggplot(data=plotdata, aes(x=elev, y=estimate))+
ggtitle("Effects of elevation on species co-occurance")+
geom_point()+
geom_ribbon(aes(ymin=lcl, ymax=ucl), alpha=0.2)+
facet_wrap(~parameter, ncol=2)
# 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.
RMark::model.average(model.set, param="pA")
##
## Model 3dropped from the model averaging because one or more beta variances are not positive
##
## Model 5dropped from the model averaging because one or more beta variances are not positive
## par.index estimate se fixed note group age
## pA g1 a0 t1 4 0.9999997 0.0003653775 1 0
## pA g1 a1 t2 5 0.9999997 0.0003653775 1 1
## pA g1 a2 t3 6 0.9999997 0.0003653775 1 2
## pA g1 a3 t4 7 0.9999997 0.0003653775 1 3
## pA g1 a4 t5 8 0.9999997 0.0003653775 1 4
## pA g1 a5 t6 9 0.9999997 0.0003653775 1 5
## pA g1 a6 t7 10 0.9999997 0.0003653775 1 6
## pA g1 a7 t8 11 0.9999997 0.0003653775 1 7
## pA g1 a8 t9 12 0.9999997 0.0003653775 1 8
## pA g1 a9 t10 13 0.9999997 0.0003653775 1 9
## pA g1 a10 t11 14 0.9999997 0.0003653775 1 10
## pA g1 a11 t12 15 0.9999997 0.0003653775 1 11
## pA g1 a12 t13 16 0.9999997 0.0003653775 1 12
## pA g1 a13 t14 17 0.9999997 0.0003653775 1 13
## pA g1 a14 t15 18 0.9999997 0.0003653775 1 14
## pA g1 a15 t16 19 0.9999997 0.0003653775 1 15
## pA g1 a16 t17 20 0.9999997 0.0003653775 1 16
## pA g1 a17 t18 21 0.9999997 0.0003653775 1 17
## time Age Time
## pA g1 a0 t1 1 0 0
## pA g1 a1 t2 2 1 1
## pA g1 a2 t3 3 2 2
## pA g1 a3 t4 4 3 3
## pA g1 a4 t5 5 4 4
## pA g1 a5 t6 6 5 5
## pA g1 a6 t7 7 6 6
## pA g1 a7 t8 8 7 7
## pA g1 a8 t9 9 8 8
## pA g1 a9 t10 10 9 9
## pA g1 a10 t11 11 10 10
## pA g1 a11 t12 12 11 11
## pA g1 a12 t13 13 12 12
## pA g1 a13 t14 14 13 13
## pA g1 a14 t15 15 14 14
## pA g1 a15 t16 16 15 15
## pA g1 a16 t17 17 16 16
## pA g1 a17 t18 18 17 17
RMark::model.average(model.set, param="pB")
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## par.index estimate se fixed note group age time
## pB g1 a0 t1 22 0.1965643 0.03004401 1 0 1
## pB g1 a1 t2 23 0.1965643 0.03004401 1 1 2
## pB g1 a2 t3 24 0.1965643 0.03004401 1 2 3
## pB g1 a3 t4 25 0.1965643 0.03004401 1 3 4
## pB g1 a4 t5 26 0.1965643 0.03004401 1 4 5
## pB g1 a5 t6 27 0.1965643 0.03004401 1 5 6
## pB g1 a6 t7 28 0.1965643 0.03004401 1 6 7
## pB g1 a7 t8 29 0.1965643 0.03004401 1 7 8
## pB g1 a8 t9 30 0.1965643 0.03004401 1 8 9
## pB g1 a9 t10 31 0.1965643 0.03004401 1 9 10
## pB g1 a10 t11 32 0.1965643 0.03004401 1 10 11
## pB g1 a11 t12 33 0.1965643 0.03004401 1 11 12
## pB g1 a12 t13 34 0.1965643 0.03004401 1 12 13
## pB g1 a13 t14 35 0.1965643 0.03004401 1 13 14
## pB g1 a14 t15 36 0.1965643 0.03004401 1 14 15
## pB g1 a15 t16 37 0.1965643 0.03004401 1 15 16
## pB g1 a16 t17 38 0.1965643 0.03004401 1 16 17
## pB g1 a17 t18 39 0.1965643 0.03004401 1 17 18
## Age Time
## pB g1 a0 t1 0 0
## pB g1 a1 t2 1 1
## pB g1 a2 t3 2 2
## pB g1 a3 t4 3 3
## pB g1 a4 t5 4 4
## pB g1 a5 t6 5 5
## pB g1 a6 t7 6 6
## pB g1 a7 t8 7 7
## pB g1 a8 t9 8 8
## pB g1 a9 t10 9 9
## pB g1 a10 t11 10 10
## pB g1 a11 t12 11 11
## pB g1 a12 t13 12 12
## pB g1 a13 t14 13 13
## pB g1 a14 t15 14 14
## pB g1 a15 t16 15 15
## pB g1 a16 t17 16 16
## pB g1 a17 t18 17 17
RMark::model.average(model.set, param="rA")
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## par.index estimate se fixed note group age time
## rA g1 a0 t1 40 0.100413 0.0104089 1 0 1
## rA g1 a1 t2 41 0.100413 0.0104089 1 1 2
## rA g1 a2 t3 42 0.100413 0.0104089 1 2 3
## rA g1 a3 t4 43 0.100413 0.0104089 1 3 4
## rA g1 a4 t5 44 0.100413 0.0104089 1 4 5
## rA g1 a5 t6 45 0.100413 0.0104089 1 5 6
## rA g1 a6 t7 46 0.100413 0.0104089 1 6 7
## rA g1 a7 t8 47 0.100413 0.0104089 1 7 8
## rA g1 a8 t9 48 0.100413 0.0104089 1 8 9
## rA g1 a9 t10 49 0.100413 0.0104089 1 9 10
## rA g1 a10 t11 50 0.100413 0.0104089 1 10 11
## rA g1 a11 t12 51 0.100413 0.0104089 1 11 12
## rA g1 a12 t13 52 0.100413 0.0104089 1 12 13
## rA g1 a13 t14 53 0.100413 0.0104089 1 13 14
## rA g1 a14 t15 54 0.100413 0.0104089 1 14 15
## rA g1 a15 t16 55 0.100413 0.0104089 1 15 16
## rA g1 a16 t17 56 0.100413 0.0104089 1 16 17
## rA g1 a17 t18 57 0.100413 0.0104089 1 17 18
## Age Time
## rA g1 a0 t1 0 0
## rA g1 a1 t2 1 1
## rA g1 a2 t3 2 2
## rA g1 a3 t4 3 3
## rA g1 a4 t5 4 4
## rA g1 a5 t6 5 5
## rA g1 a6 t7 6 6
## rA g1 a7 t8 7 7
## rA g1 a8 t9 8 8
## rA g1 a9 t10 9 9
## rA g1 a10 t11 10 10
## rA g1 a11 t12 11 11
## rA g1 a12 t13 12 12
## rA g1 a13 t14 13 13
## rA g1 a14 t15 14 14
## rA g1 a15 t16 15 15
## rA g1 a16 t17 16 16
## rA g1 a17 t18 17 17
RMark::model.average(model.set, param="rBA")
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## par.index estimate se fixed note group age
## rBA g1 a0 t1 58 0.05346544 0.006107968 1 0
## rBA g1 a1 t2 59 0.05346544 0.006107968 1 1
## rBA g1 a2 t3 60 0.05346544 0.006107968 1 2
## rBA g1 a3 t4 61 0.05346544 0.006107968 1 3
## rBA g1 a4 t5 62 0.05346544 0.006107968 1 4
## rBA g1 a5 t6 63 0.05346544 0.006107968 1 5
## rBA g1 a6 t7 64 0.05346544 0.006107968 1 6
## rBA g1 a7 t8 65 0.05346544 0.006107968 1 7
## rBA g1 a8 t9 66 0.05346544 0.006107968 1 8
## rBA g1 a9 t10 67 0.05346544 0.006107968 1 9
## rBA g1 a10 t11 68 0.05346544 0.006107968 1 10
## rBA g1 a11 t12 69 0.05346544 0.006107968 1 11
## rBA g1 a12 t13 70 0.05346544 0.006107968 1 12
## rBA g1 a13 t14 71 0.05346544 0.006107968 1 13
## rBA g1 a14 t15 72 0.05346544 0.006107968 1 14
## rBA g1 a15 t16 73 0.05346544 0.006107968 1 15
## rBA g1 a16 t17 74 0.05346544 0.006107968 1 16
## rBA g1 a17 t18 75 0.05346544 0.006107968 1 17
## time Age Time
## rBA g1 a0 t1 1 0 0
## rBA g1 a1 t2 2 1 1
## rBA g1 a2 t3 3 2 2
## rBA g1 a3 t4 4 3 3
## rBA g1 a4 t5 5 4 4
## rBA g1 a5 t6 6 5 5
## rBA g1 a6 t7 7 6 6
## rBA g1 a7 t8 8 7 7
## rBA g1 a8 t9 9 8 8
## rBA g1 a9 t10 10 9 9
## rBA g1 a10 t11 11 10 10
## rBA g1 a11 t12 12 11 11
## rBA g1 a12 t13 13 12 12
## rBA g1 a13 t14 14 13 13
## rBA g1 a14 t15 15 14 14
## rBA g1 a15 t16 16 15 15
## rBA g1 a16 t17 17 16 16
## rBA g1 a17 t18 18 17 17
RMark::model.average(model.set, param="rBa")
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## par.index estimate se fixed note group age
## rBa g1 a0 t1 76 0.05346544 0.006107959 1 0
## rBa g1 a1 t2 77 0.05346544 0.006107959 1 1
## rBa g1 a2 t3 78 0.05346544 0.006107959 1 2
## rBa g1 a3 t4 79 0.05346544 0.006107959 1 3
## rBa g1 a4 t5 80 0.05346544 0.006107959 1 4
## rBa g1 a5 t6 81 0.05346544 0.006107959 1 5
## rBa g1 a6 t7 82 0.05346544 0.006107959 1 6
## rBa g1 a7 t8 83 0.05346544 0.006107959 1 7
## rBa g1 a8 t9 84 0.05346544 0.006107959 1 8
## rBa g1 a9 t10 85 0.05346544 0.006107959 1 9
## rBa g1 a10 t11 86 0.05346544 0.006107959 1 10
## rBa g1 a11 t12 87 0.05346544 0.006107959 1 11
## rBa g1 a12 t13 88 0.05346544 0.006107959 1 12
## rBa g1 a13 t14 89 0.05346544 0.006107959 1 13
## rBa g1 a14 t15 90 0.05346544 0.006107959 1 14
## rBa g1 a15 t16 91 0.05346544 0.006107959 1 15
## rBa g1 a16 t17 92 0.05346544 0.006107959 1 16
## rBa g1 a17 t18 93 0.05346544 0.006107959 1 17
## time Age Time
## rBa g1 a0 t1 1 0 0
## rBa g1 a1 t2 2 1 1
## rBa g1 a2 t3 3 2 2
## rBa g1 a3 t4 4 3 3
## rBa g1 a4 t5 5 4 4
## rBa g1 a5 t6 6 5 5
## rBa g1 a6 t7 7 6 6
## rBa g1 a7 t8 8 7 7
## rBa g1 a8 t9 9 8 8
## rBa g1 a9 t10 10 9 9
## rBa g1 a10 t11 11 10 10
## rBa g1 a11 t12 12 11 11
## rBa g1 a12 t13 13 12 12
## rBa g1 a13 t14 14 13 13
## rBa g1 a14 t15 15 14 14
## rBa g1 a15 t16 16 15 15
## rBa g1 a16 t17 17 16 16
## rBa g1 a17 t18 18 17 17
# Compute the observed occupance of each species
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)
# 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)
# cleanup
cleanup(ask=FALSE)