### Alaska Brown Bear Abundance 2013, 2014, 2015 ### Huggins Closed Captures Data Type with Random Effects ### RMark model HugginsRE ### Apryle Craig, April 2017, apryle@uw.edu require(RMark) setwd("~/Education/PhD/MARK bear/markdata/") # 2013 bear data bears2013 = read.csv("2013 bears v2.csv", sep = ",", colClasses=c("character", "factor", "factor", "factor", "factor")) ?read.csv head(bears2013) str(bears2013) bfrun=function(){ #can use loop, then once finalized. # HugginsRE has parameters available p, c, sigmap # sigmap models the individual heterogeneity of the p's. For sigmap = 0, you obtain the same likelihood as the basic Huggins data type, so the likelihoods of the random effects data type are compatible with the basic model, and thus AIC can be used to compare models. # we would like abundance estimates by sex (n=2), stream (n=6), salmon, and neighborhood (n=2) x.proc=process.data(bears2013, model="HugginsRE", groups=c("sex", "stream", "salmon")) x.ddl=make.design.data(x.proc) # Null model, constant p & c sharing 1 parameter (one detection estimate across all) # no individual random effect # p(~1)c()sigmap(~1) # 1 parameter, p=c p.shared=list(formula=~1,share=TRUE) mysigmap=list(fixed=0) mod.1.2013=mark(x.proc, x.ddl, model.parameters=list(p=p.shared, sigmap=mysigmap), invisible=FALSE) mod.1.2013$results$derived # looks reasonable # Time model, constant p & c sharing 1 parameter (one detection estimate across all) # No individual random effect # p(~time)c()sigmap(~1) # 3 parameters, p1, p2=c2, p3=c3 p.shared=list(formula=~time,share=TRUE) mysigmap=list(fixed=0) mod.time.2013=mark(x.proc, x.ddl, model.parameters=list(p=p.shared, sigmap=mysigmap), invisible=FALSE) mod.time.2013$results$derived # looks reasonable # p(sex) with p & c sharing 1 parameter (one detection estimate across all) # no individual random effect # p(~sex)c()sigmap(~1) # 2 parameters, pmale=cmale, pfem=cfem p.shared=list(formula=~sex,share=TRUE) mysigmap=list(fixed=0) mod.sex.2013=mark(x.proc, x.ddl, model.parameters=list(p=p.shared, sigmap=mysigmap), invisible=FALSE) mod.sex.2013$results$derived # looks reasonable # p(stream) with p & c sharing 1 parameter (one detection estimate across all) # no individual random effect # p(~stream)c()sigmap(~1) # 6 parameters, pmale=cmale, pfem=cfem p.shared=list(formula=~stream,share=TRUE) mysigmap=list(fixed=0) mod.stream.2013=mark(x.proc, x.ddl, model.parameters=list(p=p.shared, sigmap=mysigmap), invisible=FALSE) mod.stream.2013$results$derived # looks reasonable return(collect.models()) } bf.out2013 <- bfrun() bf.out2013 # get abundance estimates from a model # there will be 2 x 6 x groups estimates for sex (n=2), stream (n=6) # decoding available in text file mod.1 # INPUT --- glabel(1)=sexFemale.streamEagle.NHOOD0; # INPUT --- glabel(2)=sexMale.streamEagle.NHOOD0; # INPUT --- glabel(3)=sexFemale.streamHansen.NHOOD0; # INPUT --- glabel(4)=sexMale.streamHansen.NHOOD0; # INPUT --- glabel(5)=sexFemale.streamHappy.NHOOD0; # INPUT --- glabel(6)=sexMale.streamHappy.NHOOD0; # INPUT --- glabel(7)=sexFemale.streamBear.NHOOD1; # INPUT --- glabel(8)=sexMale.streamBear.NHOOD1; # INPUT --- glabel(9)=sexFemale.streamWhitefish.NHOOD1; # INPUT --- glabel(10)=sexMale.streamWhitefish.NHOOD1; # INPUT --- glabel(11)=sexFemale.streamYako.NHOOD1; # INPUT --- glabel(12)=sexMale.streamYako.NHOOD1; mod.1.2015$results$derived #null model mod.3.2015$results$derived #top model