# Multi species single season with a list of models
# Using RMark
# Co-occurence of spotted (SO) and barred owls (BO)
# 2018-12-13 code contributed by Neil Faught
library(car)
## Loading required package: carData
library(ggplot2)
library(readxl)
library(reshape2)
library(RMark)
## This is RMark 2.2.5
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
# get the data
SO.data <- readxl::read_excel(file.path("..","SpottedAndBarredOwls.xls"),
sheet="RawData", na='-',
col_names=FALSE,
range = "B11:K161")
BO.data <- readxl::read_excel(file.path("..","SpottedAndBarredOwls.xls"),
sheet="RawData", na='-',
col_names=FALSE,
range = "M11:V161")
# 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
input.data <- 2*BO.data + SO.data
input.data
## X__1 X__2 X__3 X__4 X__5 X__6 X__7 X__8 X__9 X__10
## 1 0 1 1 1 NA NA NA NA NA NA
## 2 1 1 1 1 NA NA NA NA NA NA
## 3 0 3 2 0 NA NA NA NA NA NA
## 4 1 1 1 1 NA NA NA NA NA NA
## 5 0 0 0 1 NA NA NA NA NA NA
## 6 0 1 1 NA NA NA NA NA NA NA
## 7 1 1 1 NA NA NA NA NA NA NA
## 8 0 0 0 0 NA NA NA NA NA NA
## 9 0 2 2 2 NA NA NA NA NA NA
## 10 1 1 1 NA NA NA NA NA NA NA
## 11 1 1 1 NA NA NA NA NA NA NA
## 12 1 1 1 1 1 0 1 NA NA NA
## 13 0 0 0 0 NA NA NA NA NA NA
## 14 1 1 1 1 1 NA NA NA NA NA
## 15 1 0 0 0 0 0 0 1 1 NA
## 16 1 1 1 0 1 1 NA NA NA NA
## 17 1 0 0 0 NA NA NA NA NA NA
## 18 0 0 0 0 NA NA NA NA NA NA
## 19 1 1 1 1 1 1 1 1 1 1
## 20 0 0 0 0 NA NA NA NA NA NA
## 21 1 1 1 NA NA NA NA NA NA NA
## 22 1 1 NA NA NA NA NA NA NA NA
## 23 1 0 0 1 0 0 1 0 1 NA
## 24 0 0 0 0 NA NA NA NA NA NA
## 25 0 0 0 0 NA NA NA NA NA NA
## 26 0 1 1 1 NA NA NA NA NA NA
## 27 0 0 0 0 NA NA NA NA NA NA
## 28 3 1 0 0 0 0 NA NA NA NA
## 29 0 0 0 0 NA NA NA NA NA NA
## 30 0 1 0 1 1 NA NA NA NA NA
## 31 0 0 0 0 NA NA NA NA NA NA
## 32 0 0 1 1 1 NA NA NA NA NA
## 33 0 0 1 1 0 1 1 0 0 NA
## 34 0 0 0 0 NA NA NA NA NA NA
## 35 1 0 0 0 1 1 1 NA NA NA
## 36 1 1 0 1 1 NA NA NA NA NA
## 37 0 1 0 0 1 NA NA NA NA NA
## 38 1 1 1 NA NA NA NA NA NA NA
## 39 1 1 1 0 1 1 NA NA NA NA
## 40 0 1 0 1 1 1 NA NA NA NA
## 41 0 0 0 0 NA NA NA NA NA NA
## 42 0 0 1 0 0 NA NA NA NA NA
## 43 0 1 1 1 1 1 1 1 1 1
## 44 1 1 0 1 1 NA NA NA NA NA
## 45 2 0 0 1 0 NA NA NA NA NA
## 46 0 0 2 0 2 NA NA NA NA NA
## 47 1 1 1 1 NA NA NA NA NA NA
## 48 0 1 1 1 NA NA NA NA NA NA
## 49 1 1 3 NA NA NA NA NA NA NA
## 50 1 0 0 0 2 0 0 0 NA NA
## 51 0 0 3 NA NA NA NA NA NA NA
## 52 1 1 0 2 1 NA NA NA NA NA
## 53 1 1 1 1 1 1 NA NA NA NA
## 54 0 0 3 2 0 2 1 0 1 0
## 55 0 1 1 1 1 1 NA NA NA NA
## 56 1 0 2 0 NA NA NA NA NA NA
## 57 0 0 0 0 NA NA NA NA NA NA
## 58 1 1 1 0 0 1 1 NA NA NA
## 59 0 0 0 0 NA NA NA NA NA NA
## 60 1 0 1 1 0 0 NA NA NA NA
## 61 0 1 0 2 NA NA NA NA NA NA
## 62 0 1 1 1 NA NA NA NA NA NA
## 63 0 1 1 0 0 NA NA NA NA NA
## 64 2 0 0 0 NA NA NA NA NA NA
## 65 2 1 1 2 2 NA NA NA NA NA
## 66 0 1 1 NA NA NA NA NA NA NA
## 67 2 2 2 2 0 NA NA NA NA NA
## 68 2 0 0 2 NA NA NA NA NA NA
## 69 0 0 0 0 NA NA NA NA NA NA
## 70 1 1 1 1 1 1 1 1 NA NA
## 71 0 0 0 0 NA NA NA NA NA NA
## 72 1 0 1 1 1 1 0 1 0 1
## 73 0 2 2 2 0 NA NA NA NA NA
## 74 1 1 0 0 2 1 0 0 3 2
## 75 0 1 0 0 1 NA NA NA NA NA
## 76 2 2 2 2 NA NA NA NA NA NA
## 77 0 0 0 0 NA NA NA NA NA NA
## 78 1 0 0 0 0 0 NA NA NA NA
## 79 0 0 0 0 NA NA NA NA NA NA
## 80 0 0 0 0 NA NA NA NA NA NA
## 81 2 0 0 0 NA NA NA NA NA NA
## 82 0 0 0 0 NA NA NA NA NA NA
## 83 0 2 0 0 NA NA NA NA NA NA
## 84 0 0 2 2 NA NA NA NA NA NA
## 85 1 1 0 1 1 1 1 NA NA NA
## 86 1 0 0 0 NA NA NA NA NA NA
## 87 0 0 0 1 NA NA NA NA NA NA
## 88 1 1 1 1 1 NA NA NA NA NA
## 89 0 2 NA NA NA NA NA NA NA NA
## 90 0 0 1 0 NA NA NA NA NA NA
## 91 0 2 0 0 NA NA NA NA NA NA
## 92 1 0 1 1 1 1 NA NA NA NA
## 93 1 1 1 1 NA NA NA NA NA NA
## 94 0 0 0 0 1 1 1 NA NA NA
## 95 1 0 0 0 1 1 1 1 NA NA
## 96 0 1 0 0 NA NA NA NA NA NA
## 97 0 0 0 1 1 NA NA NA NA NA
## 98 1 1 1 NA NA NA NA NA NA NA
## 99 2 0 0 0 NA NA NA NA NA NA
## 100 2 1 1 1 1 NA NA NA NA NA
## 101 0 3 0 0 NA NA NA NA NA NA
## 102 0 0 1 1 0 NA NA NA NA NA
## 103 0 1 1 NA NA NA NA NA NA NA
## 104 1 1 1 1 NA NA NA NA NA NA
## 105 1 0 1 2 0 0 NA NA NA NA
## 106 1 0 0 1 0 1 0 0 1 1
## 107 1 1 1 NA NA NA NA NA NA NA
## 108 1 0 0 0 1 0 1 NA NA NA
## 109 0 2 0 0 NA NA NA NA NA NA
## 110 0 1 0 1 1 NA NA NA NA NA
## 111 0 0 0 0 NA NA NA NA NA NA
## 112 0 0 1 1 1 1 1 1 NA NA
## 113 2 0 0 0 NA NA NA NA NA NA
## 114 0 0 1 1 0 0 0 NA NA NA
## 115 0 2 3 2 NA NA NA NA NA NA
## 116 0 0 0 2 NA NA NA NA NA NA
## 117 0 1 1 1 NA NA NA NA NA NA
## 118 0 1 1 3 1 0 0 NA NA NA
## 119 3 1 1 1 1 1 NA NA NA NA
## 120 1 1 0 0 1 NA NA NA NA NA
## 121 0 1 1 0 0 NA NA NA NA NA
## 122 0 0 0 0 NA NA NA NA NA NA
## 123 0 1 0 0 NA NA NA NA NA NA
## 124 0 0 0 0 NA NA NA NA NA NA
## 125 0 0 0 0 NA NA NA NA NA NA
## 126 1 1 1 1 0 0 1 NA NA NA
## 127 1 0 1 1 1 NA NA NA NA NA
## 128 0 0 0 1 1 NA NA NA NA NA
## 129 1 1 1 NA NA NA NA NA NA NA
## 130 0 0 0 2 NA NA NA NA NA NA
## 131 1 0 2 3 1 NA NA NA NA NA
## 132 0 0 0 0 NA NA NA NA NA NA
## 133 1 0 0 0 0 0 1 1 0 NA
## 134 1 1 1 1 1 1 NA NA NA NA
## 135 0 2 0 0 NA NA NA NA NA NA
## 136 0 0 0 0 0 NA NA NA NA NA
## 137 0 2 2 0 NA NA NA NA NA NA
## 138 0 0 0 0 NA NA NA NA NA NA
## 139 1 0 1 1 NA NA NA NA NA NA
## 140 0 0 0 1 NA NA NA NA NA NA
## 141 0 0 0 0 NA NA NA NA NA NA
## 142 0 2 0 0 NA NA NA NA NA NA
## 143 1 1 NA NA NA NA NA NA NA NA
## 144 1 0 1 1 1 0 3 NA NA NA
## 145 0 1 0 1 0 NA NA NA NA NA
## 146 0 0 0 0 0 NA NA NA NA NA
## 147 0 0 0 0 NA NA NA NA NA NA
## 148 0 0 0 0 NA NA NA NA NA NA
## 149 0 1 1 1 1 1 1 NA NA NA
## 150 0 NA NA NA NA NA NA NA NA NA
## 151 0 1 1 0 1 1 NA NA NA NA
input.chistory <- data.frame(lapply(input.data, as.character), stringsAsFactors=FALSE)
input.chistory <- data.frame(lapply(input.chistory,
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,]
## ch freq
## 1 00101010NANANANANANA 1
## 2 10101010NANANANANANA 1
## 3 00110100NANANANANANA 1
## 4 10101010NANANANANANA 1
## 5 00000010NANANANANANA 1
## 6 001010NANANANANANANA 1
## 7 101010NANANANANANANA 1
## 8 00000000NANANANANANA 1
## 9 00010101NANANANANANA 1
## 10 101010NANANANANANANA 1
## 11 101010NANANANANANANA 1
## 12 10101010100010NANANA 1
## 13 00000000NANANANANANA 1
## 14 1010101010NANANANANA 1
## 15 100000000000001010NA 1
## 16 101010001010NANANANA 1
## 17 10000000NANANANANANA 1
## 18 00000000NANANANANANA 1
## 20 00000000NANANANANANA 1
## 21 101010NANANANANANANA 1
## 22 1010NANANANANANANANA 1
## 23 100000100000100010NA 1
## 24 00000000NANANANANANA 1
## 25 00000000NANANANANANA 1
## 26 00101010NANANANANANA 1
## 27 00000000NANANANANANA 1
## 28 111000000000NANANANA 1
## 29 00000000NANANANANANA 1
## 30 0010001010NANANANANA 1
## 31 00000000NANANANANANA 1
## 32 0000101010NANANANANA 1
## 33 000010100010100000NA 1
## 34 00000000NANANANANANA 1
## 35 10000000101010NANANA 1
## 36 1010001010NANANANANA 1
## 37 0010000010NANANANANA 1
## 38 101010NANANANANANANA 1
## 39 101010001010NANANANA 1
## 40 001000101010NANANANA 1
## 41 00000000NANANANANANA 1
## 42 0000100000NANANANANA 1
## 44 1010001010NANANANANA 1
## 45 0100001000NANANANANA 1
## 46 0000010001NANANANANA 1
## 47 10101010NANANANANANA 1
## 48 00101010NANANANANANA 1
## 49 101011NANANANANANANA 1
## 50 1000000001000000NANA 1
## 51 000011NANANANANANANA 1
## 52 1010000110NANANANANA 1
## 53 101010101010NANANANA 1
## 55 001010101010NANANANA 1
## 56 10000100NANANANANANA 1
## 57 00000000NANANANANANA 1
## 58 10101000001010NANANA 1
## 59 00000000NANANANANANA 1
## 60 100010100000NANANANA 1
## 61 00100001NANANANANANA 1
## 62 00101010NANANANANANA 1
## 63 0010100000NANANANANA 1
## 64 01000000NANANANANANA 1
## 65 0110100101NANANANANA 1
## 66 001010NANANANANANANA 1
## 67 0101010100NANANANANA 1
## 68 01000001NANANANANANA 1
## 69 00000000NANANANANANA 1
## 70 1010101010101010NANA 1
## 71 00000000NANANANANANA 1
## 73 0001010100NANANANANA 1
## 75 0010000010NANANANANA 1
## 76 01010101NANANANANANA 1
## 77 00000000NANANANANANA 1
## 78 100000000000NANANANA 1
## 79 00000000NANANANANANA 1
## 80 00000000NANANANANANA 1
## 81 01000000NANANANANANA 1
## 82 00000000NANANANANANA 1
## 83 00010000NANANANANANA 1
## 84 00000101NANANANANANA 1
## 85 10100010101010NANANA 1
## 86 10000000NANANANANANA 1
## 87 00000010NANANANANANA 1
## 88 1010101010NANANANANA 1
## 89 0001NANANANANANANANA 1
## 90 00001000NANANANANANA 1
## 91 00010000NANANANANANA 1
## 92 100010101010NANANANA 1
## 93 10101010NANANANANANA 1
## 94 00000000101010NANANA 1
## 95 1000000010101010NANA 1
## 96 00100000NANANANANANA 1
## 97 0000001010NANANANANA 1
## 98 101010NANANANANANANA 1
## 99 01000000NANANANANANA 1
## 100 0110101010NANANANANA 1
## 101 00110000NANANANANANA 1
## 102 0000101000NANANANANA 1
## 103 001010NANANANANANANA 1
## 104 10101010NANANANANANA 1
## 105 100010010000NANANANA 1
## 107 101010NANANANANANANA 1
## 108 10000000100010NANANA 1
## 109 00010000NANANANANANA 1
## 110 0010001010NANANANANA 1
## 111 00000000NANANANANANA 1
## 112 0000101010101010NANA 1
## 113 01000000NANANANANANA 1
## 114 00001010000000NANANA 1
## 115 00011101NANANANANANA 1
## 116 00000001NANANANANANA 1
## 117 00101010NANANANANANA 1
## 118 00101011100000NANANA 1
## 119 111010101010NANANANA 1
## 120 1010000010NANANANANA 1
## 121 0010100000NANANANANA 1
## 122 00000000NANANANANANA 1
## 123 00100000NANANANANANA 1
## 124 00000000NANANANANANA 1
## 125 00000000NANANANANANA 1
## 126 10101010000010NANANA 1
## 127 1000101010NANANANANA 1
## 128 0000001010NANANANANA 1
## 129 101010NANANANANANANA 1
## 130 00000001NANANANANANA 1
## 131 1000011110NANANANANA 1
## 132 00000000NANANANANANA 1
## 133 100000000000101000NA 1
## 134 101010101010NANANANA 1
## 135 00010000NANANANANANA 1
## 136 0000000000NANANANANA 1
## 137 00010100NANANANANANA 1
## 138 00000000NANANANANANA 1
## 139 10001010NANANANANANA 1
## 140 00000010NANANANANANA 1
## 141 00000000NANANANANANA 1
## 142 00010000NANANANANANA 1
## 143 1010NANANANANANANANA 1
## 144 10001010100011NANANA 1
## 145 0010001000NANANANANA 1
## 146 0000000000NANANANANA 1
## 147 00000000NANANANANANA 1
## 148 00000000NANANANANANA 1
## 149 00101010101010NANANA 1
## 150 00NANANANANANANANANA 1
## 151 001010001010NANANANA 1
input.history$ch <- gsub("NA","..", input.history$ch, fixed=TRUE)
input.history[ select,]
## ch freq
## 1 00101010............ 1
## 2 10101010............ 1
## 3 00110100............ 1
## 4 10101010............ 1
## 5 00000010............ 1
## 6 001010.............. 1
## 7 101010.............. 1
## 8 00000000............ 1
## 9 00010101............ 1
## 10 101010.............. 1
## 11 101010.............. 1
## 12 10101010100010...... 1
## 13 00000000............ 1
## 14 1010101010.......... 1
## 15 100000000000001010.. 1
## 16 101010001010........ 1
## 17 10000000............ 1
## 18 00000000............ 1
## 20 00000000............ 1
## 21 101010.............. 1
## 22 1010................ 1
## 23 100000100000100010.. 1
## 24 00000000............ 1
## 25 00000000............ 1
## 26 00101010............ 1
## 27 00000000............ 1
## 28 111000000000........ 1
## 29 00000000............ 1
## 30 0010001010.......... 1
## 31 00000000............ 1
## 32 0000101010.......... 1
## 33 000010100010100000.. 1
## 34 00000000............ 1
## 35 10000000101010...... 1
## 36 1010001010.......... 1
## 37 0010000010.......... 1
## 38 101010.............. 1
## 39 101010001010........ 1
## 40 001000101010........ 1
## 41 00000000............ 1
## 42 0000100000.......... 1
## 44 1010001010.......... 1
## 45 0100001000.......... 1
## 46 0000010001.......... 1
## 47 10101010............ 1
## 48 00101010............ 1
## 49 101011.............. 1
## 50 1000000001000000.... 1
## 51 000011.............. 1
## 52 1010000110.......... 1
## 53 101010101010........ 1
## 55 001010101010........ 1
## 56 10000100............ 1
## 57 00000000............ 1
## 58 10101000001010...... 1
## 59 00000000............ 1
## 60 100010100000........ 1
## 61 00100001............ 1
## 62 00101010............ 1
## 63 0010100000.......... 1
## 64 01000000............ 1
## 65 0110100101.......... 1
## 66 001010.............. 1
## 67 0101010100.......... 1
## 68 01000001............ 1
## 69 00000000............ 1
## 70 1010101010101010.... 1
## 71 00000000............ 1
## 73 0001010100.......... 1
## 75 0010000010.......... 1
## 76 01010101............ 1
## 77 00000000............ 1
## 78 100000000000........ 1
## 79 00000000............ 1
## 80 00000000............ 1
## 81 01000000............ 1
## 82 00000000............ 1
## 83 00010000............ 1
## 84 00000101............ 1
## 85 10100010101010...... 1
## 86 10000000............ 1
## 87 00000010............ 1
## 88 1010101010.......... 1
## 89 0001................ 1
## 90 00001000............ 1
## 91 00010000............ 1
## 92 100010101010........ 1
## 93 10101010............ 1
## 94 00000000101010...... 1
## 95 1000000010101010.... 1
## 96 00100000............ 1
## 97 0000001010.......... 1
## 98 101010.............. 1
## 99 01000000............ 1
## 100 0110101010.......... 1
## 101 00110000............ 1
## 102 0000101000.......... 1
## 103 001010.............. 1
## 104 10101010............ 1
## 105 100010010000........ 1
## 107 101010.............. 1
## 108 10000000100010...... 1
## 109 00010000............ 1
## 110 0010001010.......... 1
## 111 00000000............ 1
## 112 0000101010101010.... 1
## 113 01000000............ 1
## 114 00001010000000...... 1
## 115 00011101............ 1
## 116 00000001............ 1
## 117 00101010............ 1
## 118 00101011100000...... 1
## 119 111010101010........ 1
## 120 1010000010.......... 1
## 121 0010100000.......... 1
## 122 00000000............ 1
## 123 00100000............ 1
## 124 00000000............ 1
## 125 00000000............ 1
## 126 10101010000010...... 1
## 127 1000101010.......... 1
## 128 0000001010.......... 1
## 129 101010.............. 1
## 130 00000001............ 1
## 131 1000011110.......... 1
## 132 00000000............ 1
## 133 100000000000101000.. 1
## 134 101010101010........ 1
## 135 00010000............ 1
## 136 0000000000.......... 1
## 137 00010100............ 1
## 138 00000000............ 1
## 139 10001010............ 1
## 140 00000010............ 1
## 141 00000000............ 1
## 142 00010000............ 1
## 143 1010................ 1
## 144 10001010100011...... 1
## 145 0010001000.......... 1
## 146 0000000000.......... 1
## 147 00000000............ 1
## 148 00000000............ 1
## 149 00101010101010...... 1
## 150 00.................. 1
## 151 001010001010........ 1
# Read in survey level covariates
night = readxl::read_excel(file.path("..","SpottedAndBarredOwls.xls"),
sheet="Nite Covariate", na='-',
col_names=FALSE,
range = "B11:K161")
head(night)
## # A tibble: 6 x 10
## X__1 X__2 X__3 X__4 X__5 X__6 X__7 X__8 X__9 X__10
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0 0 NA NA NA NA NA NA
## 2 0 0 1 0 NA NA NA NA NA NA
## 3 0 0 0 0 NA NA NA NA NA NA
## 4 0 0 1 1 NA NA NA NA NA NA
## 5 0 0 0 0 NA NA NA NA NA NA
## 6 0 0 0 NA NA NA NA NA NA NA
night[is.na(night)] = 0 # all missing values never used, but RMarks expects a non-missing value when making time varying factor
# recode night as N/D to make your life easier in the long run rather than a simple 0/1 variable
# make each column a factor variable.
night[] <- lapply(night, car::recode, "0='d'; 1='n'")
night[] <- lapply(night, as.factor)
colnames(night) = paste("night", 1:10, sep="")
# tibbles (from read_excel) make the make.time.factor routine upset.
night <- data.frame(night)
head(night)
## night1 night2 night3 night4 night5 night6 night7 night8 night9 night10
## 1 d d d d d d d d d d
## 2 d d n d d d d d d d
## 3 d d d d d d d d d d
## 4 d d n n d d d d d d
## 5 d d d d d d d d d d
## 6 d d d d d d d d d d
str(night)
## 'data.frame': 151 obs. of 10 variables:
## $ night1 : Factor w/ 2 levels "d","n": 1 1 1 1 1 1 1 1 1 1 ...
## $ night2 : Factor w/ 2 levels "d","n": 1 1 1 1 1 1 1 1 1 1 ...
## $ night3 : Factor w/ 2 levels "d","n": 1 2 1 2 1 1 1 1 1 2 ...
## $ night4 : Factor w/ 2 levels "d","n": 1 1 1 2 1 1 1 1 1 1 ...
## $ night5 : Factor w/ 2 levels "d","n": 1 1 1 1 1 1 1 1 1 1 ...
## $ night6 : Factor w/ 2 levels "d","n": 1 1 1 1 1 1 1 1 1 1 ...
## $ night7 : Factor w/ 2 levels "d","n": 1 1 1 1 1 1 1 1 1 1 ...
## $ night8 : Factor w/ 2 levels "d","n": 1 1 1 1 1 1 1 1 1 1 ...
## $ night9 : Factor w/ 2 levels "d","n": 1 1 1 1 1 1 1 1 1 1 ...
## $ night10: Factor w/ 2 levels "d","n": 1 1 1 1 1 1 1 1 1 1 ...
# does every night variable have at least two levels?
lapply(night, table)
## $night1
##
## d n
## 146 5
##
## $night2
##
## d n
## 139 12
##
## $night3
##
## d n
## 131 20
##
## $night4
##
## d n
## 136 15
##
## $night5
##
## d n
## 136 15
##
## $night6
##
## d n
## 144 7
##
## $night7
##
## d n
## 144 7
##
## $night8
##
## d n
## 148 3
##
## $night9
##
## d n
## 149 2
##
## $night10
##
## d n
## 149 2
# paste the time-varying covariate (night) with the input data
night.new <- make.time.factor(night, "night", 1:10, intercept="d",delete=TRUE)
input.history = cbind(input.history, night.new)
head(input.history)
## ch freq nightn1 nightn2 nightn3 nightn4 nightn5
## 1 00101010............ 1 0 0 0 0 0
## 2 10101010............ 1 0 0 1 0 0
## 3 00110100............ 1 0 0 0 0 0
## 4 10101010............ 1 0 0 1 1 0
## 5 00000010............ 1 0 0 0 0 0
## 6 001010.............. 1 0 0 0 0 0
## nightn6 nightn7 nightn8 nightn9 nightn10
## 1 0 0 0 0 0
## 2 0 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 0
## 5 0 0 0 0 0
## 6 0 0 0 0 0
owl.data <- process.data(data=input.history,
model="2SpecConOccup") # this parameterization is more stable
summary(owl.data)
## Length Class Mode
## data 12 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 151 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 10 -none- numeric
## begin.time 1 -none- numeric
## age.unit 1 -none- numeric
## initial.ages 1 -none- numeric
## group.covariates 0 -none- NULL
## nstrata 1 -none- numeric
## strata.labels 0 -none- NULL
## counts 0 -none- NULL
## reverse 1 -none- logical
## areas 0 -none- NULL
## events 0 -none- NULL
owl.ddl <- RMark::make.design.data(owl.data) # need for covariate predictions
owl.ddl
## $PsiA
## par.index model.index group age time Age Time
## 1 1 1 1 0 1 0 0
##
## $PsiBA
## par.index model.index group age time Age Time
## 1 1 2 1 0 1 0 0
##
## $PsiBa
## par.index model.index group age time Age Time
## 1 1 3 1 0 1 0 0
##
## $pA
## par.index model.index group age time Age Time
## 1 1 4 1 0 1 0 0
## 2 2 5 1 1 2 1 1
## 3 3 6 1 2 3 2 2
## 4 4 7 1 3 4 3 3
## 5 5 8 1 4 5 4 4
## 6 6 9 1 5 6 5 5
## 7 7 10 1 6 7 6 6
## 8 8 11 1 7 8 7 7
## 9 9 12 1 8 9 8 8
## 10 10 13 1 9 10 9 9
##
## $pB
## par.index model.index group age time Age Time
## 1 1 14 1 0 1 0 0
## 2 2 15 1 1 2 1 1
## 3 3 16 1 2 3 2 2
## 4 4 17 1 3 4 3 3
## 5 5 18 1 4 5 4 4
## 6 6 19 1 5 6 5 5
## 7 7 20 1 6 7 6 6
## 8 8 21 1 7 8 7 7
## 9 9 22 1 8 9 8 8
## 10 10 23 1 9 10 9 9
##
## $rA
## par.index model.index group age time Age Time
## 1 1 24 1 0 1 0 0
## 2 2 25 1 1 2 1 1
## 3 3 26 1 2 3 2 2
## 4 4 27 1 3 4 3 3
## 5 5 28 1 4 5 4 4
## 6 6 29 1 5 6 5 5
## 7 7 30 1 6 7 6 6
## 8 8 31 1 7 8 7 7
## 9 9 32 1 8 9 8 8
## 10 10 33 1 9 10 9 9
##
## $rBA
## par.index model.index group age time Age Time
## 1 1 34 1 0 1 0 0
## 2 2 35 1 1 2 1 1
## 3 3 36 1 2 3 2 2
## 4 4 37 1 3 4 3 3
## 5 5 38 1 4 5 4 4
## 6 6 39 1 5 6 5 5
## 7 7 40 1 6 7 6 6
## 8 8 41 1 7 8 7 7
## 9 9 42 1 8 9 8 8
## 10 10 43 1 9 10 9 9
##
## $rBa
## par.index model.index group age time Age Time
## 1 1 44 1 0 1 0 0
## 2 2 45 1 1 2 1 1
## 3 3 46 1 2 3 2 2
## 4 4 47 1 3 4 3 3
## 5 5 48 1 4 5 4 4
## 6 6 49 1 5 6 5 5
## 7 7 50 1 6 7 6 6
## 8 8 51 1 7 8 7 7
## 9 9 52 1 8 9 8 8
## 10 10 53 1 9 10 9 9
##
## $pimtypes
## $pimtypes$PsiA
## $pimtypes$PsiA$pim.type
## [1] "all"
##
##
## $pimtypes$PsiBA
## $pimtypes$PsiBA$pim.type
## [1] "all"
##
##
## $pimtypes$PsiBa
## $pimtypes$PsiBa$pim.type
## [1] "all"
##
##
## $pimtypes$pA
## $pimtypes$pA$pim.type
## [1] "all"
##
##
## $pimtypes$pB
## $pimtypes$pB$pim.type
## [1] "all"
##
##
## $pimtypes$rA
## $pimtypes$rA$pim.type
## [1] "all"
##
##
## $pimtypes$rBA
## $pimtypes$rBA$pim.type
## [1] "all"
##
##
## $pimtypes$rBa
## $pimtypes$rBa$pim.type
## [1] "all"
# 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
# Notice that we need to know the levels of night that are created in the input.history from
# the make.time.factor() function.
model.list.csv <- textConnection("
PsiA, PsiBA, PsiBa, pA, pB, rA, rBA, rBa
~1, ~1, ~1, ~nightn, ~1, ~1, ~1, ~SHARE
~1, ~1, ~1, ~1, ~1, ~1, ~1, ~1
~1, ~1, ~1, ~1, ~1, ~1, ~1, ~SHARE
~1, ~1, ~SHARE, ~1, ~1, ~1, ~1, ~1
~1, ~1, ~SHARE, ~1, ~1, ~1, ~1, ~SHARE")
model.list <- read.csv(model.list.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
model.list$model.number <- 1:nrow(model.list)
model.list
## PsiA PsiBA PsiBa pA pB rA rBA rBa model.number
## 1 ~1 ~1 ~1 ~nightn ~1 ~1 ~1 ~SHARE 1
## 2 ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1 2
## 3 ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~SHARE 3
## 4 ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~1 4
## 5 ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~SHARE 5
# fit the models
myobs <- ls()
myobs <- myobs[ grepl("m...",myobs,fixed=TRUE)]
cat("Removing ", myobs, "\n")
## Removing
rm(list=myobs)
model.fits <- plyr::dlply(model.list, "model.number", function(x,input.data, input.ddl){
cat("\n\n***** Starting ", unlist(x), "\n")
#browser()
model.parameters=list(
PsiA =list(formula=as.formula(eval(x$PsiA ))),
PsiBA =list(formula=as.formula(eval(x$PsiBA))),
PsiBa =list(formula=as.formula(eval(x$PsiBa))),
pA =list(formula=as.formula(eval(x$pA ))),
pB =list(formula=as.formula(eval(x$pB ))),
rA =list(formula=as.formula(eval(x$rA ))),
rBA =list(formula=as.formula(eval(x$rBA))),
rBa =list(formula=as.formula(eval(x$rBa)))
)
if(x$rBa == '~SHARE'){
model.parameters$rBA =list(formula=as.formula(eval(x$rBA)), share=TRUE)
model.parameters$rBa = NULL
}
if(x$PsiBa == '~SHARE'){
model.parameters$PsiBA =list(formula=as.formula(eval(x$PsiBA)), share=TRUE)
model.parameters$PsiBa = NULL
}
fit <- RMark::mark(input.data, ddl=input.ddl,
model="2SpecConOccup",
model.parameters=model.parameters
#,brief=TRUE,output=FALSE, delete=TRUE
#,invisible=TRUE,output=TRUE # set for debugging
)
mnumber <- paste("m...",formatC(x$model.number, width = 3, format = "d", flag = "0"),sep="")
assign( mnumber, fit, envir=.GlobalEnv)
#browser()
fit
},input.data=owl.data, input.ddl=owl.ddl)
##
##
## ***** Starting ~1 ~1 ~1 ~nightn ~1 ~1 ~1 ~SHARE 1
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~nightn)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 8
## -2lnL: 1206.005
## AICc : 1223.019
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.8991130 0.2099278 0.4876545 1.3105714
## PsiBA:(Intercept) -0.9209423 0.2991821 -1.5073393 -0.3345454
## PsiBa:(Intercept) -0.1194094 0.4362809 -0.9745199 0.7357011
## pA:(Intercept) 0.2857420 0.1370603 0.0171037 0.5543803
## pA:nightn 4.2607453 1.1932627 1.9219505 6.5995402
## pB:(Intercept) -0.5442612 0.3259371 -1.1830979 0.0945755
## rA:(Intercept) -0.4342847 0.2299352 -0.8849577 0.0163884
## rBA:(Intercept) -1.3087226 0.2573974 -1.8132215 -0.8042237
##
##
## Real Parameter PsiA
## 1
## 0.7107672
##
##
## Real Parameter PsiBA
## 1
## 0.2847659
##
##
## Real Parameter PsiBa
## 1
## 0.4701831
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.6051156 0.6512059 0.7005837 0.670179 0.670179 0.6185174 0.6185174
## 8 9 10
## 0.5915538 0.5847189 0.5847189
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.3671969 0.3671969 0.3671969 0.3671969 0.3671969 0.3671969 0.3671969
## 8 9 10
## 0.3671969 0.3671969 0.3671969
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.3931037 0.3931037 0.3931037 0.3931037 0.3931037 0.3931037 0.3931037
## 8 9 10
## 0.3931037 0.3931037 0.3931037
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.2127007 0.2127007 0.2127007 0.2127007 0.2127007 0.2127007 0.2127007
## 8 9 10
## 0.2127007 0.2127007 0.2127007
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.2127007 0.2127007 0.2127007 0.2127007 0.2127007 0.2127007 0.2127007
## 8 9 10
## 0.2127007 0.2127007 0.2127007
##
##
## ***** Starting ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1 2
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
##
## Npar : 8
## -2lnL: 1257.98
## AICc : 1274.995
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.9777190 0.2305016 0.5259359 1.4295021
## PsiBA:(Intercept) 0.3335506 0.2979541 -0.2504395 0.9175407
## PsiBa:(Intercept) -0.1738988 0.4496160 -1.0551461 0.7073486
## pA:(Intercept) 1.5714494 0.2874767 1.0079951 2.1349037
## pB:(Intercept) -0.4505972 0.3392304 -1.1154888 0.2142944
## rA:(Intercept) -0.4444510 0.1799102 -0.7970750 -0.0918271
## rBA:(Intercept) -2.2980619 0.3464465 -2.9770969 -1.6190268
## rBa:(Intercept) -2.1102367 0.2595061 -2.6188686 -1.6016048
##
##
## Real Parameter PsiA
## 1
## 0.7266554
##
##
## Real Parameter PsiBA
## 1
## 0.582623
##
##
## Real Parameter PsiBa
## 1
## 0.4566345
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.8279901 0.8279901 0.8279901 0.8279901 0.8279901 0.8279901 0.8279901
## 8 9 10
## 0.8279901 0.8279901 0.8279901
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.3892188 0.3892188 0.3892188 0.3892188 0.3892188 0.3892188 0.3892188
## 8 9 10
## 0.3892188 0.3892188 0.3892188
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.3906809 0.3906809 0.3906809 0.3906809 0.3906809 0.3906809 0.3906809
## 8 9 10
## 0.3906809 0.3906809 0.3906809
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.0912836 0.0912836 0.0912836 0.0912836 0.0912836 0.0912836 0.0912836
## 8 9 10
## 0.0912836 0.0912836 0.0912836
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.1081058 0.1081058 0.1081058 0.1081058 0.1081058 0.1081058 0.1081058
## 8 9 10
## 0.1081058 0.1081058 0.1081058
##
##
## ***** Starting ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~SHARE 3
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 7
## -2lnL: 1258.201
## AICc : 1272.985
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.9806891 0.2311542 0.5276269 1.4337513
## PsiBA:(Intercept) 0.3120326 0.2922453 -0.2607682 0.8848334
## PsiBa:(Intercept) -0.1524833 0.4475507 -1.0296827 0.7247160
## pA:(Intercept) 1.5552715 0.2823452 1.0018749 2.1086680
## pB:(Intercept) -0.4562723 0.3367154 -1.1162345 0.2036899
## rA:(Intercept) -0.4570133 0.1773573 -0.8046337 -0.1093929
## rBA:(Intercept) -2.1762873 0.2238375 -2.6150087 -1.7375658
##
##
## Real Parameter PsiA
## 1
## 0.7272449
##
##
## Real Parameter PsiBA
## 1
## 0.5773813
##
##
## Real Parameter PsiBa
## 1
## 0.4619529
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.8256738 0.8256738 0.8256738 0.8256738 0.8256738 0.8256738 0.8256738
## 8 9 10
## 0.8256738 0.8256738 0.8256738
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.3878705 0.3878705 0.3878705 0.3878705 0.3878705 0.3878705 0.3878705
## 8 9 10
## 0.3878705 0.3878705 0.3878705
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.3876946 0.3876946 0.3876946 0.3876946 0.3876946 0.3876946 0.3876946
## 8 9 10
## 0.3876946 0.3876946 0.3876946
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.1019002 0.1019002 0.1019002 0.1019002 0.1019002 0.1019002 0.1019002
## 8 9 10
## 0.1019002 0.1019002 0.1019002
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.1019002 0.1019002 0.1019002 0.1019002 0.1019002 0.1019002 0.1019002
## 8 9 10
## 0.1019002 0.1019002 0.1019002
##
##
## ***** Starting ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~1 4
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1)
##
## Npar : 7
## -2lnL: 1258.847
## AICc : 1273.63
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.9629361 0.2262398 0.5195060 1.4063662
## PsiBA:(Intercept) 0.1941007 0.2662335 -0.3277169 0.7159183
## pA:(Intercept) 1.5054945 0.2800929 0.9565124 2.0544766
## pB:(Intercept) -0.5941336 0.3126364 -1.2069010 0.0186337
## rA:(Intercept) -0.4620421 0.1739669 -0.8030172 -0.1210670
## rBA:(Intercept) -2.2471509 0.3447311 -2.9228240 -1.5714779
## rBa:(Intercept) -2.1183133 0.2655486 -2.6387886 -1.5978381
##
##
## Real Parameter PsiA
## 1
## 0.7237093
##
##
## Real Parameter PsiBA
## 1
## 0.5483734
##
##
## Real Parameter PsiBa
## 1
## 0.5483734
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.8183925 0.8183925 0.8183925 0.8183925 0.8183925 0.8183925 0.8183925
## 8 9 10
## 0.8183925 0.8183925 0.8183925
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7 8
## 0.355687 0.355687 0.355687 0.355687 0.355687 0.355687 0.355687 0.355687
## 9 10
## 0.355687 0.355687
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7
## 0.3865015 0.3865015 0.3865015 0.3865015 0.3865015 0.3865015 0.3865015
## 8 9 10
## 0.3865015 0.3865015 0.3865015
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.0955955 0.0955955 0.0955955 0.0955955 0.0955955 0.0955955 0.0955955
## 8 9 10
## 0.0955955 0.0955955 0.0955955
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.1073296 0.1073296 0.1073296 0.1073296 0.1073296 0.1073296 0.1073296
## 8 9 10
## 0.1073296 0.1073296 0.1073296
##
##
## ***** Starting ~1 ~1 ~SHARE ~1 ~1 ~1 ~1 ~SHARE 5
##
## Output summary for 2SpecConOccup model
## Name : PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa()
##
## Npar : 6
## -2lnL: 1258.952
## AICc : 1271.535
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.9667101 0.2270528 0.5216866 1.4117336
## PsiBA:(Intercept) 0.1879055 0.2644258 -0.3303691 0.7061801
## pA:(Intercept) 1.4983802 0.2774537 0.9545709 2.0421895
## pB:(Intercept) -0.5882535 0.3107718 -1.1973662 0.0208592
## rA:(Intercept) -0.4700729 0.1723474 -0.8078737 -0.1322720
## rBA:(Intercept) -2.1645028 0.2275984 -2.6105956 -1.7184099
##
##
## Real Parameter PsiA
## 1
## 0.7244633
##
##
## Real Parameter PsiBA
## 1
## 0.5468386
##
##
## Real Parameter PsiBa
## 1
## 0.5468386
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.8173328 0.8173328 0.8173328 0.8173328 0.8173328 0.8173328 0.8173328
## 8 9 10
## 0.8173328 0.8173328 0.8173328
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.3570357 0.3570357 0.3570357 0.3570357 0.3570357 0.3570357 0.3570357
## 8 9 10
## 0.3570357 0.3570357 0.3570357
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7 8
## 0.384599 0.384599 0.384599 0.384599 0.384599 0.384599 0.384599 0.384599
## 9 10
## 0.384599 0.384599
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837
## 8 9 10
## 0.1029837 0.1029837 0.1029837
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837
## 8 9 10
## 0.1029837 0.1029837 0.1029837
# examine individual model results
model.number <-5
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: 1258.952
## AICc : 1271.535
##
## Beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.9667101 0.2270528 0.5216866 1.4117336
## PsiBA:(Intercept) 0.1879055 0.2644258 -0.3303691 0.7061801
## pA:(Intercept) 1.4983802 0.2774537 0.9545709 2.0421895
## pB:(Intercept) -0.5882535 0.3107718 -1.1973662 0.0208592
## rA:(Intercept) -0.4700729 0.1723474 -0.8078737 -0.1322720
## rBA:(Intercept) -2.1645028 0.2275984 -2.6105956 -1.7184099
##
##
## Real Parameter PsiA
## 1
## 0.7244633
##
##
## Real Parameter PsiBA
## 1
## 0.5468386
##
##
## Real Parameter PsiBa
## 1
## 0.5468386
##
##
## Real Parameter pA
## 1 2 3 4 5 6 7
## 0.8173328 0.8173328 0.8173328 0.8173328 0.8173328 0.8173328 0.8173328
## 8 9 10
## 0.8173328 0.8173328 0.8173328
##
##
## Real Parameter pB
## 1 2 3 4 5 6 7
## 0.3570357 0.3570357 0.3570357 0.3570357 0.3570357 0.3570357 0.3570357
## 8 9 10
## 0.3570357 0.3570357 0.3570357
##
##
## Real Parameter rA
## 1 2 3 4 5 6 7 8
## 0.384599 0.384599 0.384599 0.384599 0.384599 0.384599 0.384599 0.384599
## 9 10
## 0.384599 0.384599
##
##
## Real Parameter rBA
## 1 2 3 4 5 6 7
## 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837
## 8 9 10
## 0.1029837 0.1029837 0.1029837
##
##
## Real Parameter rBa
## 1 2 3 4 5 6 7
## 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837 0.1029837
## 8 9 10
## 0.1029837 0.1029837 0.1029837
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## PsiA g1 a0 t1 0.7244633 0.0453234 0.6275421 0.8040392
## PsiBA g1 a0 t1 0.5468386 0.0655263 0.4181508 0.6695566
## pA g1 a0 t1 0.8173328 0.0414238 0.7220335 0.8851560
## pB g1 a0 t1 0.3570357 0.0713411 0.2319441 0.5052146
## rA g1 a0 t1 0.3845990 0.0407916 0.3083438 0.4669801
## rBA g1 a0 t1 0.1029837 0.0210251 0.0684596 0.1520761
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## PsiA:(Intercept) 0.9667101 0.2270528 0.5216866 1.4117336
## PsiBA:(Intercept) 0.1879055 0.2644258 -0.3303691 0.7061801
## pA:(Intercept) 1.4983802 0.2774537 0.9545709 2.0421895
## pB:(Intercept) -0.5882535 0.3107718 -1.1973662 0.0208592
## rA:(Intercept) -0.4700729 0.1723474 -0.8078737 -0.1322720
## rBA:(Intercept) -2.1645028 0.2275984 -2.6105956 -1.7184099
model.fits[[model.number]]$results$derived
## $`Species Interaction Factor`
## estimate se lcl ucl
## 1 1 0 1 1
##
## $`Occupancy of Species B`
## estimate se lcl ucl
## 1 0.5468386 0.06552634 0.4181508 0.6695566
##
## $`Occupancy of Both Species`
## estimate se lcl ucl
## 1 0.3961645 0.05535094 0.2942187 0.5080088
get.real(model.fits[[model.number]], "PsiA", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiA g1 a0 t1 1 1 0.7244633 0.0453234 0.6275421
## ucl fixed note group age time Age Time
## PsiA g1 a0 t1 0.8040392 1 0 1 0 0
get.real(model.fits[[model.number]], "PsiBA", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiBA g1 a0 t1 2 2 0.5468386 0.0655263 0.4181508
## ucl fixed note group age time Age Time
## PsiBA g1 a0 t1 0.6695566 1 0 1 0 0
get.real(model.fits[[model.number]], "PsiBa", se=TRUE)
## all.diff.index par.index estimate se lcl
## PsiBa g1 a0 t1 3 2 0.5468386 0.0655263 0.4181508
## ucl fixed note group age time Age Time
## PsiBa g1 a0 t1 0.6695566 1 0 1 0 0
get.real(model.fits[[model.number]], "pA", se=TRUE)
## all.diff.index par.index estimate se lcl
## pA g1 a0 t1 4 3 0.8173328 0.0414238 0.7220335
## pA g1 a1 t2 5 3 0.8173328 0.0414238 0.7220335
## pA g1 a2 t3 6 3 0.8173328 0.0414238 0.7220335
## pA g1 a3 t4 7 3 0.8173328 0.0414238 0.7220335
## pA g1 a4 t5 8 3 0.8173328 0.0414238 0.7220335
## pA g1 a5 t6 9 3 0.8173328 0.0414238 0.7220335
## pA g1 a6 t7 10 3 0.8173328 0.0414238 0.7220335
## pA g1 a7 t8 11 3 0.8173328 0.0414238 0.7220335
## pA g1 a8 t9 12 3 0.8173328 0.0414238 0.7220335
## pA g1 a9 t10 13 3 0.8173328 0.0414238 0.7220335
## ucl fixed note group age time Age Time
## pA g1 a0 t1 0.885156 1 0 1 0 0
## pA g1 a1 t2 0.885156 1 1 2 1 1
## pA g1 a2 t3 0.885156 1 2 3 2 2
## pA g1 a3 t4 0.885156 1 3 4 3 3
## pA g1 a4 t5 0.885156 1 4 5 4 4
## pA g1 a5 t6 0.885156 1 5 6 5 5
## pA g1 a6 t7 0.885156 1 6 7 6 6
## pA g1 a7 t8 0.885156 1 7 8 7 7
## pA g1 a8 t9 0.885156 1 8 9 8 8
## pA g1 a9 t10 0.885156 1 9 10 9 9
get.real(model.fits[[model.number]], "pB", se=TRUE)
## all.diff.index par.index estimate se lcl
## pB g1 a0 t1 14 4 0.3570357 0.0713411 0.2319441
## pB g1 a1 t2 15 4 0.3570357 0.0713411 0.2319441
## pB g1 a2 t3 16 4 0.3570357 0.0713411 0.2319441
## pB g1 a3 t4 17 4 0.3570357 0.0713411 0.2319441
## pB g1 a4 t5 18 4 0.3570357 0.0713411 0.2319441
## pB g1 a5 t6 19 4 0.3570357 0.0713411 0.2319441
## pB g1 a6 t7 20 4 0.3570357 0.0713411 0.2319441
## pB g1 a7 t8 21 4 0.3570357 0.0713411 0.2319441
## pB g1 a8 t9 22 4 0.3570357 0.0713411 0.2319441
## pB g1 a9 t10 23 4 0.3570357 0.0713411 0.2319441
## ucl fixed note group age time Age Time
## pB g1 a0 t1 0.5052146 1 0 1 0 0
## pB g1 a1 t2 0.5052146 1 1 2 1 1
## pB g1 a2 t3 0.5052146 1 2 3 2 2
## pB g1 a3 t4 0.5052146 1 3 4 3 3
## pB g1 a4 t5 0.5052146 1 4 5 4 4
## pB g1 a5 t6 0.5052146 1 5 6 5 5
## pB g1 a6 t7 0.5052146 1 6 7 6 6
## pB g1 a7 t8 0.5052146 1 7 8 7 7
## pB g1 a8 t9 0.5052146 1 8 9 8 8
## pB g1 a9 t10 0.5052146 1 9 10 9 9
get.real(model.fits[[model.number]], "rA", se=TRUE)
## all.diff.index par.index estimate se lcl
## rA g1 a0 t1 24 5 0.384599 0.0407916 0.3083438
## rA g1 a1 t2 25 5 0.384599 0.0407916 0.3083438
## rA g1 a2 t3 26 5 0.384599 0.0407916 0.3083438
## rA g1 a3 t4 27 5 0.384599 0.0407916 0.3083438
## rA g1 a4 t5 28 5 0.384599 0.0407916 0.3083438
## rA g1 a5 t6 29 5 0.384599 0.0407916 0.3083438
## rA g1 a6 t7 30 5 0.384599 0.0407916 0.3083438
## rA g1 a7 t8 31 5 0.384599 0.0407916 0.3083438
## rA g1 a8 t9 32 5 0.384599 0.0407916 0.3083438
## rA g1 a9 t10 33 5 0.384599 0.0407916 0.3083438
## ucl fixed note group age time Age Time
## rA g1 a0 t1 0.4669801 1 0 1 0 0
## rA g1 a1 t2 0.4669801 1 1 2 1 1
## rA g1 a2 t3 0.4669801 1 2 3 2 2
## rA g1 a3 t4 0.4669801 1 3 4 3 3
## rA g1 a4 t5 0.4669801 1 4 5 4 4
## rA g1 a5 t6 0.4669801 1 5 6 5 5
## rA g1 a6 t7 0.4669801 1 6 7 6 6
## rA g1 a7 t8 0.4669801 1 7 8 7 7
## rA g1 a8 t9 0.4669801 1 8 9 8 8
## rA g1 a9 t10 0.4669801 1 9 10 9 9
get.real(model.fits[[model.number]], "rBA", se=TRUE)
## all.diff.index par.index estimate se lcl
## rBA g1 a0 t1 34 6 0.1029837 0.0210251 0.0684596
## rBA g1 a1 t2 35 6 0.1029837 0.0210251 0.0684596
## rBA g1 a2 t3 36 6 0.1029837 0.0210251 0.0684596
## rBA g1 a3 t4 37 6 0.1029837 0.0210251 0.0684596
## rBA g1 a4 t5 38 6 0.1029837 0.0210251 0.0684596
## rBA g1 a5 t6 39 6 0.1029837 0.0210251 0.0684596
## rBA g1 a6 t7 40 6 0.1029837 0.0210251 0.0684596
## rBA g1 a7 t8 41 6 0.1029837 0.0210251 0.0684596
## rBA g1 a8 t9 42 6 0.1029837 0.0210251 0.0684596
## rBA g1 a9 t10 43 6 0.1029837 0.0210251 0.0684596
## ucl fixed note group age time Age Time
## rBA g1 a0 t1 0.1520761 1 0 1 0 0
## rBA g1 a1 t2 0.1520761 1 1 2 1 1
## rBA g1 a2 t3 0.1520761 1 2 3 2 2
## rBA g1 a3 t4 0.1520761 1 3 4 3 3
## rBA g1 a4 t5 0.1520761 1 4 5 4 4
## rBA g1 a5 t6 0.1520761 1 5 6 5 5
## rBA g1 a6 t7 0.1520761 1 6 7 6 6
## rBA g1 a7 t8 0.1520761 1 7 8 7 7
## rBA g1 a8 t9 0.1520761 1 8 9 8 8
## rBA g1 a9 t10 0.1520761 1 9 10 9 9
get.real(model.fits[[model.number]], "rBa", se=TRUE)
## all.diff.index par.index estimate se lcl
## rBa g1 a0 t1 44 6 0.1029837 0.0210251 0.0684596
## rBa g1 a1 t2 45 6 0.1029837 0.0210251 0.0684596
## rBa g1 a2 t3 46 6 0.1029837 0.0210251 0.0684596
## rBa g1 a3 t4 47 6 0.1029837 0.0210251 0.0684596
## rBa g1 a4 t5 48 6 0.1029837 0.0210251 0.0684596
## rBa g1 a5 t6 49 6 0.1029837 0.0210251 0.0684596
## rBa g1 a6 t7 50 6 0.1029837 0.0210251 0.0684596
## rBa g1 a7 t8 51 6 0.1029837 0.0210251 0.0684596
## rBa g1 a8 t9 52 6 0.1029837 0.0210251 0.0684596
## rBa g1 a9 t10 53 6 0.1029837 0.0210251 0.0684596
## ucl fixed note group age time Age Time
## rBa g1 a0 t1 0.1520761 1 0 1 0 0
## rBa g1 a1 t2 0.1520761 1 1 2 1 1
## rBa g1 a2 t3 0.1520761 1 2 3 2 2
## rBa g1 a3 t4 0.1520761 1 3 4 3 3
## rBa g1 a4 t5 0.1520761 1 4 5 4 4
## rBa g1 a5 t6 0.1520761 1 5 6 5 5
## rBa g1 a6 t7 0.1520761 1 6 7 6 6
## rBa g1 a7 t8 0.1520761 1 7 8 7 7
## rBa g1 a8 t9 0.1520761 1 8 9 8 8
## rBa g1 a9 t10 0.1520761 1 9 10 9 9
# collect models and make AICc table
model.set <- RMark::collect.models( type="2SpecConOccup")
model.set
## model npar
## 1 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~nightn)pB(~1)rA(~1)rBA(~1)rBa() 8
## 5 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 6
## 3 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 7
## 4 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 7
## 2 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 8
## AICc DeltaAICc weight Deviance
## 1 1223.019 0.00000 1.000000e+00 1206.005
## 5 1271.535 48.51625 2.916288e-11 1258.952
## 3 1272.985 49.96583 1.412724e-11 1258.201
## 4 1273.630 50.61113 1.023133e-11 1258.847
## 2 1274.994 51.97560 5.171802e-12 1257.980
names(model.set)
## [1] "m...001" "m...002" "m...003" "m...004" "m...005"
## [6] "model.table"
model.set$model.table
## PsiA PsiBA PsiBa pA pB rA rBA rBa
## 1 ~1 ~1 ~1 ~nightn ~1 ~1 ~1
## 5 ~1 ~1 ~1 ~1 ~1 ~1
## 3 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## 4 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## 2 ~1 ~1 ~1 ~1 ~1 ~1 ~1 ~1
## model npar
## 1 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~nightn)pB(~1)rA(~1)rBA(~1)rBa() 8
## 5 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 6
## 3 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa() 7
## 4 PsiA(~1)PsiBA(~1)PsiBa()pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 7
## 2 PsiA(~1)PsiBA(~1)PsiBa(~1)pA(~1)pB(~1)rA(~1)rBA(~1)rBa(~1) 8
## AICc DeltaAICc weight Deviance
## 1 1223.019 0.00000 1.000000e+00 1206.005
## 5 1271.535 48.51625 2.916288e-11 1258.952
## 3 1272.985 49.96583 1.412724e-11 1258.201
## 4 1273.630 50.61113 1.023133e-11 1258.847
## 2 1274.994 51.97560 5.171802e-12 1257.980
# No need to do model averaging, as 100% of weight is on one model
# cleanup
cleanup(ask=FALSE)