cjs RESULTS between rMark and library(FSA) not matching

[marna]

Hi all,
I had run a CJS model using data of "dipper" with RMARK [Package RMark] and library (FSA)[ based on the Introductory Fisheries Analysis with R, written by Derek H. Ogle; 2015, Chapman and Hall], but I got different results. Would you give me some suggestions which results are correct.

I used the "dipper data" but for Female (the r code is given).
For example: period 1 and 2 the survival rate is : 0.742 using FSA, whereas Rmark gives "0.718"
I am not sure which results are the right . Furthermore, Rmark provides the survival rate and detection rate of the last detection time (period 6 to 7: 0.7655); I think it should have the value "NA". The FSA provides "NA" for the last detection time ( See below). I have attached the code for your reference.

> cut.op$survivals
estimate stderr
period 1 -> 2 0.7428571 0.23722150
period 2 -> 3 0.4468410 0.09828285
period 3 -> 4 0.4538126 0.08153503
period 4 -> 5 0.6404243 0.08323359
period 5 -> 6 0.6280454 0.08105155
period 6 -> 7 NA NA

Whereas RMark gives
1980 1981 1982 1983 1984 1985
1980 0.718 0.435 0.478 0.626 0.599 0.766
1981 0.435 0.478 0.626 0.599 0.766
1982 0.478 0.626 0.599 0.766
1983 0.626 0.599 0.766
1984 0.599 0.766
1985 0.766

### following code I used for the above result
library(FSA)
library(RMark)
rm(list=ls(all=TRUE))

data(dipper)
head(dipper)

########################
## based on RMark
#########################
dipper.process=process.data(dipper,model="CJS",begin.time=1980,groups="sex")
dipper.ddl=make.design.data(dipper.process)

Phi.dot=list(formula=~1)
p.dot=list(formula=~1)
Phi.time=list(formula=~time)
Phi.sex=list(formula=~sex)
p.time=list(formula=~time)
p.Time=list(formula=~Time)

dipper.phi.time.p.time= mark(dipper.process,dipper.ddl,model.parameters=list(Phi=Phi.time,p=p.time))
summary(dipper.phi.time.p.time)

#################################
###### this is based on FSA, see book (introductory fisheries Analysis with R)
##########################################
DIPPER<-separate(dipper, ch, into = c("a1", "a2", "a3", "a4", "a5", "a6", "a7"), sep = c(1,2,3,4,5,6), remove = F)
head(DIPPER)
DIPPER_Female<-subset(DIPPER, DIPPER$sex=="Female")
head(DIPPER_Female)
Test1<-data.frame(subset(DIPPER_Female, select=c(-9)))
head(Test1)
summary(Test1)

# convert factor to numeric
Test1[,2:8] <- lapply(Test1[,2:8], function(x) as.numeric(as.character(x)))
summary(Test1)

#################################
###### this is based on library (FSA), see book (introductory fisheries Analysis with R)
##########################################
cut.ch1<-capHistSum(Test1,cols2use=-1)
cut.ch1$methodB.top
data.frame(cut.ch1$caphist)

cut.op<-openp(Test1[,-1])
cut.op$survivals
cut.op$capture.prob

[jlaake]

I don't know FSA so I can't help there but rather than doing incomplete comparisons I suggest that you spend more time reading Cooch and White book. If you did you would see that for phi(time)p(time) the parameters for the last occasion/interval are not separately estimable. RMark/MARK give estimates but there are an infinite set you could construct.

It looks to me like you used all of the data for RMark run but only Female data with FSA. Can you see the problem with that?

[marna]

Yes you are right, I had used both groups (Male and Female) in RMark but only one group in library (FSA). After your response, I run the model using only Female data, the results are now consistent between Rmark and FSA. But I am not sure why the results of the same group (female) are different when we are using the two groups (male and female) in the model in comparison with the model using only one group? does the survival rate vary if we use several groups in the model?
Similarly why still I am getting the survival probability of the last point, is it not supposed to be "NA"?
Thanks,

[marna]

Hi jlaake,
I realized that you mentioned in your previous response about why we got NA in the output. But still wondering why different results using only one group vs. two/three groups.

[jlaake]

You used all of the data and told RMark to estimate survival by time - only. So it uses all of the data (both males and females) to estimate survival. So why would you expect it to be different for Females and Males? Now if you used sex+time then you would get a difference but it would only match FSA if you fit the interaction model of sex*time. Adding groups argument to the process.data only allows you to fit different models by the grouping variable (eg sex) but does not require all models to fit the sexes separately. That would be something you would want to evaluate with model selection.

Now I'm trying to avoid being rude but I have told you in previous messages that you need to read Cooch and White and it doesn't seem like you have. So until you can tell me honestly that you have read at least chapters 1-7 and have gone through the examples, I will not answer any more of your questions.