#Analysis of butterfly data set
#3 seasons, 2 sampling occasions per season
#Interested in occupancy of CLLESSU
#Using RMark
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
# In this case none
bfly.ddl <- make.design.data(bfly.data)
bfly.ddl
## $Psi
## par.index model.index group age time Age Time Treatment
## 1 1 1 burn 0 1 0 0 burn
## 2 2 2 burncontrol 0 1 0 0 burncontrol
## 3 3 3 disk 0 1 0 0 disk
## 4 4 4 diskcontrol 0 1 0 0 diskcontrol
## 5 5 5 fieldcontrol 0 1 0 0 fieldcontrol
##
## $Epsilon
## par.index model.index group age time Age Time Treatment
## 1 1 6 burn 0 1 0 0 burn
## 2 2 7 burn 1 2 1 1 burn
## 3 3 8 burncontrol 0 1 0 0 burncontrol
## 4 4 9 burncontrol 1 2 1 1 burncontrol
## 5 5 10 disk 0 1 0 0 disk
## 6 6 11 disk 1 2 1 1 disk
## 7 7 12 diskcontrol 0 1 0 0 diskcontrol
## 8 8 13 diskcontrol 1 2 1 1 diskcontrol
## 9 9 14 fieldcontrol 0 1 0 0 fieldcontrol
## 10 10 15 fieldcontrol 1 2 1 1 fieldcontrol
##
## $Gamma
## par.index model.index group age time Age Time Treatment
## 1 1 16 burn 0 1 0 0 burn
## 2 2 17 burn 1 2 1 1 burn
## 3 3 18 burncontrol 0 1 0 0 burncontrol
## 4 4 19 burncontrol 1 2 1 1 burncontrol
## 5 5 20 disk 0 1 0 0 disk
## 6 6 21 disk 1 2 1 1 disk
## 7 7 22 diskcontrol 0 1 0 0 diskcontrol
## 8 8 23 diskcontrol 1 2 1 1 diskcontrol
## 9 9 24 fieldcontrol 0 1 0 0 fieldcontrol
## 10 10 25 fieldcontrol 1 2 1 1 fieldcontrol
##
## $p
## par.index model.index group time session Time Treatment
## 1 1 26 burn 1 1 0 burn
## 2 2 27 burn 2 1 1 burn
## 3 3 28 burn 1 2 0 burn
## 4 4 29 burn 2 2 1 burn
## 5 5 30 burn 1 3 0 burn
## 6 6 31 burn 2 3 1 burn
## 7 7 32 burncontrol 1 1 0 burncontrol
## 8 8 33 burncontrol 2 1 1 burncontrol
## 9 9 34 burncontrol 1 2 0 burncontrol
## 10 10 35 burncontrol 2 2 1 burncontrol
## 11 11 36 burncontrol 1 3 0 burncontrol
## 12 12 37 burncontrol 2 3 1 burncontrol
## 13 13 38 disk 1 1 0 disk
## 14 14 39 disk 2 1 1 disk
## 15 15 40 disk 1 2 0 disk
## 16 16 41 disk 2 2 1 disk
## 17 17 42 disk 1 3 0 disk
## 18 18 43 disk 2 3 1 disk
## 19 19 44 diskcontrol 1 1 0 diskcontrol
## 20 20 45 diskcontrol 2 1 1 diskcontrol
## 21 21 46 diskcontrol 1 2 0 diskcontrol
## 22 22 47 diskcontrol 2 2 1 diskcontrol
## 23 23 48 diskcontrol 1 3 0 diskcontrol
## 24 24 49 diskcontrol 2 3 1 diskcontrol
## 25 25 50 fieldcontrol 1 1 0 fieldcontrol
## 26 26 51 fieldcontrol 2 1 1 fieldcontrol
## 27 27 52 fieldcontrol 1 2 0 fieldcontrol
## 28 28 53 fieldcontrol 2 2 1 fieldcontrol
## 29 29 54 fieldcontrol 1 3 0 fieldcontrol
## 30 30 55 fieldcontrol 2 3 1 fieldcontrol
##
## $pimtypes
## $pimtypes$Psi
## $pimtypes$Psi$pim.type
## [1] "all"
##
##
## $pimtypes$Epsilon
## $pimtypes$Epsilon$pim.type
## [1] "all"
##
##
## $pimtypes$Gamma
## $pimtypes$Gamma$pim.type
## [1] "all"
##
##
## $pimtypes$p
## $pimtypes$p$pim.type
## [1] "all"
#Fit Null Model
mod.fit <- RMark::mark(bfly.data,
ddl = bfly.ddl,
model="RDOccupEG",
model.parameters=list(
Psi =list(formula=~1),
p =list(formula=~1),
Epsilon = list(formula=~1),
Gamma = list(formula=~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.2167038 0.3441069 0.5422542 1.8911533
## Epsilon:(Intercept) -2.3035658 0.6402834 -3.5585213 -1.0486103
## Gamma:(Intercept) -0.2540664 0.4468943 -1.1299794 0.6218465
## p:(Intercept) -0.0117223 0.1318397 -0.2701281 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.0908281 0.0908281
##
## Group:Treatmentburncontrol
## 1 2
## 0.0908281 0.0908281
##
## Group:Treatmentdisk
## 1 2
## 0.0908281 0.0908281
##
## Group:Treatmentdiskcontrol
## 1 2
## 0.0908281 0.0908281
##
## Group:Treatmentfieldcontrol
## 1 2
## 0.0908281 0.0908281
##
##
## 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
# Estimates of real initial occupancy, detection, local exctinction,
# and local colonization probabilities
mod.fit$results$real
## estimate se lcl ucl fixed note
## Psi gburn a0 t1 0.7714830 0.0606650 0.6323367 0.8688870
## Epsilon gburn a0 t1 0.0908281 0.0528735 0.0276922 0.2594920
## Gamma gburn a0 t1 0.4368229 0.1099399 0.2441649 0.6506384
## p gburn s1 t1 0.4970694 0.0329588 0.4328757 0.5613600
# Estimates of regression coefficients
mod.fit$results$beta
## estimate se lcl ucl
## Psi:(Intercept) 1.2167038 0.3441069 0.5422542 1.8911533
## Epsilon:(Intercept) -2.3035658 0.6402834 -3.5585213 -1.0486103
## Gamma:(Intercept) -0.2540664 0.4468943 -1.1299794 0.6218465
## p:(Intercept) -0.0117223 0.1318397 -0.2701281 0.2466834
# Estimates of various derived parameters for estimated occpancy in each year,
# estimated ratio of in occupancy from year to year, and estimated log
# of the ratio of odds in occupancy from year to year
mod.fit$results$derived
## $`psi Probability Occupied`
## estimate se lcl ucl
## 1 0.7714830 0.06066501 0.6525795 0.8903864
## 2 0.8012321 0.04891390 0.7053609 0.8971034
## 3 0.8152841 0.05846486 0.7006930 0.9298752
## 4 0.7714830 0.06066501 0.6525795 0.8903864
## 5 0.8012321 0.04891390 0.7053609 0.8971034
## 6 0.8152841 0.05846486 0.7006930 0.9298752
## 7 0.7714830 0.06066501 0.6525795 0.8903864
## 8 0.8012321 0.04891390 0.7053609 0.8971034
## 9 0.8152841 0.05846486 0.7006930 0.9298752
## 10 0.7714830 0.06066501 0.6525795 0.8903864
## 11 0.8012321 0.04891390 0.7053609 0.8971034
## 12 0.8152841 0.05846486 0.7006930 0.9298752
## 13 0.7714830 0.06066501 0.6525795 0.8903864
## 14 0.8012321 0.04891390 0.7053609 0.8971034
## 15 0.8152841 0.05846486 0.7006930 0.9298752
##
## $`lambda Rate of Change`
## estimate se lcl ucl
## 1 1.038561 0.06224831 0.9165543 1.160568
## 2 1.017538 0.02719558 0.9642346 1.070841
## 3 1.038561 0.06224831 0.9165543 1.160568
## 4 1.017538 0.02719558 0.9642346 1.070841
## 5 1.038561 0.06224831 0.9165543 1.160568
## 6 1.017538 0.02719558 0.9642346 1.070841
## 7 1.038561 0.06224831 0.9165543 1.160568
## 8 1.017538 0.02719558 0.9642346 1.070841
## 9 1.038561 0.06224831 0.9165543 1.160568
## 10 1.017538 0.02719558 0.9642346 1.070841
##
## $`log odds lambda`
## estimate se lcl ucl
## 1 1.194000 0.3240320 0.5588974 1.829103
## 2 1.094946 0.1658202 0.7699380 1.419953
## 3 1.194000 0.3240320 0.5588974 1.829103
## 4 1.094946 0.1658202 0.7699380 1.419953
## 5 1.194000 0.3240320 0.5588974 1.829103
## 6 1.094946 0.1658202 0.7699380 1.419953
## 7 1.194000 0.3240320 0.5588974 1.829103
## 8 1.094946 0.1658202 0.7699380 1.419953
## 9 1.194000 0.3240320 0.5588974 1.829103
## 10 1.094946 0.1658202 0.7699380 1.419953
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))
mod.fit$results$derived$"psi Probability Occupied"$Treatment = trt
mod.fit$results$derived$"psi Probability Occupied"
## estimate se lcl ucl Treatment
## 1 0.7714830 0.06066501 0.6525795 0.8903864 burn
## 2 0.8012321 0.04891390 0.7053609 0.8971034 burn
## 3 0.8152841 0.05846486 0.7006930 0.9298752 burn
## 4 0.7714830 0.06066501 0.6525795 0.8903864 burncontrol
## 5 0.8012321 0.04891390 0.7053609 0.8971034 burncontrol
## 6 0.8152841 0.05846486 0.7006930 0.9298752 burncontrol
## 7 0.7714830 0.06066501 0.6525795 0.8903864 disk
## 8 0.8012321 0.04891390 0.7053609 0.8971034 disk
## 9 0.8152841 0.05846486 0.7006930 0.9298752 disk
## 10 0.7714830 0.06066501 0.6525795 0.8903864 diskcontrol
## 11 0.8012321 0.04891390 0.7053609 0.8971034 diskcontrol
## 12 0.8152841 0.05846486 0.7006930 0.9298752 diskcontrol
## 13 0.7714830 0.06066501 0.6525795 0.8903864 fieldcontrol
## 14 0.8012321 0.04891390 0.7053609 0.8971034 fieldcontrol
## 15 0.8152841 0.05846486 0.7006930 0.9298752 fieldcontrol
mod.fit$results$derived$"lambda Rate of Change"$Treatment = trt2
mod.fit$results$derived$"lambda Rate of Change"
## estimate se lcl ucl Treatment
## 1 1.038561 0.06224831 0.9165543 1.160568 burn
## 2 1.017538 0.02719558 0.9642346 1.070841 burn
## 3 1.038561 0.06224831 0.9165543 1.160568 burncontrol
## 4 1.017538 0.02719558 0.9642346 1.070841 burncontrol
## 5 1.038561 0.06224831 0.9165543 1.160568 disk
## 6 1.017538 0.02719558 0.9642346 1.070841 disk
## 7 1.038561 0.06224831 0.9165543 1.160568 diskcontrol
## 8 1.017538 0.02719558 0.9642346 1.070841 diskcontrol
## 9 1.038561 0.06224831 0.9165543 1.160568 fieldcontrol
## 10 1.017538 0.02719558 0.9642346 1.070841 fieldcontrol
mod.fit$results$derived$"log odds lambda"$Treatment = trt2
mod.fit$results$derived$"log odds lambda"
## estimate se lcl ucl Treatment
## 1 1.194000 0.3240320 0.5588974 1.829103 burn
## 2 1.094946 0.1658202 0.7699380 1.419953 burn
## 3 1.194000 0.3240320 0.5588974 1.829103 burncontrol
## 4 1.094946 0.1658202 0.7699380 1.419953 burncontrol
## 5 1.194000 0.3240320 0.5588974 1.829103 disk
## 6 1.094946 0.1658202 0.7699380 1.419953 disk
## 7 1.194000 0.3240320 0.5588974 1.829103 diskcontrol
## 8 1.094946 0.1658202 0.7699380 1.419953 diskcontrol
## 9 1.194000 0.3240320 0.5588974 1.829103 fieldcontrol
## 10 1.094946 0.1658202 0.7699380 1.419953 fieldcontrol
# Plot derived occupancy probabilities
plotdata = mod.fit$results$derived$"psi Probability Occupied"
plotdata$Season = rep(c(1:3),5)
head(plotdata)
## estimate se lcl ucl Treatment Season
## 1 0.7714830 0.06066501 0.6525795 0.8903864 burn 1
## 2 0.8012321 0.04891390 0.7053609 0.8971034 burn 2
## 3 0.8152841 0.05846486 0.7006930 0.9298752 burn 3
## 4 0.7714830 0.06066501 0.6525795 0.8903864 burncontrol 1
## 5 0.8012321 0.04891390 0.7053609 0.8971034 burncontrol 2
## 6 0.8152841 0.05846486 0.7006930 0.9298752 burncontrol 3
ggplot(data=plotdata, aes(x=Season,y=estimate))+
ggtitle("Estimated occupancy 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:3)+
facet_wrap(~Treatment)
cleanup(ask=FALSE)