www.phidot.org

Hi all,

I was goofing around with the release.gof command and was getting some unexpected answers when datasets included removals (-1's). I may be doing something wrong and would appreciate any insight. What appears to be a "dropping of removals" affects all closed form estimates and the chi squared test for GOF. Some code below should reproduce what I was seeing and should linearly produce:

1) simple mod of dipper dataset (no group variable to simplify)
2) sprinkle some removals in dipper dataset for testing
3) run both (run A and B)
4) export dataset with removals for run through Mark (run C)
5) results (in comments) from all three runs

Run C properly shows removals while Run A and B are identical. Am I doing something wrong with my RMark code?

thanks for the help,

Jack

Code: Select all



#software OS used
#RMark:2.14.2
#R:    2.14.2
#Mark: Window NT (x64) (version 6.1, Build 7601)
#Windows 7

#DATA---------------------------------------------------------------------------
library(RMark)

#slightly modified dipper dataset (dropping group variable)
data(dipper)
dipper.processed=process.data(data.frame('ch' = dipper$ch, stringsAsFactors = F)
                              ,model="CJS")
str(dipper.processed)

#primitive random "addition of removals" for version two
#generate random removals
#i <- which(rbinom(n=nrow(dipper.processed$data), size=1, prob=0.10) == 1)
i <- c(5,13,20,39,40,42,57,58,60,64,93,95,97,99,100,121,137,138,139,158,164,
       168,176,181,193,197,200,204,217,219,221,223,239,241,242,243,248,262,
       272,273,274,278,287,290,294)
dipper.processed2 <- dipper.processed
dipper.processed2$data$freq[i] <- -1

#RUNS---------------------------------------------------------------------------
#run both versions
original        <- release.gof(dipper.processed,  view = T)
random.removals <- release.gof(dipper.processed2, view = T)

#run A
#original run (through RMark)===================================================
#            Observed Recaptures for Group 1
#                        Group 1
#  i   R(i)                     m(i,j)                    r(i)
#            j=  2      3      4      5      6      7
# 1     22      11      2      0      0      0      0      13
# 2     60             24      1      0      0      0      25
# 3     78                    34      2      0      0      36
# 4     80                           45      1      2      48
# 5     88                                  51      0      51
# 6     98                                         52      52
#
#Chi.square df      P
#TEST2     9.4797  4 0.0502
#TEST3     3.7547  8 0.8786
#Total    13.2343 12 0.3522

#run B
#random.removals run (through RMark)============================================
#                         Observed Recaptures for Group 1
#
#                        Group 1
#  i   R(i)                     m(i,j)                    r(i)
#            j=  2      3      4      5      6      7
# 1     22      11      2      0      0      0      0      13
# 2     60             24      1      0      0      0      25
# 3     78                    34      2      0      0      36
# 4     80                           45      1      2      48
# 5     88                                  51      0      51
# 6     98                                         52      52
#
#      Chi.square df      P
#TEST2     9.4797  4 0.0502
#TEST3     3.7547  8 0.8786
#Total    13.2343 12 0.3522

#run C
#random.removals run (through Mark**)===========================================
#export input file for Mark run
write.table(data.frame('ch' =
                      paste(dipper.processed2$data$ch,
                      dipper.processed2$data$freq, ";")),
                      file.choose(), row.names = F, col.names = F, quote = F)

#       Observed Recaptures for Group 1
#                   Group 1
# R(i)                     m(i,j)                    r(i)
#       j=  2      3      4      5      6      7
#  19*      11      2      0      0      0      0      13
#  54*             24      1      0      0      0      25
#  72*                    34      2      0      0      36
#  76*                           45      1      2      48
#  81*                                  51      0      51
#  89*                                         52      52
#
#Goodness of Fit Results (TEST 2 + TEST 3) by Group
#
#        Group  Chi-square   df   P-level
#       -----  ----------  ----  -------
#         1      **14.1653**    12    **0.2903**






Yes. freq is stored under dipper.processed. dipper.processed$freq and not dipper.processed$data$freq. You created one there but it was ignored and used dipper.processed$freq

It would be more prudent to add freq to dipper and then process because that will work with groups as well. When there are groups, freq is a matrix.

Try this:

Code: Select all

> data(dipper)
> dipper.processed=process.data(data.frame('ch' = dipper$ch, stringsAsFactors = F)
+                               ,model="CJS")
>
> freq=rep(1,nrow(dipper))
> i <- c(5,13,20,39,40,42,57,58,60,64,93,95,97,99,100,121,137,138,139,158,164,
+        168,176,181,193,197,200,204,217,219,221,223,239,241,242,243,248,262,
+        272,273,274,278,287,290,294)
> freq[i] <- -1
>   
> dipper.processed2=process.data(data.frame('ch' = dipper$ch, freq=freq, stringsAsFactors = F)
+                               ,model="CJS")
>
> release.gof(dipper.processed)
RELEASE NORMAL TERMINATION
      Chi.square df      P
TEST2     9.4797  4 0.0502
TEST3     3.7547  8 0.8786
Total    13.2343 12 0.3522
>
> release.gof(dipper.processed2)
RELEASE NORMAL TERMINATION
      Chi.square df      P
TEST2     9.4797  4 0.0502
TEST3     4.6856  8 0.7906
Total    14.1653 12 0.2903
>
>


Your posting did make me stumble onto a problem in the current version of RMark which actually attempts to run mark.exe instead of rel_32.exe unless you have MarkPath specified. If anyone has been having a problem with release.gof and needs a new version let me know. Unfortunately github no longer allows posting in the download section, so I can't post temp versions there between patches on CRAN.

I see, so the data must be run through process.data and, in order to include removals in accounting, the data must include the field 'freq' or else defaults to '1' for all entries. Otherwise, dipper.processed$freq must be manipulated directly.

Thanks for the quick response, Jeff!

Jack

Hi all,

I've been using RMark to run CJS models and after performing the goodness-of-fit test with RELEASE (release.gof), I discovered that one of my data sets has some overdispersion issues. I would like to estimate c-hat using the parametric bootstrap GOF test. I've been referencing Jeff Laake's excellent guide (RMark Workshop Notes) thus far and the instructions it gives regarding this topic are to export the mark data and model and do the test in MARK. Just wondering if anyone has written code to do the test in R?

Thanks!
Rebecca

I have not. I'll let others respond. Glad the workshop notes are useful for you. --jeff

For CJS models we've found that the bootstrap results are often quite similar to median c-hat, and both provide much less evidence of lack of fit than RELEASE. Our data sets have many thousands of encounter histories with many occasions, so your experience may differ. But exporting to MARK and running the median c-hat procedure is straightforward and would be a good thing to try if you are concerned about what RELEASE is telling you.

- Dave