# 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
# Create the *.pao file
grossbeak.pao <- RPresence::createPao(input.history,
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] 1
##
## $unitcov
## TEMP
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6
## 7 7
## 8 8
## 9 9
## 10 10
## 11 11
## 12 12
## 13 13
## 14 14
## 15 15
## 16 16
## 17 17
## 18 18
## 19 19
## 20 20
## 21 21
## 22 22
## 23 23
## 24 24
## 25 25
## 26 26
## 27 27
## 28 28
## 29 29
## 30 30
## 31 31
## 32 32
## 33 33
## 34 34
## 35 35
## 36 36
## 37 37
## 38 38
## 39 39
## 40 40
## 41 41
##
## $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
# 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