N-hat differs w/ or w/o groups defined in process.data

posts related to the RMark library, which may not be of general interest to users of 'classic' MARK

N-hat differs w/ or w/o groups defined in process.data

Postby WiPhi » Mon Jan 23, 2017 3:45 pm

Hello All,
Great forum and interface here to complement the help and Mark book. I have a question that I don't think has been addressed and may not be specific to RMark, so I apologize in advance if this is not the case.

I am fitting Robust design models to estimate survival and population size with 8 years (8 secondary sessions in each year) of capture-recapture data in a long-lived mammal population (n=385 unique individuals), starting with simple dot models to error check (and no emigration).

I find very different population estimation results depending on whether I define age and sex groups in the process.data() stage, while estimates on p and S are stable and reasonable. To confirm this was not an error arising from the data, I coded a simulation to mimic my data that may be reinventing the wheel, but perhaps others might like to use it (also allows simulation of population trend).

With this simulated data, I process data once defining sex and age groups in the process.data and once without defining groups and then fit the same model using the different processed data and designed data lists. The model fit has no group effects on any parameter. The simulation shows the estimates of S and p are consistent with simulated values but N-hats are quite different. When groups are defined in the process data (but not modeled), the SE in the f0 parameters are tiny and the subsequent N-hats appear to be simply the number of marked individuals in each year.

Code: Select all
#set start pop size, pop growth rate, survival, study design parameters
nprimary<-8
nsecondary<-8
Sreal=0.75
lambda=1.025
Npopt0=150
p=0.25

# create container for all individuals in population,
# first will simulate actual population then account for imperfect detection
Nmatrix<-matrix(0,ncol=nprimary*nsecondary, nrow=Npopt0)
## all ones for first year with starting populatoin
Nmatrix[,1:nsecondary]<-1

# build a contain for ages at t0 (negative ages will arise as pop grows after t0)
ageatt0<-vector(length=Npopt0)
#  estimate number of recruits from t-1 and assign these age 0
Nyoung<-floor(ceiling(Npopt0/lambda)*(lambda-Sreal))
ageatt0[1:Nyoung]<-0
# assign remainder of starting pop random ages from 1 to max age
set.seed(6)
ageatt0[(1+Nyoung):Npopt0]<-floor(runif(Npopt0-Nyoung,1,14))

set.seed(6)
for (i in seq(1+nsecondary,(nprimary*nsecondary),by=nsecondary)){
   nrows<-nrow(Nmatrix)
   nind<-sum(Nmatrix[,i-nsecondary])
   for (j in 1:nrows){
      if(Nmatrix[j,i-nsecondary]==0){
         Nmatrix[j,i:(i+nsecondary-1)]<-0
      }
      else{
         Nmatrix[j,i:(i+nsecondary-1)]<-rbinom(1,1,Sreal)
      }
   }
   recruits<-floor(nind*(lambda-Sreal))
   Nmatrix<-rbind(Nmatrix,matrix(0,nrow=recruits,ncol=nprimary*nsecondary))
   Nmatrix[(1+nrows-recruits):(nrows+recruits),i:(i+nsecondary-1)]<-1
   
   ageatt0<-c(ageatt0,rep(-((i-1)/8),recruits))
}

pmatrix<-matrix(NA,ncol=nsecondary*nprimary,nrow=nrow(Nmatrix))
for (i in 1:length(pmatrix)){
   pmatrix[i]<-rbinom(1,1,p)
}

## now multiply the detection matrix by the live matrix
chmatrix<-pmatrix*Nmatrix

## as a result of imperfect detection, some animals will never by detected and they must me removed from capture histories
datamatrix<-chmatrix[rowSums(chmatrix)!=0,]
ageatt0<-ageatt0[rowSums(chmatrix)!=0]

## age at first detection
#container for age at first detection
yearoffirst<-function(x) (floor(min(which(x==1))/8))
ageatfirst<-ageatt0+apply(datamatrix,1,yearoffirst)

