# NOT RUN {
# }
# NOT RUN {
data(lfbcvi)
xy=cut(lfbcvi$Pred, c(-0.0001, .1, .2, .3, .4, .5, .6, .7, .8, .9, 1),
labels=c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"))
x=data.frame(lfbcvi, New=xy)
# Note that I scaled the individual covariate of day-helps with
# convergence issues
bird.data <- data.frame(object=x$ObjectID, observer=x$Observer,
detected=x$Detected, distance=x$Distance,
Region.Label=x$New, Sample.Label=x$PointID,
Day=(x$Day/max(x$Day)))
# make observer a factor variable
bird.data$observer=factor(bird.data$observer)
# Jeff Laake suggested this snippet to quickly create distance medians
# which adds bin information to the bird.data dataframe
bird.data$distbegin=0
bird.data$distend=100
bird.data$distend[bird.data$distance==12.5]=25
bird.data$distbegin[bird.data$distance==37.5]=25
bird.data$distend[bird.data$distance==37.5]=50
bird.data$distbegin[bird.data$distance==62.5]=50
bird.data$distend[bird.data$distance==62.5]=75
bird.data$distbegin[bird.data$distance==87.5]=75
bird.data$distend[bird.data$distance==87.5]=100
# Removed all survey points with distance=NA for a survey event;
# hence no observations for use in ddf() but needed later
bird.data=bird.data[complete.cases(bird.data),]
# Manipulations on full dataset for various data.frame creation for
# use in density estimation using dht()
#Samples dataframe
xx=x
x=data.frame(PointID=x$PointID, Species=x$Species,
Category=x$New, Effort=x$Effort)
x=x[!duplicated(x$PointID),]
point.num=table(x$Category)
samples=data.frame(PointID=x$PointID, Region.Label=x$Category,
Effort=x$Effort)
final.samples=data.frame(Sample.Label=samples$PointID,
Region.Label=samples$Region.Label,
Effort=samples$Effort)
#obs dataframe
obs=data.frame(ObjectID=xx$ObjectID, PointID=xx$PointID)
#used to get Region and Sample assigned to ObjectID
obs=merge(obs, samples, by=c("PointID", "PointID"))
obs=obs[!duplicated(obs$ObjectID),]
obs=data.frame(object=obs$ObjectID, Region.Label=obs$Region.Label,
Sample.Label=obs$PointID)
region.data=data.frame(Region.Label=c(1, 2, 3,4,5,6,7,8,9, 10),
Area=c(point.num[1]*3.14, point.num[2]*3.14,
point.num[3]*3.14, point.num[4]*3.14,
point.num[5]*3.14, point.num[6]*3.14,
point.num[7]*3.14, point.num[8]*3.14,
point.num[9]*3.14, point.num[10]*3.14))
# Candidate Models
BV1=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance),
data=bird.data,
method="io",
meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
BV1FI=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance),
data=bird.data,
method="io.fi",
meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
BV2=ddf(
dsmodel=~mcds(key="hr",formula=~1),
mrmodel=~glm(~distance),
data=bird.data,
method="io",
meta.data=list(binned=TRUE,point=TRUE,width=100,breaks=c(0,50,100)))
BV3=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance+observer),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV3FI=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance+observer),
data=bird.data,
method="io.fi",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV4=ddf(
dsmodel=~mcds(key="hr",formula=~1),
mrmodel=~glm(~distance+observer),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV5=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance*observer),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV5FI=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance*observer),
data=bird.data,
method="io.fi",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV6=ddf(
dsmodel=~mcds(key="hr",formula=~1),
mrmodel=~glm(~distance*observer),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV7=ddf(
dsmodel=~cds(key="hn",formula=~1),
mrmodel=~glm(~distance*Day),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV7FI=ddf(
dsmodel=~cds(key="hn",formula=~1),
mrmodel=~glm(~distance*Day),
data=bird.data,
method="io.fi",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV8=ddf(
dsmodel=~cds(key="hr",formula=~1),
mrmodel=~glm(~distance*Day),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV9=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance*observer*Day),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV9FI=ddf(
dsmodel=~mcds(key="hn",formula=~1),
mrmodel=~glm(~distance*observer*Day),
data=bird.data,
method="io.fi",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
BV10=ddf(
dsmodel=~mcds(key="hr",formula=~1),
mrmodel=~glm(~distance*observer*Day),
data=bird.data,
method="io",
meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
#BV.DS=ddf(
# dsmodel=~mcds(key="hn",formula=~1),
# data=bird.data,
# method="ds",
# meta.data=list(binned=TRUE, point=TRUE, width=100,breaks=c(0,50,100)))
#AIC table building code.
AIC = c(BV1$criterion, BV1FI$criterion, BV2$criterion, BV3$criterion,
BV3FI$criterion, BV4$criterion, BV5$criterion, BV5FI$criterion,
BV6$criterion, BV7$criterion, BV7FI$criterion, BV8$criterion,
BV9$criterion, BV9FI$criterion, BV10$criterion)
#creates a set of row names for me to check my grep() call below
rn = c("BV1", "BV1FI", "BV2", "BV3", "BV3FI", "BV4", "BV5", "BV5FI",
"BV6", "BV7", "BV7FI", "BV8", "BV9", "BV9FI", "BV10")
#Number parameters
k = c(length(BV1$par), length(BV1FI$par), length(BV2$par),
length(BV3$par), length(BV3FI$par), length(BV4$par),
length(BV5$par),length(BV5FI$par), length(BV6$par),
length(BV7$par), length(BV7FI$par), length(BV8$par),
#build AIC table
AIC.table=data.frame(AIC = AIC, rn=rn, k=k, dAIC = abs(min(AIC)-AIC) ,
likg=exp(-.5*(abs(min(AIC)-AIC))))
#row.names(AIC.table)=grep("BV", ls(), value=TRUE)
AIC.table=AIC.table[with(AIC.table, order(-likg, -dAIC, AIC, k)),]
AIC.table=data.frame(AIC.table, wi=AIC.table$likg/sum(AIC.table$likg))
AIC.table
# Model average N_hat_covered estimates
# not very clean, but I wanted to show full process, need to use
# collect.models and model.table here later on
estimate <- c(BV1$Nhat, BV1FI$Nhat, BV2$Nhat, BV3$Nhat, BV3FI$Nhat,
BV4$Nhat, BV5$Nhat, BV5FI$Nhat, BV6$Nhat, BV7$Nhat,
BV7FI$Nhat, BV8$Nhat, BV9$Nhat, BV9FI$Nhat, BV10$Nhat)
AIC.values=AIC
# had to use str() to extract here as Nhat.se is calculated in
# mrds:::summary.io, not in ddf(), so it takes a bit
std.err <- c(summary(BV1)$Nhat.se, summary(BV1FI)$Nhat.se,
summary(BV2)$Nhat.se, summary(BV3)$Nhat.se,
summary(BV3FI)$Nhat.se, summary(BV4)$Nhat.se,
summary(BV5)$Nhat.se, summary(BV5FI)$Nhat.se,
summary(BV6)$Nhat.se, summary(BV7)$Nhat.se,
summary(BV7FI)$Nhat.se,summary(BV8)$Nhat.se,
summary(BV9)$Nhat.se, summary(BV9FI)$Nhat.se,
summary(BV10)$Nhat.se)
# }
# NOT RUN {
# }
# NOT RUN {
#Not Run
#requires RMark
library(RMark)
#uses model.average structure to model average real abundance estimates for
#covered area of the surveys
mmi.list=list(estimate=estimate, AIC=AIC.values, se=std.err)
model.average(mmi.list, revised=TRUE)
#Not Run
#Summary for the top 2 models
#summary(BV5, se=TRUE)
#summary(BV5FI, se=TRUE)
#Not Run
#Best Model
#best.model=AIC.table[1,]
#Not Run
#GOF for models
#ddf.gof(BV5, breaks=c(0, 25, 50, 75, 100))
#Not Run
#Density estimation across occupancy categories
#out.BV=dht(BV5, region.data, final.samples, obs, se=TRUE,
# options=list(convert.units=.01))
#Plot--Not Run
#Composite Detection Function
#plot(BV5, which=3, showpoints=FALSE, angle=0, density=0, col="black", lwd=3,
# main="Black-capped Vireo",xlab="Distance (m)", las=1, cex.axis=1.25,
# cex.lab=1.25)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab