Time interval problems

[Cmichel]

Hi,
I have some experience using Program Mark, and I am finally moving on to using RMark for some more complex modeling. However, I have already hit a roadblock and I have spent way too much time already trying to solve it.

I have attempted to add time intervals to some more complex models with group variables, but from the real survival estimates, I can tell that the time intervals were not applied correctly. Therefore I simplified the models and removed the group effect to just try and figure out how to add time intervals correctly.

I used the information I found in the RMark conference notes, the Mark Book and on the phidot website, and put together this script:


## load the CH file

lfc <- convert.inp("All_years_river_reaches_short_no_group.inp")

## set up the basic model structure. Time intervals are kilometers (I am looking at survival through space rather than through time)

lfc.process <- process.data(lfc, model="CJS", begin.time=0, time.intervals=c(14.057, 11.176, 36.414, 34.502, 9.628, 22.531, 26.612, 37.367, 16.677, 83.101))

lfc.ddl <- make.design.data(lfc.process)

## set up the various model specifications. Null model

Phi.dot <- list(formula=~1)
p.dot <- list(formula=~1)
lfc.Phi.dot.p.dot <- mark(lfc.process, lfc.ddl, model.parameters=list(Phi=Phi.dot, p=p.dot))

# Phi and p time dependent

Phi.t <- list(formula= ~time)
p.t <- list(formula= ~time, fixed = list(index=11,value=0.77))
lfc.Phi.t.p.t <- mark(lfc.process, lfc.ddl, model.parameters=list(Phi=Phi.t,p=p.t))

--------------------------------------------------

When run, the null model seems to run fine, but the hiccup comes when RMark tries to run the time dependent model. The error message says:


