RN occupancy model in RMark - convergence problem solved

[matobler]

I had an offline exchange with Jeff about some problems I ran into with the OccupRNPoisson model and groups in RMark. Several models would not converge when groups were defined for Lambda. The same models worked fine for regular Occupancy models and worked when using the RN model in PRESENCE. The solution to the problem was to run a simple model that would converge and use that to define the initial values for the more complex models. That way all models converged and the output from Mark and PRESENCE was identical. Bellow the R code.

Note: Groups won't work for occupancy models with the Oct 1 version of mark.exe. I used the July 08 version and everything worked fine.

Code: Select all

#  create an initial model
initial=mark(occdata,model="OccupRNPoisson", groups=c("Habitat","Site","Survey"),
   model.parameters=list(r=list(formula=~1), Lambda=list(formula=~1)))

run.rnocc=function(){
   occdata.process=process.data(occdata,model="OccupRNPoisson", groups=c("Habitat","Site","Survey"))
   occdata.ddl=make.design.data(occdata.process)
# Define r models
   r.dot=list(formula=~1)
   r.site=list(formula=~Site)
   r.survey=list(formula=~Survey)
# Define Lambda models
   Lambda.dot=list(formula=~1)
   Lambda.habitat=list(formula=~Habitat)
   Lambda.site=list(formula=~Site)
   Lambda.sitehabitatint=list(formula=~Site*Habitat)
   Lambda.sitehabitat=list(formula=~Site+Habitat)
   Lambda.survey=list(formula=~Survey)
   Lambda.surveyhabitatint=list(formula=~Survey*Habitat)
   Lambda.surveyhabitat=list(formula=~Survey+Habitat)
# Create model list
   cml=create.model.list("OccupRNPoisson")
# Run and return marklist of models
  return(mark.wrapper(cml,data=occdata.process,ddl=occdata.ddl, initial=initial))
}

rnocc.results<-run.rnocc()
print(rnocc.results)

[jlaake]

As a follow-on to the message from Mathias, I'm going to jump up on the soapbox for a moment. I run into this issue routinely with both MARK and DISTANCE.

It is important to keep in mind that the algorithms used for optimization are not compleley bullet-proof. If the optimization in MARK (or any other optimization code) begins at poorly defined initial values for the parameters it is always possible that it will NOT converge to the MLEs. If you suspect that a model is not converging to the correct values then use another model (that did converge) to provide starting values for the other models. This can be done with the MARK interface or with RMark as Mathias showed.

In Mathias' example it was obvious that the model had not converged proprerly because the estimates were obviously bogus. Sometimes it is less obvious but another way to examine this is to compare nested models. Adding parameters should always improve (make it bigger) the log likelihood or conversely the negative log likelihood should get smaller. But this comparison only works with nested models which means if you have 2 models A and B, A is nested in B if B contains all the parameters in A and more. This comparison does not hold if the models are not nested. It is quite possible for a model with say 10 parameters to have a poorer fit (based on likelihood) than one with say 5 parameters if the models are not nested.

The above should not be confused with comparisons of AIC which assumes that all of the models converged and penalizes the comparison based on the number of parameters in the model. It should be obvious why that is necessary because as I said adding parameters to a model will always provide an improvement in the likelihood value (just like an R-sq in regression).

As always make sure your results make sense like Mathias did. Don't treat the software as a black box and relegate the responsibility for correct results to the software.