## container for gender, random gender assignment based on proportion male in pop. sexmale is an indicator 1=male, 0=female
propmale<-0.5
set.seed(6)
sexmale<-rbinom(nrow(datamatrix),1,propmale)

#build data.frame, ch's must be character         
chdf<-data.frame(ch=apply(datamatrix,1,paste,sep="",collapse=""),ageatfirst,sexmale)
chdf$ch<-as.character(chdf$ch)

library(RMark)
#define time intervals of robust design
time.intervals=rep(c(rep(0,nsecondary-1),1),nprimary)[-nprimary*nsecondary]

## process data with groups
sim.process.rd<-process.data(chdf,
                        model="Robust",
                        begin.time=2008,
                        time.intervals=time.intervals,
                        age.var=1,
                        groups=c("ageatfirst","sexmale"),
                        initial.ages=sort(unique(chdf$ageatfirst))
)
## process data without groups
sim.process.rd.no.grps<-process.data(chdf,
                        model="Robust",
                        begin.time=2008,
                        time.intervals=time.intervals
)

## design.data with groups
sim.ddl.rd<-make.design.data(sim.process.rd)
## design.data without groups
sim.ddl.rd.no.grps<-make.design.data(sim.process.rd.no.grps)

## model parameters
S.dot=list(formula=~1)
p.dot=list(formula=~1,share=T)
GammaDoublePrime.no.dot=list(formula=~1,fixed=0)
GammaPrime.no.dot=list(formula=~1,fixed=1)


#Robust model, groups in design data
test.model<-make.mark.model(parameters=list(S=S.dot,p=p.dot, GammaDoublePrime=GammaDoublePrime.no.dot,GammaPrime=GammaPrime.no.dot),data=sim.process.rd,ddl=sim.ddl.rd)
test.results<-run.mark.model(test.model)
#Same  model design, no groups in design data
model.no.grps<-make.mark.model(parameters=list(S=S.dot,p=p.dot, GammaDoublePrime=GammaDoublePrime.no.dot,GammaPrime=GammaPrime.no.dot),data=sim.process.rd.no.grps,ddl=sim.ddl.rd.no.grps)
no.grps.results<-run.mark.model(model.no.grps)


# N-Hat by session, groups defined
lapply(lapply(test.results$results$derived$'N Population Size',matrix,nrow=nprimary),rowSums)
#N-hat by session, no groups
lapply(lapply(no.grps.results$results$derived$'N Population Size',matrix,nrow=nprimary),rowSums)

WiPhi
 
Posts: 27
Joined: Thu Mar 31, 2016 6:15 pm

Re: N-hat differs w/ or w/o groups defined in process.data

Postby jlaake » Mon Jan 23, 2017 7:00 pm

Any time you define groups, you should have f0 defined by group*session otherwise you are saying the number missed is constant across group which is not likely unless your groups are the same size and p is constant across groups. Also it is a lot of groups given the size of your data set. What is ageatfirst?
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: N-hat differs w/ or w/o groups defined in process.data

Postby WiPhi » Mon Jan 30, 2017 4:07 pm

Thanks for the good point about modelling group*session effects on f0.

It still is not clear to how to fit a null "Robust" model without group effects when group effects are defined in the process.data() list. A first step in my modelling approach is to test whether there is evidence for group effects on S. For example, ageatfirst is the age (in years) of first detection and I can explore age effects by binning ageatfirst using add.design.data() to explore the effects of age classes on S. However, I may discover that there is no age (or sex) effect and defining groups is not important. If this turns out to be the case, I was under the assumption that fitting a null model with no group effects would produce the desired result, but the different results I get when fitting a null model using a process.data() list with groups defined versus fitting a null model to a process.data(list) without groups defined has me perplexed. Perhaps I do not understand how this model fundamentally works?

Thanks
WiPhi
 
Posts: 27
Joined: Thu Mar 31, 2016 6:15 pm


Return to RMark

Who is online

Users browsing this forum: No registered users and 10 guests