# Redbreast sunfish occupancy of streams segments over 8 seasons = 4 years * summer*spring
# 3 quadrats were sampled per segment via electrofishing total 24 sampling occasions
# 15 covariates: LN(link magnitude) and standardized Q (streamflow) for the sampling intervals;
# minimum discharge for intervals 1-7, maximum discharge intervals 1-7
# Fitting several models using the RMark package
# Single Species Multiple Seasons
# 2018-11-26 code submitted by Neil Faught
library(readxl)
library(RMark)
## This is RMark 2.2.5
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
library(ggplot2)
# Get the RMark additional functions
source(file.path("..","..","..","AdditionalFunctions","RMark.additional.functions.R"))
# 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 <- readxl::read_excel(file.path("..","sunfish.xls"),
skip=1, col_names=TRUE)#No NAs in this data set, make sure to use "na =" statement if there are
head(input.data)
## # A tibble: 6 x 39
## S11 S12 S13 S21 S22 S23 S31 S32 S33 S41 S42 S43
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 1 1 1 1 1 1 1 1 1 1
## 2 0 1 1 1 0 1 0 0 0 0 1 0
## 3 0 0 0 1 1 1 0 0 0 1 1 1
## 4 1 1 1 1 1 1 1 1 1 1 1 1
## 5 1 1 1 0 0 0 1 1 1 0 0 0
## 6 1 1 1 0 0 0 0 0 0 1 1 0
## # ... with 27 more variables: S51 <dbl>, S52 <dbl>, S53 <dbl>, S61 <dbl>,
## # S62 <dbl>, S63 <dbl>, S71 <dbl>, S72 <dbl>, S73 <dbl>, S81 <dbl>,
## # S82 <dbl>, S83 <dbl>, LN <dbl>, minQ1 <dbl>, minQ2 <dbl>, minQ3 <dbl>,
## # minQ4 <dbl>, minQ5 <dbl>, minQ6 <dbl>, minQ7 <dbl>, maxQ1 <dbl>,
## # maxQ2 <dbl>, maxQ3 <dbl>, maxQ4 <dbl>, maxQ5 <dbl>, maxQ6 <dbl>,
## # maxQ7 <dbl>
input.history <- input.data[, 1:24]
head(input.history)
## # A tibble: 6 x 24
## S11 S12 S13 S21 S22 S23 S31 S32 S33 S41 S42 S43
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 1 1 1 1 1 1 1 1 1 1
## 2 0 1 1 1 0 1 0 0 0 0 1 0
## 3 0 0 0 1 1 1 0 0 0 1 1 1
## 4 1 1 1 1 1 1 1 1 1 1 1 1
## 5 1 1 1 0 0 0 1 1 1 0 0 0
## 6 1 1 1 0 0 0 0 0 0 1 1 0
## # ... with 12 more variables: S51 <dbl>, S52 <dbl>, S53 <dbl>, S61 <dbl>,
## # S62 <dbl>, S63 <dbl>, S71 <dbl>, S72 <dbl>, S73 <dbl>, S81 <dbl>,
## # S82 <dbl>, S83 <dbl>
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 26
ncol(input.history)
## [1] 24
range(input.history, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.history))
## [1] 0
site.covar <- data.frame(Site=1:nrow(input.data), input.data[,25:ncol(input.data)])
head(site.covar)
## Site LN minQ1 minQ2 minQ3 minQ4 minQ5 minQ6 minQ7 maxQ1 maxQ2 maxQ3
## 1 1 1.32 0.61 0.16 0.64 0.92 1.00 0.10 0.27 5.59 7.58 7.03
## 2 2 4.97 0.45 0.06 0.93 0.34 0.20 0.67 0.77 6.63 1.47 3.49
## 3 3 1.02 0.12 0.21 0.03 0.26 0.99 0.81 0.47 3.19 6.73 6.35
## 4 4 2.20 0.32 0.57 0.67 0.37 0.40 0.87 0.30 6.40 5.92 2.39
## 5 5 2.89 0.22 0.85 0.02 0.78 0.47 0.43 0.45 6.36 6.87 4.86
## 6 6 2.87 0.55 0.42 0.45 0.99 0.85 0.50 0.64 4.97 1.53 3.12
## maxQ4 maxQ5 maxQ6 maxQ7
## 1 1.22 6.73 4.77 1.87
## 2 3.36 2.64 3.21 0.86
## 3 1.46 4.71 1.32 1.01
## 4 0.81 5.06 5.05 5.63
## 5 7.44 5.63 4.13 7.58
## 6 2.14 0.91 1.07 2.70
#Format the capture history to be used by RMark
input.history <- data.frame(freq=1,
ch=apply(input.history,1,paste, collapse=""), stringsAsFactors=FALSE)
head(input.history)
## freq ch
## 1 1 101111111111011010000111
## 2 1 011101000010110000110101
## 3 1 000111000111111111111111
## 4 1 111111111111101111111111
## 5 1 111000111000110011010111
## 6 1 111000000110101111111111
# Change any NA to . in the chapter history
input.history$ch <- gsub("NA",".", input.history$ch, fixed=TRUE)
head(input.history)
## freq ch
## 1 1 101111111111011010000111
## 2 1 011101000010110000110101
## 3 1 000111000111111111111111
## 4 1 111111111111101111111111
## 5 1 111000111000110011010111
## 6 1 111000000110101111111111
#Add site covariates to input history
input.history = cbind(input.history,site.covar)
head(input.history)
## freq ch Site LN minQ1 minQ2 minQ3 minQ4 minQ5
## 1 1 101111111111011010000111 1 1.32 0.61 0.16 0.64 0.92 1.00
## 2 1 011101000010110000110101 2 4.97 0.45 0.06 0.93 0.34 0.20
## 3 1 000111000111111111111111 3 1.02 0.12 0.21 0.03 0.26 0.99
## 4 1 111111111111101111111111 4 2.20 0.32 0.57 0.67 0.37 0.40
## 5 1 111000111000110011010111 5 2.89 0.22 0.85 0.02 0.78 0.47
## 6 1 111000000110101111111111 6 2.87 0.55 0.42 0.45 0.99 0.85
## minQ6 minQ7 maxQ1 maxQ2 maxQ3 maxQ4 maxQ5 maxQ6 maxQ7
## 1 0.10 0.27 5.59 7.58 7.03 1.22 6.73 4.77 1.87
## 2 0.67 0.77 6.63 1.47 3.49 3.36 2.64 3.21 0.86
## 3 0.81 0.47 3.19 6.73 6.35 1.46 4.71 1.32 1.01
## 4 0.87 0.30 6.40 5.92 2.39 0.81 5.06 5.05 5.63
## 5 0.43 0.45 6.36 6.87 4.86 7.44 5.63 4.13 7.58
## 6 0.50 0.64 4.97 1.53 3.12 2.14 0.91 1.07 2.70
#Create the RMark data structure
#8 Seasons, with 3 visits per season
max.visit.per.year <- 3
n.season <- 8
sunfish.data <- process.data(data=input.history,
model="RDOccupEG",
time.intervals=c( rep( c(rep(0,max.visit.per.year-1),1),n.season-1),
rep(0,max.visit.per.year-1)))
summary(sunfish.data)
## Length Class Mode
## data 18 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 26 -none- numeric
## nocc 1 -none- numeric
## nocc.secondary 8 -none- numeric
## time.intervals 23 -none- numeric
## begin.time 1 -none- numeric
## age.unit 1 -none- numeric
## initial.ages 1 -none- numeric
## group.covariates 0 -none- NULL
## nstrata 1 -none- numeric
## strata.labels 0 -none- NULL
## counts 0 -none- NULL
## reverse 1 -none- logical
## areas 0 -none- NULL
## events 0 -none- NULL
# add visit level covariates to the ddl's in the loop below
# In this case none
# define the list of models to fit
# Notice the commas between the column and the placement of the quotes
# Define the models.
# model.type can be RDOccupPE, RDOccupPG, RDOccupEG
#The random occupancy model cannot be fit here. See the other code in this directory.
model.list.csv <- textConnection("
p, Psi, Gamma, Epsilon, model.type
~session, ~LN, ~maxQ, ~minQ, RDOccupEG
~session, ~LN, ~minQ, ~minQ, RDOccupEG
~session, ~LN, ~minQ, ~maxQ, RDOccupEG
~session, ~1, ~time, ~time, RDOccupEG")
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 Gamma Epsilon model.type model.number
## 1 ~session ~LN ~maxQ ~minQ RDOccupEG 1
## 2 ~session ~LN ~minQ ~minQ RDOccupEG 2
## 3 ~session ~LN ~minQ ~maxQ RDOccupEG 3
## 4 ~session ~1 ~time ~time RDOccupEG 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.history){
cat("\n\n***** Starting ", unlist(x), "\n")
max.visit.per.year <- 3
n.season <- 8
# we need to process the data in the loop to allow for different data types
input.data <- process.data(data=input.history, model=x$model.type,
time.intervals=c( rep( c(rep(0,max.visit.per.year-1),1),n.season-1),
rep(0,max.visit.per.year-1)))
# set up the ddls as needed for time-varying covariates.
# you need to do this in the loop because different paraemterizations have different ddl structures
input.ddl <- make.design.data(input.data)
model.parameters=list(
Psi =list(formula=as.formula(eval(x$Psi))),
p =list(formula=as.formula(eval(x$p))),
Epsilon=list(formula=as.formula(eval(x$Epsilon))),
Gamma =list(formula=as.formula(eval(x$Gamma)))
)
if(x$model.type == "RDOccupPG"){ # psi, gamma, p formulation
model.parameters$Epsilon= NULL
}
if(x$model.type == "RDOccupPE"){ # psi, epsilon, p formulation
model.parameters$Gamma = NULL
}
fit <- RMark::mark(input.data, ddl=input.ddl,
model=x$model.type,
model.parameters=model.parameters
#,brief=TRUE,output=FALSE, delete=TRUE
#,invisible=TRUE,output=TRUE # set for debugging
)
mnumber <- paste("m...",formatC(x$model.number, width = 3, format = "d", flag = "0"),sep="")
assign( mnumber, fit, envir=.GlobalEnv)
#browser()
fit
},input.history=input.history)
##
##
## ***** Starting ~session ~LN ~maxQ ~minQ RDOccupEG 1
##
## Output summary for RDOccupEG model
## Name : Psi(~LN)Epsilon(~minQ)Gamma(~maxQ)p(~session)
##
## Npar : 14
## -2lnL: 583.0108
## AICc : 613.1869
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -0.7550115 0.9144145 -2.5472639 1.037241
## Psi:LN 0.6903310 0.3996348 -0.0929532 1.473615
## Epsilon:(Intercept) 2.0191530 0.5646610 0.9124175 3.125888
## Epsilon:minQ -8.4850749 1.7377734 -11.8911110 -5.079039
## Gamma:(Intercept) -5.5623692 2.2319366 -9.9369651 -1.187773
## Gamma:maxQ 2.4416164 0.9621997 0.5557050 4.327528
## p:(Intercept) 1.0996552 0.3303472 0.4521747 1.747136
## p:session2 0.5611615 0.4789286 -0.3775385 1.499862
## p:session3 1.0994643 0.5822952 -0.0418344 2.240763
## p:session4 0.2773194 0.4635618 -0.6312617 1.185901
## p:session5 0.2583313 0.4671173 -0.6572186 1.173881
## p:session6 0.1797970 0.4562922 -0.7145358 1.074130
## p:session7 0.2379745 0.4782112 -0.6993193 1.175268
## p:session8 0.3169032 0.4736571 -0.6114647 1.245271
##
##
## Real Parameter Psi
##
## 1
## 0.7549786
##
##
## Real Parameter Epsilon
##
## 1 2 3 4 5 6 7
## 0.0888453 0.1060246 0.1054076 0.0764922 0.0599699 0.0864962 0.1735181
##
##
## Real Parameter Gamma
##
## 1 2 3 4 5 6 7
## 0.9978548 0.9988205 0.9938317 0.9866239 0.9946926 0.9903347 0.9449258
##
##
## Real Parameter p
## Session:1
## 1 2 3
## 0.7501955 0.7501955 0.7501955
##
## Session:2
## 1 2 3
## 0.8403476 0.8403476 0.8403476
##
## Session:3
## 1 2 3
## 0.9001704 0.9001704 0.9001704
##
## Session:4
## 1 2 3
## 0.7985047 0.7985047 0.7985047
##
## Session:5
## 1 2 3
## 0.7954323 0.7954323 0.7954323
##
## Session:6
## 1 2 3
## 0.7823565 0.7823565 0.7823565
##
## Session:7
## 1 2 3
## 0.7920999 0.7920999 0.7920999
##
## Session:8
## 1 2 3
## 0.8047983 0.8047983 0.8047983
##
##
## ***** Starting ~session ~LN ~minQ ~minQ RDOccupEG 2
##
## Output summary for RDOccupEG model
## Name : Psi(~LN)Epsilon(~minQ)Gamma(~minQ)p(~session)
##
## Npar : 14
## -2lnL: 609.8562
## AICc : 640.0323
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -0.7504629 0.9158700 -2.5455681 1.044642
## Psi:LN 0.6853818 0.3983787 -0.0954405 1.466204
## Epsilon:(Intercept) 1.9051751 0.5388750 0.8489801 2.961370
## Epsilon:minQ -7.9538483 1.5928501 -11.0758340 -4.831862
## Gamma:(Intercept) -0.2143369 0.6873358 -1.5615152 1.132841
## Gamma:minQ 2.8339608 1.5445840 -0.1934239 5.861345
## p:(Intercept) 1.1049816 0.3294862 0.4591886 1.750775
## p:session2 0.5409549 0.4813974 -0.4025840 1.484494
## p:session3 1.0891509 0.5836477 -0.0547986 2.233100
## p:session4 0.2620279 0.4647756 -0.6489324 1.172988
## p:session5 0.2445785 0.4679377 -0.6725794 1.161736
## p:session6 0.1697262 0.4564760 -0.7249668 1.064419
## p:session7 0.3777586 0.4915918 -0.5857613 1.341279
## p:session8 0.2878425 0.4778359 -0.6487159 1.224401
##
##
## Real Parameter Psi
##
## 1
## 0.7533225
##
##
## Real Parameter Epsilon
##
## 1 2 3 4 5 6 7
## 0.1025099 0.1206709 0.1200232 0.0892684 0.0712705 0.1000045 0.1898882
##
##
## Real Parameter Gamma
##
## 1 2 3 4 5 6 7
## 0.7751326 0.7635287 0.7639221 0.7844893 0.7988659 0.7768378 0.7273813
##
##
## Real Parameter p
## Session:1
## 1 2 3
## 0.7511923 0.7511923 0.7511923
##
## Session:2
## 1 2 3
## 0.8383411 0.8383411 0.8383411
##
## Session:3
## 1 2 3
## 0.8997214 0.8997214 0.8997214
##
## Session:4
## 1 2 3
## 0.7968966 0.7968966 0.7968966
##
## Session:5
## 1 2 3
## 0.7940577 0.7940577 0.7940577
##
## Session:6
## 1 2 3
## 0.7815476 0.7815476 0.7815476
##
## Session:7
## 1 2 3
## 0.8149861 0.8149861 0.8149861
##
## Session:8
## 1 2 3
## 0.8010427 0.8010427 0.8010427
##
##
## ***** Starting ~session ~LN ~minQ ~maxQ RDOccupEG 3
##
## Output summary for RDOccupEG model
## Name : Psi(~LN)Epsilon(~maxQ)Gamma(~minQ)p(~session)
##
## Npar : 14
## -2lnL: 665.8409
## AICc : 696.017
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -0.7561081 0.9266373 -2.5723173 1.0601010
## Psi:LN 0.6989134 0.4104989 -0.1056644 1.5034913
## Epsilon:(Intercept) -1.1340184 0.4718322 -2.0588096 -0.2092273
## Epsilon:maxQ -0.0033602 0.0961648 -0.1918431 0.1851227
## Gamma:(Intercept) -0.3244735 0.7446792 -1.7840447 1.1350978
## Gamma:minQ 3.2043795 1.9080454 -0.5353896 6.9441487
## p:(Intercept) 1.0838412 0.3354316 0.4263953 1.7412870
## p:session2 0.5641885 0.4851907 -0.3867853 1.5151623
## p:session3 1.1191819 0.5838691 -0.0252015 2.2635653
## p:session4 0.2755946 0.4703650 -0.6463209 1.1975101
## p:session5 0.2818595 0.4693053 -0.6379789 1.2016978
## p:session6 0.1647473 0.4653032 -0.7472469 1.0767415
## p:session7 0.3864820 0.5032715 -0.5999301 1.3728941
## p:session8 0.2793255 0.4878488 -0.6768581 1.2355092
##
##
## Real Parameter Psi
##
## 1
## 0.7590768
##
##
## Real Parameter Epsilon
##
## 1 2 3 4 5 6 7
## 0.2404659 0.2403154 0.2407325 0.2409291 0.2406945 0.2408464 0.2412964
##
##
## Real Parameter Gamma
##
## 1 2 3 4 5 6 7
## 0.788708 0.7761219 0.7765499 0.7987946 0.8141756 0.7905505 0.7364364
##
##
## Real Parameter p
## Session:1
## 1 2 3
## 0.7472202 0.7472202 0.7472202
##
## Session:2
## 1 2 3
## 0.8386246 0.8386246 0.8386246
##
## Session:3
## 1 2 3
## 0.9005207 0.9005207 0.9005207
##
## Session:4
## 1 2 3
## 0.795668 0.795668 0.795668
##
## Session:5
## 1 2 3
## 0.7966846 0.7966846 0.7966846
##
## Session:6
## 1 2 3
## 0.7770554 0.7770554 0.7770554
##
## Session:7
## 1 2 3
## 0.8131065 0.8131065 0.8131065
##
## Session:8
## 1 2 3
## 0.7962739 0.7962739 0.7962739
##
##
## ***** Starting ~session ~1 ~time ~time RDOccupEG 4
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~time)Gamma(~time)p(~session)
##
## Npar : 23
## -2lnL: 667.5897
## AICc : 719.5897
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 0.8650257 0.4432493 -0.0037430 1.7337944
## Epsilon:(Intercept) -1.6396588 0.6492717 -2.9122314 -0.3670862
## Epsilon:time2 1.1526698 0.7914312 -0.3985354 2.7038750
## Epsilon:time3 0.4081064 0.8821987 -1.3210031 2.1372159
## Epsilon:time4 -0.1423801 0.9211349 -1.9478045 1.6630443
## Epsilon:time5 0.2022410 0.8766405 -1.5159745 1.9204565
## Epsilon:time6 0.5205969 0.8343609 -1.1147505 2.1559443
## Epsilon:time7 0.8391983 0.8211572 -0.7702700 2.4486665
## Gamma:(Intercept) 1.0944145 0.8493713 -0.5703532 2.7591822
## Gamma:time2 0.3256226 1.4442757 -2.5051578 3.1564031
## Gamma:time3 0.1701109 1.1742270 -2.1313741 2.4715960
## Gamma:time4 -1.1157363 1.2029987 -3.4736138 1.2421411
## Gamma:time5 -0.3599647 1.2458862 -2.8019017 2.0819723
## Gamma:time6 -0.3945936 1.2465674 -2.8378658 2.0486785
## Gamma:time7 0.7678158 1.4414191 -2.0573658 3.5929973
## p:(Intercept) 1.0824523 0.3359255 0.4240382 1.7408663
## p:session2 0.5572486 0.4875211 -0.3982928 1.5127900
## p:session3 1.1269296 0.5818937 -0.0135821 2.2674413
## p:session4 0.2603815 0.4743798 -0.6694030 1.1901660
## p:session5 0.2667445 0.4729221 -0.6601828 1.1936718
## p:session6 0.1497653 0.4691444 -0.7697578 1.0692883
## p:session7 0.4339862 0.4906744 -0.5277356 1.3957080
## p:session8 0.3077299 0.4832095 -0.6393609 1.2548206
##
##
## Real Parameter Psi
##
## 1
## 0.7037096
##
##
## Real Parameter Epsilon
##
## 1 2 3 4 5 6 7
## 0.1625115 0.3806031 0.2259098 0.1440515 0.1919455 0.2461853 0.309927
##
##
## Real Parameter Gamma
##
## 1 2 3 4 5 6 7
## 0.7492121 0.8053442 0.7798042 0.4946697 0.675781 0.6681481 0.8655567
##
##
## Real Parameter p
## Session:1
## 1 2 3
## 0.7469578 0.7469578 0.7469578
##
## Session:2
## 1 2 3
## 0.8374942 0.8374942 0.8374942
##
## Session:3
## 1 2 3
## 0.9010889 0.9010889 0.9010889
##
## Session:4
## 1 2 3
## 0.7929556 0.7929556 0.7929556
##
## Session:5
## 1 2 3
## 0.7939983 0.7939983 0.7939983
##
## Session:6
## 1 2 3
## 0.7742065 0.7742065 0.7742065
##
## Session:7
## 1 2 3
## 0.8200134 0.8200134 0.8200134
##
## Session:8
## 1 2 3
## 0.8006213 0.8006213 0.8006213
# examine individual model results
model.number <-1
summary(model.fits[[model.number]])
## Output summary for RDOccupEG model
## Name : Psi(~LN)Epsilon(~minQ)Gamma(~maxQ)p(~session)
##
## Npar : 14
## -2lnL: 583.0108
## AICc : 613.1869
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) -0.7550115 0.9144145 -2.5472639 1.037241
## Psi:LN 0.6903310 0.3996348 -0.0929532 1.473615
## Epsilon:(Intercept) 2.0191530 0.5646610 0.9124175 3.125888
## Epsilon:minQ -8.4850749 1.7377734 -11.8911110 -5.079039
## Gamma:(Intercept) -5.5623692 2.2319366 -9.9369651 -1.187773
## Gamma:maxQ 2.4416164 0.9621997 0.5557050 4.327528
## p:(Intercept) 1.0996552 0.3303472 0.4521747 1.747136
## p:session2 0.5611615 0.4789286 -0.3775385 1.499862
## p:session3 1.0994643 0.5822952 -0.0418344 2.240763
## p:session4 0.2773194 0.4635618 -0.6312617 1.185901
## p:session5 0.2583313 0.4671173 -0.6572186 1.173881
## p:session6 0.1797970 0.4562922 -0.7145358 1.074130
## p:session7 0.2379745 0.4782112 -0.6993193 1.175268
## p:session8 0.3169032 0.4736571 -0.6114647 1.245271
##
##
## Real Parameter Psi
##
## 1
## 0.7549786
##
##
## Real Parameter Epsilon
##
## 1 2 3 4 5 6 7
## 0.0888453 0.1060246 0.1054076 0.0764922 0.0599699 0.0864962 0.1735181
##
##
## Real Parameter Gamma
##
## 1 2 3 4 5 6 7
## 0.9978548 0.9988205 0.9938317 0.9866239 0.9946926 0.9903347 0.9449258
##
##
## Real Parameter p
## Session:1
## 1 2 3
## 0.7501955 0.7501955 0.7501955
##
## Session:2
## 1 2 3
## 0.8403476 0.8403476 0.8403476
##
## Session:3
## 1 2 3
## 0.9001704 0.9001704 0.9001704
##
## Session:4
## 1 2 3
## 0.7985047 0.7985047 0.7985047
##
## Session:5
## 1 2 3
## 0.7954323 0.7954323 0.7954323
##
## Session:6
## 1 2 3
## 0.7823565 0.7823565 0.7823565
##
## Session:7
## 1 2 3
## 0.7920999 0.7920999 0.7920999
##
## Session:8
## 1 2 3
## 0.8047983 0.8047983 0.8047983
model.fits[[model.number]]$results$real
## estimate se lcl ucl fixed note
## Psi g1 a0 t1 0.7549786 0.1052240 0.5026132 0.9038049
## Epsilon g1 a0 t1 0.0888453 0.0389031 0.0366242 0.2000636
## Epsilon g1 a1 t2 0.1060246 0.0424746 0.0469615 0.2220626
## Epsilon g1 a2 t3 0.1054076 0.0423566 0.0465778 0.2212955
## Epsilon g1 a3 t4 0.0764922 0.0359172 0.0296687 0.1832569
## Epsilon g1 a4 t5 0.0599699 0.0312492 0.0210713 0.1590123
## Epsilon g1 a5 t6 0.0864962 0.0383645 0.0352692 0.1969393
## Epsilon g1 a6 t7 0.1735181 0.0519634 0.0935460 0.2992850
## Gamma g1 a0 t1 0.9978548 0.0053861 0.7704283 0.9999845
## Gamma g1 a1 t2 0.9988205 0.0032360 0.7954000 0.9999946
## Gamma g1 a2 t3 0.9938317 0.0129407 0.7200482 0.9999009
## Gamma g1 a3 t4 0.9866239 0.0239659 0.6773181 0.9996143
## Gamma g1 a4 t5 0.9946926 0.0114476 0.7277449 0.9999239
## Gamma g1 a5 t6 0.9903347 0.0185654 0.6959166 0.9997821
## Gamma g1 a6 t7 0.9449258 0.0667096 0.5817518 0.9952972
## p g1 s1 t1 0.7501955 0.0619078 0.6111562 0.8515912
## p g1 s2 t1 0.8403476 0.0465239 0.7273357 0.9121751
## p g1 s3 t1 0.9001704 0.0430944 0.7788875 0.9584748
## p g1 s4 t1 0.7985047 0.0523173 0.6769198 0.8822900
## p g1 s5 t1 0.7954323 0.0537163 0.6706133 0.8813224
## p g1 s6 t1 0.7823565 0.0535991 0.6598181 0.8694859
## p g1 s7 t1 0.7920999 0.0571454 0.6586947 0.8826509
## p g1 s8 t1 0.8047983 0.0533294 0.6794412 0.8891326
model.fits[[model.number]]$results$beta
## estimate se lcl ucl
## Psi:(Intercept) -0.7550115 0.9144145 -2.5472639 1.037241
## Psi:LN 0.6903310 0.3996348 -0.0929532 1.473615
## Epsilon:(Intercept) 2.0191530 0.5646610 0.9124175 3.125888
## Epsilon:minQ -8.4850749 1.7377734 -11.8911110 -5.079039
## Gamma:(Intercept) -5.5623692 2.2319366 -9.9369651 -1.187773
## Gamma:maxQ 2.4416164 0.9621997 0.5557050 4.327528
## p:(Intercept) 1.0996552 0.3303472 0.4521747 1.747136
## p:session2 0.5611615 0.4789286 -0.3775385 1.499862
## p:session3 1.0994643 0.5822952 -0.0418344 2.240763
## p:session4 0.2773194 0.4635618 -0.6312617 1.185901
## p:session5 0.2583313 0.4671173 -0.6572186 1.173881
## p:session6 0.1797970 0.4562922 -0.7145358 1.074130
## p:session7 0.2379745 0.4782112 -0.6993193 1.175268
## p:session8 0.3169032 0.4736571 -0.6114647 1.245271
model.fits[[model.number]]$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.7549786 0.10522404 0.5487395 0.9612177
## 2 0.9323981 0.03087084 0.8718912 0.9929049
## 3 0.9010631 0.03655224 0.8294207 0.9727055
## 4 0.9044108 0.03467708 0.8364438 0.9723779
## 5 0.9295410 0.03058026 0.8696037 0.9894783
## 6 0.9438816 0.02746434 0.8900515 0.9977117
## 7 0.9178155 0.03421519 0.8507537 0.9848772
## 8 0.8362161 0.04444750 0.7490990 0.9233332
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 1.2349993 0.188177307 0.8661718 1.6038268
## 2 0.9663931 0.013790158 0.9393644 0.9934218
## 3 1.0037153 0.002867533 0.9980950 1.0093357
## 4 1.0277862 0.006648449 1.0147553 1.0408172
## 5 1.0154276 0.004398320 1.0068069 1.0240483
## 6 0.9723841 0.008459235 0.9558040 0.9889642
## 7 0.9110940 0.018476127 0.8748808 0.9473072
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 4.4762223 3.81442530 -3.0000515 11.9524960
## 2 0.6603202 0.11145845 0.4418616 0.8787787
## 3 1.0388677 0.02127775 0.9971633 1.0805721
## 4 1.3943607 0.10927711 1.1801775 1.6085438
## 5 1.2749114 0.07809242 1.1218502 1.4279725
## 6 0.6639769 0.05188654 0.5622793 0.7656746
## 7 0.4571747 0.07716560 0.3059301 0.6084193
model.fits[[model.number]]$results$derived$"psi Probability Occupied"
## estimate se lcl ucl
## 1 0.7549786 0.10522404 0.5487395 0.9612177
## 2 0.9323981 0.03087084 0.8718912 0.9929049
## 3 0.9010631 0.03655224 0.8294207 0.9727055
## 4 0.9044108 0.03467708 0.8364438 0.9723779
## 5 0.9295410 0.03058026 0.8696037 0.9894783
## 6 0.9438816 0.02746434 0.8900515 0.9977117
## 7 0.9178155 0.03421519 0.8507537 0.9848772
## 8 0.8362161 0.04444750 0.7490990 0.9233332
model.fits[[model.number]]$results$derived$"lambda Rate of Change"
## estimate se lcl ucl
## 1 1.2349993 0.188177307 0.8661718 1.6038268
## 2 0.9663931 0.013790158 0.9393644 0.9934218
## 3 1.0037153 0.002867533 0.9980950 1.0093357
## 4 1.0277862 0.006648449 1.0147553 1.0408172
## 5 1.0154276 0.004398320 1.0068069 1.0240483
## 6 0.9723841 0.008459235 0.9558040 0.9889642
## 7 0.9110940 0.018476127 0.8748808 0.9473072
model.fits[[model.number]]$results$derived$"log odds lambda"
## estimate se lcl ucl
## 1 4.4762223 3.81442530 -3.0000515 11.9524960
## 2 0.6603202 0.11145845 0.4418616 0.8787787
## 3 1.0388677 0.02127775 0.9971633 1.0805721
## 4 1.3943607 0.10927711 1.1801775 1.6085438
## 5 1.2749114 0.07809242 1.1218502 1.4279725
## 6 0.6639769 0.05188654 0.5622793 0.7656746
## 7 0.4571747 0.07716560 0.3059301 0.6084193
get.real(model.fits[[model.number]], "Psi", se=TRUE)
## all.diff.index par.index estimate se lcl
## Psi g1 a0 t1 1 1 0.7549786 0.105224 0.5026132
## ucl fixed note group age time Age Time
## Psi g1 a0 t1 0.9038049 1 0 1 0 0
get.real(model.fits[[model.number]], "p", se=TRUE)
## all.diff.index par.index estimate se lcl
## p g1 s1 t1 16 16 0.7501955 0.0619078 0.6111562
## p g1 s1 t2 17 16 0.7501955 0.0619078 0.6111562
## p g1 s1 t3 18 16 0.7501955 0.0619078 0.6111562
## p g1 s2 t1 19 17 0.8403476 0.0465239 0.7273357
## p g1 s2 t2 20 17 0.8403476 0.0465239 0.7273357
## p g1 s2 t3 21 17 0.8403476 0.0465239 0.7273357
## p g1 s3 t1 22 18 0.9001704 0.0430944 0.7788875
## p g1 s3 t2 23 18 0.9001704 0.0430944 0.7788875
## p g1 s3 t3 24 18 0.9001704 0.0430944 0.7788875
## p g1 s4 t1 25 19 0.7985047 0.0523173 0.6769198
## p g1 s4 t2 26 19 0.7985047 0.0523173 0.6769198
## p g1 s4 t3 27 19 0.7985047 0.0523173 0.6769198
## p g1 s5 t1 28 20 0.7954323 0.0537163 0.6706133
## p g1 s5 t2 29 20 0.7954323 0.0537163 0.6706133
## p g1 s5 t3 30 20 0.7954323 0.0537163 0.6706133
## p g1 s6 t1 31 21 0.7823565 0.0535991 0.6598181
## p g1 s6 t2 32 21 0.7823565 0.0535991 0.6598181
## p g1 s6 t3 33 21 0.7823565 0.0535991 0.6598181
## p g1 s7 t1 34 22 0.7920999 0.0571454 0.6586947
## p g1 s7 t2 35 22 0.7920999 0.0571454 0.6586947
## p g1 s7 t3 36 22 0.7920999 0.0571454 0.6586947
## p g1 s8 t1 37 23 0.8047983 0.0533294 0.6794412
## p g1 s8 t2 38 23 0.8047983 0.0533294 0.6794412
## p g1 s8 t3 39 23 0.8047983 0.0533294 0.6794412
## ucl fixed note group time session Time
## p g1 s1 t1 0.8515912 1 1 1 0
## p g1 s1 t2 0.8515912 1 2 1 1
## p g1 s1 t3 0.8515912 1 3 1 2
## p g1 s2 t1 0.9121751 1 1 2 0
## p g1 s2 t2 0.9121751 1 2 2 1
## p g1 s2 t3 0.9121751 1 3 2 2
## p g1 s3 t1 0.9584748 1 1 3 0
## p g1 s3 t2 0.9584748 1 2 3 1
## p g1 s3 t3 0.9584748 1 3 3 2
## p g1 s4 t1 0.8822900 1 1 4 0
## p g1 s4 t2 0.8822900 1 2 4 1
## p g1 s4 t3 0.8822900 1 3 4 2
## p g1 s5 t1 0.8813224 1 1 5 0
## p g1 s5 t2 0.8813224 1 2 5 1
## p g1 s5 t3 0.8813224 1 3 5 2
## p g1 s6 t1 0.8694859 1 1 6 0
## p g1 s6 t2 0.8694859 1 2 6 1
## p g1 s6 t3 0.8694859 1 3 6 2
## p g1 s7 t1 0.8826509 1 1 7 0
## p g1 s7 t2 0.8826509 1 2 7 1
## p g1 s7 t3 0.8826509 1 3 7 2
## p g1 s8 t1 0.8891326 1 1 8 0
## p g1 s8 t2 0.8891326 1 2 8 1
## p g1 s8 t3 0.8891326 1 3 8 2
# collect models and make AICc table
model.set <- RMark::collect.models( type=NULL)
model.set
## model npar AICc DeltaAICc
## 1 Psi(~LN)Epsilon(~minQ)Gamma(~maxQ)p(~session) 14 613.1869 0.00000
## 2 Psi(~LN)Epsilon(~minQ)Gamma(~minQ)p(~session) 14 640.0323 26.84542
## 3 Psi(~LN)Epsilon(~maxQ)Gamma(~minQ)p(~session) 14 696.0170 82.83011
## 4 Psi(~1)Epsilon(~time)Gamma(~time)p(~session) 23 719.5897 106.40281
## weight Deviance
## 1 9.999985e-01 583.0108
## 2 1.481121e-06 609.8562
## 3 0.000000e+00 665.8409
## 4 0.000000e+00 498.1687
#No need to do model averaging in this example as one model has essentially 100% of the weight
# cleanup
cleanup(ask=FALSE)