#Analysis of butterfly data set
#3 seasons, 2 sampling occasions per season
#Interested in occupancy of CLLESSU
#Using RMark
#2018-11-26 code contributed by Neil Faught
library(car)
library(readxl)
library(RMark)
## This is RMark 2.2.5
## Documentation available at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip
library(ggplot2)
library(reshape)
# Get the RMark additional functions
source(file.path("..","..","..","AdditionalFunctions","RMark.additional.functions.R"))
#Read in data
bfly<-read.csv("../Butterfly.csv")
#Data is in long format, need to convert to wide format for RMark
#Want one row for each transect
head(bfly)
## Field Border Transect Treatment Date Visit AMLA AMSN BLSW BUCK
## 1 181 181E 1810 diskcontrol 6/3/2008 1 0 0 0 0
## 2 181 181E 18100 diskcontrol 6/3/2008 1 0 0 0 0
## 3 181 181E 181000 diskcontrol 6/3/2008 1 0 1 0 0
## 4 181 181N 181N1 diskcontrol 6/3/2008 1 0 0 0 0
## 5 181 181N 181N2 diskcontrol 6/3/2008 1 1 0 0 0
## 6 181 181N 181N3 diskcontrol 6/3/2008 1 0 0 0 0
## CHSK CLLESSU CLSK DESK DUSK ETBL ETSW EUSK FISK GISW GPHA GRHA GUFR HAEM
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## 3 0 1 0 0 0 0 0 0 0 0 0 0 0 1
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 1 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## HODW LIYE MONA ORSU PECR PISW RBHA RSPP SACH SICH SLOR SOCL SOSK SPSW
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 1 1 0 0 0 0 0 0 0 0 0 0
## 3 0 0 1 0 1 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## SSSK SWSK TESK VAFR VICE
## 1 0 0 0 0 0
## 2 0 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 0 0 0
## 5 0 0 0 0 0
## 6 0 0 0 0 0
bfly<-subset(bfly,select=c(Field,Border,Transect,Treatment,Date,Visit,CLLESSU))
head(bfly)
## Field Border Transect Treatment Date Visit CLLESSU
## 1 181 181E 1810 diskcontrol 6/3/2008 1 0
## 2 181 181E 18100 diskcontrol 6/3/2008 1 0
## 3 181 181E 181000 diskcontrol 6/3/2008 1 1
## 4 181 181N 181N1 diskcontrol 6/3/2008 1 0
## 5 181 181N 181N2 diskcontrol 6/3/2008 1 0
## 6 181 181N 181N3 diskcontrol 6/3/2008 1 0
bfly$Transect=as.factor(bfly$Transect)
#Convert to wide format
junk<-melt(bfly,id.var=c("Transect","Visit"),measure.var="CLLESSU")
j2=cast(junk,Transect ~ Visit)
head(j2)
## Transect 1 2 3 4 5 6
## 1 1810 0 0 1 0 0 0
## 2 18100 0 0 0 1 0 0
## 3 181000 1 0 0 2 1 1
## 4 181N1 0 0 0 0 0 0
## 5 181N2 0 0 0 0 1 0
## 6 181N3 0 0 0 0 0 1
#Note how there are counts in these columns, need to convert to simply 0/1 for
#present/absent
j2$`1`=ifelse(j2$`1`>0,1,0)
j2$`2`=ifelse(j2$`2`>0,1,0)
j2$`3`=ifelse(j2$`3`>0,1,0)
j2$`4`=ifelse(j2$`4`>0,1,0)
j2$`5`=ifelse(j2$`5`>0,1,0)
j2$`6`=ifelse(j2$`6`>0,1,0)
head(j2)
## Transect 1 2 3 4 5 6
## 1 1810 0 0 1 0 0 0
## 2 18100 0 0 0 1 0 0
## 3 181000 1 0 0 1 1 1
## 4 181N1 0 0 0 0 0 0
## 5 181N2 0 0 0 0 1 0
## 6 181N3 0 0 0 0 0 1
input.history = j2[,2:7]
# do some basic checks on your data
# e.g. check number of sites; number of visits etc
nrow(input.history)
## [1] 129
ncol(input.history)
## [1] 6
range(input.history, na.rm=TRUE)
## [1] 0 1
sum(is.na(input.history))
## [1] 0
#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 001000
## 2 1 000100
## 3 1 100111
## 4 1 000000
## 5 1 000010
## 6 1 000001
# 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 001000
## 2 1 000100
## 3 1 100111
## 4 1 000000
## 5 1 000010
## 6 1 000001
#Is the treatment a site or survey level covariate (i.e. does each transect only
#have a single treatment applied to it over the course of the 6 sampling occasions?)
jcov = melt(bfly,id.var=c("Transect","Visit"),measure.var="Treatment")
jcov2=cast(jcov,Transect ~ Visit)
head(jcov2)
## Transect 1 2 3 4 5
## 1 1810 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 2 18100 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 3 181000 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 4 181N1 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 5 181N2 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 6 181N3 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 6
## 1 diskcontrol
## 2 diskcontrol
## 3 diskcontrol
## 4 diskcontrol
## 5 diskcontrol
## 6 diskcontrol
jcov2$test = ifelse(jcov2$`1`==jcov2$`2` & jcov2$`1`==jcov2$`3` &
jcov2$`1`==jcov2$`4` & jcov2$`1`==jcov2$`5` &
jcov2$`1`==jcov2$`6`,0,1)
head(jcov2)
## Transect 1 2 3 4 5
## 1 1810 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 2 18100 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 3 181000 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 4 181N1 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 5 181N2 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 6 181N3 diskcontrol diskcontrol diskcontrol diskcontrol diskcontrol
## 6 test
## 1 diskcontrol 0
## 2 diskcontrol 0
## 3 diskcontrol 0
## 4 diskcontrol 0
## 5 diskcontrol 0
## 6 diskcontrol 0
#Looks like this is a site level covariate
sum(jcov2$test)
## [1] 0
#Add site covariates to input history
input.history = cbind(input.history,jcov2$`1`)
names(input.history)[3] = "Treatment"
head(input.history)
## freq ch Treatment
## 1810 1 001000 diskcontrol
## 18100 1 000100 diskcontrol
## 181000 1 100111 diskcontrol
## 181N1 1 000000 diskcontrol
## 181N2 1 000010 diskcontrol
## 181N3 1 000001 diskcontrol
#Create the RMark data structure
#2 Seasons, with 3 visits per season
max.visit.per.year <- 2
n.season <- 3
bfly.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)),
group = "Treatment")
summary(bfly.data)
## Length Class Mode
## data 4 data.frame list
## model 1 -none- character
## mixtures 1 -none- numeric
## freq 5 data.frame list
## nocc 1 -none- numeric
## nocc.secondary 3 -none- numeric
## time.intervals 5 -none- numeric
## begin.time 1 -none- numeric
## age.unit 1 -none- numeric
## initial.ages 5 -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
# add visit level covariates (must do in the model fitting loop below)
# Get the parameter names
# What are the parameter names for Single Season Single Species models
setup.parameters("RDOccupEG", check=TRUE)
## [1] "Psi" "Epsilon" "Gamma" "p"
# 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~time
#The random occupancy model cannot be fit here. See the other code in this directory.
# The reserved keywords can be found by looking at the ddl structures earlier
model.list.csv <- textConnection("
p, Psi, Gamma, Epsilon, model.type
~1, ~1, ~1, ~1, RDOccupEG
~session, ~1, ~Treatment, ~Treatment, 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 ~1 ~1 ~1 ~1 RDOccupEG 1
## 2 ~session ~1 ~Treatment ~Treatment RDOccupEG 2
## 3 ~session ~1 ~time ~time RDOccupEG 3
# 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")
# we need to process the data in the loop to allow for different data types
# notice that max.visits.per.year and n.season are defined outside the function (bad programming practise)
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)),
group = "Treatment")
# 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 ~1 ~1 ~1 ~1 RDOccupEG 1
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~1)Gamma(~1)p(~1)
##
## Npar : 4
## -2lnL: 1016.416
## AICc : 1024.521
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 1.2167041 0.3441070 0.5422544 1.8911538
## Epsilon:(Intercept) -2.3035664 0.6402837 -3.5585225 -1.0486103
## Gamma:(Intercept) -0.2540664 0.4468944 -1.1299795 0.6218466
## p:(Intercept) -0.0117224 0.1318397 -0.2701282 0.2466834
##
##
## Real Parameter Psi
## Group:Treatmentburn
## 1
## 0.771483
##
## Group:Treatmentburncontrol
## 1
## 0.771483
##
## Group:Treatmentdisk
## 1
## 0.771483
##
## Group:Treatmentdiskcontrol
## 1
## 0.771483
##
## Group:Treatmentfieldcontrol
## 1
## 0.771483
##
##
## Real Parameter Epsilon
## Group:Treatmentburn
## 1 2
## 0.090828 0.090828
##
## Group:Treatmentburncontrol
## 1 2
## 0.090828 0.090828
##
## Group:Treatmentdisk
## 1 2
## 0.090828 0.090828
##
## Group:Treatmentdiskcontrol
## 1 2
## 0.090828 0.090828
##
## Group:Treatmentfieldcontrol
## 1 2
## 0.090828 0.090828
##
##
## Real Parameter Gamma
## Group:Treatmentburn
## 1 2
## 0.4368229 0.4368229
##
## Group:Treatmentburncontrol
## 1 2
## 0.4368229 0.4368229
##
## Group:Treatmentdisk
## 1 2
## 0.4368229 0.4368229
##
## Group:Treatmentdiskcontrol
## 1 2
## 0.4368229 0.4368229
##
## Group:Treatmentfieldcontrol
## 1 2
## 0.4368229 0.4368229
##
##
## Real Parameter p
## Session:1Group:Treatmentburn
## 1 2
## 0.4970694 0.4970694
##
## Session:2Group:Treatmentburn
## 1 2
## 0.4970694 0.4970694
##
## Session:3Group:Treatmentburn
## 1 2
## 0.4970694 0.4970694
##
## Session:1Group:Treatmentburncontrol
## 1 2
## 0.4970694 0.4970694
##
## Session:2Group:Treatmentburncontrol
## 1 2
## 0.4970694 0.4970694
##
## Session:3Group:Treatmentburncontrol
## 1 2
## 0.4970694 0.4970694
##
## Session:1Group:Treatmentdisk
## 1 2
## 0.4970694 0.4970694
##
## Session:2Group:Treatmentdisk
## 1 2
## 0.4970694 0.4970694
##
## Session:3Group:Treatmentdisk
## 1 2
## 0.4970694 0.4970694
##
## Session:1Group:Treatmentdiskcontrol
## 1 2
## 0.4970694 0.4970694
##
## Session:2Group:Treatmentdiskcontrol
## 1 2
## 0.4970694 0.4970694
##
## Session:3Group:Treatmentdiskcontrol
## 1 2
## 0.4970694 0.4970694
##
## Session:1Group:Treatmentfieldcontrol
## 1 2
## 0.4970694 0.4970694
##
## Session:2Group:Treatmentfieldcontrol
## 1 2
## 0.4970694 0.4970694
##
## Session:3Group:Treatmentfieldcontrol
## 1 2
## 0.4970694 0.4970694
##
##
## ***** Starting ~session ~1 ~Treatment ~Treatment RDOccupEG 2
##
## Note: only 12 parameters counted of 14 specified parameters
## AICc and parameter count have been adjusted upward
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~Treatment)Gamma(~Treatment)p(~session)
##
## Npar : 14 (unadjusted=12)
## -2lnL: 973.8954
## AICc : 1003.024 (unadjusted=998.72962)
##
## Beta
## estimate se lcl
## Psi:(Intercept) 1.8894414 0.4940712 0.9210619
## Epsilon:(Intercept) -1.4960126 0.7560140 -2.9778001
## Epsilon:Treatmentburncontrol -1.0074787 1.1591901 -3.2794913
## Epsilon:Treatmentdisk -0.0658155 1.0644471 -2.1521319
## Epsilon:Treatmentdiskcontrol -1.5905845 1.2452971 -4.0313668
## Epsilon:Treatmentfieldcontrol 0.1102600 0.9954028 -1.8407296
## Gamma:(Intercept) -15.2699440 471.6833100 -939.7692500
## Gamma:Treatmentburncontrol 14.4853470 471.6832600 -910.0138700
## Gamma:Treatmentdisk 14.8362790 471.6863600 -909.6690000
## Gamma:Treatmentdiskcontrol 4.0049516 0.0000000 4.0049516
## Gamma:Treatmentfieldcontrol 16.0601820 471.6841500 -908.4407600
## p:(Intercept) -0.3804102 0.1646492 -0.7031226
## p:session2 -0.1087126 0.2211102 -0.5420886
## p:session3 1.1256377 0.2623158 0.6114988
## ucl
## Psi:(Intercept) 2.8578209
## Epsilon:(Intercept) -0.0142251
## Epsilon:Treatmentburncontrol 1.2645339
## Epsilon:Treatmentdisk 2.0205009
## Epsilon:Treatmentdiskcontrol 0.8501978
## Epsilon:Treatmentfieldcontrol 2.0612496
## Gamma:(Intercept) 909.2293600
## Gamma:Treatmentburncontrol 938.9845600
## Gamma:Treatmentdisk 939.3415500
## Gamma:Treatmentdiskcontrol 4.0049516
## Gamma:Treatmentfieldcontrol 940.5611200
## p:(Intercept) -0.0576979
## p:session2 0.3246635
## p:session3 1.6397766
##
##
## Real Parameter Psi
## Group:Treatmentburn
## 1
## 0.8686918
##
## Group:Treatmentburncontrol
## 1
## 0.8686918
##
## Group:Treatmentdisk
## 1
## 0.8686918
##
## Group:Treatmentdiskcontrol
## 1
## 0.8686918
##
## Group:Treatmentfieldcontrol
## 1
## 0.8686918
##
##
## Real Parameter Epsilon
## Group:Treatmentburn
## 1 2
## 0.183021 0.183021
##
## Group:Treatmentburncontrol
## 1 2
## 0.0756138 0.0756138
##
## Group:Treatmentdisk
## 1 2
## 0.1733845 0.1733845
##
## Group:Treatmentdiskcontrol
## 1 2
## 0.0436635 0.0436635
##
## Group:Treatmentfieldcontrol
## 1 2
## 0.2000867 0.2000867
##
##
## Real Parameter Gamma
## Group:Treatmentburn
## 1 2
## 2.335325e-07 2.335325e-07
##
## Group:Treatmentburncontrol
## 1 2
## 0.3133298 0.3133298
##
## Group:Treatmentdisk
## 1 2
## 0.3932514 0.3932514
##
## Group:Treatmentdiskcontrol
## 1 2
## 1.281357e-05 1.281357e-05
##
## Group:Treatmentfieldcontrol
## 1 2
## 0.6878824 0.6878824
##
##
## Real Parameter p
## Session:1Group:Treatmentburn
## 1 2
## 0.406028 0.406028
##
## Session:2Group:Treatmentburn
## 1 2
## 0.3801002 0.3801002
##
## Session:3Group:Treatmentburn
## 1 2
## 0.6781379 0.6781379
##
## Session:1Group:Treatmentburncontrol
## 1 2
## 0.406028 0.406028
##
## Session:2Group:Treatmentburncontrol
## 1 2
## 0.3801002 0.3801002
##
## Session:3Group:Treatmentburncontrol
## 1 2
## 0.6781379 0.6781379
##
## Session:1Group:Treatmentdisk
## 1 2
## 0.406028 0.406028
##
## Session:2Group:Treatmentdisk
## 1 2
## 0.3801002 0.3801002
##
## Session:3Group:Treatmentdisk
## 1 2
## 0.6781379 0.6781379
##
## Session:1Group:Treatmentdiskcontrol
## 1 2
## 0.406028 0.406028
##
## Session:2Group:Treatmentdiskcontrol
## 1 2
## 0.3801002 0.3801002
##
## Session:3Group:Treatmentdiskcontrol
## 1 2
## 0.6781379 0.6781379
##
## Session:1Group:Treatmentfieldcontrol
## 1 2
## 0.406028 0.406028
##
## Session:2Group:Treatmentfieldcontrol
## 1 2
## 0.3801002 0.3801002
##
## Session:3Group:Treatmentfieldcontrol
## 1 2
## 0.6781379 0.6781379
##
##
## ***** Starting ~session ~1 ~time ~time RDOccupEG 3
##
## Note: only 7 parameters counted of 8 specified parameters
##
## AICc and parameter count have been adjusted upward
##
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~time)Gamma(~time)p(~session)
##
## Npar : 8 (unadjusted=7)
## -2lnL: 979.8915
## AICc : 996.2724 (unadjusted=994.18698)
##
## Beta
## estimate se lcl ucl
## Psi:(Intercept) 2.4924447 0.8018856 0.9207489 4.0641404
## Epsilon:(Intercept) -1.3047920 0.6051544 -2.4908946 -0.1186893
## Epsilon:time2 -1.4910187 1.3321238 -4.1019813 1.1199439
## Gamma:(Intercept) -17.3298980 325.4460000 -655.2040700 620.5442700
## Gamma:time2 16.7815120 325.4460100 -621.0926900 654.6557100
## p:(Intercept) -0.4815543 0.1610004 -0.7971151 -0.1659935
## p:session2 0.1427977 0.2858877 -0.4175421 0.7031375
## p:session3 1.1747008 0.2652653 0.6547809 1.6946208
##
##
## Real Parameter Psi
## Group:Treatmentburn
## 1
## 0.9236105
##
## Group:Treatmentburncontrol
## 1
## 0.9236105
##
## Group:Treatmentdisk
## 1
## 0.9236105
##
## Group:Treatmentdiskcontrol
## 1
## 0.9236105
##
## Group:Treatmentfieldcontrol
## 1
## 0.9236105
##
##
## Real Parameter Epsilon
## Group:Treatmentburn
## 1 2
## 0.2133596 0.057551
##
## Group:Treatmentburncontrol
## 1 2
## 0.2133596 0.057551
##
## Group:Treatmentdisk
## 1 2
## 0.2133596 0.057551
##
## Group:Treatmentdiskcontrol
## 1 2
## 0.2133596 0.057551
##
## Group:Treatmentfieldcontrol
## 1 2
## 0.2133596 0.057551
##
##
## Real Parameter Gamma
## Group:Treatmentburn
## 1 2
## 2.976602e-08 0.3662388
##
## Group:Treatmentburncontrol
## 1 2
## 2.976602e-08 0.3662388
##
## Group:Treatmentdisk
## 1 2
## 2.976602e-08 0.3662388
##
## Group:Treatmentdiskcontrol
## 1 2
## 2.976602e-08 0.3662388
##
## Group:Treatmentfieldcontrol
## 1 2
## 2.976602e-08 0.3662388
##
##
## Real Parameter p
## Session:1Group:Treatmentburn
## 1 2
## 0.3818852 0.3818852
##
## Session:2Group:Treatmentburn
## 1 2
## 0.4161116 0.4161116
##
## Session:3Group:Treatmentburn
## 1 2
## 0.6666665 0.6666665
##
## Session:1Group:Treatmentburncontrol
## 1 2
## 0.3818852 0.3818852
##
## Session:2Group:Treatmentburncontrol
## 1 2
## 0.4161116 0.4161116
##
## Session:3Group:Treatmentburncontrol
## 1 2
## 0.6666665 0.6666665
##
## Session:1Group:Treatmentdisk
## 1 2
## 0.3818852 0.3818852
##
## Session:2Group:Treatmentdisk
## 1 2
## 0.4161116 0.4161116
##
## Session:3Group:Treatmentdisk
## 1 2
## 0.6666665 0.6666665
##
## Session:1Group:Treatmentdiskcontrol
## 1 2
## 0.3818852 0.3818852
##
## Session:2Group:Treatmentdiskcontrol
## 1 2
## 0.4161116 0.4161116
##
## Session:3Group:Treatmentdiskcontrol
## 1 2
## 0.6666665 0.6666665
##
## Session:1Group:Treatmentfieldcontrol
## 1 2
## 0.3818852 0.3818852
##
## Session:2Group:Treatmentfieldcontrol
## 1 2
## 0.4161116 0.4161116
##
## Session:3Group:Treatmentfieldcontrol
## 1 2
## 0.6666665 0.6666665
# examine individual model results
model.number <-2
summary(model.fits[[model.number]])
## Output summary for RDOccupEG model
## Name : Psi(~1)Epsilon(~Treatment)Gamma(~Treatment)p(~session)
##
## Npar : 14 (unadjusted=12)
## -2lnL: 973.8954
## AICc : 1003.024 (unadjusted=998.72962)
##
## Beta
## estimate se lcl
## Psi:(Intercept) 1.8894414 0.4940712 0.9210619
## Epsilon:(Intercept) -1.4960126 0.7560140 -2.9778001
## Epsilon:Treatmentburncontrol -1.0074787 1.1591901 -3.2794913
## Epsilon:Treatmentdisk -0.0658155 1.0644471 -2.1521319
## Epsilon:Treatmentdiskcontrol -1.5905845 1.2452971 -4.0313668
## Epsilon:Treatmentfieldcontrol 0.1102600 0.9954028 -1.8407296
## Gamma:(Intercept) -15.2699440 471.6833100 -939.7692500
## Gamma:Treatmentburncontrol 14.4853470 471.6832600 -910.0138700
## Gamma:Treatmentdisk 14.8362790 471.6863600 -909.6690000
## Gamma:Treatmentdiskcontrol 4.0049516 0.0000000 4.0049516
## Gamma:Treatmentfieldcontrol 16.0601820 471.6841500 -908.4407600
## p:(Intercept) -0.3804102 0.1646492 -0.7031226
## p:session2 -0.1087126 0.2211102 -0.5420886
## p:session3 1.1256377 0.2623158 0.6114988
## ucl
## Psi:(Intercept) 2.8578209
## Epsilon:(Intercept) -0.0142251
## Epsilon:Treatmentburncontrol 1.2645339
## Epsilon:Treatmentdisk 2.0205009
## Epsilon:Treatmentdiskcontrol 0.8501978
## Epsilon:Treatmentfieldcontrol 2.0612496
## Gamma:(Intercept) 909.2293600
## Gamma:Treatmentburncontrol 938.9845600
## Gamma:Treatmentdisk 939.3415500
## Gamma:Treatmentdiskcontrol 4.0049516
## Gamma:Treatmentfieldcontrol 940.5611200
## p:(Intercept) -0.0576979
## p:session2 0.3246635
## p:session3 1.6397766
##
##
## Real Parameter Psi
## Group:Treatmentburn
## 1
## 0.8686918
##
## Group:Treatmentburncontrol
## 1
## 0.8686918
##
## Group:Treatmentdisk
## 1
## 0.8686918
##
## Group:Treatmentdiskcontrol
## 1
## 0.8686918
##
## Group:Treatmentfieldcontrol
## 1
## 0.8686918
##
##
## Real Parameter Epsilon
## Group:Treatmentburn
## 1 2
## 0.183021 0.183021
##
## Group:Treatmentburncontrol
## 1 2
## 0.0756138 0.0756138
##
## Group:Treatmentdisk
## 1 2
## 0.1733845 0.1733845
##
## Group:Treatmentdiskcontrol
## 1 2
## 0.0436635 0.0436635
##
## Group:Treatmentfieldcontrol
## 1 2
## 0.2000867 0.2000867
##
##
## Real Parameter Gamma
## Group:Treatmentburn
## 1 2
## 2.335325e-07 2.335325e-07
##
## Group:Treatmentburncontrol
## 1 2
## 0.3133298 0.3133298
##
## Group:Treatmentdisk
## 1 2
## 0.3932514 0.3932514
##
## Group:Treatmentdiskcontrol
## 1 2
## 1.281357e-05 1.281357e-05
##
## Group:Treatmentfieldcontrol
## 1 2
## 0.6878824 0.6878824
##
##
## Real Parameter p
## Session:1Group:Treatmentburn
## 1 2
## 0.406028 0.406028
##
## Session:2Group:Treatmentburn
## 1 2
## 0.3801002 0.3801002
##
## Session:3Group:Treatmentburn
## 1 2
## 0.6781379 0.6781379
##
## Session:1Group:Treatmentburncontrol
## 1 2
## 0.406028 0.406028
##
## Session:2Group:Treatmentburncontrol
## 1 2
## 0.3801002 0.3801002
##
## Session:3Group:Treatmentburncontrol
## 1 2
## 0.6781379 0.6781379
##
## Session:1Group:Treatmentdisk
## 1 2
## 0.406028 0.406028
##
## Session:2Group:Treatmentdisk
## 1 2
## 0.3801002 0.3801002
##
## Session:3Group:Treatmentdisk
## 1 2
## 0.6781379 0.6781379
##
## Session:1Group:Treatmentdiskcontrol
## 1 2
## 0.406028 0.406028
##
## Session:2Group:Treatmentdiskcontrol
## 1 2
## 0.3801002 0.3801002
##
## Session:3Group:Treatmentdiskcontrol
## 1 2
## 0.6781379 0.6781379
##
## Session:1Group:Treatmentfieldcontrol
## 1 2
## 0.406028 0.406028
##
## Session:2Group:Treatmentfieldcontrol
## 1 2
## 0.3801002 0.3801002
##
## Session:3Group:Treatmentfieldcontrol
## 1 2
## 0.6781379 0.6781379
model.fits[[model.number]]$results$real
## estimate se lcl
## Psi gburn a0 t1 8.686918e-01 0.0563569000 7.152584e-01
## Epsilon gburn a0 t1 1.830210e-01 0.1130425000 4.843890e-02
## Epsilon gburncontrol a0 t1 7.561380e-02 0.0622409000 1.408000e-02
## Epsilon gdisk a0 t1 1.733845e-01 0.1082567000 4.555150e-02
## Epsilon gdiskcontrol a0 t1 4.366350e-02 0.0448342000 5.535400e-03
## Epsilon gfieldcontrol a0 t1 2.000867e-01 0.1013709000 6.741190e-02
## Gamma gburn a0 t1 2.335325e-07 0.0001101534 1.299068e-315
## Gamma gburncontrol a0 t1 3.133298e-01 0.2809111000 3.410450e-02
## Gamma gdisk a0 t1 3.932514e-01 0.4780839000 1.260660e-02
## Gamma gdiskcontrol a0 t1 1.281357e-05 0.0000000000 1.281357e-05
## Gamma gfieldcontrol a0 t1 6.878824e-01 0.2670248000 1.614561e-01
## p gburn s1 t1 4.060280e-01 0.0397083000 3.311203e-01
## p gburn s2 t1 3.801002e-01 0.0382787000 3.084153e-01
## p gburn s3 t1 6.781379e-01 0.0430052000 5.888099e-01
## ucl fixed note
## Psi gburn a0 t1 9.457216e-01
## Epsilon gburn a0 t1 4.964438e-01
## Epsilon gburncontrol a0 t1 3.190455e-01
## Epsilon gdisk a0 t1 4.796695e-01
## Epsilon gdiskcontrol a0 t1 2.724640e-01
## Epsilon gfieldcontrol a0 t1 4.639719e-01
## Gamma gburn a0 t1 1.000000e+00
## Gamma gburncontrol a0 t1 8.550077e-01
## Gamma gdisk a0 t1 9.705028e-01
## Gamma gdiskcontrol a0 t1 1.281357e-05
## Gamma gfieldcontrol a0 t1 9.618712e-01
## p gburn s1 t1 4.855795e-01
## p gburn s2 t1 4.574265e-01
## p gburn s3 t1 7.560985e-01
model.fits[[model.number]]$results$beta
## estimate se lcl
## Psi:(Intercept) 1.8894414 0.4940712 0.9210619
## Epsilon:(Intercept) -1.4960126 0.7560140 -2.9778001
## Epsilon:Treatmentburncontrol -1.0074787 1.1591901 -3.2794913
## Epsilon:Treatmentdisk -0.0658155 1.0644471 -2.1521319
## Epsilon:Treatmentdiskcontrol -1.5905845 1.2452971 -4.0313668
## Epsilon:Treatmentfieldcontrol 0.1102600 0.9954028 -1.8407296
## Gamma:(Intercept) -15.2699440 471.6833100 -939.7692500
## Gamma:Treatmentburncontrol 14.4853470 471.6832600 -910.0138700
## Gamma:Treatmentdisk 14.8362790 471.6863600 -909.6690000
## Gamma:Treatmentdiskcontrol 4.0049516 0.0000000 4.0049516
## Gamma:Treatmentfieldcontrol 16.0601820 471.6841500 -908.4407600
## p:(Intercept) -0.3804102 0.1646492 -0.7031226
## p:session2 -0.1087126 0.2211102 -0.5420886
## p:session3 1.1256377 0.2623158 0.6114988
## ucl
## Psi:(Intercept) 2.8578209
## Epsilon:(Intercept) -0.0142251
## Epsilon:Treatmentburncontrol 1.2645339
## Epsilon:Treatmentdisk 2.0205009
## Epsilon:Treatmentdiskcontrol 0.8501978
## Epsilon:Treatmentfieldcontrol 2.0612496
## Gamma:(Intercept) 909.2293600
## Gamma:Treatmentburncontrol 938.9845600
## Gamma:Treatmentdisk 939.3415500
## Gamma:Treatmentdiskcontrol 4.0049516
## Gamma:Treatmentfieldcontrol 940.5611200
## p:(Intercept) -0.0576979
## p:session2 0.3246635
## p:session3 1.6397766
model.fits[[model.number]]$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.8686918 0.05635690 0.7582323 0.9791514
## 2 0.7097030 0.09509893 0.5233091 0.8960969
## 3 0.5798125 0.15339596 0.2791565 0.8804686
## 4 0.8686918 0.05635690 0.7582323 0.9791514
## 5 0.8441495 0.05709470 0.7322439 0.9560551
## 6 0.8291528 0.08526200 0.6620392 0.9962663
## 7 0.8686918 0.05635690 0.7582323 0.9791514
## 8 0.7697113 0.09379443 0.5858742 0.9535484
## 9 0.7268166 0.14561049 0.4414201 1.0122132
## 10 0.8686918 0.05635690 0.7582323 0.9791514
## 11 0.8307634 0.05888366 0.7153514 0.9461754
## 12 0.7944915 0.08036712 0.6369719 0.9520111
## 13 0.8686918 0.05635690 0.7582323 0.9791514
## 14 0.7852027 0.08037673 0.6276643 0.9427411
## 15 0.7758494 0.08440799 0.6104097 0.9412890
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 0.8169791 0.11304249 0.5954158 1.038542
## 2 0.8169791 0.11304246 0.5954159 1.038542
## 3 0.9717480 0.07065391 0.8332663 1.110230
## 4 0.9822345 0.04863852 0.8869030 1.077566
## 5 0.8860579 0.10909496 0.6722318 1.099884
## 6 0.9442718 0.09619295 0.7557336 1.132810
## 7 0.9563384 0.04482574 0.8684800 1.044197
## 8 0.9563391 0.04480800 0.8685154 1.044163
## 9 0.9038910 0.10322564 0.7015687 1.106213
## 10 0.9880880 0.03471453 0.9200475 1.056128
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 0.3695389 0.2271894 -0.07575225 0.8148301
## 2 0.5644304 0.1203182 0.32860674 0.8002540
## 3 0.8187235 0.4187230 -0.00197365 1.6394206
## 4 0.8960153 0.2433534 0.41904266 1.3729879
## 5 0.5052207 0.3176004 -0.11727601 1.1277175
## 6 0.7960044 0.2624950 0.28151418 1.3104946
## 7 0.7420087 0.2380300 0.27546984 1.2085476
## 8 0.7875470 0.1616894 0.47063577 1.1044582
## 9 0.5525595 0.3550524 -0.14334313 1.2484622
## 10 0.9468570 0.1492750 0.65427811 1.2394359
trt = c(rep("burn",3),rep("burncontrol",3),rep("disk",3),rep("diskcontrol",3),rep("fieldcontrol",3))
trt2 = c(rep("burn",2),rep("burncontrol",2),rep("disk",2),rep("diskcontrol",2),rep("fieldcontrol",2))
model.fits[[model.number]]$results$derived$"psi Probability Occupied"$Treatment = trt
model.fits[[model.number]]$results$derived$"psi Probability Occupied"
## estimate se lcl ucl Treatment
## 1 0.8686918 0.05635690 0.7582323 0.9791514 burn
## 2 0.7097030 0.09509893 0.5233091 0.8960969 burn
## 3 0.5798125 0.15339596 0.2791565 0.8804686 burn
## 4 0.8686918 0.05635690 0.7582323 0.9791514 burncontrol
## 5 0.8441495 0.05709470 0.7322439 0.9560551 burncontrol
## 6 0.8291528 0.08526200 0.6620392 0.9962663 burncontrol
## 7 0.8686918 0.05635690 0.7582323 0.9791514 disk
## 8 0.7697113 0.09379443 0.5858742 0.9535484 disk
## 9 0.7268166 0.14561049 0.4414201 1.0122132 disk
## 10 0.8686918 0.05635690 0.7582323 0.9791514 diskcontrol
## 11 0.8307634 0.05888366 0.7153514 0.9461754 diskcontrol
## 12 0.7944915 0.08036712 0.6369719 0.9520111 diskcontrol
## 13 0.8686918 0.05635690 0.7582323 0.9791514 fieldcontrol
## 14 0.7852027 0.08037673 0.6276643 0.9427411 fieldcontrol
## 15 0.7758494 0.08440799 0.6104097 0.9412890 fieldcontrol
model.fits[[model.number]]$results$derived$"lambda Rate of Change"$Treatment = trt2
model.fits[[model.number]]$results$derived$"lambda Rate of Change"
## estimate se lcl ucl Treatment
## 1 0.8169791 0.11304249 0.5954158 1.038542 burn
## 2 0.8169791 0.11304246 0.5954159 1.038542 burn
## 3 0.9717480 0.07065391 0.8332663 1.110230 burncontrol
## 4 0.9822345 0.04863852 0.8869030 1.077566 burncontrol
## 5 0.8860579 0.10909496 0.6722318 1.099884 disk
## 6 0.9442718 0.09619295 0.7557336 1.132810 disk
## 7 0.9563384 0.04482574 0.8684800 1.044197 diskcontrol
## 8 0.9563391 0.04480800 0.8685154 1.044163 diskcontrol
## 9 0.9038910 0.10322564 0.7015687 1.106213 fieldcontrol
## 10 0.9880880 0.03471453 0.9200475 1.056128 fieldcontrol
model.fits[[model.number]]$results$derived$"log odds lambda"$Treatment = trt2
model.fits[[model.number]]$results$derived$"log odds lambda"
## estimate se lcl ucl Treatment
## 1 0.3695389 0.2271894 -0.07575225 0.8148301 burn
## 2 0.5644304 0.1203182 0.32860674 0.8002540 burn
## 3 0.8187235 0.4187230 -0.00197365 1.6394206 burncontrol
## 4 0.8960153 0.2433534 0.41904266 1.3729879 burncontrol
## 5 0.5052207 0.3176004 -0.11727601 1.1277175 disk
## 6 0.7960044 0.2624950 0.28151418 1.3104946 disk
## 7 0.7420087 0.2380300 0.27546984 1.2085476 diskcontrol
## 8 0.7875470 0.1616894 0.47063577 1.1044582 diskcontrol
## 9 0.5525595 0.3550524 -0.14334313 1.2484622 fieldcontrol
## 10 0.9468570 0.1492750 0.65427811 1.2394359 fieldcontrol
get.real(model.fits[[model.number]], "Psi", se=TRUE)
## Warning in get.real(model.fits[[model.number]], "Psi", se = TRUE):
## Improper V-C matrix for beta estimates. Some variances non-positive.
## all.diff.index par.index estimate se
## Psi gburn a0 t1 1 1 0.8686918 0.0563569
## Psi gburncontrol a0 t1 2 1 0.8686918 0.0563569
## Psi gdisk a0 t1 3 1 0.8686918 0.0563569
## Psi gdiskcontrol a0 t1 4 1 0.8686918 0.0563569
## Psi gfieldcontrol a0 t1 5 1 0.8686918 0.0563569
## lcl ucl fixed note group age
## Psi gburn a0 t1 0.7152584 0.9457216 burn 0
## Psi gburncontrol a0 t1 0.7152584 0.9457216 burncontrol 0
## Psi gdisk a0 t1 0.7152584 0.9457216 disk 0
## Psi gdiskcontrol a0 t1 0.7152584 0.9457216 diskcontrol 0
## Psi gfieldcontrol a0 t1 0.7152584 0.9457216 fieldcontrol 0
## time Age Time Treatment
## Psi gburn a0 t1 1 0 0 burn
## Psi gburncontrol a0 t1 1 0 0 burncontrol
## Psi gdisk a0 t1 1 0 0 disk
## Psi gdiskcontrol a0 t1 1 0 0 diskcontrol
## Psi gfieldcontrol a0 t1 1 0 0 fieldcontrol
get.real(model.fits[[model.number]], "p", se=TRUE)
## Warning in get.real(model.fits[[model.number]], "p", se = TRUE):
## Improper V-C matrix for beta estimates. Some variances non-positive.
## all.diff.index par.index estimate se
## p gburn s1 t1 26 12 0.4060280 0.0397083
## p gburn s1 t2 27 12 0.4060280 0.0397083
## p gburn s2 t1 28 13 0.3801002 0.0382787
## p gburn s2 t2 29 13 0.3801002 0.0382787
## p gburn s3 t1 30 14 0.6781379 0.0430052
## p gburn s3 t2 31 14 0.6781379 0.0430052
## p gburncontrol s1 t1 32 12 0.4060280 0.0397083
## p gburncontrol s1 t2 33 12 0.4060280 0.0397083
## p gburncontrol s2 t1 34 13 0.3801002 0.0382787
## p gburncontrol s2 t2 35 13 0.3801002 0.0382787
## p gburncontrol s3 t1 36 14 0.6781379 0.0430052
## p gburncontrol s3 t2 37 14 0.6781379 0.0430052
## p gdisk s1 t1 38 12 0.4060280 0.0397083
## p gdisk s1 t2 39 12 0.4060280 0.0397083
## p gdisk s2 t1 40 13 0.3801002 0.0382787
## p gdisk s2 t2 41 13 0.3801002 0.0382787
## p gdisk s3 t1 42 14 0.6781379 0.0430052
## p gdisk s3 t2 43 14 0.6781379 0.0430052
## p gdiskcontrol s1 t1 44 12 0.4060280 0.0397083
## p gdiskcontrol s1 t2 45 12 0.4060280 0.0397083
## p gdiskcontrol s2 t1 46 13 0.3801002 0.0382787
## p gdiskcontrol s2 t2 47 13 0.3801002 0.0382787
## p gdiskcontrol s3 t1 48 14 0.6781379 0.0430052
## p gdiskcontrol s3 t2 49 14 0.6781379 0.0430052
## p gfieldcontrol s1 t1 50 12 0.4060280 0.0397083
## p gfieldcontrol s1 t2 51 12 0.4060280 0.0397083
## p gfieldcontrol s2 t1 52 13 0.3801002 0.0382787
## p gfieldcontrol s2 t2 53 13 0.3801002 0.0382787
## p gfieldcontrol s3 t1 54 14 0.6781379 0.0430052
## p gfieldcontrol s3 t2 55 14 0.6781379 0.0430052
## lcl ucl fixed note group time
## p gburn s1 t1 0.3311203 0.4855795 burn 1
## p gburn s1 t2 0.3311203 0.4855795 burn 2
## p gburn s2 t1 0.3084153 0.4574265 burn 1
## p gburn s2 t2 0.3084153 0.4574265 burn 2
## p gburn s3 t1 0.5888099 0.7560985 burn 1
## p gburn s3 t2 0.5888099 0.7560985 burn 2
## p gburncontrol s1 t1 0.3311203 0.4855795 burncontrol 1
## p gburncontrol s1 t2 0.3311203 0.4855795 burncontrol 2
## p gburncontrol s2 t1 0.3084153 0.4574265 burncontrol 1
## p gburncontrol s2 t2 0.3084153 0.4574265 burncontrol 2
## p gburncontrol s3 t1 0.5888099 0.7560985 burncontrol 1
## p gburncontrol s3 t2 0.5888099 0.7560985 burncontrol 2
## p gdisk s1 t1 0.3311203 0.4855795 disk 1
## p gdisk s1 t2 0.3311203 0.4855795 disk 2
## p gdisk s2 t1 0.3084153 0.4574265 disk 1
## p gdisk s2 t2 0.3084153 0.4574265 disk 2
## p gdisk s3 t1 0.5888099 0.7560985 disk 1
## p gdisk s3 t2 0.5888099 0.7560985 disk 2
## p gdiskcontrol s1 t1 0.3311203 0.4855795 diskcontrol 1
## p gdiskcontrol s1 t2 0.3311203 0.4855795 diskcontrol 2
## p gdiskcontrol s2 t1 0.3084153 0.4574265 diskcontrol 1
## p gdiskcontrol s2 t2 0.3084153 0.4574265 diskcontrol 2
## p gdiskcontrol s3 t1 0.5888099 0.7560985 diskcontrol 1
## p gdiskcontrol s3 t2 0.5888099 0.7560985 diskcontrol 2
## p gfieldcontrol s1 t1 0.3311203 0.4855795 fieldcontrol 1
## p gfieldcontrol s1 t2 0.3311203 0.4855795 fieldcontrol 2
## p gfieldcontrol s2 t1 0.3084153 0.4574265 fieldcontrol 1
## p gfieldcontrol s2 t2 0.3084153 0.4574265 fieldcontrol 2
## p gfieldcontrol s3 t1 0.5888099 0.7560985 fieldcontrol 1
## p gfieldcontrol s3 t2 0.5888099 0.7560985 fieldcontrol 2
## session Time Treatment
## p gburn s1 t1 1 0 burn
## p gburn s1 t2 1 1 burn
## p gburn s2 t1 2 0 burn
## p gburn s2 t2 2 1 burn
## p gburn s3 t1 3 0 burn
## p gburn s3 t2 3 1 burn
## p gburncontrol s1 t1 1 0 burncontrol
## p gburncontrol s1 t2 1 1 burncontrol
## p gburncontrol s2 t1 2 0 burncontrol
## p gburncontrol s2 t2 2 1 burncontrol
## p gburncontrol s3 t1 3 0 burncontrol
## p gburncontrol s3 t2 3 1 burncontrol
## p gdisk s1 t1 1 0 disk
## p gdisk s1 t2 1 1 disk
## p gdisk s2 t1 2 0 disk
## p gdisk s2 t2 2 1 disk
## p gdisk s3 t1 3 0 disk
## p gdisk s3 t2 3 1 disk
## p gdiskcontrol s1 t1 1 0 diskcontrol
## p gdiskcontrol s1 t2 1 1 diskcontrol
## p gdiskcontrol s2 t1 2 0 diskcontrol
## p gdiskcontrol s2 t2 2 1 diskcontrol
## p gdiskcontrol s3 t1 3 0 diskcontrol
## p gdiskcontrol s3 t2 3 1 diskcontrol
## p gfieldcontrol s1 t1 1 0 fieldcontrol
## p gfieldcontrol s1 t2 1 1 fieldcontrol
## p gfieldcontrol s2 t1 2 0 fieldcontrol
## p gfieldcontrol s2 t2 2 1 fieldcontrol
## p gfieldcontrol s3 t1 3 0 fieldcontrol
## p gfieldcontrol s3 t2 3 1 fieldcontrol
# collect models and make AICc table
model.set <- RMark::collect.models( type=NULL)
model.set
## model npar AICc
## 3 Psi(~1)Epsilon(~time)Gamma(~time)p(~session) 8 996.2724
## 2 Psi(~1)Epsilon(~Treatment)Gamma(~Treatment)p(~session) 14 1003.0244
## 1 Psi(~1)Epsilon(~1)Gamma(~1)p(~1) 4 1024.5205
## DeltaAICc weight Deviance
## 3 0.00000 9.669455e-01 306.6549
## 2 6.75201 3.305382e-02 300.6589
## 1 28.24809 7.102435e-07 343.1793
# model averaging in the usual way
RMark::model.average(model.set, "Psi")
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## par.index estimate se fixed note
## Psi gburn a0 t1 1 0.9217951 0.05741498
## Psi gburncontrol a0 t1 2 0.9217951 0.05741498
## Psi gdisk a0 t1 3 0.9217951 0.05741498
## Psi gdiskcontrol a0 t1 4 0.9217951 0.05741498
## Psi gfieldcontrol a0 t1 5 0.9217951 0.05741498
## group age time Age Time Treatment
## Psi gburn a0 t1 burn 0 1 0 0 burn
## Psi gburncontrol a0 t1 burncontrol 0 1 0 0 burncontrol
## Psi gdisk a0 t1 disk 0 1 0 0 disk
## Psi gdiskcontrol a0 t1 diskcontrol 0 1 0 0 diskcontrol
## Psi gfieldcontrol a0 t1 fieldcontrol 0 1 0 0 fieldcontrol
RMark::model.average(model.set, "p")
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
## par.index estimate se fixed note
## p gburn s1 t1 26 0.3826833 0.03830559
## p gburn s1 t2 27 0.3826833 0.03830559
## p gburn s2 t1 28 0.4149213 0.06335009
## p gburn s2 t2 29 0.4149213 0.06335009
## p gburn s3 t1 30 0.6670456 0.04677183
## p gburn s3 t2 31 0.6670456 0.04677183
## p gburncontrol s1 t1 32 0.3826833 0.03830559
## p gburncontrol s1 t2 33 0.3826833 0.03830559
## p gburncontrol s2 t1 34 0.4149213 0.06335009
## p gburncontrol s2 t2 35 0.4149213 0.06335009
## p gburncontrol s3 t1 36 0.6670456 0.04677183
## p gburncontrol s3 t2 37 0.6670456 0.04677183
## p gdisk s1 t1 38 0.3826833 0.03830559
## p gdisk s1 t2 39 0.3826833 0.03830559
## p gdisk s2 t1 40 0.4149213 0.06335009
## p gdisk s2 t2 41 0.4149213 0.06335009
## p gdisk s3 t1 42 0.6670456 0.04677183
## p gdisk s3 t2 43 0.6670456 0.04677183
## p gdiskcontrol s1 t1 44 0.3826833 0.03830559
## p gdiskcontrol s1 t2 45 0.3826833 0.03830559
## p gdiskcontrol s2 t1 46 0.4149213 0.06335009
## p gdiskcontrol s2 t2 47 0.4149213 0.06335009
## p gdiskcontrol s3 t1 48 0.6670456 0.04677183
## p gdiskcontrol s3 t2 49 0.6670456 0.04677183
## p gfieldcontrol s1 t1 50 0.3826833 0.03830559
## p gfieldcontrol s1 t2 51 0.3826833 0.03830559
## p gfieldcontrol s2 t1 52 0.4149213 0.06335009
## p gfieldcontrol s2 t2 53 0.4149213 0.06335009
## p gfieldcontrol s3 t1 54 0.6670456 0.04677183
## p gfieldcontrol s3 t2 55 0.6670456 0.04677183
## group time session Time Treatment
## p gburn s1 t1 burn 1 1 0 burn
## p gburn s1 t2 burn 2 1 1 burn
## p gburn s2 t1 burn 1 2 0 burn
## p gburn s2 t2 burn 2 2 1 burn
## p gburn s3 t1 burn 1 3 0 burn
## p gburn s3 t2 burn 2 3 1 burn
## p gburncontrol s1 t1 burncontrol 1 1 0 burncontrol
## p gburncontrol s1 t2 burncontrol 2 1 1 burncontrol
## p gburncontrol s2 t1 burncontrol 1 2 0 burncontrol
## p gburncontrol s2 t2 burncontrol 2 2 1 burncontrol
## p gburncontrol s3 t1 burncontrol 1 3 0 burncontrol
## p gburncontrol s3 t2 burncontrol 2 3 1 burncontrol
## p gdisk s1 t1 disk 1 1 0 disk
## p gdisk s1 t2 disk 2 1 1 disk
## p gdisk s2 t1 disk 1 2 0 disk
## p gdisk s2 t2 disk 2 2 1 disk
## p gdisk s3 t1 disk 1 3 0 disk
## p gdisk s3 t2 disk 2 3 1 disk
## p gdiskcontrol s1 t1 diskcontrol 1 1 0 diskcontrol
## p gdiskcontrol s1 t2 diskcontrol 2 1 1 diskcontrol
## p gdiskcontrol s2 t1 diskcontrol 1 2 0 diskcontrol
## p gdiskcontrol s2 t2 diskcontrol 2 2 1 diskcontrol
## p gdiskcontrol s3 t1 diskcontrol 1 3 0 diskcontrol
## p gdiskcontrol s3 t2 diskcontrol 2 3 1 diskcontrol
## p gfieldcontrol s1 t1 fieldcontrol 1 1 0 fieldcontrol
## p gfieldcontrol s1 t2 fieldcontrol 2 1 1 fieldcontrol
## p gfieldcontrol s2 t1 fieldcontrol 1 2 0 fieldcontrol
## p gfieldcontrol s2 t2 fieldcontrol 2 2 1 fieldcontrol
## p gfieldcontrol s3 t1 fieldcontrol 1 3 0 fieldcontrol
## p gfieldcontrol s3 t2 fieldcontrol 2 3 1 fieldcontrol
# Model average the derived parameters
# Because RMarks stores psi in different places, we standarize the dervied parameters.
# We need to do this here because collect.models re-extracts the output from MARK and wipes anything else
model.set[-length(model.set)] <- plyr::llply(model.set[-length(model.set)], function(x){RMark.add.derived(x)})
# NOTE: Because a group was used in the process.data step above, the first half
# of each derived parameter table will be for Prox=No, and the second
# will be for Prox = Yes
model.set[[1]]$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.7714830 0.06066502 0.6525796 0.8903864
## 2 0.8012322 0.04891390 0.7053609 0.8971034
## 3 0.8152842 0.05846487 0.7006930 0.9298753
## 4 0.7714830 0.06066502 0.6525796 0.8903864
## 5 0.8012322 0.04891390 0.7053609 0.8971034
## 6 0.8152842 0.05846487 0.7006930 0.9298753
## 7 0.7714830 0.06066502 0.6525796 0.8903864
## 8 0.8012322 0.04891390 0.7053609 0.8971034
## 9 0.8152842 0.05846487 0.7006930 0.9298753
## 10 0.7714830 0.06066502 0.6525796 0.8903864
## 11 0.8012322 0.04891390 0.7053609 0.8971034
## 12 0.8152842 0.05846487 0.7006930 0.9298753
## 13 0.7714830 0.06066502 0.6525796 0.8903864
## 14 0.8012322 0.04891390 0.7053609 0.8971034
## 15 0.8152842 0.05846487 0.7006930 0.9298753
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 1.038561 0.06224830 0.9165543 1.160568
## 2 1.017538 0.02719558 0.9642346 1.070841
## 3 1.038561 0.06224830 0.9165543 1.160568
## 4 1.017538 0.02719558 0.9642346 1.070841
## 5 1.038561 0.06224830 0.9165543 1.160568
## 6 1.017538 0.02719558 0.9642346 1.070841
## 7 1.038561 0.06224830 0.9165543 1.160568
## 8 1.017538 0.02719558 0.9642346 1.070841
## 9 1.038561 0.06224830 0.9165543 1.160568
## 10 1.017538 0.02719558 0.9642346 1.070841
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 1.194000 0.3240321 0.5588974 1.829103
## 2 1.094946 0.1658203 0.7699379 1.419954
## 3 1.194000 0.3240321 0.5588974 1.829103
## 4 1.094946 0.1658203 0.7699379 1.419954
## 5 1.194000 0.3240321 0.5588974 1.829103
## 6 1.094946 0.1658203 0.7699379 1.419954
## 7 1.194000 0.3240321 0.5588974 1.829103
## 8 1.094946 0.1658203 0.7699379 1.419954
## 9 1.194000 0.3240321 0.5588974 1.829103
## 10 1.094946 0.1658203 0.7699379 1.419954
##
## $all_psi
## estimate se lcl ucl
## 1 0.7714830 0.06066502 0.6525796 0.8903864
## 2 0.8012322 0.04891390 0.7053609 0.8971034
## 3 0.8152842 0.05846487 0.7006930 0.9298753
## 4 0.7714830 0.06066502 0.6525796 0.8903864
## 5 0.8012322 0.04891390 0.7053609 0.8971034
## 6 0.8152842 0.05846487 0.7006930 0.9298753
## 7 0.7714830 0.06066502 0.6525796 0.8903864
## 8 0.8012322 0.04891390 0.7053609 0.8971034
## 9 0.8152842 0.05846487 0.7006930 0.9298753
## 10 0.7714830 0.06066502 0.6525796 0.8903864
## 11 0.8012322 0.04891390 0.7053609 0.8971034
## 12 0.8152842 0.05846487 0.7006930 0.9298753
## 13 0.7714830 0.06066502 0.6525796 0.8903864
## 14 0.8012322 0.04891390 0.7053609 0.8971034
## 15 0.8152842 0.05846487 0.7006930 0.9298753
d_all_psi = RMark.model.average.derived(model.set, "all_psi")
## Loading required package: plyr
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:reshape':
##
## rename, round_any
## Loading required package: boot
##
## Attaching package: 'boot'
## The following object is masked from 'package:car':
##
## logit
d_all_psi$Treatment = trt
d_all_psi
## estimate se lcl ucl Treatment
## 1 0.9217951 0.05741502 0.7121773 0.9825017 burn
## 2 0.7259925 0.09968603 0.4980529 0.8761598 burn
## 3 0.7781056 0.06970129 0.6138271 0.8855319 burn
## 4 0.9217951 0.05741502 0.7121773 0.9825017 burncontrol
## 5 0.7304365 0.10089137 0.4981374 0.8809153 burncontrol
## 6 0.7863472 0.05513001 0.6592367 0.8750314 burncontrol
## 7 0.9217951 0.05741502 0.7121773 0.9825017 disk
## 8 0.7279760 0.09989810 0.4989100 0.8779453 disk
## 9 0.7829646 0.05953976 0.6448100 0.8775844 disk
## 10 0.9217951 0.05741502 0.7121773 0.9825017 diskcontrol
## 11 0.7299940 0.10045410 0.4989224 0.8801123 diskcontrol
## 12 0.7852015 0.05433999 0.6603303 0.8729962 diskcontrol
## 13 0.9217951 0.05741502 0.7121773 0.9825017 fieldcontrol
## 14 0.7284880 0.09976370 0.4995977 0.8782048 fieldcontrol
## 15 0.7845853 0.05453899 0.6592875 0.8727016 fieldcontrol
d_lamb_rc = RMark.model.average.derived(model.set, "lambda Rate of Change")
d_lamb_rc$Treatment = trt2
d_lamb_rc
## estimate se lcl ucl Treatment
## 1 0.7876434 0.1021118 0.5875080 0.9877787 burn
## 2 1.0715865 0.1599945 0.7580030 1.3851701 burn
## 3 0.7927591 0.1059960 0.5850106 1.0005075 burncontrol
## 4 1.0770489 0.1527921 0.7775819 1.3765158 burncontrol
## 5 0.7899267 0.1033649 0.5873351 0.9925182 disk
## 6 1.0757941 0.1544574 0.7730632 1.3785249 disk
## 7 0.7922497 0.1046988 0.5870439 0.9974555 diskcontrol
## 8 1.0761929 0.1533536 0.7756254 1.3767604 diskcontrol
## 9 0.7905161 0.1037623 0.5871458 0.9938865 fieldcontrol
## 10 1.0772423 0.1525498 0.7782501 1.3762346 fieldcontrol
d_lamb_lo = RMark.model.average.derived(model.set, "log odds lambda")
d_lamb_lo$Treatment = trt2
d_lamb_lo
## estimate se lcl ucl Treatment
## 1 0.2247028 0.1758887 -0.12003276 0.5694384 burn
## 2 1.3465060 0.7523149 -0.12800406 2.8210160 burn
## 3 0.2395501 0.2139524 -0.17978889 0.6588891 burncontrol
## 4 1.3574661 0.7441958 -0.10113085 2.8160631 burncontrol
## 5 0.2291876 0.1856137 -0.13460845 0.5929837 disk
## 6 1.3541604 0.7466713 -0.10928844 2.8176092 disk
## 7 0.2370144 0.1977472 -0.15056307 0.6245918 diskcontrol
## 8 1.3538808 0.7459350 -0.10812483 2.8158865 diskcontrol
## 9 0.2307524 0.1903171 -0.14226223 0.6037670 fieldcontrol
## 10 1.3591466 0.7423868 -0.09590471 2.8141980 fieldcontrol
# model averaging of derived parameters such as the occupancy at each time step
psi.ma <- RMark.model.average.derived(model.set, "all_psi")
psi.ma$Year <- rep(1:(nrow(psi.ma)/5),5)
psi.ma$parameter <- 'psi'
psi.ma$Treatment = trt
psi.ma
## estimate se lcl ucl Year parameter Treatment
## 1 0.9217951 0.05741502 0.7121773 0.9825017 1 psi burn
## 2 0.7259925 0.09968603 0.4980529 0.8761598 2 psi burn
## 3 0.7781056 0.06970129 0.6138271 0.8855319 3 psi burn
## 4 0.9217951 0.05741502 0.7121773 0.9825017 1 psi burncontrol
## 5 0.7304365 0.10089137 0.4981374 0.8809153 2 psi burncontrol
## 6 0.7863472 0.05513001 0.6592367 0.8750314 3 psi burncontrol
## 7 0.9217951 0.05741502 0.7121773 0.9825017 1 psi disk
## 8 0.7279760 0.09989810 0.4989100 0.8779453 2 psi disk
## 9 0.7829646 0.05953976 0.6448100 0.8775844 3 psi disk
## 10 0.9217951 0.05741502 0.7121773 0.9825017 1 psi diskcontrol
## 11 0.7299940 0.10045410 0.4989224 0.8801123 2 psi diskcontrol
## 12 0.7852015 0.05433999 0.6603303 0.8729962 3 psi diskcontrol
## 13 0.9217951 0.05741502 0.7121773 0.9825017 1 psi fieldcontrol
## 14 0.7284880 0.09976370 0.4995977 0.8782048 2 psi fieldcontrol
## 15 0.7845853 0.05453899 0.6592875 0.8727016 3 psi fieldcontrol
# likely more interested in colonization and extinction probabilities
epsilon.ma <- RMark::model.average(model.set, "Epsilon", vcv=TRUE)$estimates
## Warning in get.real(model, parameter, design = model$design.matrix, se = TRUE, :
## Improper V-C matrix for beta estimates. Some variances non-positive.
epsilon.ma$Year <- rep(1:(nrow(epsilon.ma)/5),5)
epsilon.ma$parameter <- 'epsilon'
epsilon.ma
## par.index estimate se lcl
## Epsilon gburn a0 t1 6 0.21235674 0.10211159 0.075346751
## Epsilon gburn a1 t2 7 0.06169826 0.06829413 0.006470785
## Epsilon gburncontrol a0 t1 8 0.20880652 0.10348649 0.071766202
## Epsilon gburncontrol a1 t2 9 0.05814805 0.06226577 0.006605751
## Epsilon gdisk a0 t1 10 0.21203822 0.10204619 0.075176099
## Epsilon gdisk a1 t2 11 0.06137974 0.06748897 0.006539718
## Epsilon gdiskcontrol a0 t1 12 0.20775044 0.10469861 0.070090174
## Epsilon gdiskcontrol a1 t2 13 0.05709197 0.06173456 0.006356043
## Epsilon gfieldcontrol a0 t1 14 0.21292083 0.10158873 0.076172842
## Epsilon gfieldcontrol a1 t2 15 0.06226235 0.06875719 0.006559642
## ucl fixed note group age time
## Epsilon gburn a0 t1 0.4714718 burn 0 1
## Epsilon gburn a1 t2 0.3989923 burn 1 2
## Epsilon gburncontrol a0 t1 0.4739241 burncontrol 0 1
## Epsilon gburncontrol a1 t2 0.3643519 burncontrol 1 2
## Epsilon gdisk a0 t1 0.4711330 disk 0 1
## Epsilon gdisk a1 t2 0.3938010 disk 1 2
## Epsilon gdiskcontrol a0 t1 0.4770729 diskcontrol 0 1
## Epsilon gdiskcontrol a1 t2 0.3643258 diskcontrol 1 2
## Epsilon gfieldcontrol a0 t1 0.4702110 fieldcontrol 0 1
## Epsilon gfieldcontrol a1 t2 0.4003545 fieldcontrol 1 2
## Age Time Treatment Year parameter
## Epsilon gburn a0 t1 0 0 burn 1 epsilon
## Epsilon gburn a1 t2 1 1 burn 2 epsilon
## Epsilon gburncontrol a0 t1 0 0 burncontrol 1 epsilon
## Epsilon gburncontrol a1 t2 1 1 burncontrol 2 epsilon
## Epsilon gdisk a0 t1 0 0 disk 1 epsilon
## Epsilon gdisk a1 t2 1 1 disk 2 epsilon
## Epsilon gdiskcontrol a0 t1 0 0 diskcontrol 1 epsilon
## Epsilon gdiskcontrol a1 t2 1 1 diskcontrol 2 epsilon
## Epsilon gfieldcontrol a0 t1 0 0 fieldcontrol 1 epsilon
## Epsilon gfieldcontrol a1 t2 1 1 fieldcontrol 2 epsilon
gamma.ma <- RMark::model.average(model.set, "Gamma", vcv=TRUE)$estimates
##
## Model 2dropped from the model averaging because one or more beta variances are not positive
gamma.ma$Year <- rep(1:(nrow(gamma.ma)/5),5)
gamma.ma$parameter <- 'gamma'
gamma.ma
## par.index estimate se lcl
## Gamma gburn a0 t1 16 3.506221e-07 0.0003861723 0.00000000
## Gamma gburn a1 t2 17 3.662390e-01 0.2216862566 0.08163424
## Gamma gburncontrol a0 t1 18 3.506221e-07 0.0003861723 0.00000000
## Gamma gburncontrol a1 t2 19 3.662390e-01 0.2216862566 0.08163424
## Gamma gdisk a0 t1 20 3.506221e-07 0.0003861723 0.00000000
## Gamma gdisk a1 t2 21 3.662390e-01 0.2216862566 0.08163424
## Gamma gdiskcontrol a0 t1 22 3.506221e-07 0.0003861723 0.00000000
## Gamma gdiskcontrol a1 t2 23 3.662390e-01 0.2216862566 0.08163424
## Gamma gfieldcontrol a0 t1 24 3.506221e-07 0.0003861723 0.00000000
## Gamma gfieldcontrol a1 t2 25 3.662390e-01 0.2216862566 0.08163424
## ucl fixed note group age time
## Gamma gburn a0 t1 1.0000000 burn 0 1
## Gamma gburn a1 t2 0.7897759 burn 1 2
## Gamma gburncontrol a0 t1 1.0000000 burncontrol 0 1
## Gamma gburncontrol a1 t2 0.7897759 burncontrol 1 2
## Gamma gdisk a0 t1 1.0000000 disk 0 1
## Gamma gdisk a1 t2 0.7897759 disk 1 2
## Gamma gdiskcontrol a0 t1 1.0000000 diskcontrol 0 1
## Gamma gdiskcontrol a1 t2 0.7897759 diskcontrol 1 2
## Gamma gfieldcontrol a0 t1 1.0000000 fieldcontrol 0 1
## Gamma gfieldcontrol a1 t2 0.7897759 fieldcontrol 1 2
## Age Time Treatment Year parameter
## Gamma gburn a0 t1 0 0 burn 1 gamma
## Gamma gburn a1 t2 1 1 burn 2 gamma
## Gamma gburncontrol a0 t1 0 0 burncontrol 1 gamma
## Gamma gburncontrol a1 t2 1 1 burncontrol 2 gamma
## Gamma gdisk a0 t1 0 0 disk 1 gamma
## Gamma gdisk a1 t2 1 1 disk 2 gamma
## Gamma gdiskcontrol a0 t1 0 0 diskcontrol 1 gamma
## Gamma gdiskcontrol a1 t2 1 1 diskcontrol 2 gamma
## Gamma gfieldcontrol a0 t1 0 0 fieldcontrol 1 gamma
## Gamma gfieldcontrol a1 t2 1 1 fieldcontrol 2 gamma
all.est <- plyr::rbind.fill(psi.ma, epsilon.ma, gamma.ma)
ggplot(data=all.est, aes(x=Year,y=estimate, color=parameter))+
ggtitle("Estimated occupancy, extinction, and colonization over time")+
geom_point(position=position_dodge(w=0.2))+
geom_line(position=position_dodge(w=0.2))+
ylim(0,1)+
geom_errorbar(aes(ymin=lcl, ymax=ucl), width=.1,position=position_dodge(w=0.2))+
scale_x_continuous(breaks=1:10)+
facet_wrap(~Treatment)+
xlab("Season")
cleanup(ask=FALSE)