Error RMark Multi Strata Analysis

[Jen]

I am getting an error when running multi-strata models in RMark. I updated to the most recent versions of R 2.13.0, RMark 2.0.1 and Mark_INT.EXE 23/04/2011, MARK.EXE 23/04/2011, which hasn’t resolved the problem. OS: Windows Vista

The data has been modelled successfully in Mark so I am happy that it is set up correctly. The .inp file consists of over 1100 capture histories, with each consisting of 76 capture events. It contains 4 stratum and sex (m/f) as a group. I am able to convert.inp, process.data and make.design.data however when I try to run any model it errors after the model has run for a few moments.
An example of one of my R Codes:
CaptureData=convert.inp("Capture.inp",group.df=data.frame(sex=c("Male","Female")),use.comments=FALSE)
run.mstrata=function()
{
mstrata.processed=process.data(CaptureData,model="Multistrata",groups="sex")
mstrata.ddl=make.design.data(mstrata.processed)
Psi.sex.plus.stratum=list(formula=~sex+stratum)
p.time=list(formula=~time)
S.sex.plus.stratum=list(formula=~sex+stratum)
model.list=create.model.list("Multistrata")
mstrata.results=mark.wrapper(model.list,data=mstrata.processed,ddl=mstrata.ddl)
return(mstrata.results)
}
mstrata.results=run.mstrata()

I get the visual Fortran run-time error:

Forrtl:severe (161) Program Exception- array bounds exceeded

Which then results in the model failing to run and this error in R:

Error in extract.mark.output(out, model, adjust, realvcv)

I’ve looked in the output files and there’s no obvious error message, some model outputs even contain AIC values etc (I can see on the forum that this error has occurred before but not in the same context).

Any suggestions would be greatly appreciated.

Thank you

Jen

[jlaake]

Jen-

Thanks for taking the time to update to make sure that was not the problem. Glad Bret got you squared away this morning. The error you are getting is from mark.exe and not from RMark. My guess is that mark is failing when it creates the real v-c matrix. You could see that if you set invisible=FALSE in mark.wrapper. A window will open and it will show you the running of mark.exe. If it fits the model and computes the beta v-c matrix and then fails, it is occurring when it is creating the real v-c matrix. Or you can look at the output file it created in the directory to see how far it got. This is an issue that Gary would have to deal with but there are some ways to avoid it. But first, is sex+stratum for Psi the model that you want to fit? More than likely you want

Psi.sex.plus.stratum=list(formula=~sex+stratum:tostratum,remove.intercept=TRUE)

Including only stratum in the model means that the probability that I transition to other states only depends on the stratum that I'm in, so if you have 4 states A,B,C,D then

Psi(AB)=Psi(AC)=Psi(AD)
Psi(BA)=Psi(BC)=Psi(BD)
Psi(CA)=Psi(CB)=Psi(CD)
Psi(DA)=Psi(DB)=Psi(DC)

If that is what you want then the model you wrote is fine with additive sex differences. The model you wrote would have 5 parameters. If what you intended was different transition probabilities from each state to every other state plus an additive sex difference (13 parameters), then you want to use the model I specified. The stratum:tostratum is used because there is not complete crossing of the stratum and tostratum factors (ie a different one is missing in each), unless you set subtract.stratum to the same value for each stratum. When you use the : operator with 2 factor variables then you don't need an intercept and -1+stratum:tostratum would work except when you are adding another factor variable like sex. In that case you have to use the remove.intercept argument.

Now none of that addresses your problem. The basic issue is that RMark is creating too many real parameters and mark.exe is failing to provide sufficient array sizes for it to handle that many for the v-c matrix of the real parametrers. You can reduce the number of real parameters in two ways:
1) set pim.type="time" for Psi

mstrata.ddl=make.design.data(mstrata.processed,parameters=list(Psi=list(pim.type="time")))

2) If that is still not a sufficient reduction you'll need to use link="logit" on the formula for Psi

Psi.sex.plus.stratum=list(formula=~sex+stratum:tostratum,remove.intercept=TRUE,link="logit")

The default is "mlogit" which doesn't need a penalty to make sure the sum of the transition probabilities =1 whereas logit does. This makes a difference because RMark will then simplify the PIMS (ie reduce the number of real parameters to the minimum needed for the formula). This can't be done with the mlogit link due to the way that mark.exe handles the summing for calculation of mlogit parameters.

Make sure that the model is converging by starting with a simple model and then using the simple model for starting values of more complex models. See the initial argument for mark function in RMark. Also, consider using simmulated annealing by using options="SIMANNEAL" in mark. All arguments of mark can be used with mark.wrapper.

Hope this helps.

--jeff

[bacollier]

Jen,
I am running some multi-strata models as we speak and I fought this battle with Jeff and Gary's help a while back, I had about 1100 individuals, 10 sampling occasions, 2 sex, 2 age, and 16 unique strata and I was estimating something like 60000 parameters initially in MARK. Start with the pim.type="time" Jeff suggested, that should reduce the number of parameters enough to get a S(.), Psi(.), p(.) model to run initially. Once I did that I did not have to do Jeff's 2 to get all the rest of mine running (although I am not trying to run a model with 76 encounter occasions so you might).

Also, I assume that all your transitions are occurring, e.g., you have individuals that transition between a couple of strata at least once? If not you can probably adjust your design data to fix that strata transitions so they are not estimated.

Also, Jeff is right as you will probably want to use simulated annealing in the future, but stay away from that right now until you get your basic models up and running. I started a simple S(.), Psi(.), p(.) model on my multistrata data using simulated annealing on Tuesday (57 hours and counting) and it is still running on my laptop; when I don't use simulated annealing the model took about 8 hour to work its way though. You will have more parameters than I had for it too weed through, so it might be Christmas before yours finished:)

[Jen]

Thank you for your clear explanations. I used the pim.type="time" code which has prevented the model errors so far.

I will definately be taking on board your other suggestions when running more complex models.

Thanks for a great forum.

Jen

[bacollier]

Thank you for your clear explanations. I used the pim.type="time" code which has prevented the model errors so far.

I will definately be taking on board your other suggestions when running more complex models.

Thanks for a great forum.

Jen

Jen,
Since we are working on the same general thing right now, I have another tip (credit to Jeff for suggesting a better way to do it than the approach I was using).

After you run your model(s), be sure you save the workspace after each run, that way you don't have to go back and re-run the models again.

So, after your mark() call, use save.image to save the workspace. I use a working directory, so mine is save.image("F:\\Bret Research\\SwaziBat\\.Rdata"). I was wrapping save.image() into a function, but Jeff suggested I use a call to save.image() after each model run, so each one is saved to the .Rdata file, which is really useful for long-running individual models like I think you have. If you do File-->Save Workspace, the location of where you want to save it will pop up in your R window appropriately formatted, which is nice.

Just trying to help you not make the mistakes I have already made.

Bret

[Eldar]

Hi Jen,
I had this problem before and found two things that affect it:
1. The sequence in each strata needs to be longer than 3.
2. In my case errors occured without simmanealing. But with simannealing "on" everything works fine (though it takes time..). Running simanneal models in parallel really helps.
Eldar