Error in if (x4 > x2) { : argument is of length zero

********Following model failed to run : Phi(~time)p(~time) **********

Error in `row.names<-.data.frame`(`*tmp*`, value = c("Phi:(Intercept)", :
duplicate 'row.names' are not allowed
In addition: Warning message:
non-unique values when setting 'row.names': ‘Phi:time105.777’, ‘Phi:time154.92’, ‘Phi:time192.287’

---------------------------------------------------

I then tried to tweak the script in many different ways, and the closest I’ve come to it running correctly is by changing the time dependent model to run with continuous time (~Time) rather than time as a factor (~time). So, same script as above, with the only modification in the time dependent model:

Phi.t <- list(formula= ~Time)
p.t <- list(formula= ~Time, fixed = list(index=11,value=0.77))

This script actually runs without error, and gives me real survival estimates that seemed to be normalized for time, but are still off from the “correct” parameters I get running this in Mark itself. Plus, if I take RMark’s estimates and raise them to the power of the time intervals, they should equal the survival estimates using the default time intervals, but in this case, they don’t.

Could anyone give me any advice on what I am doing wrong? It seems like a super simple problem, but I can’t figure it out!

Thanks!!!

[jlaake]

Are you using the most recent version of RMark? It is 1.9.6. I had some problems with unequal time intervals that I fixed not that long ago and it generated the error that you listed in your post. Let me know and I'll help you pursue this. --jeff

[cooch]

Hi,
I have some experience using Program Mark, and I am finally moving on to using RMark for some more complex modeling...

RMark doesn't let you do more complex modeling than MARK. It provides you with convenient - and powerful - alternative ways to build models in MARK. There are no models that RMark accommodates that aren't available in MARK, since MARK in fact is doing all the estimation!

[mcmelnychuk]

After you try Jeff's suggestion, if you get the same message you could try removing the intercept to see if that would give more realistic km^(-1) survival rates. (I'm not sure if it would change things.)

Phi.t <- list(formula= ~ -1 + time)

Note that if you have more than one group of smolts (I'm guessing!) released in different locations but you want to incorporate all groups in the same analysis (e.g., for estimating common p), the distances of each river segment will differ among groups so they won't have the same "time intervals". If this is the case, one option is to incorporate those segment distances you listed as a group covariate for survival in each segment:

Phi.dist <- list(formula= ~ dist)
#where dist is the vector of segment & group-specific covariate values

This has the advantage that different covariate values can be used for different groups, but has the drawback that it won't return km^(-1) survival rates (along with SE for those "instantaneous" rates). Instead, it will just return the survival probability for each segment, which are constrained as a function of the distance of each segment.

Cheers, Mike

[Cmichel]

jlaake:
I installed the newest version of RMark (1.9.6). I was using 1.9.3
I also went ahead and installed the newest version of R.

I re ran the script and, unfortunately, the Phi(~time)p(~time) model still could not run, and produced a similar but shorter error message as last time:

"Error in if (x4 > x2) { : argument is of length zero

********Following model failed to run : Phi(~time)p(~time) ********** "

However, I also ran a Phi(~time)p(~1) model, and it did run and produced the same survival estimates as Mark produced (with the exact same model). So we are on the right track in the sense that this is the first time I have produced any time normalized survival estimates that match Mark's estimates. So I think installing the newest version of Rmark helped.
But I still need to get Phi(t)p(t) to run.

But, mcmelnychuk's suggestion worked for Phi(t)p(t)! I used Phi.t <- list(formula= ~ -1 + time) which removed the intercept, and it gave me the exact same survival estimates for Phi(t)p(t) in Mark.

Guess it was a relatively simple fix afterall. Thanks for the help!!!

cooch:
You are right, I misspoke. I am using RMark to make more complex models easier for me to run and rerun, thus making the process less error-prone and time consuming.

mcmelnychuk:

Thanks for your suggestion on removing the intercept. I had tried that before, but now that I have the newest version of RMark, it actually worked!
You are right about the smolts, I am working down in California with Bruce MacFarlane. As for different groups having different "time" intervals, we purposely released fish at monitor sites, so that those releases could also count as the first detection. This, fortunately, makes so that we have no "portions" of river reaches traveled by different groups. In other words, I should be able to use the same time intervals for all release groups, as long as the capture history is coded right.

[mcmelnychuk]

Glad it worked!

This, fortunately, makes so that we have no "portions" of river reaches traveled by different groups.

You guys get points for planning ahead (our situation was a bit messier in B.C.) If you have some groups released near river stations which are part way along the downstream migration route of earlier released groups, that provides a good opportunity to assess tag-related mortality (or other forms of mortality that might occur immediately after release, wherever the release locations may be). Email me off-list if you're interested in discussing - mmel@uw.edu

[jlaake]

I'm glad that solved your problem but I'm perplexed that the Phi(~time)p(~time) did not work and the Phi(~-1+time)p(~time) worked. They do produce different design matrices, For example with 3 phis, ~time produces the DM

1 0 0
1 1 0
1 0 1

whereas ~-1+time produces the DM

1 0 0
0 1 0
0 0 1

Sometimes the DM can affect the convergence properties of the fit even though in reality they should produce the same results. Can you look in the .out file created for the run that failed? It may show that there was convergence problems. There is a section of Appendix C in Cooch and White that lists some of the common errors that you can get with RMark and what you can do to find the cause. Your most recent error is in there. I probably need to write some more sophisticated parsing code for the output to handle these cases where MARK.EXE is not converging.

[jlaake]

One other comment. Not sure what you were doing with the ~Time model and what you compared it to with MARK but ~Time with unequal time intervals will not yield the same interval specific values of Phi for ~Time with default time intervals of 1. That will only work with ~time because there is a separate estimate for each interval. The ~Time model is a linear model which will fit differently if Time values are 1,2,3,... versus 11.3,22.4,....

In reference to the comment made by Evan, it is important to always keep in mind that RMark simply creates the input file for MARK.EXE. If there is a difference between RMark and MARK interfaces with the "same" model, then you need to look at the data, links, PIMS and the DMs and see how they differ. That is the only way the results can be different because they both use MARK.EXE for the model fitting.

regards --jeff