# Multi Species, Single Season Occupancy analyais
# Cove MV, Gardner B, Simons TR, O'Connell AF (2018)
# Co-occurrence dynamics of endangered Lower Keys marsh rabbits and free-ranging domestic cats:
# Prey responses to an exotic predator removal program.
# Ecology and Evolution 8(8): 4042-4052.
# https://doi.org/10.1002/ece3.3954
# Data obtained from the Dryad data package:
#
# Cove MV, Gardner B, Simons TR, O'Connell AF (2018)
# Data from: Co-occurrence dynamics of endangered Lower Keys marsh rabbits and free-ranging domestic cats:
# prey responses to an exotic predator removal program.
# Dryad Digital Repository.
# https://doi.org/10.5061/dryad.748pd64
#
# Note that the orginal data on the Dryad repository had incorrect cat numbers for 2013/2014 and 2015 data.
# I wrote to the author and got the corrected data.
# Briefly: between 2013 and 2015, camera traps were used to survey marsh
# rabbits and free-ranging
# cats at 84 sites in the National Key Deer Refuge, Big Pine
# Key, Florida, USA. We used dynamic occupancy models to determine factors associated
# with marsh rabbit occurrence, colonization, extinction, and the co-occurrence
# of marsh rabbits and cats during a period of predator removal
# ***** Only the first "year" of data (2013/2014) prior to any manipulation of the number of cats **** is used
# 2020-01-07 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
# load libraries
library(readxl)
library(RMark)
## This is RMark 2.2.7
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
library(ggplot2)
# Get the Mark additional functions
source(file.path("..","..","..","..","AdditionalFunctions","RMark.additional.functions.R"))
# get the data read in
# Data for detections should be a data frame with rows corresponding to sites
# and columns to visits.
# The usual 1=detected; 0=not detected; NA=not visited is used.
input.data <- read.csv(file.path("..","..","Cove_etal_Ecology_Evolution_2018_data-original.csv"),
header=TRUE, as.is=TRUE, strip.white=TRUE, na.string=".")
head(input.data)
## Site UTM.N UTM.E X2015.Start.Date X2015.End.Date Rabbits.2013.2014 X
## 1 Fs003 2735114 460522.3 3/28/15 4/12/15 0 0
## 2 Fs005 2733915 460235.5 3/26/15 4/10/15 0 0
## 3 Fs006 2733787 460643.2 4/30/15 5/15/15 0 1
## 4 Fs007 2733456 460429.5 3/27/15 4/11/15 1 0
## 5 Fs008 2733476 461280.8 4/30/15 5/15/15 0 1
## 6 Fs009 2733213 460392.6 3/27/15 4/11/15 0 0
## X.1 X.2 X.3 X.4 X.5 X.6 X.7 X.8 X.9 X.10 X.11 X.12 X.13 X.14 Rabbits.2015
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 1 0 1 0 0 1 0 0 0 1 0 0
## 3 0 0 0 0 0 0 0 1 0 1 1 1 0 1 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## X.15 X.16 X.17 X.18 X.19 X.20 X.21 X.22 X.23 X.24 X.25 X.26 X.27 X.28 X.29
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 NA
## 2 0 1 0 1 1 1 1 0 1 0 1 1 0 1 NA
## 3 0 0 0 1 1 1 0 0 0 1 0 0 1 0 NA
## 4 1 1 1 1 1 1 0 0 0 0 1 0 0 0 NA
## 5 0 0 0 0 0 0 1 0 0 0 0 0 0 0 NA
## 6 0 0 1 0 0 0 1 0 1 1 0 1 0 0 NA
## X.30 Cats.2013.2014 X.31 X.32 X.33 X.34 X.35 X.36 X.37 X.38 X.39 X.40 X.41
## 1 FS003 0 0 0 0 0 0 0 0 0 0 0 0
## 2 FS005 0 0 0 0 0 0 0 0 0 0 1 0
## 3 FS006 0 0 0 0 0 0 0 0 0 0 0 0
## 4 FS007 0 0 0 0 0 0 0 0 0 0 0 0
## 5 FS008 0 0 0 0 0 0 0 0 0 0 0 0
## 6 FS009 0 0 0 0 0 0 0 0 0 0 0 0
## X.42 X.43 X.44 X.45 Cats.2015 X.46 X.47 X.48 X.49 X.50 X.51 X.52 X.53 X.54
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## X.55 X.56 X.57 X.58 X.59 X.60 year_2013 dist_cam_z Total.CatCaps_2013_14
## 1 0 0 0 0 0 NA 1 0.86 0
## 2 0 0 0 0 1 NA 1 -1.49 1
## 3 0 0 0 0 0 NA 1 NA 0
## 4 0 0 1 0 0 NA 1 -0.80 1
## 5 0 0 0 0 0 NA 1 NA 0
## 6 0 0 0 0 0 NA 1 -0.32 1
## CatIndividuals_2013_14 Total.CatCaps_2015 CatIndividuals_2015
## 1 0 0 0
## 2 1 1 1
## 3 0 0 0
## 4 1 1 1
## 5 0 0 0
## 6 1 1 1
## RabbitPresence_2013_14 HumanTrail Rabbitat Rabbitat_buff HumanTrail.1
## 1 0 0 1 0 0
## 2 1 0 1 0 0
## 3 1 0 1 0 0
## 4 1 0 1 0 0
## 5 1 1 1 1 0
## 6 1 1 0 0 0
## dist_rd_z dist_dev_z dist_cat2014_z dist_cat2015_z bin_2014cat bin_2015cat
## 1 1.51 1.05 0.51 2.79 0 0
## 2 -0.46 -0.83 -1.60 1.73 1 0
## 3 -0.08 0.12 -0.88 1.21 1 0
## 4 0.71 0.00 -0.75 1.14 1 0
## 5 0.78 2.20 0.34 0.32 0 0
## 6 1.45 1.02 -0.27 0.88 0 0
## bin_20142015 DESCRIPTION freshwater salt_button upland hectares_z
## 1 0 Buttonwood 0 1 0 -0.74
## 2 1 Freshwater Wetland 1 0 0 -0.78
## 3 1 Freshwater Wetland 1 0 0 -0.78
## 4 1 Freshwater Wetland 1 0 0 -0.78
## 5 0 Pineland 0 0 1 0.51
## 6 0 Pineland 0 0 1 1.84
## LiDar_z
## 1 -0.35
## 2 -0.88
## 3 -1.13
## 4 -0.84
## 5 -0.32
## 6 -0.07
# OOPs...# compute the difference in individual cats. These are problematic as the values are the same
# in the two sessions
input.data$CatIndividuals_2013_14
## [1] 0 1 0 1 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 2 0 0 0 0 1 2 0 0 0
## [39] 0 0 0 0 1 0 2 1 1 2 1 2 3 2 2 2 1 1 0 1 1 2 0 0 0 2 0 2 1 0 7 0 4 1 1 2 0 2
## [77] 1 0 0 1 0 0 0 2
input.data$CatIndividuals_2015
## [1] 0 1 0 1 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 2 0 0 0 0 1 2 0 0 0
## [39] 0 0 0 0 1 0 2 1 1 2 1 2 3 2 2 2 1 1 0 1 1 2 0 0 0 2 0 2 1 0 7 0 4 1 1 2 0 2
## [77] 1 0 0 1 0 0 0 2
input.data$delta.icats <- input.data$CatIndividuals_2013_14 - input.data$CatIndividuals_2015
input.data$delta.icats
## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [77] 0 0 0 0 0 0 0 0
input.data$Total.CatCaps_2013_14
## [1] 0 1 0 1 0 1 1 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [26] 0 1 0 7 0 0 0 0 1 2 0 0 0 0 0 0 0 1 0 2 1 3 3 1 3
## [51] 4 3 2 5 3 1 0 3 1 2 0 0 0 6 0 7 2 0 16 0 10 1 5 15 0
## [76] 8 1 0 0 1 0 0 0 5
input.data$Total.CatCaps_2015
## [1] 0 1 0 1 0 1 1 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [26] 0 1 0 7 0 0 0 0 1 2 0 0 0 0 0 0 0 1 0 2 1 3 3 1 3
## [51] 4 3 2 5 3 1 0 3 1 2 0 0 0 6 0 7 2 0 16 0 10 1 5 15 0
## [76] 8 1 1 0 1 0 0 0 5
# here is the updated cat data
input.data2 <- read.csv(file.path("..","..","Cove_etal_Ecology_Evolution_2018_data-updated.csv"),
header=TRUE, as.is=TRUE, strip.white=TRUE, na.string=".")
head(input.data2)
## year_2013 North dist_cam_z Catsyear1 Capturesyear1 Individualsyear1
## 1 1 1 0.8568465 0 0 0
## 2 1 1 -1.4944675 1 2 2
## 3 1 1 0.0000000 0 0 0
## 4 1 1 -0.8029046 0 0 0
## 5 1 1 0.0000000 0 0 0
## 6 1 1 -0.3188105 0 0 0
## Capturesyear2 Individualsyear2 twoyear_indchange twoyear_capschange
## 1 0 0 0 0
## 2 1 1 -1 -1
## 3 0 0 0 0
## 4 1 1 1 1
## 5 0 0 0 0
## 6 1 1 1 1
## capchange_z RabbitPresence HumanTrail Rabbitat Rabbitat_buff dist_rd_z
## 1 0.30246914 0 0 1 0 1.50897691
## 2 -0.06063907 1 0 1 0 -0.46101737
## 3 0.30246914 1 0 1 0 -0.07882644
## 4 0.66557734 1 0 1 0 0.70617147
## 5 0.30246914 1 1 1 1 0.78310282
## 6 0.66557734 1 1 0 0 1.44628482
## dist_dev_z dist_cat2014_z dist_cat2015_z bin_2014cat bin_2015cat
## 1 1.046661022 0.5141839 2.7878397 0 0
## 2 -0.834653866 -1.5959094 1.7330113 1 0
## 3 0.115757425 -0.8770415 1.2084800 1 0
## 4 -0.000375798 -0.7457188 1.1412177 1 0
## 5 2.200931956 0.3372080 0.3171738 0 0
## 6 1.020459541 -0.2725414 0.8839006 0 0
## bin_20142015 DESCRIPTION freshwater salt_button upland hectares_z
## 1 0 Buttonwood 0 1 0 -0.7414135
## 2 1 Freshwater Wetland 1 0 0 -0.7808261
## 3 1 Freshwater Wetland 1 0 0 -0.7808261
## 4 1 Freshwater Wetland 1 0 0 -0.7808261
## 5 0 Pineland 0 0 1 0.5121787
## 6 0 Pineland 0 0 1 1.8352219
## LiDar_z Sum_Cats Min_Cats Max_Cats Sum_Captur Min_Captur Max_Captur
## 1 -0.35022907 0 0 0 0 0 0
## 2 -0.88143611 2 1 1 6 2 4
## 3 -1.12883639 2 0 1 8 0 4
## 4 -0.84022781 0 0 0 0 0 0
## 5 -0.31962944 0 0 0 0 0 0
## 6 -0.06844092 0 0 0 0 0 0
## Sum_Indivi
## 1 0
## 2 4
## 3 5
## 4 0
## 5 0
## 6 0
input.data2 <- input.data2[,c("Individualsyear1","Individualsyear2","twoyear_indchange","RabbitPresence")]
input.data <- cbind(input.data, input.data2) # the second files doesn't have site ids but appears to be sorted in correct order
input.data$CatIndividuals_2013_14 <- NULL
input.data$Total.CatCaps_2015 <- NULL
input.data$Individualsyear1
## [1] 0 2 0 0 0 0 0 2 0 0 2 0 1 2 0 2 1 0 1 0 0 1 3 2 4 1 0 0 2 0 0 1 1 0 1 4 2 0
## [39] 0 2 0 1 1 1 2 2 0 2 2 1 2 2 1 4 2 4 1 1 0 0 1 3 0 2 1 1 1 1 5 3 3 0 3 3 0 3
## [77] 1 0 1 4 0 0 1 3
input.data$Individualsyear2
## [1] 0 1 0 1 0 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 2 0 0 0 0 1 2 0 0 0
## [39] 0 0 0 0 1 0 2 1 1 2 1 2 3 2 2 2 1 1 0 1 1 2 0 0 0 2 0 2 1 0 7 0 4 1 1 2 0 2
## [77] 1 0 0 1 0 0 0 2
input.data$twoyear_indchange
## [1] 0 -1 0 1 0 1 1 -2 0 0 -1 0 -1 -2 0 -2 -1 0 -1 0 0 -1 -3 -2 -4
## [26] -1 1 0 0 0 0 -1 -1 1 1 -4 -2 0 0 -2 0 -1 0 -1 0 -1 1 0 -1 1
## [51] 1 0 1 -2 -1 -3 -1 0 1 2 -1 -3 0 0 -1 1 0 -1 2 -3 1 1 -2 -1 0
## [76] -1 0 0 -1 -3 0 0 -1 -1
# Groan, covariate names must be 10 characters or less (a limitation in RMark)
# we need to rename some of the variables
input.data <- plyr::rename(input.data, c("salt_button" ="salt_butt",
"Individualsyear1" ="icats1",
"dist_cat2014_z" ="dist2014z",
"bin_2014cat" ="icat2014",
"twoyear_indchange" ="delta_cats"))
names(input.data)
## [1] "Site" "UTM.N" "UTM.E"
## [4] "X2015.Start.Date" "X2015.End.Date" "Rabbits.2013.2014"
## [7] "X" "X.1" "X.2"
## [10] "X.3" "X.4" "X.5"
## [13] "X.6" "X.7" "X.8"
## [16] "X.9" "X.10" "X.11"
## [19] "X.12" "X.13" "X.14"
## [22] "Rabbits.2015" "X.15" "X.16"
## [25] "X.17" "X.18" "X.19"
## [28] "X.20" "X.21" "X.22"
## [31] "X.23" "X.24" "X.25"
## [34] "X.26" "X.27" "X.28"
## [37] "X.29" "X.30" "Cats.2013.2014"
## [40] "X.31" "X.32" "X.33"
## [43] "X.34" "X.35" "X.36"
## [46] "X.37" "X.38" "X.39"
## [49] "X.40" "X.41" "X.42"
## [52] "X.43" "X.44" "X.45"
## [55] "Cats.2015" "X.46" "X.47"
## [58] "X.48" "X.49" "X.50"
## [61] "X.51" "X.52" "X.53"
## [64] "X.54" "X.55" "X.56"
## [67] "X.57" "X.58" "X.59"
## [70] "X.60" "year_2013" "dist_cam_z"
## [73] "Total.CatCaps_2013_14" "CatIndividuals_2015" "RabbitPresence_2013_14"
## [76] "HumanTrail" "Rabbitat" "Rabbitat_buff"
## [79] "HumanTrail.1" "dist_rd_z" "dist_dev_z"
## [82] "dist2014z" "dist_cat2015_z" "icat2014"
## [85] "bin_2015cat" "bin_20142015" "DESCRIPTION"
## [88] "freshwater" "salt_butt" "upland"
## [91] "hectares_z" "LiDar_z" "delta.icats"
## [94] "icats1" "Individualsyear2" "delta_cats"
## [97] "RabbitPresence"
# the paper collapses the different habitats into 3 categories, and creates indicator variables
unique(input.data$DESCRIPTION)
## [1] "Buttonwood" "Freshwater Wetland" "Pineland"
## [4] "Hammock" "Scrub Mangrove" "Exotic"
## [7] "Salt Marsh"
xtabs(~DESCRIPTION+salt_butt, data=input.data,exclude=NULL, na.action=na.pass)
## salt_butt
## DESCRIPTION 0 1
## Buttonwood 0 7
## Exotic 1 0
## Freshwater Wetland 27 0
## Hammock 6 0
## Pineland 41 0
## Salt Marsh 0 1
## Scrub Mangrove 0 1
xtabs(~DESCRIPTION+freshwater, data=input.data,exclude=NULL, na.action=na.pass)
## freshwater
## DESCRIPTION 0 1
## Buttonwood 7 0
## Exotic 1 0
## Freshwater Wetland 0 27
## Hammock 6 0
## Pineland 41 0
## Salt Marsh 1 0
## Scrub Mangrove 1 0
xtabs(~DESCRIPTION+upland, data=input.data,exclude=NULL, na.action=na.pass)
## upland
## DESCRIPTION 0 1
## Buttonwood 7 0
## Exotic 0 1
## Freshwater Wetland 27 0
## Hammock 0 6
## Pineland 0 41
## Salt Marsh 1 0
## Scrub Mangrove 1 0
# distance to sites with cat removed in 2014 vs binary indicator
input.data$icat2014
## [1] 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 1 1 1
## [39] 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 1 1 0 1 1 1 1 1 0 0
## [77] 0 0 1 1 1 0 0 1
input.data$dist2014z
## [1] 0.51 -1.60 -0.88 -0.75 0.34 -0.27 1.31 0.92 0.41 1.20 0.67 0.46
## [13] 0.37 0.18 0.55 0.91 0.11 0.44 1.42 -0.45 1.49 1.47 -0.90 -1.29
## [25] 0.29 1.02 1.05 1.38 1.17 0.11 -0.26 -1.14 -0.66 -0.19 -0.39 -1.57
## [37] -1.15 -1.14 -0.62 -1.52 -1.06 0.07 -0.80 -1.14 -0.62 0.61 -0.23 -0.04
## [49] 1.51 1.53 1.45 1.10 0.69 2.20 1.83 0.27 0.27 0.58 0.55 -0.32
## [61] 0.07 -1.22 -0.71 -0.36 0.19 -1.33 -0.83 -0.73 -0.49 -1.09 -1.40 -0.88
## [73] -1.09 -0.99 0.62 1.61 1.06 -0.09 -0.60 -1.12 -1.20 -0.09 -0.16 -1.18
ggplot(data=input.data, aes(x=dist2014z, y=icat2014))+
ggtitle("indicatior if cats captures within 500m buffer of camera trap")+
geom_point()
# The author defines a set of global covariates that are used
global.covar.vars <- c("HumanTrail","Rabbitat","dist_dev_z","freshwater","salt_butt","hectares_z","LiDar_z")
setdiff(global.covar.vars, names(input.data))
## character(0)
#####################################################################
#####################################################################
#####################################################################
# First model occupancy for the rabbits in the first year
# See Table 1 of the paper for the results from a dynamic occupancy
rabbit.detect.vars <- c("Rabbits.2013.2014","X", paste0("X.",1:14))
rabbit.input.history <- input.data[, rabbit.detect.vars]
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(rabbit.input.history)
## [1] 84
ncol(rabbit.input.history)
## [1] 16
range(rabbit.input.history, na.rm=TRUE)
## [1] 0 1
sum(is.na(rabbit.input.history))
## [1] 52
sum(rabbit.input.history, na.rm=TRUE)
## [1] 58
other.covar.vars <- c("RabbitPresence_2013_14", # were rabbits were detected
"icats1", # number of individual cats in year 1
"dist2014z", # distance to site with a cat in 2014
"icat2014") # indicator if cat captured within 500m of the site in 2014
setdiff(other.covar.vars, names(input.data))
## character(0)
site.covar <- input.data[, c(global.covar.vars,other.covar.vars)]
row.names(site.covar) <- NULL
head(site.covar)
## HumanTrail Rabbitat dist_dev_z freshwater salt_butt hectares_z LiDar_z
## 1 0 1 1.05 0 1 -0.74 -0.35
## 2 0 1 -0.83 1 0 -0.78 -0.88
## 3 0 1 0.12 1 0 -0.78 -1.13
## 4 0 1 0.00 1 0 -0.78 -0.84
## 5 1 1 2.20 0 0 0.51 -0.32
## 6 1 0 1.02 0 0 1.84 -0.07
## RabbitPresence_2013_14 icats1 dist2014z icat2014
## 1 0 0 0.51 0
## 2 1 2 -1.60 1
## 3 1 0 -0.88 1
## 4 1 0 -0.75 1
## 5 1 0 0.34 0
## 6 1 0 -0.27 0
#Format the capture history to be used by RMark and deal with missing values
input.history <- data.frame(freq=1,
ch=apply(rabbit.input.history,1,paste, collapse=""), stringsAsFactors=FALSE)
head(input.history,n=20)
## freq ch
## 1 1 0000000000000000
## 2 1 0000010100100010
## 3 1 0100000001011101
## 4 1 1000000000000000
## 5 1 0100000000000000
## 6 1 0010000000001000
## 7 1 0000000010000011
## 8 1 0000000000000000
## 9 1 0000000001000001
## 10 1 0000000000000000
## 11 1 0000000000000000
## 12 1 0000000000000000
## 13 1 0000000000000000
## 14 1 0000000000000000
## 15 1 1000000000000000
## 16 1 0000000000000000
## 17 1 0001100000000000
## 18 1 0000000000000000
## 19 1 0000000000000000
## 20 1 100000000001000NA
input.history$ch <- gsub("NA",".",input.history$ch, fixed=TRUE)
head(input.history,n=20)
## freq ch
## 1 1 0000000000000000
## 2 1 0000010100100010
## 3 1 0100000001011101
## 4 1 1000000000000000
## 5 1 0100000000000000
## 6 1 0010000000001000
## 7 1 0000000010000011
## 8 1 0000000000000000
## 9 1 0000000001000001
## 10 1 0000000000000000
## 11 1 0000000000000000
## 12 1 0000000000000000
## 13 1 0000000000000000
## 14 1 0000000000000000
## 15 1 1000000000000000
## 16 1 0000000000000000
## 17 1 0001100000000000
## 18 1 0000000000000000
## 19 1 0000000000000000
## 20 1 100000000001000.
#Add site covariates to input history
input.history = cbind(input.history,site.covar)
head(input.history)
## freq ch HumanTrail Rabbitat dist_dev_z freshwater salt_butt
## 1 1 0000000000000000 0 1 1.05 0 1
## 2 1 0000010100100010 0 1 -0.83 1 0
## 3 1 0100000001011101 0 1 0.12 1 0
## 4 1 1000000000000000 0 1 0.00 1 0
## 5 1 0100000000000000 1 1 2.20 0 0
## 6 1 0010000000001000 1 0 1.02 0 0
## hectares_z LiDar_z RabbitPresence_2013_14 icats1 dist2014z icat2014
## 1 -0.74 -0.35 0 0 0.51 0
## 2 -0.78 -0.88 1 2 -1.60 1
## 3 -0.78 -1.13 1 0 -0.88 1
## 4 -0.78 -0.84 1 0 -0.75 1
## 5 0.51 -0.32 1 0 0.34 0
## 6 1.84 -0.07 1 0 -0.27 0
# With 16 detections visits, look at the naive occupancy vs. the different covariates
prop.table(xtabs(~RabbitPresence_2013_14+HumanTrail, data=input.history, exclude=NULL, na.action=na.pass), margin=2)
## HumanTrail
## RabbitPresence_2013_14 0 1
## 0 0.5967742 0.5909091
## 1 0.4032258 0.4090909
prop.table(xtabs(~RabbitPresence_2013_14+Rabbitat, data=input.history, exclude=NULL, na.action=na.pass), margin=2)
## Rabbitat
## RabbitPresence_2013_14 0 1
## 0 0.8800000 0.4745763
## 1 0.1200000 0.5254237
prop.table(xtabs(~RabbitPresence_2013_14+freshwater, data=input.history, exclude=NULL, na.action=na.pass), margin=2)
## freshwater
## RabbitPresence_2013_14 0 1
## 0 0.7894737 0.1851852
## 1 0.2105263 0.8148148
prop.table(xtabs(~RabbitPresence_2013_14+salt_butt, data=input.history, exclude=NULL, na.action=na.pass), margin=2)
## salt_butt
## RabbitPresence_2013_14 0 1
## 0 0.5600000 0.8888889
## 1 0.4400000 0.1111111
prop.table(xtabs(~RabbitPresence_2013_14+icats1, data=input.history, exclude=NULL, na.action=na.pass), margin=2)
## icats1
## RabbitPresence_2013_14 0 1 2 3 4
## 0 0.4666667 0.6521739 0.5882353 0.7500000 0.8000000
## 1 0.5333333 0.3478261 0.4117647 0.2500000 0.2000000
## icats1
## RabbitPresence_2013_14 5
## 0 1.0000000
## 1 0.0000000
xtabs(~freshwater+salt_butt, data=input.history, exclude=NULL, na.action=na.pass)
## salt_butt
## freshwater 0 1
## 0 48 9
## 1 27 0
#Create the RMark data structure
# 1 years with 16 camera trapping days/year
max.visit <- 16
rabbit.data <- process.data(data=input.history,
model="Occupancy"
)
rabbit.data
## $data
## freq ch HumanTrail Rabbitat dist_dev_z freshwater salt_butt
## 1 1 0000000000000000 0 1 1.05 0 1
## 2 1 0000010100100010 0 1 -0.83 1 0
## 3 1 0100000001011101 0 1 0.12 1 0
## 4 1 1000000000000000 0 1 0.00 1 0
## 5 1 0100000000000000 1 1 2.20 0 0
## 6 1 0010000000001000 1 0 1.02 0 0
## 7 1 0000000010000011 0 1 -0.35 0 0
## 8 1 0000000000000000 0 0 1.00 0 0
## 9 1 0000000001000001 0 1 0.98 0 0
## 10 1 0000000000000000 0 1 -0.24 0 0
## 11 1 0000000000000000 0 0 2.18 0 0
## 12 1 0000000000000000 1 1 0.50 0 0
## 13 1 0000000000000000 0 1 0.40 0 0
## 14 1 0000000000000000 0 0 -0.34 0 0
## 15 1 1000000000000000 0 1 1.13 0 0
## 16 1 0000000000000000 0 0 0.29 0 0
## 17 1 0001100000000000 0 0 1.61 0 0
## 18 1 0000000000000000 0 0 0.97 0 0
## 19 1 0000000000000000 0 0 -0.38 0 0
## 20 1 100000000001000. 0 1 1.38 1 0
## 21 1 000000000010000. 0 1 0.69 1 0
## 22 1 000000000000000. 0 0 0.21 0 0
## 23 1 0000000000100100 0 1 -1.01 1 0
## 24 1 0000000000010000 0 1 -0.48 1 0
## 25 1 000000000000000. 0 0 -0.08 0 0
## 26 1 0000000000000000 0 1 -0.14 0 0
## 27 1 0000001000000000 0 1 0.06 1 0
## 28 1 0000000100000000 0 0 -0.79 0 0
## 29 1 0000000000000000 1 0 1.59 0 0
## 30 1 000000000000000. 0 1 1.68 0 0
## 31 1 000000000000000. 0 1 2.08 0 1
## 32 1 000000000000000. 0 0 -0.72 0 0
## 33 1 000000000000000. 1 1 -0.54 0 1
## 34 1 100000000000000. 0 1 0.41 1 0
## 35 1 1100101000000000 0 1 -1.13 1 0
## 36 1 000000000000000. 0 1 -0.85 1 0
## 37 1 000000000000000. 0 1 -0.93 1 0
## 38 1 0000100000000001 0 1 -0.95 1 0
## 39 1 000001000000000. 0 1 -0.76 1 0
## 40 1 000010000000001. 0 1 -0.31 1 0
## 41 1 0000000000000000 0 1 -0.01 1 0
## 42 1 0000000000000000 0 0 -1.01 0 0
## 43 1 0000000000000000 0 1 -0.54 1 0
## 44 1 0000000000000000 0 1 -0.96 0 0
## 45 1 0000000000000000 0 0 -0.85 0 0
## 46 1 000000000000000. 0 1 -0.96 0 1
## 47 1 000000000000000. 0 1 -0.49 1 0
## 48 1 000000000010000. 1 1 2.24 1 0
## 49 1 000000000000000. 0 0 -0.69 0 0
## 50 1 000000000000000. 1 0 0.06 0 0
## 51 1 000000000000000. 1 0 0.38 0 0
## 52 1 000000000000000. 1 1 1.26 0 0
## 53 1 000000010000000. 1 0 0.38 0 0
## 54 1 000000000000000. 0 1 -0.85 0 0
## 55 1 000000000000000. 0 1 -0.57 0 0
## 56 1 000000000000000. 0 1 -0.98 0 0
## 57 1 000000000000000. 1 1 1.22 0 1
## 58 1 000000000000000. 1 1 1.62 0 1
## 59 1 000000000000000. 1 1 0.71 0 0
## 60 1 000000000000000. 0 0 -0.59 0 0
## 61 1 000010000000000. 0 1 -0.82 0 0
## 62 1 000000000000000. 0 1 0.22 0 0
## 63 1 000000001000000. 0 1 -0.75 1 0
## 64 1 000000000000000. 0 1 -0.47 1 0
## 65 1 000000000000000. 1 1 0.34 0 1
## 66 1 000000000000000. 1 1 -0.40 1 0
## 67 1 000000000000000. 0 1 -0.29 1 0
## 68 1 000000000001000. 1 1 0.18 1 0
## 69 1 000000000000000. 0 0 -0.70 0 0
## 70 1 000000000000000. 0 1 -0.59 0 0
## 71 1 000000000000000. 0 0 -0.87 0 0
## 72 1 000001000000000. 0 1 -0.59 1 0
## 73 1 000000000000000. 1 1 -0.65 0 0
## 74 1 000000000000000. 0 0 -0.81 0 0
## 75 1 000000000000010. 0 1 2.88 0 0
## 76 1 000000000000000. 1 1 2.65 0 0
## 77 1 000000000000000. 1 0 0.45 0 0
## 78 1 000000000000000. 0 1 -0.26 0 0
## 79 1 000010001000110. 0 1 -0.68 1 0
## 80 1 000000000000000. 1 1 -0.75 0 0
## 81 1 000000000000000. 0 1 -0.15 1 0
## 82 1 010101000000011. 1 1 0.26 0 1
## 83 1 000000000001100. 1 1 0.03 0 1
## 84 1 000000000000000. 0 0 -0.96 0 0
## hectares_z LiDar_z RabbitPresence_2013_14 icats1 dist2014z icat2014
## 1 -0.74 -0.35 0 0 0.51 0
## 2 -0.78 -0.88 1 2 -1.60 1
## 3 -0.78 -1.13 1 0 -0.88 1
## 4 -0.78 -0.84 1 0 -0.75 1
## 5 0.51 -0.32 1 0 0.34 0
## 6 1.84 -0.07 1 0 -0.27 0
## 7 0.51 -0.33 0 0 1.31 0
## 8 0.51 0.43 1 2 0.92 0
## 9 1.84 -0.35 1 0 0.41 0
## 10 1.84 0.09 0 0 1.20 0
## 11 -0.71 -0.21 0 2 0.67 0
## 12 0.38 -0.50 1 0 0.46 0
## 13 0.38 1.81 0 1 0.37 0
## 14 1.84 1.47 0 2 0.18 0
## 15 1.84 -0.86 1 0 0.55 0
## 16 0.38 1.62 0 2 0.91 0
## 17 1.84 0.11 1 1 0.11 0
## 18 0.51 1.28 0 0 0.44 0
## 19 0.51 1.00 0 1 1.42 0
## 20 -0.78 -0.80 1 0 -0.45 0
## 21 -1.38 -0.82 1 0 1.49 0
## 22 0.51 1.58 0 1 1.47 0
## 23 -0.78 -0.86 1 3 -0.90 1
## 24 -0.78 -1.17 1 2 -1.29 1
## 25 1.84 2.23 0 4 0.29 0
## 26 0.38 0.61 0 1 1.02 0
## 27 -1.12 -0.92 0 0 1.05 0
## 28 1.84 0.58 0 0 1.38 0
## 29 -0.71 0.54 0 2 1.17 0
## 30 -1.21 -0.89 0 0 0.11 0
## 31 -1.09 -1.05 0 0 -0.26 0
## 32 -1.01 0.88 0 1 -1.14 1
## 33 -1.31 -1.76 0 1 -0.66 0
## 34 0.18 -1.12 1 0 -0.19 0
## 35 -0.03 -1.24 1 1 -0.39 0
## 36 0.18 -0.77 1 4 -1.57 1
## 37 0.18 -0.44 0 2 -1.15 1
## 38 0.18 -0.66 1 0 -1.14 1
## 39 -0.03 -0.67 0 0 -0.62 0
## 40 0.18 -1.23 1 2 -1.52 1
## 41 0.18 -1.03 1 0 -1.06 1
## 42 1.84 2.33 0 1 0.07 0
## 43 -0.03 -1.17 0 1 -0.80 1
## 44 -1.49 -0.58 0 1 -1.14 1
## 45 -1.22 -0.07 0 2 -0.62 0
## 46 -1.28 -0.42 0 2 0.61 0
## 47 -0.95 -0.47 0 0 -0.23 0
## 48 -0.58 -0.20 1 2 -0.04 0
## 49 0.54 0.69 0 2 1.51 0
## 50 0.54 1.96 0 1 1.53 0
## 51 1.95 0.59 0 2 1.45 0
## 52 1.95 -0.09 1 2 1.10 0
## 53 0.54 2.65 0 1 0.69 0
## 54 -1.18 -0.69 0 4 2.20 0
## 55 -1.18 -0.02 0 2 1.83 0
## 56 -1.28 0.13 0 4 0.27 0
## 57 -1.08 -1.22 0 1 0.27 0
## 58 -1.08 -0.65 0 1 0.58 0
## 59 -1.13 -0.38 0 0 0.55 0
## 60 1.95 2.10 0 0 -0.32 0
## 61 1.95 -0.61 1 1 0.07 0
## 62 -0.80 1.24 0 3 -1.22 1
## 63 0.14 -1.26 1 0 -0.71 1
## 64 0.14 -1.11 1 2 -0.36 0
## 65 -1.15 -0.97 0 1 0.19 0
## 66 0.36 -0.51 1 1 -1.33 1
## 67 0.36 -0.95 1 1 -0.83 1
## 68 0.14 -0.92 1 1 -0.73 1
## 69 0.19 0.99 0 5 -0.49 0
## 70 0.19 0.16 0 3 -1.09 1
## 71 0.19 0.09 0 3 -1.40 1
## 72 0.14 -0.79 1 0 -0.88 1
## 73 0.19 0.24 1 3 -1.09 1
## 74 0.19 0.32 0 3 -0.99 1
## 75 0.67 -0.35 0 0 0.62 0
## 76 1.95 1.34 0 3 1.61 0
## 77 0.54 1.59 0 1 1.06 0
## 78 -1.01 -0.26 1 0 -0.09 0
## 79 0.14 -1.05 1 1 -0.60 1
## 80 -0.80 0.79 0 4 -1.12 1
## 81 0.36 -0.43 1 0 -1.20 1
## 82 -0.98 -1.24 0 0 -0.09 0
## 83 -1.08 -0.96 1 1 -0.16 0
## 84 0.19 0.70 0 3 -1.18 1
##
## $model
## [1] "Occupancy"
##
## $mixtures
## [1] 1
##
## $freq
## group1
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1
## 7 1
## 8 1
## 9 1
## 10 1
## 11 1
## 12 1
## 13 1
## 14 1
## 15 1
## 16 1
## 17 1
## 18 1
## 19 1
## 20 1
## 21 1
## 22 1
## 23 1
## 24 1
## 25 1
## 26 1
## 27 1
## 28 1
## 29 1
## 30 1
## 31 1
## 32 1
## 33 1
## 34 1
## 35 1
## 36 1
## 37 1
## 38 1
## 39 1
## 40 1
## 41 1
## 42 1
## 43 1
## 44 1
## 45 1
## 46 1
## 47 1
## 48 1
## 49 1
## 50 1
## 51 1
## 52 1
## 53 1
## 54 1
## 55 1
## 56 1
## 57 1
## 58 1
## 59 1
## 60 1
## 61 1
## 62 1
## 63 1
## 64 1
## 65 1
## 66 1
## 67 1
## 68 1
## 69 1
## 70 1
## 71 1
## 72 1
## 73 1
## 74 1
## 75 1
## 76 1
## 77 1
## 78 1
## 79 1
## 80 1
## 81 1
## 82 1
## 83 1
## 84 1
##
## $nocc
## [1] 16
##
## $nocc.secondary
## NULL
##
## $time.intervals
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $begin.time
## [1] 1
##
## $age.unit
## [1] 1
##
## $initial.ages
## [1] 0
##
## $group.covariates
## NULL
##
## $nstrata
## [1] 1
##
## $strata.labels
## NULL
##
## $counts
## NULL
##
## $reverse
## [1] FALSE
##
## $areas
## NULL
##
## $events
## NULL
# Get the parameter names
# What are the parameter names for Single Season Single Species models
setup.parameters("Occupancy", check=TRUE)
## [1] "p" "Psi"
# The top model for the dyanmic occuanpcy is psi(global), .... p(humtrail)
# If you look at the supplemental table for model 20 you see that the "global" terms are
# this is psi(humtrail+rabbitat+devel+fresh+salt_butt+patch+LiDar)
# Notice the commas between the column and the placement of the quotes
# Define the models.
# model.type can be
#The random occupancy model cannot be fit here. See the other code in this directory.
# The reserved keywords can be found by looking at the ddl structures earlier
global <- paste0("~",paste0(global.covar.vars, collapse="+"),collapse="")
model.list.csv <- textConnection("
p, Psi, model.type, comment
~HumanTrail, ~global+icats1, Occupancy, (20) in paper
~HumanTrail, ~1, Occupancy,
~1, ~global+icats1, Occupancy,
~1, ~dist_dev_z, Occupancy,
~1, ~Rabbitat, Occupancy,
~1, ~freshwater, Occupancy,
~1, ~icats1, Occupancy,
~1, ~HumanTrail, Occupancy,
~1, ~freshwater+icats1, Occupancy,
~1, ~freshwater+icats1+Rabbitat, Occupancy")
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
## p Psi model.type comment
## 1 ~HumanTrail ~global+icats1 Occupancy (20) in paper
## 2 ~HumanTrail ~1 Occupancy
## 3 ~1 ~global+icats1 Occupancy
## 4 ~1 ~dist_dev_z Occupancy
## 5 ~1 ~Rabbitat Occupancy
## 6 ~1 ~freshwater Occupancy
## 7 ~1 ~icats1 Occupancy
## 8 ~1 ~HumanTrail Occupancy
## 9 ~1 ~freshwater+icats1 Occupancy
## 10 ~1 ~freshwater+icats1+Rabbitat Occupancy
## model.number
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6
## 7 7
## 8 8
## 9 9
## 10 10
# substitute for some of the shortcuts
model.list$Psi <- gsub("~global",global, model.list$Psi, )
model.list
## p
## 1 ~HumanTrail
## 2 ~HumanTrail
## 3 ~1
## 4 ~1
## 5 ~1
## 6 ~1
## 7 ~1
## 8 ~1
## 9 ~1
## 10 ~1
## Psi
## 1 ~HumanTrail+Rabbitat+dist_dev_z+freshwater+salt_butt+hectares_z+LiDar_z+icats1
## 2 ~1
## 3 ~HumanTrail+Rabbitat+dist_dev_z+freshwater+salt_butt+hectares_z+LiDar_z+icats1
## 4 ~dist_dev_z
## 5 ~Rabbitat
## 6 ~freshwater
## 7 ~icats1
## 8 ~HumanTrail
## 9 ~freshwater+icats1
## 10 ~freshwater+icats1+Rabbitat
## model.type comment model.number
## 1 Occupancy (20) in paper 1
## 2 Occupancy 2
## 3 Occupancy 3
## 4 Occupancy 4
## 5 Occupancy 5
## 6 Occupancy 6
## 7 Occupancy 7
## 8 Occupancy 8
## 9 Occupancy 9
## 10 Occupancy 10
# fit the models
myobs <- ls()
myobs <- myobs[ grepl("m...",myobs,fixed=TRUE)]
cat("Removing ", myobs, "\n")
## Removing
rm(list=myobs)
cleanup(ask=FALSE)
model.fits <- plyr::dlply(model.list, "model.number", function(x,input.history){
cat("\n\n***** Starting ", unlist(x), "\n")
# we need to process the data in the loop to allow for different data types
# notice that max.visits.per.year and n.season are defined outside the function (bad programming practise)
input.data<- process.data(data=input.history, model=x$model.type)
# set up the ddls as needed for time-varying covariates.
# you need to do this in the loop because different paraemterizations have different ddl structures
input.ddl <- make.design.data(input.data)
model.parameters=list(
Psi =list(formula=as.formula(eval(x$Psi))),
p =list(formula=as.formula(eval(x$p ))) )
fit <- RMark::mark(input.data, ddl=input.ddl,
model=x$model.type,
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.history=input.history)
##
##
## ***** Starting ~HumanTrail ~HumanTrail+Rabbitat+dist_dev_z+freshwater+salt_butt+hectares_z+LiDar_z+icats1 Occupancy (20) in paper 1
##
## Output summary for Occupancy model
## Name : p(~HumanTrail)Psi(~HumanTrail + Rabbitat + dist_dev_z + freshwater + salt_butt + hectares_z + LiDar_z + icats1)
##
## Npar : 11
## -2lnL: 405.1902
## AICc : 430.8569
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -2.3304590 0.1627693 -2.6494868 -2.011431
## p:HumanTrail -0.0365147 0.3686350 -0.7590394 0.686010
## Psi:(Intercept) -1.9715646 3.9313280 -9.6769676 5.733839
## Psi:HumanTrail 10.6505360 11.0306750 -10.9695860 32.270659
## Psi:Rabbitat -0.6847231 2.9409807 -6.4490454 5.079599
## Psi:dist_dev_z -0.5417020 0.7882886 -2.0867478 1.003344
## Psi:freshwater 43.8216010 43.7563750 -41.9408970 129.584100
## Psi:salt_butt 9.5082147 12.4825200 -14.9575250 33.973954
## Psi:hectares_z 7.7346818 7.6075254 -7.1760682 22.645432
## Psi:LiDar_z -0.9842652 1.3207657 -3.5729661 1.604436
## Psi:icats1 -11.4199110 10.9643300 -32.9100000 10.070177
##
##
## Real Parameter p
## 1 2 3 4 5 6 7
## 0.0878621 0.0878621 0.0878621 0.0878621 0.0878621 0.0878621 0.0878621
## 8 9 10 11 12 13 14
## 0.0878621 0.0878621 0.0878621 0.0878621 0.0878621 0.0878621 0.0878621
## 15 16
## 0.0878621 0.0878621
##
##
## Real Parameter Psi
## 1
## 0.8258169
##
##
## ***** Starting ~HumanTrail ~1 Occupancy 2
##
## Output summary for Occupancy model
## Name : p(~HumanTrail)Psi(~1)
##
## Npar : 3
## -2lnL: 448.8881
## AICc : 455.1881
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -2.1652649 0.1984234 -2.5541749 -1.776355
## p:HumanTrail -0.1424886 0.4000019 -0.9264923 0.641515
## Psi:(Intercept) -0.2110573 0.2939012 -0.7871036 0.364989
##
##
## Real Parameter p
## 1 2 3 4 5 6 7
## 0.0995187 0.0995187 0.0995187 0.0995187 0.0995187 0.0995187 0.0995187
## 8 9 10 11 12 13 14
## 0.0995187 0.0995187 0.0995187 0.0995187 0.0995187 0.0995187 0.0995187
## 15 16
## 0.0995187 0.0995187
##
##
## Real Parameter Psi
## 1
## 0.4474307
##
##
## ***** Starting ~1 ~HumanTrail+Rabbitat+dist_dev_z+freshwater+salt_butt+hectares_z+LiDar_z+icats1 Occupancy 3
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~HumanTrail + Rabbitat + dist_dev_z + freshwater + salt_butt + hectares_z + LiDar_z + icats1)
##
## Npar : 10
## -2lnL: 405.2001
## AICc : 428.2138
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -2.3381125 0.1443965 -2.621130 -2.055095
## Psi:(Intercept) -1.9945896 4.0450748 -9.922936 5.933757
## Psi:HumanTrail 10.6649800 11.4168420 -11.712031 33.041992
## Psi:Rabbitat -0.6680538 2.9061555 -6.364119 5.028011
## Psi:dist_dev_z -0.5430811 0.7892060 -2.089925 1.003763
## Psi:freshwater 44.0956310 45.6993110 -45.475019 133.666280
## Psi:salt_butt 9.6261074 13.1506540 -16.149176 35.401391
## Psi:hectares_z 7.7768542 7.9413688 -7.788229 23.341937
## Psi:LiDar_z -0.9688752 1.3047814 -3.526247 1.588496
## Psi:icats1 -11.4891070 11.4535090 -33.937984 10.959771
##
##
## Real Parameter p
## 1 2 3 4 5 6 7
## 0.0880153 0.0880153 0.0880153 0.0880153 0.0880153 0.0880153 0.0880153
## 8 9 10 11 12 13 14
## 0.0880153 0.0880153 0.0880153 0.0880153 0.0880153 0.0880153 0.0880153
## 15 16
## 0.0880153 0.0880153
##
##
## Real Parameter Psi
## 1
## 0.8268588
##
##
## ***** Starting ~1 ~dist_dev_z Occupancy 4
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~dist_dev_z)
##
## Npar : 3
## -2lnL: 448.3872
## AICc : 454.6872
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -2.1940917 0.1813106 -2.5494605 -1.8387229
## Psi:(Intercept) -0.2399991 0.2939109 -0.8160644 0.3360662
## Psi:dist_dev_z 0.2139611 0.2743336 -0.3237328 0.7516549
##
##
## Real Parameter p
## 1 2 3 4 5 6 7
## 0.1002823 0.1002823 0.1002823 0.1002823 0.1002823 0.1002823 0.1002823
## 8 9 10 11 12 13 14
## 0.1002823 0.1002823 0.1002823 0.1002823 0.1002823 0.1002823 0.1002823
## 15 16
## 0.1002823 0.1002823
##
##
## Real Parameter Psi
## 1
## 0.4449115
##
##
## ***** Starting ~1 ~Rabbitat Occupancy 5
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~Rabbitat)
##
## Npar : 3
## -2lnL: 442.3904
## AICc : 448.6904
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -2.195926 0.1815495 -2.5517628 -1.8400888
## Psi:(Intercept) -1.391517 0.5785782 -2.5255306 -0.2575041
## Psi:Rabbitat 1.596797 0.6625250 0.2982483 2.8953463
##
##
## Real Parameter p
## 1 2 3 4 5 6 7 8
## 0.100117 0.100117 0.100117 0.100117 0.100117 0.100117 0.100117 0.100117
## 9 10 11 12 13 14 15 16
## 0.100117 0.100117 0.100117 0.100117 0.100117 0.100117 0.100117 0.100117
##
##
## Real Parameter Psi
## 1
## 0.4329176
##
##
## ***** Starting ~1 ~freshwater Occupancy 6
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~freshwater)
##
## Npar : 3
## -2lnL: 432.4914
## AICc : 438.7914
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -2.195872 0.1815383 -2.5516874 -1.840057
## Psi:(Intercept) -1.032017 0.3602240 -1.7380563 -0.325978
## Psi:freshwater 2.646822 0.9479552 0.7888297 4.504814
##
##
## Real Parameter p
## 1 2 3 4 5 6 7
## 0.1001218 0.1001218 0.1001218 0.1001218 0.1001218 0.1001218 0.1001218
## 8 9 10 11 12 13 14
## 0.1001218 0.1001218 0.1001218 0.1001218 0.1001218 0.1001218 0.1001218
## 15 16
## 0.1001218 0.1001218
##
##
## Real Parameter Psi
## 1
## 0.4548104
##
##
## ***** Starting ~1 ~icats1 Occupancy 7
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~icats1)
##
## Npar : 3
## -2lnL: 434.194
## AICc : 440.494
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -2.1854604 0.1788186 -2.5359449 -1.8349758
## Psi:(Intercept) 0.8633936 0.4890935 -0.0952298 1.8220169
## Psi:icats1 -0.9956022 0.3177901 -1.6184708 -0.3727336
##
##
## Real Parameter p
## 1 2 3 4 5 6 7
## 0.1010638 0.1010638 0.1010638 0.1010638 0.1010638 0.1010638 0.1010638
## 8 9 10 11 12 13 14
## 0.1010638 0.1010638 0.1010638 0.1010638 0.1010638 0.1010638 0.1010638
## 15 16
## 0.1010638 0.1010638
##
##
## Real Parameter Psi
## 1
## 0.4030045
##
##
## ***** Starting ~1 ~HumanTrail Occupancy 8
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~HumanTrail)
##
## Npar : 3
## -2lnL: 448.8339
## AICc : 455.1339
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -2.1944736 0.1814201 -2.550057 -1.8388902
## Psi:(Intercept) -0.1522397 0.3327460 -0.804422 0.4999425
## Psi:HumanTrail -0.2579762 0.6034353 -1.440709 0.9247569
##
##
## Real Parameter p
## 1 2 3 4 5 6 7
## 0.1002479 0.1002479 0.1002479 0.1002479 0.1002479 0.1002479 0.1002479
## 8 9 10 11 12 13 14
## 0.1002479 0.1002479 0.1002479 0.1002479 0.1002479 0.1002479 0.1002479
## 15 16
## 0.1002479 0.1002479
##
##
## Real Parameter Psi
## 1
## 0.4452689
##
##
## ***** Starting ~1 ~freshwater+icats1 Occupancy 9
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~freshwater + icats1)
##
## Npar : 4
## -2lnL: 416.9018
## AICc : 425.4081
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -2.2999524 0.1670603 -2.6273906 -1.9725142
## Psi:(Intercept) 0.6619266 0.7786362 -0.8642005 2.1880536
## Psi:freshwater 5.6681064 3.2987674 -0.7974778 12.1336910
## Psi:icats1 -1.8492627 0.8176438 -3.4518445 -0.2466809
##
##
## Real Parameter p
## 1 2 3 4 5 6 7
## 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269
## 8 9 10 11 12 13 14
## 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269
## 15 16
## 0.0911269 0.0911269
##
##
## Real Parameter Psi
## 1
## 0.5374856
##
##
## ***** Starting ~1 ~freshwater+icats1+Rabbitat Occupancy 10
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~freshwater + icats1 + Rabbitat)
##
## Npar : 5
## -2lnL: 416.8612
## AICc : 427.6304
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -2.3022866 0.1657357 -2.627128 -1.9774447
## Psi:(Intercept) 0.8284251 1.1585794 -1.442390 3.0992406
## Psi:freshwater 5.9401779 3.6085631 -1.132606 13.0129620
## Psi:icats1 -1.9084839 0.8795672 -3.632436 -0.1845322
## Psi:Rabbitat -0.2063470 1.0361540 -2.237209 1.8245149
##
##
## Real Parameter p
## 1 2 3 4 5 6 7
## 0.0909338 0.0909338 0.0909338 0.0909338 0.0909338 0.0909338 0.0909338
## 8 9 10 11 12 13 14
## 0.0909338 0.0909338 0.0909338 0.0909338 0.0909338 0.0909338 0.0909338
## 15 16
## 0.0909338 0.0909338
##
##
## Real Parameter Psi
## 1
## 0.5459968
# examine individual model results
model.number <-9
summary(model.fits[[model.number]])
## Output summary for Occupancy model
## Name : p(~1)Psi(~freshwater + icats1)
##
## Npar : 4
## -2lnL: 416.9018
## AICc : 425.4081
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -2.2999524 0.1670603 -2.6273906 -1.9725142
## Psi:(Intercept) 0.6619266 0.7786362 -0.8642005 2.1880536
## Psi:freshwater 5.6681064 3.2987674 -0.7974778 12.1336910
## Psi:icats1 -1.8492627 0.8176438 -3.4518445 -0.2466809
##
##
## Real Parameter p
## 1 2 3 4 5 6 7
## 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269
## 8 9 10 11 12 13 14
## 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269 0.0911269
## 15 16
## 0.0911269 0.0911269
##
##
## Real Parameter Psi
## 1
## 0.5374856
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## p g1 a0 t1 0.0911269 0.0138364 0.0673963 0.1221191
## Psi g1 a0 t1 0.5374856 0.2080525 0.1839027 0.8569979
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## p:(Intercept) -2.2999524 0.1670603 -2.6273906 -1.9725142
## Psi:(Intercept) 0.6619266 0.7786362 -0.8642005 2.1880536
## Psi:freshwater 5.6681064 3.2987674 -0.7974778 12.1336910
## Psi:icats1 -1.8492627 0.8176438 -3.4518445 -0.2466809
model.fits[[model.number]]$results$derived
## $Occupancy
## estimate lcl ucl
## 1 0.5374856 0.1839027 0.8569979
model.fits[[model.number]]$results$derived$"psi Probability Occupied"
## NULL
model.fits[[model.number]]$results$derived$"lambda Rate of Change"
## NULL
model.fits[[model.number]]$results$derived$"log odds lambda"
## NULL
get.real(model.fits[[model.number]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi g1 a0 t1 17 2 0.5374856 0.2080525 0.1839027 0.8569979
## fixed note group age time Age Time
## Psi g1 a0 t1 1 0 1 0 0
get.real(model.fits[[model.number]], "p", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## p g1 a0 t1 1 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a1 t2 2 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a2 t3 3 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a3 t4 4 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a4 t5 5 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a5 t6 6 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a6 t7 7 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a7 t8 8 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a8 t9 9 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a9 t10 10 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a10 t11 11 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a11 t12 12 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a12 t13 13 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a13 t14 14 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a14 t15 15 1 0.0911269 0.0138364 0.0673963 0.1221191
## p g1 a15 t16 16 1 0.0911269 0.0138364 0.0673963 0.1221191
## fixed note group age time Age Time
## p g1 a0 t1 1 0 1 0 0
## p g1 a1 t2 1 1 2 1 1
## p g1 a2 t3 1 2 3 2 2
## p g1 a3 t4 1 3 4 3 3
## p g1 a4 t5 1 4 5 4 4
## p g1 a5 t6 1 5 6 5 5
## p g1 a6 t7 1 6 7 6 6
## p g1 a7 t8 1 7 8 7 7
## p g1 a8 t9 1 8 9 8 8
## p g1 a9 t10 1 9 10 9 9
## p g1 a10 t11 1 10 11 10 10
## p g1 a11 t12 1 11 12 11 11
## p g1 a12 t13 1 12 13 12 12
## p g1 a13 t14 1 13 14 13 13
## p g1 a14 t15 1 14 15 14 14
## p g1 a15 t16 1 15 16 15 15
# collect models and make AICc table
model.set <- RMark::collect.models( type=NULL)
model.set
## model
## 9 p(~1)Psi(~freshwater + icats1)
## 10 p(~1)Psi(~freshwater + icats1 + Rabbitat)
## 3 p(~1)Psi(~HumanTrail + Rabbitat + dist_dev_z + freshwater + salt_butt + hectares_z + LiDar_z + icats1)
## 1 p(~HumanTrail)Psi(~HumanTrail + Rabbitat + dist_dev_z + freshwater + salt_butt + hectares_z + LiDar_z + icats1)
## 6 p(~1)Psi(~freshwater)
## 7 p(~1)Psi(~icats1)
## 5 p(~1)Psi(~Rabbitat)
## 4 p(~1)Psi(~dist_dev_z)
## 8 p(~1)Psi(~HumanTrail)
## 2 p(~HumanTrail)Psi(~1)
## npar AICc DeltaAICc weight Deviance
## 9 4 425.4081 0.000000 6.088512e-01 416.9018
## 10 5 427.6304 2.222322 2.004196e-01 416.8612
## 3 10 428.2138 2.805710 1.497128e-01 405.2001
## 1 11 430.8569 5.448788 3.993208e-02 405.1902
## 6 3 438.7914 13.383351 7.557070e-04 432.4914
## 7 3 440.4940 15.085891 3.225905e-04 434.1940
## 5 3 448.6904 23.282351 5.355659e-06 442.3904
## 4 3 454.6872 29.279141 2.670708e-07 448.3872
## 8 3 455.1339 29.725791 2.136179e-07 448.8339
## 2 3 455.1881 29.779971 2.079087e-07 448.8881
# model averaging in the usual way
RMark::model.average(model.set, "Psi")
## par.index estimate se fixed note group age time Age
## Psi g1 a0 t1 17 0.5939216 0.2463569 1 0 1 0
## Time
## Psi g1 a0 t1 0
RMark::model.average(model.set, "p")
## par.index estimate se fixed note group age time Age
## p g1 a0 t1 1 0.09050204 0.01347597 1 0 1 0
## p g1 a1 t2 2 0.09050204 0.01347597 1 1 2 1
## p g1 a2 t3 3 0.09050204 0.01347597 1 2 3 2
## p g1 a3 t4 4 0.09050204 0.01347597 1 3 4 3
## p g1 a4 t5 5 0.09050204 0.01347597 1 4 5 4
## p g1 a5 t6 6 0.09050204 0.01347597 1 5 6 5
## p g1 a6 t7 7 0.09050204 0.01347597 1 6 7 6
## p g1 a7 t8 8 0.09050204 0.01347597 1 7 8 7
## p g1 a8 t9 9 0.09050204 0.01347597 1 8 9 8
## p g1 a9 t10 10 0.09050204 0.01347597 1 9 10 9
## p g1 a10 t11 11 0.09050204 0.01347597 1 10 11 10
## p g1 a11 t12 12 0.09050204 0.01347597 1 11 12 11
## p g1 a12 t13 13 0.09050204 0.01347597 1 12 13 12
## p g1 a13 t14 14 0.09050204 0.01347597 1 13 14 13
## p g1 a14 t15 15 0.09050204 0.01347597 1 14 15 14
## p g1 a15 t16 16 0.09050204 0.01347597 1 15 16 15
## Time
## p g1 a0 t1 0
## p g1 a1 t2 1
## p g1 a2 t3 2
## p g1 a3 t4 3
## p g1 a4 t5 4
## p g1 a5 t6 5
## p g1 a6 t7 6
## p g1 a7 t8 7
## p g1 a8 t9 8
## p g1 a9 t10 9
## p g1 a10 t11 10
## p g1 a11 t12 11
## p g1 a12 t13 12
## p g1 a13 t14 13
## p g1 a14 t15 14
## p g1 a15 t16 15
# find occupancy for combination of freshwater and icats1
all.ddl <- make.design.data(rabbit.data)
all.ddl$Psi
## par.index model.index group age time Age Time
## 1 1 17 1 0 1 0 0
Temp.df <- expand.grid(freshwater=c(0,1), icats1=0:4)
model.fits[[9]]$results$beta
## estimate se lcl ucl
## p:(Intercept) -2.2999524 0.1670603 -2.6273906 -1.9725142
## Psi:(Intercept) 0.6619266 0.7786362 -0.8642005 2.1880536
## Psi:freshwater 5.6681064 3.2987674 -0.7974778 12.1336910
## Psi:icats1 -1.8492627 0.8176438 -3.4518445 -0.2466809
pred.psi <- covariate.predictions(model.fits[[9]], indices=17, data=Temp.df)
head(pred.psi$estimates)
## vcv.index model.index par.index freshwater icats1 estimate se
## 1 1 2 17 0 0 0.65969304 0.17480238
## 2 2 2 17 1 0 0.99822119 0.00672987
## 3 3 2 17 0 1 0.23373571 0.09529901
## 4 4 2 17 1 1 0.98880213 0.03445483
## 5 5 2 17 0 2 0.04579958 0.04977129
## 6 6 2 17 1 2 0.93286203 0.15767693
## lcl ucl fixed
## 1 0.296468349 0.8991690
## 2 0.250032137 0.9999989
## 3 0.097070014 0.4639477
## 4 0.165459572 0.9999746
## 5 0.005123596 0.3090775
## 6 0.090887192 0.9994824
ggplot(data=pred.psi$estimates, aes(x=icats1, y=estimate, color=as.factor(freshwater)))+
ggtitle("Estimated rabbit occupancy as a function of cats seen at a site")+
geom_point()+
geom_line()+
geom_ribbon(aes(ymin=lcl, ymax=ucl), alpha=0.1)
#####################################################################
#####################################################################
#####################################################################
# Next model static occupancy for the cats
# See Tables S2 and S3 of the paper for dynamic occupancy model
global.covar.vars <- c("HumanTrail","Rabbitat","dist_dev_z","freshwater","salt_butt","hectares_z","LiDar_z")
setdiff(global.covar.vars, names(input.data))
## character(0)
cat.detect.vars <- c("Cats.2013.2014", paste0("X.",31:45))
cat.input.history <- input.data[, cat.detect.vars]
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(cat.input.history)
## [1] 84
ncol(cat.input.history)
## [1] 16
range(cat.input.history, na.rm=TRUE)
## [1] 0 1
sum(is.na(cat.input.history))
## [1] 52
sum(cat.input.history, na.rm=TRUE)
## [1] 163
other.covar.vars <- c("Total.CatCaps_2013_14") # number of cats detected
setdiff(other.covar.vars, names(input.data))
## character(0)
site.covar <- input.data[, c(global.covar.vars,other.covar.vars)]
row.names(site.covar) <- NULL
head(site.covar)
## HumanTrail Rabbitat dist_dev_z freshwater salt_butt hectares_z LiDar_z
## 1 0 1 1.05 0 1 -0.74 -0.35
## 2 0 1 -0.83 1 0 -0.78 -0.88
## 3 0 1 0.12 1 0 -0.78 -1.13
## 4 0 1 0.00 1 0 -0.78 -0.84
## 5 1 1 2.20 0 0 0.51 -0.32
## 6 1 0 1.02 0 0 1.84 -0.07
## Total.CatCaps_2013_14
## 1 0
## 2 1
## 3 0
## 4 1
## 5 0
## 6 1
#Format the capture history to be used by RMark and deal with missing values
input.history <- data.frame(freq=1,
ch=apply(cat.input.history,1,paste, collapse=""), stringsAsFactors=FALSE)
head(input.history,n=20)
## freq ch
## 1 1 0000000000000000
## 2 1 0000000000101000
## 3 1 0000000000000000
## 4 1 0000000000000000
## 5 1 0000000000000000
## 6 1 0000000000000000
## 7 1 0000000000000000
## 8 1 0001010001000000
## 9 1 0000000000000000
## 10 1 0000000000000000
## 11 1 0000010000000100
## 12 1 0000000000000000
## 13 1 0100001000000000
## 14 1 1010000001001010
## 15 1 0000000000000000
## 16 1 0100000000000000
## 17 1 0000000001000000
## 18 1 0000000000000000
## 19 1 0001000000000000
## 20 1 000000000000000NA
input.history$ch <- gsub("NA",".",input.history$ch, fixed=TRUE)
head(input.history,n=20)
## freq ch
## 1 1 0000000000000000
## 2 1 0000000000101000
## 3 1 0000000000000000
## 4 1 0000000000000000
## 5 1 0000000000000000
## 6 1 0000000000000000
## 7 1 0000000000000000
## 8 1 0001010001000000
## 9 1 0000000000000000
## 10 1 0000000000000000
## 11 1 0000010000000100
## 12 1 0000000000000000
## 13 1 0100001000000000
## 14 1 1010000001001010
## 15 1 0000000000000000
## 16 1 0100000000000000
## 17 1 0000000001000000
## 18 1 0000000000000000
## 19 1 0001000000000000
## 20 1 000000000000000.
#Add site covariates to input history
input.history = cbind(input.history,site.covar)
head(input.history)
## freq ch HumanTrail Rabbitat dist_dev_z freshwater salt_butt
## 1 1 0000000000000000 0 1 1.05 0 1
## 2 1 0000000000101000 0 1 -0.83 1 0
## 3 1 0000000000000000 0 1 0.12 1 0
## 4 1 0000000000000000 0 1 0.00 1 0
## 5 1 0000000000000000 1 1 2.20 0 0
## 6 1 0000000000000000 1 0 1.02 0 0
## hectares_z LiDar_z Total.CatCaps_2013_14
## 1 -0.74 -0.35 0
## 2 -0.78 -0.88 1
## 3 -0.78 -1.13 0
## 4 -0.78 -0.84 1
## 5 0.51 -0.32 0
## 6 1.84 -0.07 1
# With 16 detections visits, look at the naive occupancy vs. the different covariates
input.history$CatCaps <- as.numeric(input.history$Total.CatCaps_2013_14>0)
prop.table(xtabs(~CatCaps+HumanTrail, data=input.history, exclude=NULL, na.action=na.pass), margin=2)
## HumanTrail
## CatCaps 0 1
## 0 0.6290323 0.3636364
## 1 0.3709677 0.6363636
prop.table(xtabs(~CatCaps+Rabbitat, data=input.history, exclude=NULL, na.action=na.pass), margin=2)
## Rabbitat
## CatCaps 0 1
## 0 0.4400000 0.6101695
## 1 0.5600000 0.3898305
prop.table(xtabs(~CatCaps+freshwater, data=input.history, exclude=NULL, na.action=na.pass), margin=2)
## freshwater
## CatCaps 0 1
## 0 0.5614035 0.5555556
## 1 0.4385965 0.4444444
prop.table(xtabs(~CatCaps+salt_butt, data=input.history, exclude=NULL, na.action=na.pass), margin=2)
## salt_butt
## CatCaps 0 1
## 0 0.5333333 0.7777778
## 1 0.4666667 0.2222222
xtabs(~freshwater+salt_butt, data=input.history, exclude=NULL, na.action=na.pass)
## salt_butt
## freshwater 0 1
## 0 48 9
## 1 27 0
#Create the RMark data structure
cat.data <- process.data(data=input.history,
model="Occupancy")
cat.data
## $data
## freq ch HumanTrail Rabbitat dist_dev_z freshwater salt_butt
## 1 1 0000000000000000 0 1 1.05 0 1
## 2 1 0000000000101000 0 1 -0.83 1 0
## 3 1 0000000000000000 0 1 0.12 1 0
## 4 1 0000000000000000 0 1 0.00 1 0
## 5 1 0000000000000000 1 1 2.20 0 0
## 6 1 0000000000000000 1 0 1.02 0 0
## 7 1 0000000000000000 0 1 -0.35 0 0
## 8 1 0001010001000000 0 0 1.00 0 0
## 9 1 0000000000000000 0 1 0.98 0 0
## 10 1 0000000000000000 0 1 -0.24 0 0
## 11 1 0000010000000100 0 0 2.18 0 0
## 12 1 0000000000000000 1 1 0.50 0 0
## 13 1 0100001000000000 0 1 0.40 0 0
## 14 1 1010000001001010 0 0 -0.34 0 0
## 15 1 0000000000000000 0 1 1.13 0 0
## 16 1 0100000000000000 0 0 0.29 0 0
## 17 1 0000000001000000 0 0 1.61 0 0
## 18 1 0000000000000000 0 0 0.97 0 0
## 19 1 0001000000000000 0 0 -0.38 0 0
## 20 1 000000000000000. 0 1 1.38 1 0
## 21 1 000000000000000. 0 1 0.69 1 0
## 22 1 000100000000000. 0 0 0.21 0 0
## 23 1 0000000010100011 0 1 -1.01 1 0
## 24 1 0001000000001101 0 1 -0.48 1 0
## 25 1 000000111000010. 0 0 -0.08 0 0
## 26 1 0000001000000001 0 1 -0.14 0 0
## 27 1 0000000000000000 0 1 0.06 1 0
## 28 1 0000000000000000 0 0 -0.79 0 0
## 29 1 1000010000100001 1 0 1.59 0 0
## 30 1 000000000000000. 0 1 1.68 0 0
## 31 1 000000000000000. 0 1 2.08 0 1
## 32 1 001000000000000. 0 0 -0.72 0 0
## 33 1 100000000000000. 1 1 -0.54 0 1
## 34 1 000000000000000. 0 1 0.41 1 0
## 35 1 0000000100000000 0 1 -1.13 1 0
## 36 1 001010101100000. 0 1 -0.85 1 0
## 37 1 100000101001000. 0 1 -0.93 1 0
## 38 1 0000000000000000 0 1 -0.95 1 0
## 39 1 000000000000000. 0 1 -0.76 1 0
## 40 1 001100001001000. 0 1 -0.31 1 0
## 41 1 1000000000000000 0 1 -0.01 1 0
## 42 1 0000001000000000 0 0 -1.01 0 0
## 43 1 0010000000000000 0 1 -0.54 1 0
## 44 1 0000000000100000 0 1 -0.96 0 0
## 45 1 0000000010000001 0 0 -0.85 0 0
## 46 1 100010000000000. 0 1 -0.96 0 1
## 47 1 000000000000000. 0 1 -0.49 1 0
## 48 1 010001000010000. 1 1 2.24 1 0
## 49 1 000101010010000. 0 0 -0.69 0 0
## 50 1 000000000000001. 1 0 0.06 0 0
## 51 1 000000010000100. 1 0 0.38 0 0
## 52 1 010000001000000. 1 1 1.26 0 0
## 53 1 000000010000000. 1 0 0.38 0 0
## 54 1 011011110010111. 0 1 -0.85 0 0
## 55 1 010001100001001. 0 1 -0.57 0 0
## 56 1 010110000101010. 0 1 -0.98 0 0
## 57 1 000000000100000. 1 1 1.22 0 1
## 58 1 000000000001000. 1 1 1.62 0 1
## 59 1 000000000000000. 1 1 0.71 0 0
## 60 1 000000000000000. 0 0 -0.59 0 0
## 61 1 000000010000000. 0 1 -0.82 0 0
## 62 1 000010100100000. 0 1 0.22 0 0
## 63 1 000000000000000. 0 1 -0.75 1 0
## 64 1 000100000010100. 0 1 -0.47 1 0
## 65 1 000000000000100. 1 1 0.34 0 1
## 66 1 000001000000000. 1 1 -0.40 1 0
## 67 1 000010000000000. 0 1 -0.29 1 0
## 68 1 000010000000000. 1 1 0.18 1 0
## 69 1 010111110111011. 0 0 -0.70 0 0
## 70 1 001000010000001. 0 1 -0.59 0 0
## 71 1 001101000000100. 0 0 -0.87 0 0
## 72 1 000000000000000. 0 1 -0.59 1 0
## 73 1 100010001100001. 1 1 -0.65 0 0
## 74 1 001000101110111. 0 0 -0.81 0 0
## 75 1 000000000000000. 0 1 2.88 0 0
## 76 1 010001000101010. 1 1 2.65 0 0
## 77 1 000000001001000. 1 0 0.45 0 0
## 78 1 000000000000000. 0 1 -0.26 0 0
## 79 1 000101001000000. 0 1 -0.68 1 0
## 80 1 000011100011011. 1 1 -0.75 0 0
## 81 1 000000000000000. 0 1 -0.15 1 0
## 82 1 000000000000000. 1 1 0.26 0 1
## 83 1 000000000010000. 1 1 0.03 0 1
## 84 1 010011001100110. 0 0 -0.96 0 0
## hectares_z LiDar_z Total.CatCaps_2013_14 CatCaps
## 1 -0.74 -0.35 0 0
## 2 -0.78 -0.88 1 1
## 3 -0.78 -1.13 0 0
## 4 -0.78 -0.84 1 1
## 5 0.51 -0.32 0 0
## 6 1.84 -0.07 1 1
## 7 0.51 -0.33 1 1
## 8 0.51 0.43 0 0
## 9 1.84 -0.35 0 0
## 10 1.84 0.09 0 0
## 11 -0.71 -0.21 2 1
## 12 0.38 -0.50 0 0
## 13 0.38 1.81 0 0
## 14 1.84 1.47 0 0
## 15 1.84 -0.86 0 0
## 16 0.38 1.62 0 0
## 17 1.84 0.11 0 0
## 18 0.51 1.28 0 0
## 19 0.51 1.00 0 0
## 20 -0.78 -0.80 0 0
## 21 -1.38 -0.82 0 0
## 22 0.51 1.58 0 0
## 23 -0.78 -0.86 0 0
## 24 -0.78 -1.17 0 0
## 25 1.84 2.23 0 0
## 26 0.38 0.61 0 0
## 27 -1.12 -0.92 1 1
## 28 1.84 0.58 0 0
## 29 -0.71 0.54 7 1
## 30 -1.21 -0.89 0 0
## 31 -1.09 -1.05 0 0
## 32 -1.01 0.88 0 0
## 33 -1.31 -1.76 0 0
## 34 0.18 -1.12 1 1
## 35 -0.03 -1.24 2 1
## 36 0.18 -0.77 0 0
## 37 0.18 -0.44 0 0
## 38 0.18 -0.66 0 0
## 39 -0.03 -0.67 0 0
## 40 0.18 -1.23 0 0
## 41 0.18 -1.03 0 0
## 42 1.84 2.33 0 0
## 43 -0.03 -1.17 1 1
## 44 -1.49 -0.58 0 0
## 45 -1.22 -0.07 2 1
## 46 -1.28 -0.42 1 1
## 47 -0.95 -0.47 3 1
## 48 -0.58 -0.20 3 1
## 49 0.54 0.69 1 1
## 50 0.54 1.96 3 1
## 51 1.95 0.59 4 1
## 52 1.95 -0.09 3 1
## 53 0.54 2.65 2 1
## 54 -1.18 -0.69 5 1
## 55 -1.18 -0.02 3 1
## 56 -1.28 0.13 1 1
## 57 -1.08 -1.22 0 0
## 58 -1.08 -0.65 3 1
## 59 -1.13 -0.38 1 1
## 60 1.95 2.10 2 1
## 61 1.95 -0.61 0 0
## 62 -0.80 1.24 0 0
## 63 0.14 -1.26 0 0
## 64 0.14 -1.11 6 1
## 65 -1.15 -0.97 0 0
## 66 0.36 -0.51 7 1
## 67 0.36 -0.95 2 1
## 68 0.14 -0.92 0 0
## 69 0.19 0.99 16 1
## 70 0.19 0.16 0 0
## 71 0.19 0.09 10 1
## 72 0.14 -0.79 1 1
## 73 0.19 0.24 5 1
## 74 0.19 0.32 15 1
## 75 0.67 -0.35 0 0
## 76 1.95 1.34 8 1
## 77 0.54 1.59 1 1
## 78 -1.01 -0.26 0 0
## 79 0.14 -1.05 0 0
## 80 -0.80 0.79 1 1
## 81 0.36 -0.43 0 0
## 82 -0.98 -1.24 0 0
## 83 -1.08 -0.96 0 0
## 84 0.19 0.70 5 1
##
## $model
## [1] "Occupancy"
##
## $mixtures
## [1] 1
##
## $freq
## group1
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1
## 7 1
## 8 1
## 9 1
## 10 1
## 11 1
## 12 1
## 13 1
## 14 1
## 15 1
## 16 1
## 17 1
## 18 1
## 19 1
## 20 1
## 21 1
## 22 1
## 23 1
## 24 1
## 25 1
## 26 1
## 27 1
## 28 1
## 29 1
## 30 1
## 31 1
## 32 1
## 33 1
## 34 1
## 35 1
## 36 1
## 37 1
## 38 1
## 39 1
## 40 1
## 41 1
## 42 1
## 43 1
## 44 1
## 45 1
## 46 1
## 47 1
## 48 1
## 49 1
## 50 1
## 51 1
## 52 1
## 53 1
## 54 1
## 55 1
## 56 1
## 57 1
## 58 1
## 59 1
## 60 1
## 61 1
## 62 1
## 63 1
## 64 1
## 65 1
## 66 1
## 67 1
## 68 1
## 69 1
## 70 1
## 71 1
## 72 1
## 73 1
## 74 1
## 75 1
## 76 1
## 77 1
## 78 1
## 79 1
## 80 1
## 81 1
## 82 1
## 83 1
## 84 1
##
## $nocc
## [1] 16
##
## $nocc.secondary
## NULL
##
## $time.intervals
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $begin.time
## [1] 1
##
## $age.unit
## [1] 1
##
## $initial.ages
## [1] 0
##
## $group.covariates
## NULL
##
## $nstrata
## [1] 1
##
## $strata.labels
## NULL
##
## $counts
## NULL
##
## $reverse
## [1] FALSE
##
## $areas
## NULL
##
## $events
## NULL
# all time-specific covariates must be added to the ddl in the loop below
# Get the parameter names
# What are the parameter names for Single Season Single Species models
setup.parameters("Occupancy", check=TRUE)
## [1] "p" "Psi"
# The top model is psi(global), gamma(ind change in cats) epsilon (dot) p(humtrail)
# If you look at the supplemental table for model 20 you see that
# this is psi(humtrail+catat+devel+fresh+salt_butt+patch+LiDar)
# Notice the commas between the column and the placement of the quotes
global <- paste0("~",paste0(global.covar.vars, collapse="+"),collapse="")
model.list.csv <- textConnection("
p, Psi, model.type, comment
~1, ~global, Occupancy, (15) Table S2
~1, ~dist_dev_z, Occupancy,
~1, ~Rabbitat, Occupancy,
~1, ~freshwater, Occupancy,
~1, ~HumanTrail, Occupancy,
~1, ~freshwater+Rabbitat, Occupancy,
~1, ~1, Occupancy,
~1, ~hectares_z+dist_dev_z, Occupancy")
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
## p Psi model.type comment model.number
## 1 ~1 ~global Occupancy (15) Table S2 1
## 2 ~1 ~dist_dev_z Occupancy 2
## 3 ~1 ~Rabbitat Occupancy 3
## 4 ~1 ~freshwater Occupancy 4
## 5 ~1 ~HumanTrail Occupancy 5
## 6 ~1 ~freshwater+Rabbitat Occupancy 6
## 7 ~1 ~1 Occupancy 7
## 8 ~1 ~hectares_z+dist_dev_z Occupancy 8
# substitute for some of the shortcuts
model.list$Psi <- gsub("~global",global, model.list$Psi, )
model.list
## p Psi
## 1 ~1 ~HumanTrail+Rabbitat+dist_dev_z+freshwater+salt_butt+hectares_z+LiDar_z
## 2 ~1 ~dist_dev_z
## 3 ~1 ~Rabbitat
## 4 ~1 ~freshwater
## 5 ~1 ~HumanTrail
## 6 ~1 ~freshwater+Rabbitat
## 7 ~1 ~1
## 8 ~1 ~hectares_z+dist_dev_z
## model.type comment model.number
## 1 Occupancy (15) Table S2 1
## 2 Occupancy 2
## 3 Occupancy 3
## 4 Occupancy 4
## 5 Occupancy 5
## 6 Occupancy 6
## 7 Occupancy 7
## 8 Occupancy 8
# fit the models
myobs <- ls()
myobs <- myobs[ grepl("m...",myobs,fixed=TRUE)]
cat("Removing ", myobs, "\n")
## Removing m...001 m...002 m...003 m...004 m...005 m...006 m...007 m...008 m...009 m...010
rm(list=myobs)
cleanup(ask=FALSE)
model.fits <- plyr::dlply(model.list, "model.number", function(x,input.history){
cat("\n\n***** Starting ", unlist(x), "\n")
# we need to process the data in the loop to allow for different data types
# notice that max.visits.per.year and n.season are defined outside the function (bad programming practise)
input.data<- process.data(data=input.history, model=x$model.type)
# set up the ddls as needed for time-varying covariates.
# you need to do this in the loop because different paraemterizations have different ddl structures
input.ddl <- make.design.data(input.data)
model.parameters=list(
Psi =list(formula=as.formula(eval(x$Psi))),
p =list(formula=as.formula(eval(x$p)))
)
fit <- RMark::mark(input.data, ddl=input.ddl,
model=x$model.type,
model.parameters=model.parameters,
se=TRUE
#,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.history=input.history)
##
##
## ***** Starting ~1 ~HumanTrail+Rabbitat+dist_dev_z+freshwater+salt_butt+hectares_z+LiDar_z Occupancy (15) Table S2 1
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~HumanTrail + Rabbitat + dist_dev_z + freshwater + salt_butt + hectares_z + LiDar_z)
##
## Npar : 9
## -2lnL: 910.8596
## AICc : 931.2921
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.4852457 0.0966141 -1.6746093 -1.2958821
## Psi:(Intercept) 1.6277101 1.0614847 -0.4528000 3.7082202
## Psi:HumanTrail 1.7847479 0.9309654 -0.0399443 3.6094402
## Psi:Rabbitat -1.4161672 1.2061910 -3.7803017 0.9479673
## Psi:dist_dev_z -0.9676143 0.3831503 -1.7185889 -0.2166398
## Psi:freshwater 0.4237291 0.9246481 -1.3885812 2.2360394
## Psi:salt_butt 0.6479538 1.4202754 -2.1357861 3.4316937
## Psi:hectares_z -0.2622384 0.3811695 -1.0093305 0.4848538
## Psi:LiDar_z 0.8824997 0.6559389 -0.4031406 2.1681400
##
##
## Real Parameter p
## all.diff.index par.index estimate se lcl ucl
## p g1 a0 t1 1 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a1 t2 2 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a2 t3 3 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a3 t4 4 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a4 t5 5 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a5 t6 6 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a6 t7 7 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a7 t8 8 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a8 t9 9 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a9 t10 10 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a10 t11 11 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a11 t12 12 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a12 t13 13 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a13 t14 14 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a14 t15 15 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a15 t16 16 1 0.1846364 0.0145448 0.1578106 0.2148589
## fixed
## p g1 a0 t1
## p g1 a1 t2
## p g1 a2 t3
## p g1 a3 t4
## p g1 a4 t5
## p g1 a5 t6
## p g1 a6 t7
## p g1 a7 t8
## p g1 a8 t9
## p g1 a9 t10
## p g1 a10 t11
## p g1 a11 t12
## p g1 a12 t13
## p g1 a13 t14
## p g1 a14 t15
## p g1 a15 t16
##
##
## Real Parameter Psi
## all.diff.index par.index estimate se lcl ucl
## Psi g1 a0 t1 17 2 0.7569719 0.0837883 0.5605646 0.8837925
## fixed
## Psi g1 a0 t1
##
##
## ***** Starting ~1 ~dist_dev_z Occupancy 2
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~dist_dev_z)
##
## Npar : 3
## -2lnL: 926.8677
## AICc : 933.1677
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.4800707 0.0941193 -1.664544 -1.2955969
## Psi:(Intercept) 0.8579453 0.2700573 0.328633 1.3872576
## Psi:dist_dev_z -0.5564771 0.2555065 -1.057270 -0.0556844
##
##
## Real Parameter p
## all.diff.index par.index estimate se lcl ucl
## p g1 a0 t1 1 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a1 t2 2 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a2 t3 3 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a3 t4 4 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a4 t5 5 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a5 t6 6 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a6 t7 7 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a7 t8 8 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a8 t9 9 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a9 t10 10 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a10 t11 11 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a11 t12 12 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a12 t13 13 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a13 t14 14 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a14 t15 15 1 0.1854167 0.0142155 0.1591529 0.214907
## p g1 a15 t16 16 1 0.1854167 0.0142155 0.1591529 0.214907
## fixed
## p g1 a0 t1
## p g1 a1 t2
## p g1 a2 t3
## p g1 a3 t4
## p g1 a4 t5
## p g1 a5 t6
## p g1 a6 t7
## p g1 a7 t8
## p g1 a8 t9
## p g1 a9 t10
## p g1 a10 t11
## p g1 a11 t12
## p g1 a12 t13
## p g1 a13 t14
## p g1 a14 t15
## p g1 a15 t16
##
##
## Real Parameter Psi
## all.diff.index par.index estimate se lcl ucl
## Psi g1 a0 t1 17 2 0.6919363 0.0566101 0.571668 0.7907934
## fixed
## Psi g1 a0 t1
##
##
## ***** Starting ~1 ~Rabbitat Occupancy 3
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~Rabbitat)
##
## Npar : 3
## -2lnL: 926.0183
## AICc : 932.3183
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.481175 0.0943624 -1.6661251 -1.296224
## Psi:(Intercept) 1.944027 0.7028725 0.5663969 3.321657
## Psi:Rabbitat -1.530903 0.7537195 -3.0081933 -0.053613
##
##
## Real Parameter p
## all.diff.index par.index estimate se lcl ucl
## p g1 a0 t1 1 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a1 t2 2 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a2 t3 3 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a3 t4 4 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a4 t5 5 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a5 t6 6 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a6 t7 7 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a7 t8 8 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a8 t9 9 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a9 t10 10 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a10 t11 11 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a11 t12 12 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a12 t13 13 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a13 t14 14 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a14 t15 15 1 0.18525 0.0142424 0.1589415 0.2148011
## p g1 a15 t16 16 1 0.18525 0.0142424 0.1589415 0.2148011
## fixed
## p g1 a0 t1
## p g1 a1 t2
## p g1 a2 t3
## p g1 a3 t4
## p g1 a4 t5
## p g1 a5 t6
## p g1 a6 t7
## p g1 a7 t8
## p g1 a8 t9
## p g1 a9 t10
## p g1 a10 t11
## p g1 a11 t12
## p g1 a12 t13
## p g1 a13 t14
## p g1 a14 t15
## p g1 a15 t16
##
##
## Real Parameter Psi
## all.diff.index par.index estimate se lcl ucl
## Psi g1 a0 t1 17 2 0.7044855 0.060381 0.5745115 0.8080234
## fixed
## Psi g1 a0 t1
##
##
## ***** Starting ~1 ~freshwater Occupancy 4
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~freshwater)
##
## Npar : 3
## -2lnL: 930.1465
## AICc : 936.4465
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.4815439 0.0943966 -1.6665612 -1.2965265
## Psi:(Intercept) 1.0044638 0.3250277 0.3674095 1.6415180
## Psi:freshwater -0.6784594 0.5220632 -1.7017033 0.3447845
##
##
## Real Parameter p
## all.diff.index par.index estimate se lcl ucl
## p g1 a0 t1 1 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a1 t2 2 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a2 t3 3 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a3 t4 4 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a4 t5 5 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a5 t6 6 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a6 t7 7 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a7 t8 8 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a8 t9 9 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a9 t10 10 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a10 t11 11 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a11 t12 12 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a12 t13 13 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a13 t14 14 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a14 t15 15 1 0.1851943 0.0142442 0.1588832 0.2147502
## p g1 a15 t16 16 1 0.1851943 0.0142442 0.1588832 0.2147502
## fixed
## p g1 a0 t1
## p g1 a1 t2
## p g1 a2 t3
## p g1 a3 t4
## p g1 a4 t5
## p g1 a5 t6
## p g1 a6 t7
## p g1 a7 t8
## p g1 a8 t9
## p g1 a9 t10
## p g1 a10 t11
## p g1 a11 t12
## p g1 a12 t13
## p g1 a13 t14
## p g1 a14 t15
## p g1 a15 t16
##
##
## Real Parameter Psi
## all.diff.index par.index estimate se lcl ucl
## Psi g1 a0 t1 17 2 0.6870551 0.055529 0.5695938 0.7845837
## fixed
## Psi g1 a0 t1
##
##
## ***** Starting ~1 ~HumanTrail Occupancy 5
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~HumanTrail)
##
## Npar : 3
## -2lnL: 929.9099
## AICc : 936.2099
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.4816925 0.0944098 -1.6667357 -1.296649
## Psi:(Intercept) 0.5754842 0.2823177 0.0221415 1.128827
## Psi:HumanTrail 0.8479414 0.6587810 -0.4432694 2.139152
##
##
## Real Parameter p
## all.diff.index par.index estimate se lcl ucl
## p g1 a0 t1 1 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a1 t2 2 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a2 t3 3 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a3 t4 4 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a4 t5 5 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a5 t6 6 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a6 t7 7 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a7 t8 8 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a8 t9 9 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a9 t10 10 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a10 t11 11 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a11 t12 12 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a12 t13 13 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a13 t14 14 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a14 t15 15 1 0.1851719 0.0142449 0.1588599 0.2147295
## p g1 a15 t16 16 1 0.1851719 0.0142449 0.1588599 0.2147295
## fixed
## p g1 a0 t1
## p g1 a1 t2
## p g1 a2 t3
## p g1 a3 t4
## p g1 a4 t5
## p g1 a5 t6
## p g1 a6 t7
## p g1 a7 t8
## p g1 a8 t9
## p g1 a9 t10
## p g1 a10 t11
## p g1 a11 t12
## p g1 a12 t13
## p g1 a13 t14
## p g1 a14 t15
## p g1 a15 t16
##
##
## Real Parameter Psi
## all.diff.index par.index estimate se lcl ucl
## Psi g1 a0 t1 17 2 0.6894532 0.0561453 0.5704272 0.7877697
## fixed
## Psi g1 a0 t1
##
##
## ***** Starting ~1 ~freshwater+Rabbitat Occupancy 6
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~freshwater + Rabbitat)
##
## Npar : 4
## -2lnL: 925.9356
## AICc : 934.442
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.4811307 0.0943585 -1.6660733 -1.2961880
## Psi:(Intercept) 1.9439842 0.7028450 0.5664081 3.3215604
## Psi:freshwater -0.1617186 0.5623827 -1.2639887 0.9405514
## Psi:Rabbitat -1.4563877 0.7981478 -3.0207575 0.1079820
##
##
## Real Parameter p
## all.diff.index par.index estimate se lcl ucl
## p g1 a0 t1 1 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a1 t2 2 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a2 t3 3 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a3 t4 4 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a4 t5 5 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a5 t6 6 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a6 t7 7 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a7 t8 8 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a8 t9 9 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a9 t10 10 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a10 t11 11 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a11 t12 12 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a12 t13 13 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a13 t14 14 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a14 t15 15 1 0.1852567 0.0142422 0.1589484 0.2148073
## p g1 a15 t16 16 1 0.1852567 0.0142422 0.1589484 0.2148073
## fixed
## p g1 a0 t1
## p g1 a1 t2
## p g1 a2 t3
## p g1 a3 t4
## p g1 a4 t5
## p g1 a5 t6
## p g1 a6 t7
## p g1 a7 t8
## p g1 a8 t9
## p g1 a9 t10
## p g1 a10 t11
## p g1 a11 t12
## p g1 a12 t13
## p g1 a13 t14
## p g1 a14 t15
## p g1 a15 t16
##
##
## Real Parameter Psi
## all.diff.index par.index estimate se lcl ucl
## Psi g1 a0 t1 17 2 0.704551 0.0603981 0.5745312 0.8081084
## fixed
## Psi g1 a0 t1
##
##
## ***** Starting ~1 ~1 Occupancy 7
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~1)
##
## Npar : 2
## -2lnL: 931.8143
## AICc : 935.9624
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.4818738 0.0944267 -1.6669500 -1.296798
## Psi:(Intercept) 0.7703676 0.2531895 0.2741162 1.266619
##
##
## Real Parameter p
## all.diff.index par.index estimate se lcl ucl
## p g1 a0 t1 1 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a1 t2 2 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a2 t3 3 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a3 t4 4 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a4 t5 5 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a5 t6 6 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a6 t7 7 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a7 t8 8 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a8 t9 9 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a9 t10 10 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a10 t11 11 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a11 t12 12 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a12 t13 13 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a13 t14 14 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a14 t15 15 1 0.1851446 0.0142458 0.1588312 0.2147045
## p g1 a15 t16 16 1 0.1851446 0.0142458 0.1588312 0.2147045
## fixed
## p g1 a0 t1
## p g1 a1 t2
## p g1 a2 t3
## p g1 a3 t4
## p g1 a4 t5
## p g1 a5 t6
## p g1 a6 t7
## p g1 a7 t8
## p g1 a8 t9
## p g1 a9 t10
## p g1 a10 t11
## p g1 a11 t12
## p g1 a12 t13
## p g1 a13 t14
## p g1 a14 t15
## p g1 a15 t16
##
##
## Real Parameter Psi
## all.diff.index par.index estimate se lcl ucl
## Psi g1 a0 t1 17 2 0.6836004 0.0547626 0.5681031 0.7801634
## fixed
## Psi g1 a0 t1
##
##
## ***** Starting ~1 ~hectares_z+dist_dev_z Occupancy 8
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~hectares_z + dist_dev_z)
##
## Npar : 4
## -2lnL: 926.8675
## AICc : 935.3738
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.4800881 0.0941291 -1.6645812 -1.2955950
## Psi:(Intercept) 0.8577813 0.2702743 0.3280437 1.3875189
## Psi:hectares_z 0.0039090 0.2529017 -0.4917783 0.4995963
## Psi:dist_dev_z -0.5568887 0.2568993 -1.0604114 -0.0533660
##
##
## Real Parameter p
## all.diff.index par.index estimate se lcl ucl
## p g1 a0 t1 1 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a1 t2 2 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a2 t3 3 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a3 t4 4 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a4 t5 5 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a5 t6 6 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a6 t7 7 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a7 t8 8 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a8 t9 9 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a9 t10 10 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a10 t11 11 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a11 t12 12 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a12 t13 13 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a13 t14 14 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a14 t15 15 1 0.1854141 0.0142169 0.159148 0.2149073
## p g1 a15 t16 16 1 0.1854141 0.0142169 0.159148 0.2149073
## fixed
## p g1 a0 t1
## p g1 a1 t2
## p g1 a2 t3
## p g1 a3 t4
## p g1 a4 t5
## p g1 a5 t6
## p g1 a6 t7
## p g1 a7 t8
## p g1 a8 t9
## p g1 a9 t10
## p g1 a10 t11
## p g1 a11 t12
## p g1 a12 t13
## p g1 a13 t14
## p g1 a14 t15
## p g1 a15 t16
##
##
## Real Parameter Psi
## all.diff.index par.index estimate se lcl ucl
## Psi g1 a0 t1 17 2 0.6919472 0.0566182 0.57166 0.7908158
## fixed
## Psi g1 a0 t1
# examine individula model results
model.number <-1
summary(model.fits[[model.number]])
## Output summary for Occupancy model
## Name : p(~1)Psi(~HumanTrail + Rabbitat + dist_dev_z + freshwater + salt_butt + hectares_z + LiDar_z)
##
## Npar : 9
## -2lnL: 910.8596
## AICc : 931.2921
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -1.4852457 0.0966141 -1.6746093 -1.2958821
## Psi:(Intercept) 1.6277101 1.0614847 -0.4528000 3.7082202
## Psi:HumanTrail 1.7847479 0.9309654 -0.0399443 3.6094402
## Psi:Rabbitat -1.4161672 1.2061910 -3.7803017 0.9479673
## Psi:dist_dev_z -0.9676143 0.3831503 -1.7185889 -0.2166398
## Psi:freshwater 0.4237291 0.9246481 -1.3885812 2.2360394
## Psi:salt_butt 0.6479538 1.4202754 -2.1357861 3.4316937
## Psi:hectares_z -0.2622384 0.3811695 -1.0093305 0.4848538
## Psi:LiDar_z 0.8824997 0.6559389 -0.4031406 2.1681400
##
##
## Real Parameter p
## 1 2 3 4 5 6 7
## 0.1846364 0.1846364 0.1846364 0.1846364 0.1846364 0.1846364 0.1846364
## 8 9 10 11 12 13 14
## 0.1846364 0.1846364 0.1846364 0.1846364 0.1846364 0.1846364 0.1846364
## 15 16
## 0.1846364 0.1846364
##
##
## Real Parameter Psi
## 1
## 0.7569719
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## p g1 a0 t1 0.1846364 0.0145448 0.1578106 0.2148589
## Psi g1 a0 t1 0.7569719 0.0837883 0.5605646 0.8837925
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## p:(Intercept) -1.4852457 0.0966141 -1.6746093 -1.2958821
## Psi:(Intercept) 1.6277101 1.0614847 -0.4528000 3.7082202
## Psi:HumanTrail 1.7847479 0.9309654 -0.0399443 3.6094402
## Psi:Rabbitat -1.4161672 1.2061910 -3.7803017 0.9479673
## Psi:dist_dev_z -0.9676143 0.3831503 -1.7185889 -0.2166398
## Psi:freshwater 0.4237291 0.9246481 -1.3885812 2.2360394
## Psi:salt_butt 0.6479538 1.4202754 -2.1357861 3.4316937
## Psi:hectares_z -0.2622384 0.3811695 -1.0093305 0.4848538
## Psi:LiDar_z 0.8824997 0.6559389 -0.4031406 2.1681400
model.fits[[model.number]]$results$derived
## $Occupancy
## estimate lcl ucl
## 1 0.7569719 0.5605646 0.8837925
model.fits[[model.number]]$results$derived$"psi Probability Occupied"
## NULL
model.fits[[model.number]]$results$derived$"lambda Rate of Change"
## NULL
model.fits[[model.number]]$results$derived$"log odds lambda"
## NULL
get.real(model.fits[[model.number]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi g1 a0 t1 17 2 0.7569719 0.0837883 0.5605646 0.8837925
## fixed note group age time Age Time
## Psi g1 a0 t1 1 0 1 0 0
get.real(model.fits[[model.number]], "p", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## p g1 a0 t1 1 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a1 t2 2 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a2 t3 3 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a3 t4 4 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a4 t5 5 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a5 t6 6 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a6 t7 7 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a7 t8 8 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a8 t9 9 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a9 t10 10 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a10 t11 11 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a11 t12 12 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a12 t13 13 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a13 t14 14 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a14 t15 15 1 0.1846364 0.0145448 0.1578106 0.2148589
## p g1 a15 t16 16 1 0.1846364 0.0145448 0.1578106 0.2148589
## fixed note group age time Age Time
## p g1 a0 t1 1 0 1 0 0
## p g1 a1 t2 1 1 2 1 1
## p g1 a2 t3 1 2 3 2 2
## p g1 a3 t4 1 3 4 3 3
## p g1 a4 t5 1 4 5 4 4
## p g1 a5 t6 1 5 6 5 5
## p g1 a6 t7 1 6 7 6 6
## p g1 a7 t8 1 7 8 7 7
## p g1 a8 t9 1 8 9 8 8
## p g1 a9 t10 1 9 10 9 9
## p g1 a10 t11 1 10 11 10 10
## p g1 a11 t12 1 11 12 11 11
## p g1 a12 t13 1 12 13 12 12
## p g1 a13 t14 1 13 14 13 13
## p g1 a14 t15 1 14 15 14 14
## p g1 a15 t16 1 15 16 15 15
# collect models and make AICc table
model.set <- RMark::collect.models( type=NULL)
model.set
## model
## 1 p(~1)Psi(~HumanTrail + Rabbitat + dist_dev_z + freshwater + salt_butt + hectares_z + LiDar_z)
## 3 p(~1)Psi(~Rabbitat)
## 2 p(~1)Psi(~dist_dev_z)
## 6 p(~1)Psi(~freshwater + Rabbitat)
## 8 p(~1)Psi(~hectares_z + dist_dev_z)
## 7 p(~1)Psi(~1)
## 5 p(~1)Psi(~HumanTrail)
## 4 p(~1)Psi(~freshwater)
## npar AICc DeltaAICc weight Deviance
## 1 9 931.2921 0.000000 0.38679584 910.8596
## 3 3 932.3183 1.026228 0.23154708 926.0183
## 2 3 933.1677 1.875658 0.15142163 926.8677
## 6 4 934.4420 3.149897 0.08007380 925.9356
## 8 4 935.3738 4.081747 0.05025066 926.8675
## 7 2 935.9624 4.670376 0.03743886 348.1174
## 5 3 936.2099 4.917868 0.03308114 929.9099
## 4 3 936.4465 5.154418 0.02939099 930.1465
#####################################################################
#####################################################################
#####################################################################
# Finally model MultiSpecies, Single Season using the 2013/2014 data
# prior to any manipulation of cats
# This was NOT done in the paper and so there is nothing to compare
# the results with.
# let us look at the raw data of naive occupancy
xtabs(~RabbitPresence+(icats1>0), data=input.data, exclude=NULL, na.action=na.pass)
## icats1 > 0
## RabbitPresence FALSE TRUE
## 0 14 36
## 1 16 18
prop.table(xtabs(~RabbitPresence+(icats1>0), data=input.data, exclude=NULL, na.action=na.pass), margin=1)
## icats1 > 0
## RabbitPresence FALSE TRUE
## 0 0.2800000 0.7200000
## 1 0.4705882 0.5294118
prop.table(xtabs(~RabbitPresence+(icats1>0), data=input.data, exclude=NULL, na.action=na.pass), margin=2)
## icats1 > 0
## RabbitPresence FALSE TRUE
## 0 0.4666667 0.6666667
## 1 0.5333333 0.3333333
# MARK wants a 2 "digit" code xy for each visit where
# x = 0/1 if species A is not detected/detected
# and y = 0/1 if species B is not detected/detected
# It is not clear if we should make the rabbits or the cats
# the "dominant species", so we should try if both ways.
# The 2*species code is the dominant species
input.history2 <- 1*rabbit.input.history + 2*cat.input.history
colnames(input.history2) <- NULL
input.history2[1:2,]
##
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 1 0 1 0 0 3 0 2 0 1 0
input.chistory <- data.frame(lapply(input.history2,
car::recode, " 0='00'; 1='10';2='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
## 20 100000000000000000000010000000NA 1
## 21 000000000000000000001000000000NA 1
## 22 000000010000000000000000000000NA 1
## 25 000000000000010101000000000100NA 1
## 30 000000000000000000000000000000NA 1
input.history$ch <- gsub("NA","..", input.history$ch, fixed=TRUE)
input.history[ select,][1:5,]
## ch freq
## 20 100000000000000000000010000000.. 1
## 21 000000000000000000001000000000.. 1
## 22 000000010000000000000000000000.. 1
## 25 000000000000010101000000000100.. 1
## 30 000000000000000000000000000000.. 1
other.covar.vars <- NULL #
setdiff(other.covar.vars, names(input.data))
## NULL
site.covar <- input.data[, c(global.covar.vars,other.covar.vars)]
row.names(site.covar) <- NULL
head(site.covar)
## HumanTrail Rabbitat dist_dev_z freshwater salt_butt hectares_z LiDar_z
## 1 0 1 1.05 0 1 -0.74 -0.35
## 2 0 1 -0.83 1 0 -0.78 -0.88
## 3 0 1 0.12 1 0 -0.78 -1.13
## 4 0 1 0.00 1 0 -0.78 -0.84
## 5 1 1 2.20 0 0 0.51 -0.32
## 6 1 0 1.02 0 0 1.84 -0.07
#Add site covariates to input history
input.history = cbind(input.history,site.covar)
head(input.history)
## ch freq HumanTrail Rabbitat dist_dev_z
## 1 00000000000000000000000000000000 1 0 1 1.05
## 2 00000000001000100000110001001000 1 0 1 -0.83
## 3 00100000000000000010001010100010 1 0 1 0.12
## 4 10000000000000000000000000000000 1 0 1 0.00
## 5 00100000000000000000000000000000 1 1 1 2.20
## 6 00001000000000000000000010000000 1 1 0 1.02
## freshwater salt_butt hectares_z LiDar_z
## 1 0 1 -0.74 -0.35
## 2 1 0 -0.78 -0.88
## 3 1 0 -0.78 -1.13
## 4 1 0 -0.78 -0.84
## 5 0 0 0.51 -0.32
## 6 0 0 1.84 -0.07
catrabbit.data <- process.data(data=input.history,
model="2SpecConOccup") # this parameterization is more stable
# are there any time-specific covariates? If so add them to the ddl's
catrabbit.ddl <- RMark::make.design.data(catrabbit.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
~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 ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1 1
## 2 ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~SHARE 2
## 3 ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~1 3
## 4 ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~SHARE 4
# fit the models
myobs <- ls()
myobs <- myobs[ grepl("m...",myobs,fixed=TRUE)]
cat("Removing ", myobs, "\n")
## Removing m...001 m...002 m...003 m...004 m...005 m...006 m...007 m...008
rm(list=myobs)
cleanup(ask=FALSE)
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=catrabbit.data, input.ddl=catrabbit.ddl)
##
##
## ***** Starting ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1 1
##
## Note: only 7 parameters counted of 8 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(~1)
##
## Npar : 8 (unadjusted=7)
## -2lnL: 1339.072
## AICc : 1356.992 (unadjusted=1354.5455)
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 1.2277412 0.3236003 0.5934846 1.8619978
## PsiBA:(Intercept) 3.0012987 0.9114373 1.2148816 4.7877159
## PsiBa:(Intercept) 11.8791190 3649.8832000 -7141.8921000 7165.6503000
## pA:(Intercept) -0.9503081 0.6110370 -2.1479407 0.2473245
## pB:(Intercept) -0.6199148 0.1721175 -0.9572651 -0.2825645
## rA:(Intercept) -3.0165145 0.1758224 -3.3611264 -2.6719026
## rBA:(Intercept) -2.0696979 0.4862678 -3.0227829 -1.1166129
## rBa:(Intercept) -2.7045988 0.1710577 -3.0398719 -2.3693257
##
##
## Real Parameter PsiA
## 1
## 0.773423
##
##
## Real Parameter PsiBA
## 1
## 0.9526328
##
##
## Real Parameter PsiBa
## 1
## 0.9999931
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.2788229 0.2788229 0.2788229 0.2788229 0.2788229 0.2788229 0.2788229
## 8 9 10 11 12 13 14
## 0.2788229 0.2788229 0.2788229 0.2788229 0.2788229 0.2788229 0.2788229
## 15 16
## 0.2788229 0.2788229
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.3498008 0.3498008 0.3498008 0.3498008 0.3498008 0.3498008 0.3498008
## 8 9 10 11 12 13 14
## 0.3498008 0.3498008 0.3498008 0.3498008 0.3498008 0.3498008 0.3498008
## 15 16
## 0.3498008 0.3498008
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.0466854 0.0466854 0.0466854 0.0466854 0.0466854 0.0466854 0.0466854
## 8 9 10 11 12 13 14
## 0.0466854 0.0466854 0.0466854 0.0466854 0.0466854 0.0466854 0.0466854
## 15 16
## 0.0466854 0.0466854
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.1120771 0.1120771 0.1120771 0.1120771 0.1120771 0.1120771 0.1120771
## 8 9 10 11 12 13 14
## 0.1120771 0.1120771 0.1120771 0.1120771 0.1120771 0.1120771 0.1120771
## 15 16
## 0.1120771 0.1120771
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.0627025 0.0627025 0.0627025 0.0627025 0.0627025 0.0627025 0.0627025
## 8 9 10 11 12 13 14
## 0.0627025 0.0627025 0.0627025 0.0627025 0.0627025 0.0627025 0.0627025
## 15 16
## 0.0627025 0.0627025
##
##
## ***** Starting ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~SHARE 2
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 7
## -2lnL: 1340.432
## AICc : 1355.906
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 1.2034609 0.4319049 0.3569271 2.0499946
## PsiBA:(Intercept) 3.0527099 0.8863354 1.3154924 4.7899273
## PsiBa:(Intercept) 3.1648743 5.7647788 -8.1340924 14.4638410
## pA:(Intercept) -0.9120670 0.5857430 -2.0601233 0.2359894
## pB:(Intercept) -0.6040225 0.1729063 -0.9429189 -0.2651261
## rA:(Intercept) -3.0046578 0.1863118 -3.3698289 -2.6394868
## rBA:(Intercept) -2.6412230 0.1816779 -2.9973117 -2.2851342
##
##
## Real Parameter PsiA
## 1
## 0.7691399
##
##
## Real Parameter PsiBA
## 1
## 0.9548994
##
##
## Real Parameter PsiBa
## 1
## 0.9594908
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.2865771 0.2865771 0.2865771 0.2865771 0.2865771 0.2865771 0.2865771
## 8 9 10 11 12 13 14
## 0.2865771 0.2865771 0.2865771 0.2865771 0.2865771 0.2865771 0.2865771
## 15 16
## 0.2865771 0.2865771
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.3534239 0.3534239 0.3534239 0.3534239 0.3534239 0.3534239 0.3534239
## 8 9 10 11 12 13 14
## 0.3534239 0.3534239 0.3534239 0.3534239 0.3534239 0.3534239 0.3534239
## 15 16
## 0.3534239 0.3534239
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.0472159 0.0472159 0.0472159 0.0472159 0.0472159 0.0472159 0.0472159
## 8 9 10 11 12 13 14
## 0.0472159 0.0472159 0.0472159 0.0472159 0.0472159 0.0472159 0.0472159
## 15 16
## 0.0472159 0.0472159
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7 8
## 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532
## 9 10 11 12 13 14 15 16
## 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7 8
## 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532
## 9 10 11 12 13 14 15 16
## 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532 0.066532
##
##
## ***** Starting ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~1 3
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
##
## Npar : 7
## -2lnL: 1339.114
## AICc : 1354.588
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 1.1698897 0.3284455 0.5261364 1.8136429
## PsiBA:(Intercept) 3.0093018 0.9163917 1.2131740 4.8054296
## pA:(Intercept) -0.9451849 0.6101729 -2.1411239 0.2507541
## pB:(Intercept) -0.6192793 0.1730546 -0.9584664 -0.2800922
## rA:(Intercept) -3.0009521 0.1739264 -3.3418477 -2.6600564
## rBA:(Intercept) -2.0710323 0.4860247 -3.0236408 -1.1184238
## rBa:(Intercept) -2.6871116 0.1747149 -3.0295528 -2.3446704
##
##
## Real Parameter PsiA
## 1
## 0.7631251
##
##
## Real Parameter PsiBA
## 1
## 0.9529926
##
##
## Real Parameter PsiBa
## 1
## 0.9529926
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.2798542 0.2798542 0.2798542 0.2798542 0.2798542 0.2798542 0.2798542
## 8 9 10 11 12 13 14
## 0.2798542 0.2798542 0.2798542 0.2798542 0.2798542 0.2798542 0.2798542
## 15 16
## 0.2798542 0.2798542
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.3499454 0.3499454 0.3499454 0.3499454 0.3499454 0.3499454 0.3499454
## 8 9 10 11 12 13 14
## 0.3499454 0.3499454 0.3499454 0.3499454 0.3499454 0.3499454 0.3499454
## 15 16
## 0.3499454 0.3499454
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.0473829 0.0473829 0.0473829 0.0473829 0.0473829 0.0473829 0.0473829
## 8 9 10 11 12 13 14
## 0.0473829 0.0473829 0.0473829 0.0473829 0.0473829 0.0473829 0.0473829
## 15 16
## 0.0473829 0.0473829
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.1119444 0.1119444 0.1119444 0.1119444 0.1119444 0.1119444 0.1119444
## 8 9 10 11 12 13 14
## 0.1119444 0.1119444 0.1119444 0.1119444 0.1119444 0.1119444 0.1119444
## 15 16
## 0.1119444 0.1119444
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.0637382 0.0637382 0.0637382 0.0637382 0.0637382 0.0637382 0.0637382
## 8 9 10 11 12 13 14
## 0.0637382 0.0637382 0.0637382 0.0637382 0.0637382 0.0637382 0.0637382
## 15 16
## 0.0637382 0.0637382
##
##
## ***** Starting ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~SHARE 4
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 6
## -2lnL: 1340.432
## AICc : 1353.523
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 1.197714 0.3300233 0.5508683 1.8445595
## PsiBA:(Intercept) 3.053326 0.8855898 1.3175700 4.7890821
## pA:(Intercept) -0.911727 0.5854105 -2.0591317 0.2356776
## pB:(Intercept) -0.603978 0.1729881 -0.9430347 -0.2649214
## rA:(Intercept) -3.003141 0.1713816 -3.3390484 -2.6672327
## rBA:(Intercept) -2.639586 0.1636411 -2.9603225 -2.3188495
##
##
## Real Parameter PsiA
## 1
## 0.7681179
##
##
## Real Parameter PsiBA
## 1
## 0.9549259
##
##
## Real Parameter PsiBa
## 1
## 0.9549259
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466
## 8 9 10 11 12 13 14
## 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466
## 15 16
## 0.2866466 0.2866466
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341
## 8 9 10 11 12 13 14
## 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341
## 15 16
## 0.3534341 0.3534341
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842
## 8 9 10 11 12 13 14
## 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842
## 15 16
## 0.0472842 0.0472842
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338
## 8 9 10 11 12 13 14
## 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338
## 15 16
## 0.0666338 0.0666338
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338
## 8 9 10 11 12 13 14
## 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338
## 15 16
## 0.0666338 0.0666338
# examine individual model results
model.number <-4
summary(model.fits[[model.number]])
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 6
## -2lnL: 1340.432
## AICc : 1353.523
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 1.197714 0.3300233 0.5508683 1.8445595
## PsiBA:(Intercept) 3.053326 0.8855898 1.3175700 4.7890821
## pA:(Intercept) -0.911727 0.5854105 -2.0591317 0.2356776
## pB:(Intercept) -0.603978 0.1729881 -0.9430347 -0.2649214
## rA:(Intercept) -3.003141 0.1713816 -3.3390484 -2.6672327
## rBA:(Intercept) -2.639586 0.1636411 -2.9603225 -2.3188495
##
##
## Real Parameter PsiA
## 1
## 0.7681179
##
##
## Real Parameter PsiBA
## 1
## 0.9549259
##
##
## Real Parameter PsiBa
## 1
## 0.9549259
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466
## 8 9 10 11 12 13 14
## 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466 0.2866466
## 15 16
## 0.2866466 0.2866466
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341
## 8 9 10 11 12 13 14
## 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341 0.3534341
## 15 16
## 0.3534341 0.3534341
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842
## 8 9 10 11 12 13 14
## 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842 0.0472842
## 15 16
## 0.0472842 0.0472842
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338
## 8 9 10 11 12 13 14
## 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338
## 15 16
## 0.0666338 0.0666338
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338
## 8 9 10 11 12 13 14
## 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338 0.0666338
## 15 16
## 0.0666338 0.0666338
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## PsiA g1 a0 t1 0.7681179 0.0587814 0.6343370 0.8634871
## PsiBA g1 a0 t1 0.9549259 0.0381179 0.7887771 0.9917486
## pA g1 a0 t1 0.2866466 0.1197049 0.1131329 0.5586482
## pB g1 a0 t1 0.3534341 0.0395310 0.2802878 0.4341543
## rA g1 a0 t1 0.0472842 0.0077205 0.0342556 0.0649348
## rBA g1 a0 t1 0.0666338 0.0101774 0.0492509 0.0895738
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## PsiA:(Intercept) 1.197714 0.3300233 0.5508683 1.8445595
## PsiBA:(Intercept) 3.053326 0.8855898 1.3175700 4.7890821
## pA:(Intercept) -0.911727 0.5854105 -2.0591317 0.2356776
## pB:(Intercept) -0.603978 0.1729881 -0.9430347 -0.2649214
## rA:(Intercept) -3.003141 0.1713816 -3.3390484 -2.6672327
## rBA:(Intercept) -2.639586 0.1636411 -2.9603225 -2.3188495
model.fits[[model.number]]$results$derived
## $`Species Interaction Factor`
## estimate lcl ucl
## 1 1 1 1
##
## $`Occupancy of Species B`
## estimate lcl ucl
## 1 0.9549259 0.7887771 0.9917486
##
## $`Occupancy of Both Species`
## estimate lcl ucl
## 1 0.7334956 0.5866224 0.8422204
get.real(model.fits[[model.number]], "PsiA", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## PsiA g1 a0 t1 1 1 0.7681179 0.0587814 0.634337 0.8634871
## fixed note group age time Age Time
## PsiA g1 a0 t1 1 0 1 0 0
get.real(model.fits[[model.number]], "PsiBA", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## PsiBA g1 a0 t1 2 2 0.9549259 0.0381179 0.7887771 0.9917486
## fixed note group age time Age Time
## PsiBA g1 a0 t1 1 0 1 0 0
get.real(model.fits[[model.number]], "PsiBa", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## PsiBa g1 a0 t1 3 2 0.9549259 0.0381179 0.7887771 0.9917486
## fixed note group age time Age Time
## PsiBa g1 a0 t1 1 0 1 0 0
get.real(model.fits[[model.number]], "pA", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## pA g1 a0 t1 4 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a1 t2 5 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a2 t3 6 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a3 t4 7 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a4 t5 8 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a5 t6 9 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a6 t7 10 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a7 t8 11 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a8 t9 12 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a9 t10 13 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a10 t11 14 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a11 t12 15 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a12 t13 16 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a13 t14 17 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a14 t15 18 3 0.2866466 0.1197049 0.1131329 0.5586482
## pA g1 a15 t16 19 3 0.2866466 0.1197049 0.1131329 0.5586482
## fixed note group age time Age Time
## pA g1 a0 t1 1 0 1 0 0
## pA g1 a1 t2 1 1 2 1 1
## pA g1 a2 t3 1 2 3 2 2
## pA g1 a3 t4 1 3 4 3 3
## pA g1 a4 t5 1 4 5 4 4
## pA g1 a5 t6 1 5 6 5 5
## pA g1 a6 t7 1 6 7 6 6
## pA g1 a7 t8 1 7 8 7 7
## pA g1 a8 t9 1 8 9 8 8
## pA g1 a9 t10 1 9 10 9 9
## pA g1 a10 t11 1 10 11 10 10
## pA g1 a11 t12 1 11 12 11 11
## pA g1 a12 t13 1 12 13 12 12
## pA g1 a13 t14 1 13 14 13 13
## pA g1 a14 t15 1 14 15 14 14
## pA g1 a15 t16 1 15 16 15 15
get.real(model.fits[[model.number]], "pB", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## pB g1 a0 t1 20 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a1 t2 21 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a2 t3 22 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a3 t4 23 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a4 t5 24 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a5 t6 25 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a6 t7 26 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a7 t8 27 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a8 t9 28 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a9 t10 29 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a10 t11 30 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a11 t12 31 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a12 t13 32 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a13 t14 33 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a14 t15 34 4 0.3534341 0.039531 0.2802878 0.4341543
## pB g1 a15 t16 35 4 0.3534341 0.039531 0.2802878 0.4341543
## fixed note group age time Age Time
## pB g1 a0 t1 1 0 1 0 0
## pB g1 a1 t2 1 1 2 1 1
## pB g1 a2 t3 1 2 3 2 2
## pB g1 a3 t4 1 3 4 3 3
## pB g1 a4 t5 1 4 5 4 4
## pB g1 a5 t6 1 5 6 5 5
## pB g1 a6 t7 1 6 7 6 6
## pB g1 a7 t8 1 7 8 7 7
## pB g1 a8 t9 1 8 9 8 8
## pB g1 a9 t10 1 9 10 9 9
## pB g1 a10 t11 1 10 11 10 10
## pB g1 a11 t12 1 11 12 11 11
## pB g1 a12 t13 1 12 13 12 12
## pB g1 a13 t14 1 13 14 13 13
## pB g1 a14 t15 1 14 15 14 14
## pB g1 a15 t16 1 15 16 15 15
get.real(model.fits[[model.number]], "rA", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## rA g1 a0 t1 36 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a1 t2 37 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a2 t3 38 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a3 t4 39 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a4 t5 40 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a5 t6 41 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a6 t7 42 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a7 t8 43 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a8 t9 44 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a9 t10 45 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a10 t11 46 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a11 t12 47 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a12 t13 48 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a13 t14 49 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a14 t15 50 5 0.0472842 0.0077205 0.0342556 0.0649348
## rA g1 a15 t16 51 5 0.0472842 0.0077205 0.0342556 0.0649348
## fixed note group age time Age Time
## rA g1 a0 t1 1 0 1 0 0
## rA g1 a1 t2 1 1 2 1 1
## rA g1 a2 t3 1 2 3 2 2
## rA g1 a3 t4 1 3 4 3 3
## rA g1 a4 t5 1 4 5 4 4
## rA g1 a5 t6 1 5 6 5 5
## rA g1 a6 t7 1 6 7 6 6
## rA g1 a7 t8 1 7 8 7 7
## rA g1 a8 t9 1 8 9 8 8
## rA g1 a9 t10 1 9 10 9 9
## rA g1 a10 t11 1 10 11 10 10
## rA g1 a11 t12 1 11 12 11 11
## rA g1 a12 t13 1 12 13 12 12
## rA g1 a13 t14 1 13 14 13 13
## rA g1 a14 t15 1 14 15 14 14
## rA g1 a15 t16 1 15 16 15 15
get.real(model.fits[[model.number]], "rBA", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## rBA g1 a0 t1 52 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a1 t2 53 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a2 t3 54 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a3 t4 55 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a4 t5 56 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a5 t6 57 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a6 t7 58 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a7 t8 59 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a8 t9 60 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a9 t10 61 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a10 t11 62 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a11 t12 63 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a12 t13 64 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a13 t14 65 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a14 t15 66 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBA g1 a15 t16 67 6 0.0666338 0.0101774 0.0492509 0.0895738
## fixed note group age time Age Time
## rBA g1 a0 t1 1 0 1 0 0
## rBA g1 a1 t2 1 1 2 1 1
## rBA g1 a2 t3 1 2 3 2 2
## rBA g1 a3 t4 1 3 4 3 3
## rBA g1 a4 t5 1 4 5 4 4
## rBA g1 a5 t6 1 5 6 5 5
## rBA g1 a6 t7 1 6 7 6 6
## rBA g1 a7 t8 1 7 8 7 7
## rBA g1 a8 t9 1 8 9 8 8
## rBA g1 a9 t10 1 9 10 9 9
## rBA g1 a10 t11 1 10 11 10 10
## rBA g1 a11 t12 1 11 12 11 11
## rBA g1 a12 t13 1 12 13 12 12
## rBA g1 a13 t14 1 13 14 13 13
## rBA g1 a14 t15 1 14 15 14 14
## rBA g1 a15 t16 1 15 16 15 15
get.real(model.fits[[model.number]], "rBa", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## rBa g1 a0 t1 68 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a1 t2 69 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a2 t3 70 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a3 t4 71 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a4 t5 72 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a5 t6 73 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a6 t7 74 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a7 t8 75 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a8 t9 76 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a9 t10 77 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a10 t11 78 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a11 t12 79 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a12 t13 80 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a13 t14 81 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a14 t15 82 6 0.0666338 0.0101774 0.0492509 0.0895738
## rBa g1 a15 t16 83 6 0.0666338 0.0101774 0.0492509 0.0895738
## fixed note group age time Age Time
## rBa g1 a0 t1 1 0 1 0 0
## rBa g1 a1 t2 1 1 2 1 1
## rBa g1 a2 t3 1 2 3 2 2
## rBa g1 a3 t4 1 3 4 3 3
## rBa g1 a4 t5 1 4 5 4 4
## rBa g1 a5 t6 1 5 6 5 5
## rBa g1 a6 t7 1 6 7 6 6
## rBa g1 a7 t8 1 7 8 7 7
## rBa g1 a8 t9 1 8 9 8 8
## rBa g1 a9 t10 1 9 10 9 9
## rBa g1 a10 t11 1 10 11 10 10
## rBa g1 a11 t12 1 11 12 11 11
## rBa g1 a12 t13 1 12 13 12 12
## rBa g1 a13 t14 1 13 14 13 13
## rBa g1 a14 t15 1 14 15 14 14
## rBa g1 a15 t16 1 15 16 15 15
# collect models and make AICc table
model.set <- RMark::collect.models( type="2SpecConOccup")
model.set
## model npar AICc
## 4 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 6 1353.523
## 3 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 7 1354.587
## 2 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 7 1355.906
## 1 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 8 1356.992
## DeltaAICc weight Deviance
## 4 0.000000 0.48361460 1340.432
## 3 1.064175 0.28406434 1339.114
## 2 2.382375 0.14695123 1340.432
## 1 3.468591 0.08536984 1339.072
names(model.set)
## [1] "m...001" "m...002" "m...003" "m...004" "model.table"
model.set$model.table
## PsiA PsiBA PsiBa pA pB rA rBA rBa
## 4 ~1 ~1 ~1 ~1 ~1 ~1
## 3 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## 2 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## 1 ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## model npar AICc
## 4 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 6 1353.523
## 3 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 7 1354.587
## 2 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 7 1355.906
## 1 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 8 1356.992
## DeltaAICc weight Deviance
## 4 0.000000 0.48361460 1340.432
## 3 1.064175 0.28406434 1339.114
## 2 2.382375 0.14695123 1340.432
## 1 3.468591 0.08536984 1339.072
# model averaged values at average covariate value
mean(input.data$RabbitPresence)
## [1] 0.4047619
PsiA.ma <- RMark::model.average(model.set, param="PsiA")
PsiA.ma
## par.index estimate se fixed note group age time Age
## PsiA g1 a0 t1 1 0.7673027 0.06180711 1 0 1 0
## Time
## PsiA g1 a0 t1 0
PsiBA.ma <- RMark::model.average(model.set, param="PsiBA")
PsiBA.ma
## par.index estimate se fixed note group age time Age
## PsiBA g1 a0 t1 2 0.9541771 0.039254 1 0 1 0
## Time
## PsiBA g1 a0 t1 0
PsiBa.ma <- RMark::model.average(model.set, param="PsiBa")
PsiBa.ma
## par.index estimate se fixed note group age time Age
## PsiBa g1 a0 t1 3 0.9588949 0.09367842 1 0 1 0
## Time
## PsiBa g1 a0 t1 0
# look at detection probabiities
RMark::model.average(model.set, param="pA")
## par.index estimate se fixed note group age time Age
## pA g1 a0 t1 4 0.284039 0.1209681 1 0 1 0
## pA g1 a1 t2 5 0.284039 0.1209681 1 1 2 1
## pA g1 a2 t3 6 0.284039 0.1209681 1 2 3 2
## pA g1 a3 t4 7 0.284039 0.1209681 1 3 4 3
## pA g1 a4 t5 8 0.284039 0.1209681 1 4 5 4
## pA g1 a5 t6 9 0.284039 0.1209681 1 5 6 5
## pA g1 a6 t7 10 0.284039 0.1209681 1 6 7 6
## pA g1 a7 t8 11 0.284039 0.1209681 1 7 8 7
## pA g1 a8 t9 12 0.284039 0.1209681 1 8 9 8
## pA g1 a9 t10 13 0.284039 0.1209681 1 9 10 9
## pA g1 a10 t11 14 0.284039 0.1209681 1 10 11 10
## pA g1 a11 t12 15 0.284039 0.1209681 1 11 12 11
## pA g1 a12 t13 16 0.284039 0.1209681 1 12 13 12
## pA g1 a13 t14 17 0.284039 0.1209681 1 13 14 13
## pA g1 a14 t15 18 0.284039 0.1209681 1 14 15 14
## pA g1 a15 t16 19 0.284039 0.1209681 1 15 16 15
## Time
## pA g1 a0 t1 0
## pA g1 a1 t2 1
## pA g1 a2 t3 2
## pA g1 a3 t4 3
## pA g1 a4 t5 4
## pA g1 a5 t6 5
## pA g1 a6 t7 6
## pA g1 a7 t8 7
## pA g1 a8 t9 8
## pA g1 a9 t10 9
## pA g1 a10 t11 10
## pA g1 a11 t12 11
## pA g1 a12 t13 12
## pA g1 a13 t14 13
## pA g1 a14 t15 14
## pA g1 a15 t16 15
RMark::model.average(model.set, param="pB")
## par.index estimate se fixed note group age time Age
## pB g1 a0 t1 20 0.3521314 0.03948552 1 0 1 0
## pB g1 a1 t2 21 0.3521314 0.03948552 1 1 2 1
## pB g1 a2 t3 22 0.3521314 0.03948552 1 2 3 2
## pB g1 a3 t4 23 0.3521314 0.03948552 1 3 4 3
## pB g1 a4 t5 24 0.3521314 0.03948552 1 4 5 4
## pB g1 a5 t6 25 0.3521314 0.03948552 1 5 6 5
## pB g1 a6 t7 26 0.3521314 0.03948552 1 6 7 6
## pB g1 a7 t8 27 0.3521314 0.03948552 1 7 8 7
## pB g1 a8 t9 28 0.3521314 0.03948552 1 8 9 8
## pB g1 a9 t10 29 0.3521314 0.03948552 1 9 10 9
## pB g1 a10 t11 30 0.3521314 0.03948552 1 10 11 10
## pB g1 a11 t12 31 0.3521314 0.03948552 1 11 12 11
## pB g1 a12 t13 32 0.3521314 0.03948552 1 12 13 12
## pB g1 a13 t14 33 0.3521314 0.03948552 1 13 14 13
## pB g1 a14 t15 34 0.3521314 0.03948552 1 14 15 14
## pB g1 a15 t16 35 0.3521314 0.03948552 1 15 16 15
## Time
## pB g1 a0 t1 0
## pB g1 a1 t2 1
## pB g1 a2 t3 2
## pB g1 a3 t4 3
## pB g1 a4 t5 4
## pB g1 a5 t6 5
## pB g1 a6 t7 6
## pB g1 a7 t8 7
## pB g1 a8 t9 8
## pB g1 a9 t10 9
## pB g1 a10 t11 10
## pB g1 a11 t12 11
## pB g1 a12 t13 12
## pB g1 a13 t14 13
## pB g1 a14 t15 14
## pB g1 a15 t16 15
RMark::model.average(model.set, param="rA")
## par.index estimate se fixed note group age time Age
## rA g1 a0 t1 36 0.04725107 0.007868764 1 0 1 0
## rA g1 a1 t2 37 0.04725107 0.007868764 1 1 2 1
## rA g1 a2 t3 38 0.04725107 0.007868764 1 2 3 2
## rA g1 a3 t4 39 0.04725107 0.007868764 1 3 4 3
## rA g1 a4 t5 40 0.04725107 0.007868764 1 4 5 4
## rA g1 a5 t6 41 0.04725107 0.007868764 1 5 6 5
## rA g1 a6 t7 42 0.04725107 0.007868764 1 6 7 6
## rA g1 a7 t8 43 0.04725107 0.007868764 1 7 8 7
## rA g1 a8 t9 44 0.04725107 0.007868764 1 8 9 8
## rA g1 a9 t10 45 0.04725107 0.007868764 1 9 10 9
## rA g1 a10 t11 46 0.04725107 0.007868764 1 10 11 10
## rA g1 a11 t12 47 0.04725107 0.007868764 1 11 12 11
## rA g1 a12 t13 48 0.04725107 0.007868764 1 12 13 12
## rA g1 a13 t14 49 0.04725107 0.007868764 1 13 14 13
## rA g1 a14 t15 50 0.04725107 0.007868764 1 14 15 14
## rA g1 a15 t16 51 0.04725107 0.007868764 1 15 16 15
## Time
## rA g1 a0 t1 0
## rA g1 a1 t2 1
## rA g1 a2 t3 2
## rA g1 a3 t4 3
## rA g1 a4 t5 4
## rA g1 a5 t6 5
## rA g1 a6 t7 6
## rA g1 a7 t8 7
## rA g1 a8 t9 8
## rA g1 a9 t10 9
## rA g1 a10 t11 10
## rA g1 a11 t12 11
## rA g1 a12 t13 12
## rA g1 a13 t14 13
## rA g1 a14 t15 14
## rA g1 a15 t16 15
RMark::model.average(model.set, param="rBA")
## par.index estimate se fixed note group age time Age
## rBA g1 a0 t1 52 0.08336944 0.03756712 1 0 1 0
## rBA g1 a1 t2 53 0.08336944 0.03756712 1 1 2 1
## rBA g1 a2 t3 54 0.08336944 0.03756712 1 2 3 2
## rBA g1 a3 t4 55 0.08336944 0.03756712 1 3 4 3
## rBA g1 a4 t5 56 0.08336944 0.03756712 1 4 5 4
## rBA g1 a5 t6 57 0.08336944 0.03756712 1 5 6 5
## rBA g1 a6 t7 58 0.08336944 0.03756712 1 6 7 6
## rBA g1 a7 t8 59 0.08336944 0.03756712 1 7 8 7
## rBA g1 a8 t9 60 0.08336944 0.03756712 1 8 9 8
## rBA g1 a9 t10 61 0.08336944 0.03756712 1 9 10 9
## rBA g1 a10 t11 62 0.08336944 0.03756712 1 10 11 10
## rBA g1 a11 t12 63 0.08336944 0.03756712 1 11 12 11
## rBA g1 a12 t13 64 0.08336944 0.03756712 1 12 13 12
## rBA g1 a13 t14 65 0.08336944 0.03756712 1 13 14 13
## rBA g1 a14 t15 66 0.08336944 0.03756712 1 14 15 14
## rBA g1 a15 t16 67 0.08336944 0.03756712 1 15 16 15
## Time
## rBA g1 a0 t1 0
## rBA g1 a1 t2 1
## rBA g1 a2 t3 2
## rBA g1 a3 t4 3
## rBA g1 a4 t5 4
## rBA g1 a5 t6 5
## rBA g1 a6 t7 6
## rBA g1 a7 t8 7
## rBA g1 a8 t9 8
## rBA g1 a9 t10 9
## rBA g1 a10 t11 10
## rBA g1 a11 t12 11
## rBA g1 a12 t13 12
## rBA g1 a13 t14 13
## rBA g1 a14 t15 14
## rBA g1 a15 t16 15
RMark::model.average(model.set, param="rBa")
## par.index estimate se fixed note group age time Age
## rBa g1 a0 t1 68 0.06546068 0.01051834 1 0 1 0
## rBa g1 a1 t2 69 0.06546068 0.01051834 1 1 2 1
## rBa g1 a2 t3 70 0.06546068 0.01051834 1 2 3 2
## rBa g1 a3 t4 71 0.06546068 0.01051834 1 3 4 3
## rBa g1 a4 t5 72 0.06546068 0.01051834 1 4 5 4
## rBa g1 a5 t6 73 0.06546068 0.01051834 1 5 6 5
## rBa g1 a6 t7 74 0.06546068 0.01051834 1 6 7 6
## rBa g1 a7 t8 75 0.06546068 0.01051834 1 7 8 7
## rBa g1 a8 t9 76 0.06546068 0.01051834 1 8 9 8
## rBa g1 a9 t10 77 0.06546068 0.01051834 1 9 10 9
## rBa g1 a10 t11 78 0.06546068 0.01051834 1 10 11 10
## rBa g1 a11 t12 79 0.06546068 0.01051834 1 11 12 11
## rBa g1 a12 t13 80 0.06546068 0.01051834 1 12 13 12
## rBa g1 a13 t14 81 0.06546068 0.01051834 1 13 14 13
## rBa g1 a14 t15 82 0.06546068 0.01051834 1 14 15 14
## rBa g1 a15 t16 83 0.06546068 0.01051834 1 15 16 15
## Time
## rBa g1 a0 t1 0
## rBa g1 a1 t2 1
## rBa g1 a2 t3 2
## rBa g1 a3 t4 3
## rBa g1 a4 t5 4
## rBa g1 a5 t6 5
## rBa g1 a6 t7 6
## rBa g1 a7 t8 7
## rBa g1 a8 t9 8
## rBa g1 a9 t10 9
## rBa g1 a10 t11 10
## rBa g1 a11 t12 11
## rBa g1 a12 t13 12
## rBa g1 a13 t14 13
## rBa g1 a14 t15 14
## rBa g1 a15 t16 15
cleanup(ask=FALSE)