# Example using model averaging on the pronghorn dataset
# 256 sites, two sampling occasions per site
# Covariates used in this example:
# 1. sagebrush (continuos) - Sagebrush density
# 2. aspect (Categorical) - Compass direction slope faces (N,S,E,W)
# Code contributed by Neil Faught - 26/11/2018
# 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("..","pronghorn.csv"), header=TRUE, as.is=TRUE, strip.white=TRUE)
head(input.data)
## Plot Survey.1 Survey.2 sagebrush slope DW aspect
## 1 1 0 0 9.0 0 25 W
## 2 2 0 1 18.0 5 150 S
## 3 3 0 0 8.4 45 150 W
## 4 4 0 0 3.2 65 375 E
## 5 5 0 1 12.0 5 375 S
## 6 6 1 1 7.8 5 150 S
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.data)
## [1] 256
range(input.data[, c("Survey.1","Survey.2")], na.rm=TRUE)
## [1] 0 1
sum(is.na(input.data[, c("Survey.1","Survey.2")]))
## [1] 0
input.history <- input.data[, c("Survey.1","Survey.2")]
head(input.history)
## Survey.1 Survey.2
## 1 0 0
## 2 0 1
## 3 0 0
## 4 0 0
## 5 0 1
## 6 1 1
site.covar <- input.data[, c("sagebrush","aspect")]
head(site.covar)
## sagebrush aspect
## 1 9.0 W
## 2 18.0 S
## 3 8.4 W
## 4 3.2 E
## 5 12.0 S
## 6 7.8 S
#Convert aspect to a factor
site.covar$aspect = as.factor(site.covar$aspect)
# 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 00
## 2 1 01
## 3 1 00
## 4 1 00
## 5 1 01
## 6 1 11
# 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 00
## 2 1 01
## 3 1 00
## 4 1 00
## 5 1 01
## 6 1 11
#Add aspect and sagebrush data to the capture history
input.history = cbind(input.history, site.covar)
head(input.history)
## freq ch sagebrush aspect
## 1 1 00 9.0 W
## 2 1 01 18.0 S
## 3 1 00 8.4 W
## 4 1 00 3.2 E
## 5 1 01 12.0 S
## 6 1 11 7.8 S
#Create the data structure
pronghorn.data <- process.data(data = input.history,
group="aspect",
model = "Occupancy")
summary(pronghorn.data)
## Length Class Mode
## data 5 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 4 data.frame list
## nocc 1 -none- numeric
## nocc.secondary 0 -none- NULL
## time.intervals 2 -none- numeric
## begin.time 1 -none- numeric
## age.unit 1 -none- numeric
## initial.ages 4 -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)
pronghorn.ddl <- make.design.data(pronghorn.data)
pronghorn.ddl
## $p
## par.index model.index group age time Age Time aspect
## 1 1 1 E 0 1 0 0 E
## 2 2 2 E 1 2 1 1 E
## 3 3 3 N 0 1 0 0 N
## 4 4 4 N 1 2 1 1 N
## 5 5 5 S 0 1 0 0 S
## 6 6 6 S 1 2 1 1 S
## 7 7 7 W 0 1 0 0 W
## 8 8 8 W 1 2 1 1 W
##
## $Psi
## par.index model.index group age time Age Time aspect
## 1 1 9 E 0 1 0 0 E
## 2 2 10 N 0 1 0 0 N
## 3 3 11 S 0 1 0 0 S
## 4 4 12 W 0 1 0 0 W
##
## $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, ~aspect
~1, ~sagebrush")
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 ~aspect 3
## 4 ~1 ~sagebrush 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=pronghorn.data, input.ddl=pronghorn.ddl)
##
##
## ***** Starting ~1 ~1 1
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~1)
##
## Npar : 2
## -2lnL: 635.3553
## AICc : 639.4027
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.1690764 0.1934604 -0.5482587 0.2101059
## Psi:(Intercept) 0.8864741 0.3312481 0.2372279 1.5357203
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4578313 0.4578313
## Group:aspectN 0.4578313 0.4578313
## Group:aspectS 0.4578313 0.4578313
## Group:aspectW 0.4578313 0.4578313
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.708162
## Group:aspectN 0.708162
## Group:aspectS 0.708162
## Group:aspectW 0.708162
##
##
## ***** Starting ~time ~1 2
##
## Output summary for Occupancy model
## Name : p(~time)Psi(~1)
##
## Npar : 3
## -2lnL: 629.9226
## AICc : 636.0179
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.3877655 0.2101736 -0.7997058 0.0241747
## p:time2 0.4989912 0.2174142 0.0728594 0.9251229
## Psi:(Intercept) 0.8270163 0.3137117 0.2121413 1.4418914
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4042553 0.5277778
## Group:aspectN 0.4042553 0.5277778
## Group:aspectS 0.4042553 0.5277778
## Group:aspectW 0.4042553 0.5277778
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.6957237
## Group:aspectN 0.6957237
## Group:aspectS 0.6957237
## Group:aspectW 0.6957237
##
##
## ***** Starting ~1 ~aspect 3
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~aspect)
##
## Npar : 5
## -2lnL: 631.3945
## AICc : 641.6345
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.1690763 0.1934605 -0.5482588 0.2101062
## Psi:(Intercept) 0.1936200 0.5262422 -0.8378148 1.2250547
## Psi:aspectN 2.1260634 1.7786193 -1.3600305 5.6121573
## Psi:aspectS 0.7702649 0.7347789 -0.6699018 2.2104317
## Psi:aspectW 0.6481796 0.5756880 -0.4801689 1.7765281
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4578313 0.4578313
## Group:aspectN 0.4578313 0.4578313
## Group:aspectS 0.4578313 0.4578313
## Group:aspectW 0.4578313 0.4578313
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.5482543
## Group:aspectN 0.9104941
## Group:aspectS 0.7238990
## Group:aspectW 0.6988441
##
##
## ***** Starting ~1 ~sagebrush 4
##
## Output summary for Occupancy model
## Name : p(~1)Psi(~sagebrush)
##
## Npar : 3
## -2lnL: 634.5521
## AICc : 640.6473
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.1631302 0.1920302 -0.5395094 0.2132490
## Psi:(Intercept) 0.6837763 0.3725475 -0.0464167 1.4139694
## Psi:sagebrush 0.0476779 0.0554558 -0.0610155 0.1563713
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4593077 0.4593077
## Group:aspectN 0.4593077 0.4593077
## Group:aspectS 0.4593077 0.4593077
## Group:aspectW 0.4593077 0.4593077
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7072814
## Group:aspectN 0.7072814
## Group:aspectS 0.7072814
## Group:aspectW 0.7072814
# examine individula model results
model.number <-4
summary(model.fits[[model.number]])
## Output summary for Occupancy model
## Name : p(~1)Psi(~sagebrush)
##
## Npar : 3
## -2lnL: 634.5521
## AICc : 640.6473
##
## Beta
## estimate se lcl ucl
## p:(Intercept) -0.1631302 0.1920302 -0.5395094 0.2132490
## Psi:(Intercept) 0.6837763 0.3725475 -0.0464167 1.4139694
## Psi:sagebrush 0.0476779 0.0554558 -0.0610155 0.1563713
##
##
## Real Parameter p
## 1 2
## Group:aspectE 0.4593077 0.4593077
## Group:aspectN 0.4593077 0.4593077
## Group:aspectS 0.4593077 0.4593077
## Group:aspectW 0.4593077 0.4593077
##
##
## Real Parameter Psi
## 1
## Group:aspectE 0.7072814
## Group:aspectN 0.7072814
## Group:aspectS 0.7072814
## Group:aspectW 0.7072814
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## p gE a0 t1 0.4593077 0.0476896 0.3683017 0.5531111
## Psi gE a0 t1 0.7072814 0.0681243 0.5590438 0.8215892
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## p:(Intercept) -0.1631302 0.1920302 -0.5395094 0.2132490
## Psi:(Intercept) 0.6837763 0.3725475 -0.0464167 1.4139694
## Psi:sagebrush 0.0476779 0.0554558 -0.0610155 0.1563713
model.fits[[model.number]]$results$derived
## $Occupancy
## estimate lcl ucl
## 1 0.7072814 0.5590438 0.8215892
## 2 0.7072814 0.5590438 0.8215892
## 3 0.7072814 0.5590438 0.8215892
## 4 0.7072814 0.5590438 0.8215892
get.real(model.fits[[model.number]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gE a0 t1 9 2 0.7072814 0.0681243 0.5590438 0.8215892
## Psi gN a0 t1 10 2 0.7072814 0.0681243 0.5590438 0.8215892
## Psi gS a0 t1 11 2 0.7072814 0.0681243 0.5590438 0.8215892
## Psi gW a0 t1 12 2 0.7072814 0.0681243 0.5590438 0.8215892
## fixed note group age time Age Time aspect
## Psi gE a0 t1 E 0 1 0 0 E
## Psi gN a0 t1 N 0 1 0 0 N
## Psi gS a0 t1 S 0 1 0 0 S
## Psi gW a0 t1 W 0 1 0 0 W
get.real(model.fits[[model.number]], "p", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## p gE a0 t1 1 1 0.4593077 0.0476896 0.3683017 0.5531111
## p gE a1 t2 2 1 0.4593077 0.0476896 0.3683017 0.5531111
## p gN a0 t1 3 1 0.4593077 0.0476896 0.3683017 0.5531111
## p gN a1 t2 4 1 0.4593077 0.0476896 0.3683017 0.5531111
## p gS a0 t1 5 1 0.4593077 0.0476896 0.3683017 0.5531111
## p gS a1 t2 6 1 0.4593077 0.0476896 0.3683017 0.5531111
## p gW a0 t1 7 1 0.4593077 0.0476896 0.3683017 0.5531111
## p gW a1 t2 8 1 0.4593077 0.0476896 0.3683017 0.5531111
## fixed note group age time Age Time aspect
## p gE a0 t1 E 0 1 0 0 E
## p gE a1 t2 E 1 2 1 1 E
## p gN a0 t1 N 0 1 0 0 N
## p gN a1 t2 N 1 2 1 1 N
## p gS a0 t1 S 0 1 0 0 S
## p gS a1 t2 S 1 2 1 1 S
## p gW a0 t1 W 0 1 0 0 W
## p gW a1 t2 W 1 2 1 1 W
# collect models and make AICc table
model.set <- RMark::collect.models( type="Occupancy")
model.set
## model npar AICc DeltaAICc weight Deviance
## 2 p(~time)Psi(~1) 3 636.0179 0.000000 0.74450664 11.97795
## 1 p(~1)Psi(~1) 2 639.4027 3.384853 0.13704309 17.41061
## 4 p(~1)Psi(~sagebrush) 3 640.6473 4.629460 0.07355194 634.55211
## 3 p(~1)Psi(~aspect) 5 641.6345 5.616642 0.04489834 13.44983
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
## 2 ~time ~1 p(~time)Psi(~1) 3 636.0179 0.000000 0.74450664
## 1 ~1 ~1 p(~1)Psi(~1) 2 639.4027 3.384853 0.13704309
## 4 ~1 ~sagebrush p(~1)Psi(~sagebrush) 3 640.6473 4.629460 0.07355194
## 3 ~1 ~aspect p(~1)Psi(~aspect) 5 641.6345 5.616642 0.04489834
## Deviance
## 2 11.97795
## 1 17.41061
## 4 634.55211
## 3 13.44983
# model averaged values
get.real(model.set[[1]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gE a0 t1 9 2 0.708162 0.0684586 0.5590304 0.8228417
## Psi gN a0 t1 10 2 0.708162 0.0684586 0.5590304 0.8228417
## Psi gS a0 t1 11 2 0.708162 0.0684586 0.5590304 0.8228417
## Psi gW a0 t1 12 2 0.708162 0.0684586 0.5590304 0.8228417
## fixed note group age time Age Time aspect
## Psi gE a0 t1 E 0 1 0 0 E
## Psi gN a0 t1 N 0 1 0 0 N
## Psi gS a0 t1 S 0 1 0 0 S
## Psi gW a0 t1 W 0 1 0 0 W
get.real(model.set[[2]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gE a0 t1 9 3 0.6957237 0.0664103 0.5528373 0.8087474
## Psi gN a0 t1 10 3 0.6957237 0.0664103 0.5528373 0.8087474
## Psi gS a0 t1 11 3 0.6957237 0.0664103 0.5528373 0.8087474
## Psi gW a0 t1 12 3 0.6957237 0.0664103 0.5528373 0.8087474
## fixed note group age time Age Time aspect
## Psi gE a0 t1 E 0 1 0 0 E
## Psi gN a0 t1 N 0 1 0 0 N
## Psi gS a0 t1 S 0 1 0 0 S
## Psi gW a0 t1 W 0 1 0 0 W
get.real(model.set[[3]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gE a0 t1 9 2 0.5482543 0.1303352 0.3019952 0.7729519
## Psi gN a0 t1 10 3 0.9104941 0.1447667 0.2382940 0.9969859
## Psi gS a0 t1 11 4 0.7238990 0.1182751 0.4511638 0.8931892
## Psi gW a0 t1 12 5 0.6988441 0.0771689 0.5307447 0.8264201
## fixed note group age time Age Time aspect
## Psi gE a0 t1 E 0 1 0 0 E
## Psi gN a0 t1 N 0 1 0 0 N
## Psi gS a0 t1 S 0 1 0 0 S
## Psi gW a0 t1 W 0 1 0 0 W
get.real(model.set[[4]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl ucl
## Psi gE a0 t1 9 2 0.7072814 0.0681243 0.5590438 0.8215892
## Psi gN a0 t1 10 2 0.7072814 0.0681243 0.5590438 0.8215892
## Psi gS a0 t1 11 2 0.7072814 0.0681243 0.5590438 0.8215892
## Psi gW a0 t1 12 2 0.7072814 0.0681243 0.5590438 0.8215892
## fixed note group age time Age Time aspect
## Psi gE a0 t1 E 0 1 0 0 E
## Psi gN a0 t1 N 0 1 0 0 N
## Psi gS a0 t1 S 0 1 0 0 S
## Psi gW a0 t1 W 0 1 0 0 W
# covariate predictions for continuous covariates
# such as sagebrush
# need to find the index number of the parameter
pronghorn.ddl$Psi
## par.index model.index group age time Age Time aspect
## 1 1 9 E 0 1 0 0 E
## 2 2 10 N 0 1 0 0 N
## 3 3 11 S 0 1 0 0 S
## 4 4 12 W 0 1 0 0 W
new.sagebrush <- data.frame(sagebrush=seq(min(site.covar$sagebrush),
max(site.covar$sagebrush), .05))
# Note that this just makes predictions for observations from the east aspect
# Change the indices= line to view plots for different aspects. To view plots
# from all aspects at once, use indices=9:12
psi.ma <- covariate.predictions(model.set,
indices=9,
data=new.sagebrush)
psi.ma$estimates
## vcv.index model.index par.index covdata estimate se lcl
## 1 1 2 9 0.00 0.6885165 0.07882101 0.5182039
## 2 2 2 9 0.05 0.6885556 0.07878108 0.5183281
## 3 3 2 9 0.10 0.6885946 0.07874180 0.5184509
## 4 4 2 9 0.15 0.6886336 0.07870318 0.5185721
## 5 5 2 9 0.20 0.6886726 0.07866520 0.5186918
## 6 6 2 9 0.25 0.6887116 0.07862788 0.5188101
## 7 7 2 9 0.30 0.6887505 0.07859121 0.5189268
## 8 8 2 9 0.35 0.6887894 0.07855518 0.5190420
## 9 9 2 9 0.40 0.6888282 0.07851980 0.5191558
## 10 10 2 9 0.45 0.6888671 0.07848507 0.5192680
## 11 11 2 9 0.50 0.6889058 0.07845098 0.5193788
## 12 12 2 9 0.55 0.6889446 0.07841753 0.5194880
## 13 13 2 9 0.60 0.6889833 0.07838473 0.5195958
## 14 14 2 9 0.65 0.6890220 0.07835257 0.5197021
## 15 15 2 9 0.70 0.6890607 0.07832104 0.5198069
## 16 16 2 9 0.75 0.6890993 0.07829016 0.5199103
## 17 17 2 9 0.80 0.6891379 0.07825991 0.5200121
## 18 18 2 9 0.85 0.6891765 0.07823029 0.5201125
## 19 19 2 9 0.90 0.6892150 0.07820131 0.5202114
## 20 20 2 9 0.95 0.6892535 0.07817297 0.5203089
## 21 21 2 9 1.00 0.6892920 0.07814525 0.5204049
## 22 22 2 9 1.05 0.6893305 0.07811816 0.5204994
## 23 23 2 9 1.10 0.6893689 0.07809170 0.5205925
## 24 24 2 9 1.15 0.6894072 0.07806587 0.5206842
## 25 25 2 9 1.20 0.6894456 0.07804066 0.5207743
## 26 26 2 9 1.25 0.6894839 0.07801608 0.5208631
## 27 27 2 9 1.30 0.6895222 0.07799212 0.5209504
## 28 28 2 9 1.35 0.6895604 0.07796877 0.5210362
## 29 29 2 9 1.40 0.6895986 0.07794605 0.5211207
## 30 30 2 9 1.45 0.6896368 0.07792394 0.5212036
## 31 31 2 9 1.50 0.6896750 0.07790245 0.5212852
## 32 32 2 9 1.55 0.6897131 0.07788157 0.5213653
## 33 33 2 9 1.60 0.6897512 0.07786131 0.5214440
## 34 34 2 9 1.65 0.6897892 0.07784165 0.5215213
## 35 35 2 9 1.70 0.6898272 0.07782260 0.5215972
## 36 36 2 9 1.75 0.6898652 0.07780416 0.5216717
## 37 37 2 9 1.80 0.6899032 0.07778632 0.5217447
## 38 38 2 9 1.85 0.6899411 0.07776908 0.5218164
## 39 39 2 9 1.90 0.6899790 0.07775245 0.5218867
## 40 40 2 9 1.95 0.6900168 0.07773641 0.5219555
## 41 41 2 9 2.00 0.6900547 0.07772097 0.5220230
## 42 42 2 9 2.05 0.6900924 0.07770613 0.5220891
## 43 43 2 9 2.10 0.6901302 0.07769187 0.5221538
## 44 44 2 9 2.15 0.6901679 0.07767821 0.5222171
## 45 45 2 9 2.20 0.6902056 0.07766514 0.5222791
## 46 46 2 9 2.25 0.6902432 0.07765266 0.5223397
## 47 47 2 9 2.30 0.6902809 0.07764076 0.5223989
## 48 48 2 9 2.35 0.6903184 0.07762944 0.5224568
## 49 49 2 9 2.40 0.6903560 0.07761871 0.5225133
## 50 50 2 9 2.45 0.6903935 0.07760855 0.5225684
## 51 51 2 9 2.50 0.6904310 0.07759897 0.5226222
## 52 52 2 9 2.55 0.6904684 0.07758996 0.5226747
## 53 53 2 9 2.60 0.6905058 0.07758153 0.5227259
## 54 54 2 9 2.65 0.6905432 0.07757366 0.5227757
## 55 55 2 9 2.70 0.6905806 0.07756637 0.5228241
## 56 56 2 9 2.75 0.6906179 0.07755964 0.5228713
## 57 57 2 9 2.80 0.6906552 0.07755347 0.5229171
## 58 58 2 9 2.85 0.6906924 0.07754786 0.5229616
## 59 59 2 9 2.90 0.6907296 0.07754282 0.5230048
## 60 60 2 9 2.95 0.6907668 0.07753833 0.5230467
## 61 61 2 9 3.00 0.6908039 0.07753439 0.5230874
## 62 62 2 9 3.05 0.6908410 0.07753101 0.5231267
## 63 63 2 9 3.10 0.6908781 0.07752817 0.5231647
## 64 64 2 9 3.15 0.6909151 0.07752589 0.5232014
## 65 65 2 9 3.20 0.6909521 0.07752414 0.5232369
## 66 66 2 9 3.25 0.6909891 0.07752294 0.5232711
## 67 67 2 9 3.30 0.6910260 0.07752229 0.5233040
## 68 68 2 9 3.35 0.6910629 0.07752217 0.5233357
## 69 69 2 9 3.40 0.6910998 0.07752258 0.5233661
## 70 70 2 9 3.45 0.6911366 0.07752353 0.5233952
## 71 71 2 9 3.50 0.6911734 0.07752501 0.5234232
## 72 72 2 9 3.55 0.6912101 0.07752701 0.5234498
## 73 73 2 9 3.60 0.6912469 0.07752955 0.5234752
## 74 74 2 9 3.65 0.6912835 0.07753260 0.5234994
## 75 75 2 9 3.70 0.6913202 0.07753618 0.5235224
## 76 76 2 9 3.75 0.6913568 0.07754027 0.5235442
## 77 77 2 9 3.80 0.6913934 0.07754488 0.5235647
## 78 78 2 9 3.85 0.6914299 0.07755000 0.5235841
## 79 79 2 9 3.90 0.6914664 0.07755564 0.5236022
## 80 80 2 9 3.95 0.6915029 0.07756178 0.5236191
## 81 81 2 9 4.00 0.6915393 0.07756843 0.5236349
## 82 82 2 9 4.05 0.6915757 0.07757558 0.5236494
## 83 83 2 9 4.10 0.6916121 0.07758323 0.5236628
## 84 84 2 9 4.15 0.6916484 0.07759137 0.5236750
## 85 85 2 9 4.20 0.6916847 0.07760002 0.5236860
## 86 86 2 9 4.25 0.6917210 0.07760915 0.5236959
## 87 87 2 9 4.30 0.6917572 0.07761878 0.5237046
## 88 88 2 9 4.35 0.6917934 0.07762889 0.5237122
## 89 89 2 9 4.40 0.6918295 0.07763949 0.5237186
## 90 90 2 9 4.45 0.6918656 0.07765057 0.5237238
## 91 91 2 9 4.50 0.6919017 0.07766212 0.5237280
## 92 92 2 9 4.55 0.6919378 0.07767416 0.5237310
## 93 93 2 9 4.60 0.6919738 0.07768667 0.5237329
## 94 94 2 9 4.65 0.6920097 0.07769965 0.5237337
## 95 95 2 9 4.70 0.6920457 0.07771309 0.5237333
## 96 96 2 9 4.75 0.6920816 0.07772701 0.5237319
## 97 97 2 9 4.80 0.6921174 0.07774138 0.5237293
## 98 98 2 9 4.85 0.6921532 0.07775622 0.5237257
## 99 99 2 9 4.90 0.6921890 0.07777152 0.5237210
## 100 100 2 9 4.95 0.6922248 0.07778726 0.5237152
## 101 101 2 9 5.00 0.6922605 0.07780347 0.5237083
## 102 102 2 9 5.05 0.6922962 0.07782012 0.5237004
## 103 103 2 9 5.10 0.6923318 0.07783721 0.5236914
## 104 104 2 9 5.15 0.6923674 0.07785476 0.5236813
## 105 105 2 9 5.20 0.6924030 0.07787274 0.5236702
## 106 106 2 9 5.25 0.6924385 0.07789116 0.5236581
## 107 107 2 9 5.30 0.6924740 0.07791002 0.5236449
## 108 108 2 9 5.35 0.6925094 0.07792931 0.5236306
## 109 109 2 9 5.40 0.6925449 0.07794902 0.5236154
## 110 110 2 9 5.45 0.6925802 0.07796917 0.5235991
## 111 111 2 9 5.50 0.6926156 0.07798974 0.5235818
## 112 112 2 9 5.55 0.6926509 0.07801073 0.5235636
## 113 113 2 9 5.60 0.6926862 0.07803215 0.5235443
## 114 114 2 9 5.65 0.6927214 0.07805397 0.5235240
## 115 115 2 9 5.70 0.6927566 0.07807621 0.5235027
## 116 116 2 9 5.75 0.6927917 0.07809887 0.5234805
## 117 117 2 9 5.80 0.6928269 0.07812192 0.5234572
## 118 118 2 9 5.85 0.6928619 0.07814539 0.5234331
## 119 119 2 9 5.90 0.6928970 0.07816925 0.5234079
## 120 120 2 9 5.95 0.6929320 0.07819352 0.5233818
## 121 121 2 9 6.00 0.6929670 0.07821818 0.5233547
## 122 122 2 9 6.05 0.6930019 0.07824324 0.5233267
## 123 123 2 9 6.10 0.6930368 0.07826868 0.5232977
## 124 124 2 9 6.15 0.6930716 0.07829452 0.5232679
## 125 125 2 9 6.20 0.6931065 0.07832074 0.5232371
## 126 126 2 9 6.25 0.6931412 0.07834734 0.5232053
## 127 127 2 9 6.30 0.6931760 0.07837432 0.5231727
## 128 128 2 9 6.35 0.6932107 0.07840168 0.5231391
## 129 129 2 9 6.40 0.6932454 0.07842941 0.5231047
## 130 130 2 9 6.45 0.6932800 0.07845751 0.5230693
## 131 131 2 9 6.50 0.6933146 0.07848598 0.5230331
## 132 132 2 9 6.55 0.6933491 0.07851482 0.5229960
## 133 133 2 9 6.60 0.6933837 0.07854402 0.5229580
## 134 134 2 9 6.65 0.6934181 0.07857358 0.5229192
## 135 135 2 9 6.70 0.6934526 0.07860349 0.5228794
## 136 136 2 9 6.75 0.6934870 0.07863376 0.5228389
## 137 137 2 9 6.80 0.6935214 0.07866438 0.5227974
## 138 138 2 9 6.85 0.6935557 0.07869535 0.5227552
## 139 139 2 9 6.90 0.6935900 0.07872667 0.5227120
## 140 140 2 9 6.95 0.6936242 0.07875833 0.5226681
## 141 141 2 9 7.00 0.6936584 0.07879033 0.5226233
## 142 142 2 9 7.05 0.6936926 0.07882267 0.5225777
## 143 143 2 9 7.10 0.6937268 0.07885534 0.5225313
## 144 144 2 9 7.15 0.6937608 0.07888834 0.5224841
## 145 145 2 9 7.20 0.6937949 0.07892167 0.5224361
## 146 146 2 9 7.25 0.6938289 0.07895533 0.5223873
## 147 147 2 9 7.30 0.6938629 0.07898931 0.5223377
## 148 148 2 9 7.35 0.6938969 0.07902362 0.5222874
## 149 149 2 9 7.40 0.6939308 0.07905824 0.5222362
## 150 150 2 9 7.45 0.6939646 0.07909317 0.5221843
## 151 151 2 9 7.50 0.6939985 0.07912842 0.5221316
## 152 152 2 9 7.55 0.6940323 0.07916398 0.5220782
## 153 153 2 9 7.60 0.6940660 0.07919985 0.5220240
## 154 154 2 9 7.65 0.6940997 0.07923602 0.5219691
## 155 155 2 9 7.70 0.6941334 0.07927249 0.5219134
## 156 156 2 9 7.75 0.6941670 0.07930926 0.5218570
## 157 157 2 9 7.80 0.6942006 0.07934632 0.5217999
## 158 158 2 9 7.85 0.6942342 0.07938368 0.5217420
## 159 159 2 9 7.90 0.6942677 0.07942133 0.5216834
## 160 160 2 9 7.95 0.6943012 0.07945927 0.5216242
## 161 161 2 9 8.00 0.6943346 0.07949749 0.5215642
## 162 162 2 9 8.05 0.6943680 0.07953600 0.5215035
## 163 163 2 9 8.10 0.6944014 0.07957478 0.5214422
## 164 164 2 9 8.15 0.6944347 0.07961384 0.5213801
## 165 165 2 9 8.20 0.6944680 0.07965318 0.5213174
## 166 166 2 9 8.25 0.6945013 0.07969279 0.5212540
## 167 167 2 9 8.30 0.6945345 0.07973266 0.5211899
## 168 168 2 9 8.35 0.6945676 0.07977281 0.5211252
## 169 169 2 9 8.40 0.6946008 0.07981321 0.5210599
## 170 170 2 9 8.45 0.6946339 0.07985388 0.5209938
## 171 171 2 9 8.50 0.6946669 0.07989481 0.5209272
## 172 172 2 9 8.55 0.6946999 0.07993599 0.5208599
## 173 173 2 9 8.60 0.6947329 0.07997742 0.5207919
## 174 174 2 9 8.65 0.6947658 0.08001911 0.5207234
## 175 175 2 9 8.70 0.6947987 0.08006104 0.5206542
## 176 176 2 9 8.75 0.6948316 0.08010322 0.5205844
## 177 177 2 9 8.80 0.6948644 0.08014564 0.5205141
## 178 178 2 9 8.85 0.6948972 0.08018830 0.5204431
## 179 179 2 9 8.90 0.6949299 0.08023120 0.5203715
## 180 180 2 9 8.95 0.6949626 0.08027434 0.5202993
## 181 181 2 9 9.00 0.6949953 0.08031770 0.5202266
## 182 182 2 9 9.05 0.6950279 0.08036130 0.5201532
## 183 183 2 9 9.10 0.6950605 0.08040512 0.5200793
## 184 184 2 9 9.15 0.6950930 0.08044917 0.5200049
## 185 185 2 9 9.20 0.6951255 0.08049344 0.5199298
## 186 186 2 9 9.25 0.6951580 0.08053793 0.5198543
## 187 187 2 9 9.30 0.6951904 0.08058264 0.5197781
## 188 188 2 9 9.35 0.6952228 0.08062756 0.5197015
## 189 189 2 9 9.40 0.6952552 0.08067269 0.5196243
## 190 190 2 9 9.45 0.6952875 0.08071804 0.5195465
## 191 191 2 9 9.50 0.6953197 0.08076359 0.5194683
## 192 192 2 9 9.55 0.6953520 0.08080934 0.5193895
## 193 193 2 9 9.60 0.6953841 0.08085530 0.5193102
## 194 194 2 9 9.65 0.6954163 0.08090146 0.5192304
## 195 195 2 9 9.70 0.6954484 0.08094781 0.5191500
## 196 196 2 9 9.75 0.6954805 0.08099436 0.5190692
## 197 197 2 9 9.80 0.6955125 0.08104110 0.5189879
## 198 198 2 9 9.85 0.6955445 0.08108804 0.5189061
## 199 199 2 9 9.90 0.6955764 0.08113516 0.5188239
## 200 200 2 9 9.95 0.6956084 0.08118246 0.5187411
## 201 201 2 9 10.00 0.6956402 0.08122995 0.5186579
## 202 202 2 9 10.05 0.6956721 0.08127762 0.5185742
## 203 203 2 9 10.10 0.6957038 0.08132547 0.5184901
## 204 204 2 9 10.15 0.6957356 0.08137349 0.5184055
## 205 205 2 9 10.20 0.6957673 0.08142169 0.5183204
## 206 206 2 9 10.25 0.6957990 0.08147005 0.5182349
## 207 207 2 9 10.30 0.6958306 0.08151859 0.5181490
## 208 208 2 9 10.35 0.6958622 0.08156729 0.5180627
## 209 209 2 9 10.40 0.6958938 0.08161616 0.5179759
## 210 210 2 9 10.45 0.6959253 0.08166519 0.5178887
## 211 211 2 9 10.50 0.6959567 0.08171438 0.5178010
## 212 212 2 9 10.55 0.6959882 0.08176372 0.5177130
## 213 213 2 9 10.60 0.6960196 0.08181322 0.5176245
## 214 214 2 9 10.65 0.6960509 0.08186288 0.5175357
## 215 215 2 9 10.70 0.6960822 0.08191268 0.5174464
## 216 216 2 9 10.75 0.6961135 0.08196264 0.5173568
## 217 217 2 9 10.80 0.6961447 0.08201274 0.5172668
## 218 218 2 9 10.85 0.6961759 0.08206298 0.5171764
## 219 219 2 9 10.90 0.6962071 0.08211336 0.5170856
## 220 220 2 9 10.95 0.6962382 0.08216389 0.5169944
## 221 221 2 9 11.00 0.6962693 0.08221455 0.5169029
## 222 222 2 9 11.05 0.6963003 0.08226535 0.5168110
## 223 223 2 9 11.10 0.6963313 0.08231628 0.5167188
## 224 224 2 9 11.15 0.6963623 0.08236734 0.5166262
## 225 225 2 9 11.20 0.6963932 0.08241853 0.5165332
## 226 226 2 9 11.25 0.6964241 0.08246985 0.5164400
## 227 227 2 9 11.30 0.6964549 0.08252129 0.5163463
## 228 228 2 9 11.35 0.6964857 0.08257285 0.5162524
## 229 229 2 9 11.40 0.6965164 0.08262454 0.5161581
## 230 230 2 9 11.45 0.6965472 0.08267634 0.5160635
## 231 231 2 9 11.50 0.6965778 0.08272826 0.5159686
## 232 232 2 9 11.55 0.6966085 0.08278029 0.5158733
## 233 233 2 9 11.60 0.6966391 0.08283243 0.5157778
## 234 234 2 9 11.65 0.6966696 0.08288469 0.5156820
## 235 235 2 9 11.70 0.6967001 0.08293705 0.5155858
## 236 236 2 9 11.75 0.6967306 0.08298952 0.5154894
## 237 237 2 9 11.80 0.6967610 0.08304210 0.5153926
## 238 238 2 9 11.85 0.6967914 0.08309477 0.5152956
## 239 239 2 9 11.90 0.6968218 0.08314755 0.5151983
## 240 240 2 9 11.95 0.6968521 0.08320042 0.5151008
## 241 241 2 9 12.00 0.6968824 0.08325339 0.5150029
## 242 242 2 9 12.05 0.6969126 0.08330646 0.5149048
## 243 243 2 9 12.10 0.6969428 0.08335961 0.5148064
## 244 244 2 9 12.15 0.6969729 0.08341286 0.5147078
## 245 245 2 9 12.20 0.6970030 0.08346620 0.5146089
## 246 246 2 9 12.25 0.6970331 0.08351962 0.5145098
## 247 247 2 9 12.30 0.6970632 0.08357313 0.5144104
## 248 248 2 9 12.35 0.6970931 0.08362672 0.5143108
## 249 249 2 9 12.40 0.6971231 0.08368039 0.5142109
## 250 250 2 9 12.45 0.6971530 0.08373414 0.5141108
## 251 251 2 9 12.50 0.6971829 0.08378797 0.5140105
## 252 252 2 9 12.55 0.6972127 0.08384187 0.5139099
## 253 253 2 9 12.60 0.6972425 0.08389584 0.5138092
## 254 254 2 9 12.65 0.6972722 0.08394989 0.5137082
## 255 255 2 9 12.70 0.6973020 0.08400401 0.5136070
## 256 256 2 9 12.75 0.6973316 0.08405820 0.5135056
## 257 257 2 9 12.80 0.6973613 0.08411245 0.5134040
## 258 258 2 9 12.85 0.6973908 0.08416677 0.5133022
## 259 259 2 9 12.90 0.6974204 0.08422115 0.5132002
## 260 260 2 9 12.95 0.6974499 0.08427559 0.5130980
## 261 261 2 9 13.00 0.6974794 0.08433009 0.5129956
## 262 262 2 9 13.05 0.6975088 0.08438465 0.5128931
## 263 263 2 9 13.10 0.6975382 0.08443927 0.5127904
## 264 264 2 9 13.15 0.6975675 0.08449394 0.5126874
## 265 265 2 9 13.20 0.6975969 0.08454866 0.5125844
## 266 266 2 9 13.25 0.6976261 0.08460343 0.5124811
## 267 267 2 9 13.30 0.6976553 0.08465826 0.5123777
## 268 268 2 9 13.35 0.6976845 0.08471313 0.5122741
## 269 269 2 9 13.40 0.6977137 0.08476804 0.5121704
## 270 270 2 9 13.45 0.6977428 0.08482300 0.5120665
## 271 271 2 9 13.50 0.6977719 0.08487801 0.5119625
## 272 272 2 9 13.55 0.6978009 0.08493305 0.5118583
## 273 273 2 9 13.60 0.6978299 0.08498814 0.5117540
## 274 274 2 9 13.65 0.6978588 0.08504326 0.5116496
## 275 275 2 9 13.70 0.6978877 0.08509842 0.5115450
## 276 276 2 9 13.75 0.6979166 0.08515361 0.5114403
## 277 277 2 9 13.80 0.6979454 0.08520884 0.5113354
## 278 278 2 9 13.85 0.6979742 0.08526410 0.5112305
## 279 279 2 9 13.90 0.6980029 0.08531939 0.5111254
## 280 280 2 9 13.95 0.6980316 0.08537471 0.5110202
## 281 281 2 9 14.00 0.6980603 0.08543005 0.5109149
## 282 282 2 9 14.05 0.6980889 0.08548542 0.5108095
## 283 283 2 9 14.10 0.6981175 0.08554082 0.5107040
## 284 284 2 9 14.15 0.6981461 0.08559623 0.5105984
## 285 285 2 9 14.20 0.6981746 0.08565167 0.5104926
## 286 286 2 9 14.25 0.6982030 0.08570713 0.5103868
## 287 287 2 9 14.30 0.6982315 0.08576261 0.5102809
## 288 288 2 9 14.35 0.6982598 0.08581810 0.5101749
## 289 289 2 9 14.40 0.6982882 0.08587361 0.5100689
## 290 290 2 9 14.45 0.6983165 0.08592914 0.5099627
## 291 291 2 9 14.50 0.6983447 0.08598467 0.5098565
## 292 292 2 9 14.55 0.6983730 0.08604022 0.5097502
## 293 293 2 9 14.60 0.6984012 0.08609578 0.5096438
## 294 294 2 9 14.65 0.6984293 0.08615134 0.5095374
## 295 295 2 9 14.70 0.6984574 0.08620692 0.5094309
## 296 296 2 9 14.75 0.6984855 0.08626250 0.5093243
## 297 297 2 9 14.80 0.6985135 0.08631808 0.5092177
## 298 298 2 9 14.85 0.6985415 0.08637367 0.5091110
## 299 299 2 9 14.90 0.6985694 0.08642926 0.5090042
## 300 300 2 9 14.95 0.6985973 0.08648485 0.5088974
## 301 301 2 9 15.00 0.6986252 0.08654043 0.5087906
## 302 302 2 9 15.05 0.6986530 0.08659602 0.5086837
## 303 303 2 9 15.10 0.6986808 0.08665160 0.5085768
## 304 304 2 9 15.15 0.6987085 0.08670718 0.5084699
## 305 305 2 9 15.20 0.6987362 0.08676275 0.5083629
## 306 306 2 9 15.25 0.6987639 0.08681832 0.5082559
## 307 307 2 9 15.30 0.6987915 0.08687388 0.5081488
## 308 308 2 9 15.35 0.6988191 0.08692943 0.5080418
## 309 309 2 9 15.40 0.6988466 0.08698496 0.5079347
## 310 310 2 9 15.45 0.6988741 0.08704049 0.5078276
## 311 311 2 9 15.50 0.6989016 0.08709600 0.5077204
## 312 312 2 9 15.55 0.6989290 0.08715150 0.5076133
## 313 313 2 9 15.60 0.6989564 0.08720698 0.5075062
## 314 314 2 9 15.65 0.6989837 0.08726244 0.5073990
## 315 315 2 9 15.70 0.6990110 0.08731789 0.5072918
## 316 316 2 9 15.75 0.6990383 0.08737332 0.5071847
## 317 317 2 9 15.80 0.6990655 0.08742873 0.5070775
## 318 318 2 9 15.85 0.6990927 0.08748411 0.5069703
## 319 319 2 9 15.90 0.6991198 0.08753948 0.5068632
## 320 320 2 9 15.95 0.6991469 0.08759482 0.5067560
## 321 321 2 9 16.00 0.6991740 0.08765013 0.5066489
## 322 322 2 9 16.05 0.6992010 0.08770542 0.5065418
## 323 323 2 9 16.10 0.6992280 0.08776069 0.5064347
## 324 324 2 9 16.15 0.6992549 0.08781592 0.5063276
## 325 325 2 9 16.20 0.6992818 0.08787113 0.5062205
## 326 326 2 9 16.25 0.6993087 0.08792630 0.5061134
## 327 327 2 9 16.30 0.6993355 0.08798145 0.5060064
## 328 328 2 9 16.35 0.6993623 0.08803656 0.5058994
## 329 329 2 9 16.40 0.6993891 0.08809164 0.5057925
## 330 330 2 9 16.45 0.6994158 0.08814669 0.5056855
## 331 331 2 9 16.50 0.6994424 0.08820170 0.5055786
## 332 332 2 9 16.55 0.6994690 0.08825667 0.5054718
## 333 333 2 9 16.60 0.6994956 0.08831161 0.5053650
## 334 334 2 9 16.65 0.6995222 0.08836651 0.5052582
## 335 335 2 9 16.70 0.6995487 0.08842137 0.5051515
## 336 336 2 9 16.75 0.6995751 0.08847619 0.5050448
## 337 337 2 9 16.80 0.6996016 0.08853097 0.5049382
## 338 338 2 9 16.85 0.6996279 0.08858571 0.5048316
## 339 339 2 9 16.90 0.6996543 0.08864040 0.5047251
## 340 340 2 9 16.95 0.6996806 0.08869505 0.5046186
## 341 341 2 9 17.00 0.6997069 0.08874966 0.5045122
## 342 342 2 9 17.05 0.6997331 0.08880422 0.5044059
## 343 343 2 9 17.10 0.6997593 0.08885873 0.5042996
## 344 344 2 9 17.15 0.6997854 0.08891320 0.5041934
## 345 345 2 9 17.20 0.6998115 0.08896762 0.5040872
## 346 346 2 9 17.25 0.6998376 0.08902199 0.5039811
## 347 347 2 9 17.30 0.6998636 0.08907631 0.5038751
## 348 348 2 9 17.35 0.6998896 0.08913057 0.5037692
## 349 349 2 9 17.40 0.6999156 0.08918479 0.5036633
## 350 350 2 9 17.45 0.6999415 0.08923895 0.5035575
## 351 351 2 9 17.50 0.6999673 0.08929306 0.5034518
## 352 352 2 9 17.55 0.6999932 0.08934712 0.5033462
## 353 353 2 9 17.60 0.7000190 0.08940112 0.5032407
## 354 354 2 9 17.65 0.7000447 0.08945506 0.5031352
## 355 355 2 9 17.70 0.7000704 0.08950895 0.5030299
## 356 356 2 9 17.75 0.7000961 0.08956278 0.5029246
## 357 357 2 9 17.80 0.7001217 0.08961655 0.5028194
## 358 358 2 9 17.85 0.7001473 0.08967027 0.5027143
## 359 359 2 9 17.90 0.7001729 0.08972392 0.5026093
## 360 360 2 9 17.95 0.7001984 0.08977751 0.5025044
## 361 361 2 9 18.00 0.7002239 0.08983104 0.5023996
## ucl fixed
## 1 0.8195849
## 2 0.8195653
## 3 0.8195465
## 4 0.8195285
## 5 0.8195113
## 6 0.8194950
## 7 0.8194795
## 8 0.8194648
## 9 0.8194511
## 10 0.8194382
## 11 0.8194261
## 12 0.8194148
## 13 0.8194045
## 14 0.8193949
## 15 0.8193861
## 16 0.8193783
## 17 0.8193712
## 18 0.8193650
## 19 0.8193595
## 20 0.8193549
## 21 0.8193512
## 22 0.8193482
## 23 0.8193462
## 24 0.8193448
## 25 0.8193444
## 26 0.8193447
## 27 0.8193459
## 28 0.8193478
## 29 0.8193506
## 30 0.8193542
## 31 0.8193586
## 32 0.8193638
## 33 0.8193699
## 34 0.8193767
## 35 0.8193842
## 36 0.8193926
## 37 0.8194018
## 38 0.8194117
## 39 0.8194224
## 40 0.8194339
## 41 0.8194464
## 42 0.8194594
## 43 0.8194733
## 44 0.8194879
## 45 0.8195033
## 46 0.8195194
## 47 0.8195364
## 48 0.8195540
## 49 0.8195725
## 50 0.8195917
## 51 0.8196117
## 52 0.8196324
## 53 0.8196538
## 54 0.8196762
## 55 0.8196991
## 56 0.8197227
## 57 0.8197471
## 58 0.8197724
## 59 0.8197982
## 60 0.8198248
## 61 0.8198521
## 62 0.8198800
## 63 0.8199088
## 64 0.8199384
## 65 0.8199685
## 66 0.8199993
## 67 0.8200309
## 68 0.8200632
## 69 0.8200961
## 70 0.8201298
## 71 0.8201641
## 72 0.8201992
## 73 0.8202349
## 74 0.8202714
## 75 0.8203084
## 76 0.8203461
## 77 0.8203844
## 78 0.8204234
## 79 0.8204632
## 80 0.8205036
## 81 0.8205447
## 82 0.8205863
## 83 0.8206286
## 84 0.8206716
## 85 0.8207151
## 86 0.8207593
## 87 0.8208041
## 88 0.8208496
## 89 0.8208956
## 90 0.8209423
## 91 0.8209897
## 92 0.8210376
## 93 0.8210861
## 94 0.8211353
## 95 0.8211849
## 96 0.8212353
## 97 0.8212862
## 98 0.8213376
## 99 0.8213897
## 100 0.8214423
## 101 0.8214956
## 102 0.8215494
## 103 0.8216036
## 104 0.8216586
## 105 0.8217139
## 106 0.8217700
## 107 0.8218266
## 108 0.8218837
## 109 0.8219413
## 110 0.8219994
## 111 0.8220582
## 112 0.8221173
## 113 0.8221771
## 114 0.8222374
## 115 0.8222982
## 116 0.8223595
## 117 0.8224213
## 118 0.8224837
## 119 0.8225464
## 120 0.8226098
## 121 0.8226735
## 122 0.8227379
## 123 0.8228026
## 124 0.8228678
## 125 0.8229335
## 126 0.8229996
## 127 0.8230663
## 128 0.8231335
## 129 0.8232010
## 130 0.8232690
## 131 0.8233375
## 132 0.8234064
## 133 0.8234758
## 134 0.8235455
## 135 0.8236157
## 136 0.8236865
## 137 0.8237575
## 138 0.8238289
## 139 0.8239008
## 140 0.8239731
## 141 0.8240459
## 142 0.8241190
## 143 0.8241926
## 144 0.8242665
## 145 0.8243408
## 146 0.8244155
## 147 0.8244904
## 148 0.8245660
## 149 0.8246417
## 150 0.8247179
## 151 0.8247945
## 152 0.8248714
## 153 0.8249487
## 154 0.8250264
## 155 0.8251043
## 156 0.8251826
## 157 0.8252614
## 158 0.8253404
## 159 0.8254197
## 160 0.8254993
## 161 0.8255794
## 162 0.8256598
## 163 0.8257404
## 164 0.8258214
## 165 0.8259026
## 166 0.8259841
## 167 0.8260660
## 168 0.8261482
## 169 0.8262306
## 170 0.8263135
## 171 0.8263965
## 172 0.8264799
## 173 0.8265635
## 174 0.8266473
## 175 0.8267316
## 176 0.8268159
## 177 0.8269007
## 178 0.8269855
## 179 0.8270708
## 180 0.8271563
## 181 0.8272419
## 182 0.8273279
## 183 0.8274141
## 184 0.8275006
## 185 0.8275872
## 186 0.8276741
## 187 0.8277612
## 188 0.8278486
## 189 0.8279362
## 190 0.8280240
## 191 0.8281120
## 192 0.8282002
## 193 0.8282887
## 194 0.8283773
## 195 0.8284661
## 196 0.8285552
## 197 0.8286444
## 198 0.8287337
## 199 0.8288234
## 200 0.8289131
## 201 0.8290030
## 202 0.8290932
## 203 0.8291834
## 204 0.8292739
## 205 0.8293645
## 206 0.8294553
## 207 0.8295463
## 208 0.8296375
## 209 0.8297286
## 210 0.8298200
## 211 0.8299116
## 212 0.8300033
## 213 0.8300952
## 214 0.8301871
## 215 0.8302792
## 216 0.8303714
## 217 0.8304638
## 218 0.8305563
## 219 0.8306489
## 220 0.8307416
## 221 0.8308345
## 222 0.8309274
## 223 0.8310205
## 224 0.8311136
## 225 0.8312069
## 226 0.8313002
## 227 0.8313937
## 228 0.8314872
## 229 0.8315809
## 230 0.8316746
## 231 0.8317685
## 232 0.8318624
## 233 0.8319563
## 234 0.8320504
## 235 0.8321446
## 236 0.8322387
## 237 0.8323330
## 238 0.8324274
## 239 0.8325218
## 240 0.8326161
## 241 0.8327108
## 242 0.8328053
## 243 0.8328998
## 244 0.8329946
## 245 0.8330893
## 246 0.8331841
## 247 0.8332789
## 248 0.8333738
## 249 0.8334686
## 250 0.8335635
## 251 0.8336585
## 252 0.8337535
## 253 0.8338484
## 254 0.8339436
## 255 0.8340385
## 256 0.8341336
## 257 0.8342288
## 258 0.8343239
## 259 0.8344190
## 260 0.8345141
## 261 0.8346093
## 262 0.8347044
## 263 0.8347995
## 264 0.8348947
## 265 0.8349899
## 266 0.8350850
## 267 0.8351801
## 268 0.8352753
## 269 0.8353704
## 270 0.8354655
## 271 0.8355606
## 272 0.8356557
## 273 0.8357507
## 274 0.8358458
## 275 0.8359409
## 276 0.8360359
## 277 0.8361309
## 278 0.8362258
## 279 0.8363207
## 280 0.8364155
## 281 0.8365105
## 282 0.8366052
## 283 0.8367001
## 284 0.8367947
## 285 0.8368895
## 286 0.8369842
## 287 0.8370787
## 288 0.8371733
## 289 0.8372678
## 290 0.8373623
## 291 0.8374567
## 292 0.8375510
## 293 0.8376454
## 294 0.8377396
## 295 0.8378338
## 296 0.8379279
## 297 0.8380220
## 298 0.8381159
## 299 0.8382099
## 300 0.8383038
## 301 0.8383975
## 302 0.8384913
## 303 0.8385850
## 304 0.8386785
## 305 0.8387719
## 306 0.8388654
## 307 0.8389588
## 308 0.8390521
## 309 0.8391452
## 310 0.8392384
## 311 0.8393313
## 312 0.8394243
## 313 0.8395170
## 314 0.8396098
## 315 0.8397025
## 316 0.8397951
## 317 0.8398876
## 318 0.8399800
## 319 0.8400723
## 320 0.8401644
## 321 0.8402565
## 322 0.8403486
## 323 0.8404405
## 324 0.8405322
## 325 0.8406240
## 326 0.8407154
## 327 0.8408070
## 328 0.8408983
## 329 0.8409896
## 330 0.8410807
## 331 0.8411718
## 332 0.8412627
## 333 0.8413536
## 334 0.8414442
## 335 0.8415349
## 336 0.8416253
## 337 0.8417156
## 338 0.8418059
## 339 0.8418959
## 340 0.8419860
## 341 0.8420758
## 342 0.8421656
## 343 0.8422553
## 344 0.8423448
## 345 0.8424341
## 346 0.8425233
## 347 0.8426125
## 348 0.8427016
## 349 0.8427904
## 350 0.8428791
## 351 0.8429677
## 352 0.8430562
## 353 0.8431445
## 354 0.8432328
## 355 0.8433208
## 356 0.8434088
## 357 0.8434965
## 358 0.8435842
## 359 0.8436717
## 360 0.8437591
## 361 0.8438464
ggplot(data=psi.ma$estimates, aes(x=covdata, y=estimate))+
geom_point()+
geom_ribbon(aes(ymin=lcl, ymax=ucl), alpha=0.2)+
labs(x = "Sagebrush", y = "Estimated Occupancy")
# covariate predictions for categorical covariates
# such as aspect
pronghorn.ddl$Psi
## par.index model.index group age time Age Time aspect
## 1 1 9 E 0 1 0 0 E
## 2 2 10 N 0 1 0 0 N
## 3 3 11 S 0 1 0 0 S
## 4 4 12 W 0 1 0 0 W
psi.ma <- model.average(model.set, param="Psi", vcv=TRUE)
names(psi.ma)
## [1] "estimates" "vcv.real"
psi.ma$estimates
## par.index estimate se lcl ucl fixed note
## Psi gE a0 t1 9 0.6916572 0.07759342 0.5236778 0.8206820
## Psi gN a0 t1 10 0.7079212 0.08462465 0.5207758 0.8438908
## Psi gS a0 t1 11 0.6995434 0.07033702 0.5471662 0.8177267
## Psi gW a0 t1 12 0.6984185 0.06751821 0.5526800 0.8127610
## group age time Age Time aspect
## Psi gE a0 t1 E 0 1 0 0 E
## Psi gN a0 t1 N 0 1 0 0 N
## Psi gS a0 t1 S 0 1 0 0 S
## Psi gW a0 t1 W 0 1 0 0 W
ggplot(data=psi.ma$estimates, aes(x=aspect, y=estimate))+
geom_point()+
geom_errorbar(aes(ymin=lcl,
ymax=ucl), width=0.2)+
ylim(0,1)+
labs(x = "Aspect", y = "Estimated Occupancy")
# or use covariate predictions
psi.ma2 <- covariate.predictions(model.set, indices=9:12)
psi.ma2$estimates <- cbind(psi.ma2$estimates, pronghorn.ddl$Psi[,"aspect", drop=FALSE])
psi.ma2$estimates
## vcv.index model.index par.index index estimate se lcl
## 1 2 2 9 9 0.6916572 0.07759342 0.5236778
## 2 2 2 10 10 0.7079212 0.08462465 0.5207758
## 3 2 2 11 11 0.6995434 0.07033702 0.5471662
## 4 2 2 12 12 0.6984185 0.06751821 0.5526800
## ucl fixed aspect
## 1 0.8206820 E
## 2 0.8438908 N
## 3 0.8177267 S
## 4 0.8127610 W
ggplot(data=psi.ma2$estimates, aes(x=aspect, y=estimate))+
geom_point()+
geom_errorbar(aes(ymin=lcl,
ymax=ucl), width=0.2)+
ylim(0,1)+
labs(x = "Aspect", y = "Estimated Occupancy")
# 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 gE a0 t1 9 0.6916572 0.07759342 0.5236778 0.8206820
## Psi gN a0 t1 10 0.7079212 0.08462465 0.5207758 0.8438908
## Psi gS a0 t1 11 0.6995434 0.07033702 0.5471662 0.8177267
## Psi gW a0 t1 12 0.6984185 0.06751821 0.5526800 0.8127610
## group age time Age Time aspect
## Psi gE a0 t1 E 0 1 0 0 E
## Psi gN a0 t1 N 0 1 0 0 N
## Psi gS a0 t1 S 0 1 0 0 S
## Psi gW a0 t1 W 0 1 0 0 W
##
## $vcv.real
## 9 10 11 12
## 9 0.006020739 0.006313922 0.005246962 0.005052484
## 10 0.006313922 0.007161332 0.005740953 0.005536641
## 11 0.005246962 0.005740953 0.004947297 0.004600086
## 12 0.005052484 0.005536641 0.004600086 0.004558708
psi.ma2$estimates
## vcv.index model.index par.index index estimate se lcl
## 1 2 2 9 9 0.6916572 0.07759342 0.5236778
## 2 2 2 10 10 0.7079212 0.08462465 0.5207758
## 3 2 2 11 11 0.6995434 0.07033702 0.5471662
## 4 2 2 12 12 0.6984185 0.06751821 0.5526800
## ucl fixed aspect
## 1 0.8206820 E
## 2 0.8438908 N
## 3 0.8177267 S
## 4 0.8127610 W
# 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 measure 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
pronghorn.ddl$Psi
## par.index model.index group age time Age Time aspect
## 1 1 9 E 0 1 0 0 E
## 2 2 10 N 0 1 0 0 N
## 3 3 11 S 0 1 0 0 S
## 4 4 12 W 0 1 0 0 W
site.covariates <- merge(pronghorn.ddl$Psi, site.covariates)
head(site.covariates)
## aspect par.index model.index group age time Age Time sagebrush
## 1 E 1 9 E 0 1 0 0 6.1
## 2 E 1 9 E 0 1 0 0 4.8
## 3 E 1 9 E 0 1 0 0 3.0
## 4 E 1 9 E 0 1 0 0 3.2
## 5 E 1 9 E 0 1 0 0 3.2
## 6 E 1 9 E 0 1 0 0 1.8
# 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=sagebrush, 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(~aspect, ncol=1)
#cleanup
cleanup(ask=FALSE)