# Single Species Single Season Occupancy models using RPresence
# Blue Gross Beaks.
#Downloaded from https://sites.google.com/site/asrworkshop/home/schedule/r-occupancy-1
#An occupancy study was made on Blue Grosbeaks (Guiraca caerulea)
# on 41 old fields planted to longleaf pines (Pinus palustris)
# in southern Georgia, USA.
# Surveys were 500 m transects across each field
# and were completed three times during the breeding season in 2001.
# Columns in the file are:
# field - field number
# v1, v2, v3 - detection histories for each site on each of 3 visit during the 2001 breeding season.
# field.size - size of the files
# bqi - unknown
# crop.hist - crop history
# crop1, crop2 - indicator variables for the crop history
# count1, count2, count3 - are actual counts of birds detected in each visit
# RPresence package
library(readxl)
library(RPresence)
library(ggplot2)
# 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("..","blgr.csv"),
header=TRUE, as.is=TRUE, strip.white=TRUE)
head(input.data)
## field v1 v2 v3 field.size bqi crop.hist crop1 crop2 count1 count2 count3
## 1 1 1 1 1 14.0 1 crop 1 0 1 2 2
## 2 2 1 1 0 12.7 1 crop 1 0 2 2 0
## 3 3 0 0 0 15.7 0 grass 0 1 0 0 0
## 4 4 0 1 0 19.5 0 grass 0 1 0 2 0
## 5 5 1 0 1 13.5 0 crop 1 0 1 0 1
## 6 6 0 0 1 9.6 0 mixed 0 1 0 0 2
## X logFS
## 1 NA 1.1461280
## 2 NA 1.1038037
## 3 NA 1.1958997
## 4 NA 1.2900346
## 5 NA 1.1303338
## 6 NA 0.9822712
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.data)
## [1] 41
range(input.data[, c("v1","v2","v3")], na.rm=TRUE)
## [1] 0 1
sum(is.na(input.data[, c("v1","v2","v3")]))
## [1] 0
input.history <- input.data[, c("v1","v2","v3")]
head(input.history)
## v1 v2 v3
## 1 1 1 1
## 2 1 1 0
## 3 0 0 0
## 4 0 1 0
## 5 1 0 1
## 6 0 0 1
site.covar <- input.data[, c("field","field.size")]
site.covar$logFS <- log(site.covar$field.size)
head(site.covar)
## field field.size logFS
## 1 1 14.0 2.639057
## 2 2 12.7 2.541602
## 3 3 15.7 2.753661
## 4 4 19.5 2.970414
## 5 5 13.5 2.602690
## 6 6 9.6 2.261763
# Create the *.pao file
grossbeak.pao <- RPresence::createPao(input.history,
unitcov=site.covar,
title='Grossbeak SSSS')
grossbeak.pao
## $nunits
## [1] 41
##
## $nsurveys
## [1] 3
##
## $nseasons
## [1] 1
##
## $nmethods
## [1] 1
##
## $det.data
## v1 v2 v3
## 1 1 1 1
## 2 1 1 0
## 3 0 0 0
## 4 0 1 0
## 5 1 0 1
## 6 0 0 1
## 7 0 0 1
## 8 1 1 1
## 9 1 1 0
## 10 1 1 1
## 11 1 1 0
## 12 0 0 0
## 13 0 0 0
## 14 0 0 1
## 15 1 1 1
## 16 0 0 1
## 17 0 0 1
## 18 0 0 0
## 19 0 1 1
## 20 0 0 0
## 21 1 0 0
## 22 0 1 0
## 23 1 0 0
## 24 1 1 1
## 25 1 1 1
## 26 0 1 1
## 27 0 0 1
## 28 0 1 0
## 29 1 1 0
## 30 0 1 1
## 31 0 0 0
## 32 1 1 1
## 33 1 0 0
## 34 1 1 0
## 35 0 0 0
## 36 0 0 0
## 37 0 1 0
## 38 0 1 1
## 39 1 1 1
## 40 1 0 1
## 41 0 1 0
##
## $nunitcov
## [1] 3
##
## $unitcov
## field field.size logFS
## 1 1 14.0 2.639057
## 2 2 12.7 2.541602
## 3 3 15.7 2.753661
## 4 4 19.5 2.970414
## 5 5 13.5 2.602690
## 6 6 9.6 2.261763
## 7 7 44.0 3.784190
## 8 8 9.4 2.240710
## 9 9 19.6 2.975530
## 10 10 7.0 1.945910
## 11 11 20.9 3.039749
## 12 12 17.8 2.879198
## 13 13 14.4 2.667228
## 14 14 24.9 3.214868
## 15 15 6.4 1.856298
## 16 16 53.9 3.987130
## 17 17 10.3 2.332144
## 18 18 12.0 2.484907
## 19 19 6.8 1.916923
## 20 20 17.6 2.867899
## 21 21 18.8 2.933857
## 22 22 28.0 3.332205
## 23 23 7.9 2.066863
## 24 24 38.3 3.645450
## 25 25 15.9 2.766319
## 26 26 11.3 2.424803
## 27 27 10.5 2.351375
## 28 28 10.3 2.332144
## 29 29 18.2 2.901422
## 30 30 23.0 3.135494
## 31 31 30.2 3.407842
## 32 32 6.7 1.902108
## 33 33 10.0 2.302585
## 34 34 10.0 2.302585
## 35 35 18.1 2.895912
## 36 36 6.8 1.916923
## 37 37 9.2 2.219203
## 38 38 10.7 2.370244
## 39 39 9.0 2.197225
## 40 40 9.1 2.208274
## 41 41 10.4 2.341806
##
## $nsurvcov
## [1] 1
##
## $survcov
## SURVEY
## 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 2
## 43 2
## 44 2
## 45 2
## 46 2
## 47 2
## 48 2
## 49 2
## 50 2
## 51 2
## 52 2
## 53 2
## 54 2
## 55 2
## 56 2
## 57 2
## 58 2
## 59 2
## 60 2
## 61 2
## 62 2
## 63 2
## 64 2
## 65 2
## 66 2
## 67 2
## 68 2
## 69 2
## 70 2
## 71 2
## 72 2
## 73 2
## 74 2
## 75 2
## 76 2
## 77 2
## 78 2
## 79 2
## 80 2
## 81 2
## 82 2
## 83 3
## 84 3
## 85 3
## 86 3
## 87 3
## 88 3
## 89 3
## 90 3
## 91 3
## 92 3
## 93 3
## 94 3
## 95 3
## 96 3
## 97 3
## 98 3
## 99 3
## 100 3
## 101 3
## 102 3
## 103 3
## 104 3
## 105 3
## 106 3
## 107 3
## 108 3
## 109 3
## 110 3
## 111 3
## 112 3
## 113 3
## 114 3
## 115 3
## 116 3
## 117 3
## 118 3
## 119 3
## 120 3
## 121 3
## 122 3
## 123 3
##
## $nsurveyseason
## [1] 3
##
## $title
## [1] "Grossbeak SSSS"
##
## $unitnames
## [1] "unit1" "unit2" "unit3" "unit4" "unit5" "unit6" "unit7"
## [8] "unit8" "unit9" "unit10" "unit11" "unit12" "unit13" "unit14"
## [15] "unit15" "unit16" "unit17" "unit18" "unit19" "unit20" "unit21"
## [22] "unit22" "unit23" "unit24" "unit25" "unit26" "unit27" "unit28"
## [29] "unit29" "unit30" "unit31" "unit32" "unit33" "unit34" "unit35"
## [36] "unit36" "unit37" "unit38" "unit39" "unit40" "unit41"
##
## $surveynames
## [1] "1-1" "1-2" "1-3"
##
## $paoname
## [1] "pres.pao"
##
## $frq
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1
##
## attr(,"class")
## [1] "pao"
# Fit some models.
# Note that formula DO NOT HAVE AN = SIGN
mod.pdot <- RPresence::occMod(model=list(psi~1, p~1),
type="so", data=grossbeak.pao)
## PRESENCE Version 2.12.18.
summary(mod.pdot)
## Model name=psi()p()
## AIC=172.1898
## -2*log-likelihood=168.1898
## num. par=2
names(mod.pdot)
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "npar" "aic" "beta" "real"
## [11] "derived" "gof" "warnings" "version"
mod.pdot$beta$psi
## est se
## A1_psi 2.038689 0.75139
# look at estimated occupancy probability. RPresence gives for EACH site in case it depends on covariates
head(mod.pdot$real$psi)
## est se lower_0.95 upper_0.95
## psi_unit1 0.8847997 0.07658864 0.6378374 0.9710101
## psi_unit2 0.8847997 0.07658864 0.6378374 0.9710101
## psi_unit3 0.8847997 0.07658864 0.6378374 0.9710101
## psi_unit4 0.8847997 0.07658864 0.6378374 0.9710101
## psi_unit5 0.8847997 0.07658864 0.6378374 0.9710101
## psi_unit6 0.8847997 0.07658864 0.6378374 0.9710101
mod.pdot.psi <-mod.pdot$real$psi[1,] # occupancy probability
mod.pdot.psi
## est se lower_0.95 upper_0.95
## psi_unit1 0.8847997 0.07658864 0.6378374 0.9710101
# look at the estimated probability of detection. It gives an estimate for every site at very visit
head(mod.pdot$real$p)
## est se lower_0.95 upper_0.95
## p1_unit1 0.5513167 0.05987549 0.4332949 0.663829
## p1_unit2 0.5513167 0.05987549 0.4332949 0.663829
## p1_unit3 0.5513167 0.05987549 0.4332949 0.663829
## p1_unit4 0.5513167 0.05987549 0.4332949 0.663829
## p1_unit5 0.5513167 0.05987549 0.4332949 0.663829
## p1_unit6 0.5513167 0.05987549 0.4332949 0.663829
mod.pdot.p <- mod.pdot$real$p[seq(1, by=nrow(input.history), length.out=ncol(input.history)),]
mod.pdot.p
## est se lower_0.95 upper_0.95
## p1_unit1 0.5513167 0.05987549 0.4332949 0.663829
## p2_unit1 0.5513167 0.05987549 0.4332949 0.663829
## p3_unit1 0.5513167 0.05987549 0.4332949 0.663829
# Look at the posterior probability of detection
names(mod.pdot$derived)
## [1] "psi_c"
mod.pdot$derived$psi_c
## est se lower_0.95 upper_0.95
## unit1 1.0000000 0.0000000 1.00000000 1.0000000
## unit2 1.0000000 0.0000000 1.00000000 1.0000000
## unit3 0.4095987 0.2419634 0.08893652 0.8313801
## unit4 1.0000000 0.0000000 1.00000000 1.0000000
## unit5 1.0000000 0.0000000 1.00000000 1.0000000
## unit6 1.0000000 0.0000000 1.00000000 1.0000000
## unit7 1.0000000 0.0000000 1.00000000 1.0000000
## unit8 1.0000000 0.0000000 1.00000000 1.0000000
## unit9 1.0000000 0.0000000 1.00000000 1.0000000
## unit10 1.0000000 0.0000000 1.00000000 1.0000000
## unit11 1.0000000 0.0000000 1.00000000 1.0000000
## unit12 0.4095987 0.2419634 0.08893652 0.8313801
## unit13 0.4095987 0.2419634 0.08893652 0.8313801
## unit14 1.0000000 0.0000000 1.00000000 1.0000000
## unit15 1.0000000 0.0000000 1.00000000 1.0000000
## unit16 1.0000000 0.0000000 1.00000000 1.0000000
## unit17 1.0000000 0.0000000 1.00000000 1.0000000
## unit18 0.4095987 0.2419634 0.08893652 0.8313801
## unit19 1.0000000 0.0000000 1.00000000 1.0000000
## unit20 0.4095987 0.2419634 0.08893652 0.8313801
## unit21 1.0000000 0.0000000 1.00000000 1.0000000
## unit22 1.0000000 0.0000000 1.00000000 1.0000000
## unit23 1.0000000 0.0000000 1.00000000 1.0000000
## unit24 1.0000000 0.0000000 1.00000000 1.0000000
## unit25 1.0000000 0.0000000 1.00000000 1.0000000
## unit26 1.0000000 0.0000000 1.00000000 1.0000000
## unit27 1.0000000 0.0000000 1.00000000 1.0000000
## unit28 1.0000000 0.0000000 1.00000000 1.0000000
## unit29 1.0000000 0.0000000 1.00000000 1.0000000
## unit30 1.0000000 0.0000000 1.00000000 1.0000000
## unit31 0.4095987 0.2419634 0.08893652 0.8313801
## unit32 1.0000000 0.0000000 1.00000000 1.0000000
## unit33 1.0000000 0.0000000 1.00000000 1.0000000
## unit34 1.0000000 0.0000000 1.00000000 1.0000000
## unit35 0.4095987 0.2419634 0.08893652 0.8313801
## unit36 0.4095987 0.2419634 0.08893652 0.8313801
## unit37 1.0000000 0.0000000 1.00000000 1.0000000
## unit38 1.0000000 0.0000000 1.00000000 1.0000000
## unit39 1.0000000 0.0000000 1.00000000 1.0000000
## unit40 1.0000000 0.0000000 1.00000000 1.0000000
## unit41 1.0000000 0.0000000 1.00000000 1.0000000
# alternatively
RPresence::print_one_site_estimates(mod.pdot, site = 1)
## psi()p()
## est se lower_0.95 upper_0.95
## psi_psi_unit1 0.8847997 0.07658864 0.6378374 0.9710101
## p_p1_unit1 0.5513167 0.05987549 0.4332949 0.6638290
# Fit some models.
# Note that formula DO NOT HAVE AN = SIGN
mod.psilogFS.pdot <- RPresence::occMod(model=list(psi~logFS, p~1),
type="so", data=grossbeak.pao)
## PRESENCE Version 2.12.18.
summary(mod.psilogFS.pdot)
## Model name=psi(logFS)p()
## AIC=173.9232
## -2*log-likelihood=167.9232
## num. par=3
names(mod.psilogFS.pdot)
## [1] "modname" "model" "dmat" "data" "outfile"
## [6] "neg2loglike" "npar" "aic" "beta" "real"
## [11] "derived" "gof" "warnings" "version"
mod.psilogFS.pdot$beta$psi
## est se
## A1_psi 3.444160 2.931248
## A2_psi.psi.logFS -0.532543 1.020910
# look at estimated occupancy probability. RPresence gives for EACH site in case it depends on covariates
head(mod.psilogFS.pdot$real$psi)
## est se lower_0.95 upper_0.95
## psi_unit1 0.8848058 0.07585133 0.6411044 0.9706119
## psi_unit2 0.8899909 0.07536329 0.6415229 0.9733853
## psi_unit3 0.8784376 0.07800863 0.6331892 0.9680005
## psi_unit4 0.8655632 0.08788248 0.5943654 0.9658594
## psi_unit5 0.8867651 0.07553952 0.6419078 0.9716006
## psi_unit6 0.9037561 0.07845782 0.6158014 0.9821472
mod.psilogFS.pdot.psi <-mod.psilogFS.pdot$real$psi[1,] # occupancy probability
mod.psilogFS.pdot.psi
## est se lower_0.95 upper_0.95
## psi_unit1 0.8848058 0.07585133 0.6411044 0.9706119
# plot of psi vs logfs #1
# individual psi
mod.psilogFS.pdot$real$psi
## est se lower_0.95 upper_0.95
## psi_unit1 0.8848058 0.07585133 0.6411044 0.9706119
## psi_unit2 0.8899909 0.07536329 0.6415229 0.9733853
## psi_unit3 0.8784376 0.07800863 0.6331892 0.9680005
## psi_unit4 0.8655632 0.08788248 0.5943654 0.9658594
## psi_unit5 0.8867651 0.07553952 0.6419078 0.9716006
## psi_unit6 0.9037561 0.07845782 0.6158014 0.9821472
## psi_unit7 0.8067324 0.19612731 0.2618178 0.9800500
## psi_unit8 0.9047269 0.07885235 0.6125177 0.9827724
## psi_unit9 0.8652460 0.08821567 0.5930674 0.9658572
## psi_unit10 0.9174262 0.08501169 0.5519207 0.9901201
## psi_unit11 0.8612084 0.09280304 0.5753293 0.9660096
## psi_unit12 0.8711161 0.08273203 0.6145781 0.9662724
## psi_unit13 0.8832678 0.07620954 0.6399298 0.9698931
## psi_unit14 0.8496814 0.10907555 0.5145690 0.9678889
## psi_unit15 0.9209701 0.08684277 0.5292319 0.9917899
## psi_unit16 0.7893217 0.23808083 0.1846388 0.9841235
## psi_unit17 0.9004463 0.07725138 0.6255461 0.9799884
## psi_unit18 0.8929123 0.07554074 0.6393279 0.9751379
## psi_unit19 0.9185881 0.08561457 0.5447552 0.9906883
## psi_unit20 0.8717902 0.08219707 0.6166906 0.9663735
## psi_unit21 0.8678126 0.08563900 0.6031375 0.9659392
## psi_unit22 0.8415249 0.12294166 0.4657488 0.9700104
## psi_unit23 0.9124135 0.08243521 0.5797373 0.9874479
## psi_unit24 0.8179919 0.17054526 0.3224878 0.9769769
## psi_unit25 0.8777159 0.07836703 0.6317931 0.9677681
## psi_unit26 0.8959346 0.07602844 0.6352044 0.9770470
## psi_unit27 0.8995244 0.07695854 0.6278600 0.9793840
## psi_unit28 0.9004463 0.07725138 0.6255461 0.9799884
## psi_unit29 0.8697815 0.08384979 0.6101705 0.9661057
## psi_unit30 0.8550008 0.10102641 0.5441857 0.9668033
## psi_unit31 0.8360789 0.13307070 0.4319979 0.9715952
## psi_unit32 0.9191762 0.08591912 0.5410266 0.9909682
## psi_unit33 0.9018485 0.07773443 0.6216919 0.9809066
## psi_unit34 0.9018485 0.07773443 0.6216919 0.9809066
## psi_unit35 0.8701135 0.08356452 0.6112945 0.9661431
## psi_unit36 0.9185881 0.08561457 0.5447552 0.9906883
## psi_unit37 0.9057095 0.07926800 0.6089978 0.9833994
## psi_unit38 0.8986126 0.07668950 0.6299769 0.9787868
## psi_unit39 0.9067044 0.07970412 0.6052334 0.9840272
## psi_unit40 0.9062054 0.07948354 0.6071467 0.9837133
## psi_unit41 0.8999841 0.07710201 0.6267281 0.9796854
# covariate values
site.covar
## field field.size logFS
## 1 1 14.0 2.639057
## 2 2 12.7 2.541602
## 3 3 15.7 2.753661
## 4 4 19.5 2.970414
## 5 5 13.5 2.602690
## 6 6 9.6 2.261763
## 7 7 44.0 3.784190
## 8 8 9.4 2.240710
## 9 9 19.6 2.975530
## 10 10 7.0 1.945910
## 11 11 20.9 3.039749
## 12 12 17.8 2.879198
## 13 13 14.4 2.667228
## 14 14 24.9 3.214868
## 15 15 6.4 1.856298
## 16 16 53.9 3.987130
## 17 17 10.3 2.332144
## 18 18 12.0 2.484907
## 19 19 6.8 1.916923
## 20 20 17.6 2.867899
## 21 21 18.8 2.933857
## 22 22 28.0 3.332205
## 23 23 7.9 2.066863
## 24 24 38.3 3.645450
## 25 25 15.9 2.766319
## 26 26 11.3 2.424803
## 27 27 10.5 2.351375
## 28 28 10.3 2.332144
## 29 29 18.2 2.901422
## 30 30 23.0 3.135494
## 31 31 30.2 3.407842
## 32 32 6.7 1.902108
## 33 33 10.0 2.302585
## 34 34 10.0 2.302585
## 35 35 18.1 2.895912
## 36 36 6.8 1.916923
## 37 37 9.2 2.219203
## 38 38 10.7 2.370244
## 39 39 9.0 2.197225
## 40 40 9.1 2.208274
## 41 41 10.4 2.341806
both <- cbind(mod.psilogFS.pdot$real$psi, site.covar)
ggplot(data=both, aes(x=logFS, y=est))+
geom_point()+
geom_ribbon(aes(ymin=lower_0.95 , ymax=upper_0.95 ), alpha=0.2)
# plot #2
# what are beta values
mod.psilogFS.pdot$beta$psi
## est se
## A1_psi 3.444160 2.931248
## A2_psi.psi.logFS -0.532543 1.020910
plotdata <- data.frame(logFS=seq(1,5,.1))
plotdata$logitpsi <-mod.psilogFS.pdot$beta$psi[1,1]+
mod.psilogFS.pdot$beta$psi[2,1]*plotdata$logFS
plotdata$psi <- 1/(1+exp(-plotdata$logitpsi))
ggplot(data=plotdata, aes(x=logFS, y=psi))+
geom_point()
# look at the estimated probability of detection. It gives an estimate for every site at very visit
head(mod.psilogFS.pdot$real$p)
## est se lower_0.95 upper_0.95
## p1_unit1 0.553009 0.05932818 0.4359597 0.6644632
## p1_unit2 0.553009 0.05932818 0.4359597 0.6644632
## p1_unit3 0.553009 0.05932818 0.4359597 0.6644632
## p1_unit4 0.553009 0.05932818 0.4359597 0.6644632
## p1_unit5 0.553009 0.05932818 0.4359597 0.6644632
## p1_unit6 0.553009 0.05932818 0.4359597 0.6644632
mod.psilogFS.pdot.p <- mod.psilogFS.pdot$real$p[seq(1, by=nrow(input.history), length.out=ncol(input.history)),]
mod.psilogFS.pdot.p
## est se lower_0.95 upper_0.95
## p1_unit1 0.553009 0.05932818 0.4359597 0.6644632
## p2_unit1 0.553009 0.05932818 0.4359597 0.6644632
## p3_unit1 0.553009 0.05932818 0.4359597 0.6644632
# Look at the posterior probability of detection
names(mod.psilogFS.pdot$derived)
## [1] "psi_c"
mod.psilogFS.pdot$derived$psi_c
## est se lower_0.95 upper_0.95
## unit1 1.0000000 0.0000000 1.00000000 1.0000000
## unit2 1.0000000 0.0000000 1.00000000 1.0000000
## unit3 0.3922335 0.2323411 0.08720936 0.8134109
## unit4 1.0000000 0.0000000 1.00000000 1.0000000
## unit5 1.0000000 0.0000000 1.00000000 1.0000000
## unit6 1.0000000 0.0000000 1.00000000 1.0000000
## unit7 1.0000000 0.0000000 1.00000000 1.0000000
## unit8 1.0000000 0.0000000 1.00000000 1.0000000
## unit9 1.0000000 0.0000000 1.00000000 1.0000000
## unit10 1.0000000 0.0000000 1.00000000 1.0000000
## unit11 1.0000000 0.0000000 1.00000000 1.0000000
## unit12 0.3764162 0.2294300 0.08162113 0.8039152
## unit13 0.4032589 0.2364252 0.08968146 0.8225502
## unit14 1.0000000 0.0000000 1.00000000 1.0000000
## unit15 1.0000000 0.0000000 1.00000000 1.0000000
## unit16 1.0000000 0.0000000 1.00000000 1.0000000
## unit17 1.0000000 0.0000000 1.00000000 1.0000000
## unit18 0.4268266 0.2505321 0.09096295 0.8471360
## unit19 1.0000000 0.0000000 1.00000000 1.0000000
## unit20 0.3778297 0.2295487 0.08221218 0.8045726
## unit21 1.0000000 0.0000000 1.00000000 1.0000000
## unit22 1.0000000 0.0000000 1.00000000 1.0000000
## unit23 1.0000000 0.0000000 1.00000000 1.0000000
## unit24 1.0000000 0.0000000 1.00000000 1.0000000
## unit25 1.0000000 0.0000000 1.00000000 1.0000000
## unit26 1.0000000 0.0000000 1.00000000 1.0000000
## unit27 1.0000000 0.0000000 1.00000000 1.0000000
## unit28 1.0000000 0.0000000 1.00000000 1.0000000
## unit29 1.0000000 0.0000000 1.00000000 1.0000000
## unit30 1.0000000 0.0000000 1.00000000 1.0000000
## unit31 0.3129609 0.2487853 0.04504100 0.8147944
## unit32 1.0000000 0.0000000 1.00000000 1.0000000
## unit33 1.0000000 0.0000000 1.00000000 1.0000000
## unit34 1.0000000 0.0000000 1.00000000 1.0000000
## unit35 0.3743293 0.2293051 0.08071796 0.8030165
## unit36 0.5019164 0.3309777 0.06996740 0.9310237
## unit37 1.0000000 0.0000000 1.00000000 1.0000000
## unit38 1.0000000 0.0000000 1.00000000 1.0000000
## unit39 1.0000000 0.0000000 1.00000000 1.0000000
## unit40 1.0000000 0.0000000 1.00000000 1.0000000
## unit41 1.0000000 0.0000000 1.00000000 1.0000000
# alternatively
RPresence::print_one_site_estimates(mod.psilogFS.pdot, site = 1)
## psi(logFS)p()
## est se lower_0.95 upper_0.95
## psi_psi_unit1 0.8848058 0.07585133 0.6411044 0.9706119
## p_p1_unit1 0.5530090 0.05932818 0.4359597 0.6644632
# Model where p(t) varies across survey occasions
#
mod.pt <- RPresence::occMod(model=list(psi~1, p~SURVEY), type="so", data=grossbeak.pao)
## PRESENCE Version 2.12.18.
summary(mod.pt)
## Model name=psi()p(SURVEY)
## AIC=175.295
## -2*log-likelihood=167.295
## num. par=4
mod.pt$real$psi[1,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.882712 0.07618293 0.640179 0.9695455
mod.pt$real$p[seq(1, by=nrow(input.history), length.out=ncol(input.history)),]
## est se lower_0.95 upper_0.95
## p1_unit1 0.4973585 0.08915117 0.3297054 0.6656078
## p2_unit1 0.6078827 0.09022755 0.4247047 0.7650074
## p3_unit1 0.5526207 0.09009054 0.3768495 0.7161558
print_one_site_estimates(mod.pt, site = 1)
## psi()p(SURVEY)
## est se lower_0.95 upper_0.95
## psi_psi_unit1 0.8827120 0.07618293 0.6401790 0.9695455
## p_p1_unit1 0.4973585 0.08915117 0.3297054 0.6656078
fitted(mod.pt, param="psi")
## est se lower_0.95 upper_0.95
## psi_unit1 0.882712 0.07618293 0.640179 0.9695455
## psi_unit2 0.882712 0.07618293 0.640179 0.9695455
## psi_unit3 0.882712 0.07618293 0.640179 0.9695455
## psi_unit4 0.882712 0.07618293 0.640179 0.9695455
## psi_unit5 0.882712 0.07618293 0.640179 0.9695455
## psi_unit6 0.882712 0.07618293 0.640179 0.9695455
## psi_unit7 0.882712 0.07618293 0.640179 0.9695455
## psi_unit8 0.882712 0.07618293 0.640179 0.9695455
## psi_unit9 0.882712 0.07618293 0.640179 0.9695455
## psi_unit10 0.882712 0.07618293 0.640179 0.9695455
## psi_unit11 0.882712 0.07618293 0.640179 0.9695455
## psi_unit12 0.882712 0.07618293 0.640179 0.9695455
## psi_unit13 0.882712 0.07618293 0.640179 0.9695455
## psi_unit14 0.882712 0.07618293 0.640179 0.9695455
## psi_unit15 0.882712 0.07618293 0.640179 0.9695455
## psi_unit16 0.882712 0.07618293 0.640179 0.9695455
## psi_unit17 0.882712 0.07618293 0.640179 0.9695455
## psi_unit18 0.882712 0.07618293 0.640179 0.9695455
## psi_unit19 0.882712 0.07618293 0.640179 0.9695455
## psi_unit20 0.882712 0.07618293 0.640179 0.9695455
## psi_unit21 0.882712 0.07618293 0.640179 0.9695455
## psi_unit22 0.882712 0.07618293 0.640179 0.9695455
## psi_unit23 0.882712 0.07618293 0.640179 0.9695455
## psi_unit24 0.882712 0.07618293 0.640179 0.9695455
## psi_unit25 0.882712 0.07618293 0.640179 0.9695455
## psi_unit26 0.882712 0.07618293 0.640179 0.9695455
## psi_unit27 0.882712 0.07618293 0.640179 0.9695455
## psi_unit28 0.882712 0.07618293 0.640179 0.9695455
## psi_unit29 0.882712 0.07618293 0.640179 0.9695455
## psi_unit30 0.882712 0.07618293 0.640179 0.9695455
## psi_unit31 0.882712 0.07618293 0.640179 0.9695455
## psi_unit32 0.882712 0.07618293 0.640179 0.9695455
## psi_unit33 0.882712 0.07618293 0.640179 0.9695455
## psi_unit34 0.882712 0.07618293 0.640179 0.9695455
## psi_unit35 0.882712 0.07618293 0.640179 0.9695455
## psi_unit36 0.882712 0.07618293 0.640179 0.9695455
## psi_unit37 0.882712 0.07618293 0.640179 0.9695455
## psi_unit38 0.882712 0.07618293 0.640179 0.9695455
## psi_unit39 0.882712 0.07618293 0.640179 0.9695455
## psi_unit40 0.882712 0.07618293 0.640179 0.9695455
## psi_unit41 0.882712 0.07618293 0.640179 0.9695455
RPresence::print_one_site_estimates(mod.pt, site = 1)
## psi()p(SURVEY)
## est se lower_0.95 upper_0.95
## psi_psi_unit1 0.8827120 0.07618293 0.6401790 0.9695455
## p_p1_unit1 0.4973585 0.08915117 0.3297054 0.6656078
# Model averaging
models<-list(mod.pdot,mod.pt)
results<-RPresence::createAicTable(models)
summary(results)
## Model DAIC wgt npar neg2ll warn.conv warn.VC
## 1 psi()p() 0.00 0.83 2 168.19 0 0
## 2 psi()p(SURVEY) 3.11 0.17 4 167.29 0 0
RPresence::modAvg(results, param="psi")[1,]
## est se lower_0.95 upper_0.95
## psi_unit1 0.884435 0.07652202 0.6382408 0.9707587