www.phidot.org

Hello,

I'm fitting multistate CMR models in RMark to estimate seasonal mortality of marked individuals based on live-recaptures and dead recoveries following the approach of Gauthier and Lebreton (2008). This method fixes survival parameters so that transition probabilities from live states to dead states become mortality probabilities. I have 2114 individuals in the analysis that covers 20 intervals (10 years; each year has a breeding and non-breeding season). I'm fitting 420 models comprising combinations of detection and psi formulas using the approach in the RMark appendix (i.e., create.model.list followed by mark.wrapper). Everything is working at this point, except that my likelihoods commonly get 'hung-up' in a local minima when fitted. I've started using simulated annealing and that helps immensely, but I still end up having to fit some models more than once to get beyond the local minima. I would like to subset the mark.list results from the model fitting function based on deviances so I can select models that need to be re-fit based on being in a local minima and then (possibly) use 'initial' to provide starting values from the initial model fitting to help the likelihood move beyond the local minima. Long story short, I can't figure out how to successfully subset mark.list in a way that allows me to maintain the subset as a mark.list (and not just a mark.list$model.table) so I can fit the subset of models and provide the initial subset results as starting values.

A bit of script to help explain:
Code: Select all

Code: Select all

lesc.models<-function()
{...formula statements for S [fixed], Psi, and p here...
 cml<-create.model.list("Multistrata")
  results<-mark.wrapper(cml, data=lesc.process, ddl=lesc.ddl, parm.specific=T, output=F, options="SIMANNEAL")
  return(results)
}

#Run the function to create the models and store the results in lesc.model.results, which is a marklist
lesc.models.results <- lesc.models()

#Identify the models in the mark.list model.table based on deviance values
re.fit <- as.numeric(row.names(lesc.models.results$model.table[lesc.models.results$model.table$Deviance>7000,]))


From here I'd like to subset the mark.list 'lesc.models.results' based on the models identified in the 're.fit' object and run the model fitting function again based on the subsetted list of models using 'initial' to provide starting values.

Any help would be greatly appreciated.

You can subset a list with [] which returns a list. Don't use [[]] as that is used to extract a single list element. It is not returned as a list and can thus only return a single element.

Code: Select all

newlist=lesc.models.results[re.fit]
newlist$model.table=model.table(newlist)


You may also have to set the class as a marklist. Alternatively look at ?remove.mark

Jeff,

Thanks, that helped. I used 'remove.mark' to drop models from the marklist that were not caught in a local minima based on deviance values. I noted in your help file for create.model.list you mention 'creat[ing] all the combinations and then delete the ones from the dataframe that do not make sense...and re-run the function'. Therefore, I next tried to run the model fitting function by calling the subsetted marklist (i.e., replaced 'cml' with 'lesc.models.results.2' in the above code). Based on help files for create.model.list and mark.wrapper I thought I was on the right track of calling the model.table from lesc.models.results.2 as it provides a data.frame containing columns for each formula specification for S, p, and Psi. I tried this a number of ways (i.e., subsetting the model.table for the first 3 and 4 columns, calling the entire model.table, and calling the marklist 'lesc.models.results.2' directly) but didn't succeed (received a different error for each attempt). Alternatively, maybe I'm painting myself into a corner by trying to fit the subsetted model list via the model fitting function; perhaps there is a more direct way to specify a marklist for fitting? Any ideas/insights are greatly appreciated.

Thanks again.

Jeff

The documentation was saying you could do something like the following simple example to exclude some models. It removes the Phi(time)p(time) model.

Code: Select all

function()
{
   Phi.1=list(formula=~1)
   Phi.2=list(formula=~time)
   p.1=list(formula=~1)
   p.2=list(formula=~time)
   cml=create.model.list("CJS")
   cml=cml[-4,]
   dp=process.data(dipper)
   ddl=make.design.data(dp)
   mark.wrapper(cml,data=dp,dd=ddl)
}


The mark.wrapper code is expecting a dataframe with character names of variables containing the formulas for each parameter. (e.g., "Phi.1") and that is then used to extract the parameter specification. mark.wrapper does't have any facility to pull out the formula combinations from an existing model list although it could certainly be done with added code.

Jeff,

That makes sense, and I was able to follow your example to drop the models in 'cml' in the function to only fit the subset of models I selected.

Thanks!

Jeff