# Single Species Single Season Occupancy models using RMark
# Blue Gross Beaks.
# Downloaded from https://sites.google.com/site/asrworkshop/home/schedule/r-occupancy-1
# Changed BQI to Y/N from 0/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.
# bqi - bobwhite quality initiative applied to this field
# field.size - size of the files
# crop.hist - crop history
# crop1, crop2 - indicator variables for the crop history
# count1, count2, count3 - are actual counts of birds detected in each visit
# 2018-08-15 Code contributed by Carl James Schwarz (cschwarz.stat.sfu.cs@gmail.com)
# RMark package
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 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 X
## 1 1 1 1 1 14.0 y crop 1 0 1 2 2 NA
## 2 2 1 1 0 12.7 y crop 1 0 2 2 0 NA
## 3 3 0 0 0 15.7 n mixed 0 1 0 0 0 NA
## 4 4 0 1 0 19.5 n mixed 0 1 0 2 0 NA
## 5 5 1 0 1 13.5 n crop 1 0 1 0 1 NA
## 6 6 0 0 1 9.6 n mixed 0 1 0 0 2 NA
## logFS
## 1 1.1461280
## 2 1.1038037
## 3 1.1958997
## 4 1.2900346
## 5 1.1303338
## 6 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","bqi")]
site.covar$logFS <- log(site.covar$field.size)
head(site.covar)
## field field.size bqi logFS
## 1 1 14.0 y 2.639057
## 2 2 12.7 y 2.541602
## 3 3 15.7 n 2.753661
## 4 4 19.5 n 2.970414
## 5 5 13.5 n 2.602690
## 6 6 9.6 n 2.261763
# Extract the history records and create a capture history
input.history <- data.frame(freq=1,
ch=apply(input.history,1,paste, collapse=""),
stringsAsFactors=FALSE)
head(input.history)
## freq ch
## 1 1 111
## 2 1 110
## 3 1 000
## 4 1 010
## 5 1 101
## 6 1 001
# Change any NA to . in the capture history
input.history$ch <- gsub("NA",".", input.history$ch, fixed=TRUE)
head(input.history)
## freq ch
## 1 1 111
## 2 1 110
## 3 1 000
## 4 1 010
## 5 1 101
## 6 1 001
#Add field size data to the capture history
input.history = cbind(input.history, site.covar)
head(input.history)
## freq ch field field.size bqi logFS
## 1 1 111 1 14.0 y 2.639057
## 2 1 110 2 12.7 y 2.541602
## 3 1 000 3 15.7 n 2.753661
## 4 1 010 4 19.5 n 2.970414
## 5 1 101 5 13.5 n 2.602690
## 6 1 001 6 9.6 n 2.261763
#Create the data structure
grossbeak.data <- process.data(data = input.history,
group="bqi",
model = "Occupancy")
## Warning in process.data(data = input.history, group = "bqi", model = "Occupancy"):
## bqi is not a factor variable. Coercing to factor.
summary(grossbeak.data)
## Length Class Mode
## data 7 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 2 data.frame list
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 3 -none- numeric
## begin.time 1 -none- numeric
## age.unit 1 -none- numeric
## initial.ages 2 -none- numeric
## group.covariates 1 data.frame list
## 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
# modify the ddl if needed (e.g. for site level covariates)
grossbeak.ddl <- make.design.data(grossbeak.data)
grossbeak.ddl
## $p
## par.index model.index group age time Age Time bqi
## 1 1 1 n 0 1 0 0 n
## 2 2 2 n 1 2 1 1 n
## 3 3 3 n 2 3 2 2 n
## 4 4 4 y 0 1 0 0 y
## 5 5 5 y 1 2 1 1 y
## 6 6 6 y 2 3 2 2 y
##
## $Psi
## par.index model.index group age time Age Time bqi
## 1 1 7 n 0 1 0 0 n
## 2 2 8 y 0 1 0 0 y
##
## $pimtypes
## $pimtypes$p
## $pimtypes$p$pim.type
## [1] "all"
##
##
## $pimtypes$Psi
## $pimtypes$Psi$pim.type
## [1] "all"
# What are the parameter names for Single Season Single Species models
setup.parameters("Occupancy", check=TRUE)
## [1] "p" "Psi"
# Get the list of models. NOtice NO equal signs here
model.list.csv <- textConnection("
p, Psi
~1, ~1
~time, ~1
~1, ~logFS
~1, ~bqi")
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.number
## 1 ~1 ~1 1
## 2 ~time ~1 2
## 3 ~1 ~logFS 3
## 4 ~1 ~bqi 4
# 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")
fit <- RMark::mark(input.data, ddl=input.ddl,
model="Occupancy",
model.parameters=list(
Psi =list(formula=as.formula(eval(x$Psi))),
p =list(formula=as.formula(eval(x$p)))
)
#,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=grossbeak.data, input.ddl=grossbeak.ddl)
##
##
## ***** Starting ~1 ~1 1
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~1)
##
## Npar : 2
## -2lnL: 168.1898
## AICc : 172.5055
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.2059921 0.2420584 -0.2684423 0.6804266
## Psi:(Intercept) 2.0386894 0.7514126 0.5659207 3.5114580
##
##
## Real Parameter p
## 1 2 3
## Group:bqin 0.5513167 0.5513167 0.5513167
## Group:bqiy 0.5513167 0.5513167 0.5513167
##
##
## Real Parameter Psi
## 1
## Group:bqin 0.8847997
## Group:bqiy 0.8847997
##
##
## ***** Starting ~time ~1 2
##
## Output summary for Occupancy model
## Name : p(~time)Psi(~1)
##
## Npar : 4
## -2lnL: 167.295
## AICc : 176.4061
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.0105659 0.3566281 -0.7095571 0.6884252
## p:time2 0.4489867 0.4772539 -0.4864310 1.3844044
## p:time3 0.2218307 0.4717450 -0.7027895 1.1464509
## Psi:(Intercept) 2.0183671 0.7358609 0.5760797 3.4606545
##
##
## Real Parameter p
## 1 2 3
## Group:bqin 0.4973585 0.6078827 0.5526206
## Group:bqiy 0.4973585 0.6078827 0.5526206
##
##
## Real Parameter Psi
## 1
## Group:bqin 0.8827121
## Group:bqiy 0.8827121
##
##
## ***** Starting ~1 ~logFS 3
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~logFS)
##
## Npar : 3
## -2lnL: 167.9232
## AICc : 174.5719
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.2128360 0.240008 -0.2575797 0.6832518
## Psi:(Intercept) 3.4441599 2.901345 -2.2424764 9.1307962
## Psi:logFS -0.5325425 1.010487 -2.5130972 1.4480121
##
##
## Real Parameter p
## 1 2 3
## Group:bqin 0.5530091 0.5530091 0.5530091
## Group:bqiy 0.5530091 0.5530091 0.5530091
##
##
## Real Parameter Psi
## 1
## Group:bqin 0.8851824
## Group:bqiy 0.8851824
##
##
## ***** Starting ~1 ~bqi 4
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~bqi)
##
## Npar : 3
## -2lnL: 167.1211
## AICc : 173.7697
##
## Beta
## estimate se lcl ucl
## p:(Intercept) 0.2059921 0.2420585 -0.2684424 0.6804267
## Psi:(Intercept) 2.6900661 1.4090681 -0.0717074 5.4518395
## Psi:bqiy -1.3937646 1.5516219 -4.4349436 1.6474144
##
##
## Real Parameter p
## 1 2 3
## Group:bqin 0.5513167 0.5513167 0.5513167
## Group:bqiy 0.5513167 0.5513167 0.5513167
##
##
## Real Parameter Psi
## 1
## Group:bqin 0.9364379
## Group:bqiy 0.7852119
# examine individula model results
model.number <-2
summary(model.fits[[model.number]])
## Output summary for Occupancy model
## Name : p(~time)Psi(~1)
##
## Npar : 4
## -2lnL: 167.295
## AICc : 176.4061
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.0105659 0.3566281 -0.7095571 0.6884252
## p:time2 0.4489867 0.4772539 -0.4864310 1.3844044
## p:time3 0.2218307 0.4717450 -0.7027895 1.1464509
## Psi:(Intercept) 2.0183671 0.7358609 0.5760797 3.4606545
##
##
## Real Parameter p
## 1 2 3
## Group:bqin 0.4973585 0.6078827 0.5526206
## Group:bqiy 0.4973585 0.6078827 0.5526206
##
##
## Real Parameter Psi
## 1
## Group:bqin 0.8827121
## Group:bqiy 0.8827121
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## p gn a0 t1 0.4973585 0.0891545 0.3296967 0.6656165
## p gn a1 t2 0.6078827 0.0902286 0.4246992 0.7650113
## p gn a2 t3 0.5526206 0.0900910 0.3768454 0.7161592
## Psi gn a0 t1 0.8827121 0.0761848 0.6401649 0.9695473
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## p:(Intercept) -0.0105659 0.3566281 -0.7095571 0.6884252
## p:time2 0.4489867 0.4772539 -0.4864310 1.3844044
## p:time3 0.2218307 0.4717450 -0.7027895 1.1464509
## Psi:(Intercept) 2.0183671 0.7358609 0.5760797 3.4606545
model.fits[[model.number]]$results$derived
## $Occupancy
## estimate lcl ucl
## 1 0.8827121 0.6401649 0.9695473
## 2 0.8827121 0.6401649 0.9695473
get.real(model.fits[[model.number]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gn a0 t1 7 4 0.8827121 0.0761848 0.6401649 0.9695473
## Psi gy a0 t1 8 4 0.8827121 0.0761848 0.6401649 0.9695473
## fixed note group age time Age Time bqi
## Psi gn a0 t1 n 0 1 0 0 n
## Psi gy a0 t1 y 0 1 0 0 y
get.real(model.fits[[model.number]], "p", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## p gn a0 t1 1 1 0.4973585 0.0891545 0.3296967 0.6656165
## p gn a1 t2 2 2 0.6078827 0.0902286 0.4246992 0.7650113
## p gn a2 t3 3 3 0.5526206 0.0900910 0.3768454 0.7161592
## p gy a0 t1 4 1 0.4973585 0.0891545 0.3296967 0.6656165
## p gy a1 t2 5 2 0.6078827 0.0902286 0.4246992 0.7650113
## p gy a2 t3 6 3 0.5526206 0.0900910 0.3768454 0.7161592
## fixed note group age time Age Time bqi
## p gn a0 t1 n 0 1 0 0 n
## p gn a1 t2 n 1 2 1 1 n
## p gn a2 t3 n 2 3 2 2 n
## p gy a0 t1 y 0 1 0 0 y
## p gy a1 t2 y 1 2 1 1 y
## p gy a2 t3 y 2 3 2 2 y
# collect models and make AICc table
model.set <- RMark::collect.models( type="Occupancy")
model.set
## model npar AICc DeltaAICc weight Deviance
## 1 p(~1)Psi(~1) 2 172.5055 0.000000 0.49270963 13.55421
## 4 p(~1)Psi(~bqi) 3 173.7697 1.264169 0.26186666 12.48551
## 3 p(~1)Psi(~logFS) 3 174.5719 2.066329 0.17534499 167.92322
## 2 p(~time)Psi(~1) 4 176.4061 3.900602 0.07007872 12.65949
names(model.set)
## [1] "m...001" "m...002" "m...003" "m...004" "model.table"
model.set$model.table
## p Psi model npar AICc DeltaAICc weight Deviance
## 1 ~1 ~1 p(~1)Psi(~1) 2 172.5055 0.000000 0.49270963 13.55421
## 4 ~1 ~bqi p(~1)Psi(~bqi) 3 173.7697 1.264169 0.26186666 12.48551
## 3 ~1 ~logFS p(~1)Psi(~logFS) 3 174.5719 2.066329 0.17534499 167.92322
## 2 ~time ~1 p(~time)Psi(~1) 4 176.4061 3.900602 0.07007872 12.65949
# model averaged values
get.real(model.set[[1]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gn a0 t1 7 2 0.8847997 0.0765908 0.6378214 0.971012
## Psi gy a0 t1 8 2 0.8847997 0.0765908 0.6378214 0.971012
## fixed note group age time Age Time bqi
## Psi gn a0 t1 n 0 1 0 0 n
## Psi gy a0 t1 y 0 1 0 0 y
get.real(model.set[[2]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gn a0 t1 7 4 0.8827121 0.0761848 0.6401649 0.9695473
## Psi gy a0 t1 8 4 0.8827121 0.0761848 0.6401649 0.9695473
## fixed note group age time Age Time bqi
## Psi gn a0 t1 n 0 1 0 0 n
## Psi gy a0 t1 y 0 1 0 0 y
get.real(model.set[[3]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gn a0 t1 7 2 0.8851824 0.0757111 0.6416142 0.9707594
## Psi gy a0 t1 8 2 0.8851824 0.0757111 0.6416142 0.9707594
## fixed note group age time Age Time bqi
## Psi gn a0 t1 n 0 1 0 0 n
## Psi gy a0 t1 y 0 1 0 0 y
# covariate predictions for continuous covariates
# such as logFS
# need to find the index number of the parameter
grossbeak.ddl$Psi
## par.index model.index group age time Age Time bqi
## 1 1 7 n 0 1 0 0 n
## 2 2 8 y 0 1 0 0 y
new.logFS <- data.frame(logFS=seq(min(site.covar$logFS),
max(site.covar$logFS), .05))
psi.ma <- covariate.predictions(model.set,
indices=7,
data=new.logFS)
psi.ma$estimates
## vcv.index model.index par.index covdata estimate se lcl
## 1 1 2 7 1.856298 0.9045180 0.08354983 0.5871974
## 2 2 2 7 1.906298 0.9041744 0.08330891 0.5890285
## 3 3 2 7 1.956298 0.9038230 0.08306975 0.5908425
## 4 4 2 7 2.006298 0.9034637 0.08283458 0.5926272
## 5 5 2 7 2.056298 0.9030963 0.08260588 0.5943689
## 6 6 2 7 2.106298 0.9027209 0.08238642 0.5960529
## 7 7 2 7 2.156298 0.9023371 0.08217928 0.5976630
## 8 8 2 7 2.206298 0.9019448 0.08198787 0.5991814
## 9 9 2 7 2.256298 0.9015440 0.08181591 0.6005893
## 10 10 2 7 2.306298 0.9011345 0.08166749 0.6018661
## 11 11 2 7 2.356298 0.9007161 0.08154707 0.6029899
## 12 12 2 7 2.406298 0.9002887 0.08145946 0.6039372
## 13 13 2 7 2.456298 0.8998522 0.08140986 0.6046829
## 14 14 2 7 2.506298 0.8994064 0.08140382 0.6052007
## 15 15 2 7 2.556298 0.8989512 0.08144726 0.6054626
## 16 16 2 7 2.606298 0.8984865 0.08154644 0.6054395
## 17 17 2 7 2.656298 0.8980121 0.08170793 0.6051010
## 18 18 2 7 2.706298 0.8975280 0.08193860 0.6044155
## 19 19 2 7 2.756298 0.8970339 0.08224555 0.6033504
## 20 20 2 7 2.806298 0.8965297 0.08263609 0.6018726
## 21 21 2 7 2.856298 0.8960154 0.08311765 0.5999484
## 22 22 2 7 2.906298 0.8954908 0.08369774 0.5975438
## 23 23 2 7 2.956298 0.8949557 0.08438383 0.5946251
## 24 24 2 7 3.006298 0.8944101 0.08518336 0.5911589
## 25 25 2 7 3.056298 0.8938538 0.08610352 0.5871131
## 26 26 2 7 3.106298 0.8932868 0.08715130 0.5824567
## 27 27 2 7 3.156298 0.8927089 0.08833335 0.5771605
## 28 28 2 7 3.206298 0.8921200 0.08965590 0.5711981
## 29 29 2 7 3.256298 0.8915200 0.09112472 0.5645457
## 30 30 2 7 3.306298 0.8909088 0.09274503 0.5571833
## 31 31 2 7 3.356298 0.8902864 0.09452147 0.5490950
## 32 32 2 7 3.406298 0.8896526 0.09645812 0.5402696
## 33 33 2 7 3.456298 0.8890073 0.09855831 0.5307017
## 34 34 2 7 3.506298 0.8883505 0.10082484 0.5203915
## 35 35 2 7 3.556298 0.8876821 0.10325984 0.5093458
## 36 36 2 7 3.606298 0.8870021 0.10586483 0.4975786
## 37 37 2 7 3.656298 0.8863103 0.10864072 0.4851113
## 38 38 2 7 3.706298 0.8856067 0.11158787 0.4719731
## 39 39 2 7 3.756298 0.8848913 0.11470611 0.4582013
## 40 40 2 7 3.806298 0.8841640 0.11799484 0.4438407
## 41 41 2 7 3.856298 0.8834248 0.12145284 0.4289450
## 42 42 2 7 3.906298 0.8826737 0.12507867 0.4135745
## 43 43 2 7 3.956298 0.8819107 0.12887042 0.3977967
## ucl fixed
## 1 0.9843966
## 2 0.9841566
## 3 0.9839116
## 4 0.9836621
## 5 0.9834091
## 6 0.9831533
## 7 0.9828961
## 8 0.9826386
## 9 0.9823823
## 10 0.9821289
## 11 0.9818804
## 12 0.9816386
## 13 0.9814060
## 14 0.9811851
## 15 0.9809783
## 16 0.9807887
## 17 0.9806191
## 18 0.9804725
## 19 0.9803521
## 20 0.9802610
## 21 0.9802020
## 22 0.9801782
## 23 0.9801921
## 24 0.9802460
## 25 0.9803418
## 26 0.9804811
## 27 0.9806648
## 28 0.9808932
## 29 0.9811663
## 30 0.9814830
## 31 0.9818421
## 32 0.9822414
## 33 0.9826783
## 34 0.9831497
## 35 0.9836521
## 36 0.9841817
## 37 0.9847343
## 38 0.9853056
## 39 0.9858914
## 40 0.9864873
## 41 0.9870892
## 42 0.9876930
## 43 0.9882949
ggplot(data=psi.ma$estimates, aes(x=covdata, y=estimate))+
geom_point()+
geom_ribbon(aes(ymin=lcl, ymax=ucl), alpha=0.2)
# covariate predictions for categorical covariates
# such a bqi
grossbeak.ddl$Psi
## par.index model.index group age time Age Time bqi
## 1 1 7 n 0 1 0 0 n
## 2 2 8 y 0 1 0 0 y
psi.ma <- model.average(model.set, param="Psi", vcv=TRUE)
psi.ma$estimates
## par.index estimate se lcl ucl fixed note
## Psi gn a0 t1 7 0.8982429 0.08162158 0.6053062 0.9806984
## Psi gy a0 t1 8 0.8586418 0.10527853 0.5259844 0.9708036
## group age time Age Time bqi
## Psi gn a0 t1 n 0 1 0 0 n
## Psi gy a0 t1 y 0 1 0 0 y
ggplot(data=psi.ma$estimates, aes(x=bqi, y=estimate))+
geom_point()+
geom_errorbar(aes(ymin=lcl,
ymax=ucl), width=0.2)+
ylim(0,1)
# or use covariate predictions
psi.ma2 <- covariate.predictions(model.set, indices=7:8)
psi.ma2$estimates <- cbind(psi.ma2$estimates, grossbeak.ddl$Psi[,"bqi", drop=FALSE])
psi.ma2$estimates
## vcv.index model.index par.index index estimate se lcl ucl
## 1 2 2 7 7 0.8982428 0.0816216 0.6053061 0.9806984
## 2 2 2 8 8 0.8586418 0.1052785 0.5259845 0.9708036
## fixed bqi
## 1 n
## 2 y
ggplot(data=psi.ma2$estimates, aes(x=bqi, y=estimate))+
geom_point()+
geom_errorbar(aes(ymin=lcl,
ymax=ucl), width=0.2)+
ylim(0,1)
# caution when you have both continous and categorical covariates
# as it is not clear what value of the continuous covariate is used
psi.ma
## $estimates
## par.index estimate se lcl ucl fixed note
## Psi gn a0 t1 7 0.8982429 0.08162158 0.6053062 0.9806984
## Psi gy a0 t1 8 0.8586418 0.10527853 0.5259844 0.9708036
## group age time Age Time bqi
## Psi gn a0 t1 n 0 1 0 0 n
## Psi gy a0 t1 y 0 1 0 0 y
##
## $vcv.real
## 7 8
## 7 0.006662083 0.006571441
## 8 0.006571441 0.011083568
psi.ma2$estimates
## vcv.index model.index par.index index estimate se lcl ucl
## 1 2 2 7 7 0.8982428 0.0816216 0.6053061 0.9806984
## 2 2 2 8 8 0.8586418 0.1052785 0.5259845 0.9708036
## fixed bqi
## 1 n
## 2 y
# It is often convenient to estimate occupancy for each site using the
# values of the covariates specific for that site.
# This is similar to RPresence which gives estimates at each site.
# We create a data frame with the covariates etc as measured on the input.history dataframe.
# But we also need to add the appropriate index number for psi to the covariate data frame to account
# for the groups definition
site.covariates <- input.history
site.covariates$freq <- NULL
site.covariates$ch <- NULL
grossbeak.ddl$Psi
## par.index model.index group age time Age Time bqi
## 1 1 7 n 0 1 0 0 n
## 2 2 8 y 0 1 0 0 y
site.covariates <- merge(grossbeak.ddl$Psi, site.covariates)
head(site.covariates)
## bqi par.index model.index group age time Age Time field field.size logFS
## 1 n 1 7 n 0 1 0 0 27 10.5 2.351375
## 2 n 1 7 n 0 1 0 0 6 9.6 2.261763
## 3 n 1 7 n 0 1 0 0 3 15.7 2.753661
## 4 n 1 7 n 0 1 0 0 4 19.5 2.970414
## 5 n 1 7 n 0 1 0 0 5 13.5 2.602690
## 6 n 1 7 n 0 1 0 0 10 7.0 1.945910
# we only want a prediction for that particular model.index.
# we need to create a variable "index" that has the model.index
site.covariates$index <- site.covariates$model.index
site.pred <- covariate.predictions(model.set,
data=site.covariates)
ggplot(site.pred$estimates, aes(x=logFS, y=estimate))+
ggtitle("Model averaged predictions of occupancy for each site")+
geom_point()+
geom_ribbon(aes(ymin=lcl, ymax=ucl), alpha=0.2)+
facet_wrap(~bqi, ncol=1)
#cleanup
cleanup(ask=FALSE)