www.phidot.org

In fitting a model using CJSMixture (a CJS model with Pledger mixtures), I "discovered" that the correct number of pi parameters (C-1) will not be estimated with a default pi model when the number of mixtures is greater than 2. To get an estimate for pi for each class (or rather all but the last class, because the pi parameters must sum to 1), specify a model for pi that fits the mixture parameter:

Code: Select all

pi.mix=<-list(formula=~mixture)


This works because the default design data has a mixture column with C-1 values (thanks Jeff).

I am certain this has probably been spelled out before somewhere or is common knowledge to others but I thought it might deserve it's own topic for this forum. I noticed it when i saw my model (assuming a mixture of 3 classes) fitting 3 estimates for p but only one estimate for pi (rather than 2). (As an aside, I am not sure how to interpret the single pi estimate, but my guess is that the single estimate suggest equal probability of assignment to the first two classes and the remainder (1-2*pi) is the probability of assignment for the third class?)

The mixture column in the default design data for pi would also allow for models to test for heterogeneity within groups (such as sex) using a mixture*group pi model. Pledger et al. 2003 discuss this briefly but don't illustrate an example and I would be very keen to see an example of this being implement if anyone knows of one.

You are correct that the default formula for pi is ~1 and as you found you can use ~mixture to get separate estimates. You can change these values yourself by modifying the values in parameters.txt which is in the RMark subdirectory of your R library.

The values for CJSMixture are

CJSMixture pi 0 1 0 Square all logit ~1 TRUE -1
CJSMixture Phi 0 -1 1 STriang all logit ~mixture TRUE 0
CJSMixture p 1 -1 0 STriang all logit ~mixture TRUE 0

You would just need to change ~1 to ~mixture for pi. However, anytime you update RMark you'll have to change it. Creating default formula wasn't probably the best thing to do because now that I have them, I can't really change them without possibly messing up someone's code that assumes a particular default formula.

--jeff

Good solution, Jeff. The parameters.txt is a nice interface for controlling model behavior.

I would also not recommend changing the file (or making this a permanent edit in future releases) for another reason. If the default was ~mixture and a mixture of only 2 classes is fit (perhaps the most common mixture strucutre fit?), a "contrasts can be applied only to factors with 2 or more levels" error will result in model fitting. This would then require a =~1 pi model specification for the 2 class mixture model rather than simply assuming the default pi model will work.

Had forgotten but that is the reason I chose to use ~1 because mixture=2 is the default.

On a related topic, In trying to compare models that do or do not include a mixture. I can fix pi at 0 via:

Code: Select all

pi.fix=<-list(formula=~1,fixed=0)


which produces synonymous results as a CJS model without mixture (despite being fit using a call to CJSMixture...a very nice capability).

Using this "fixed at 0" pi model (along with the default, ~1,that estimates pi), a suite of mixture and no mixture models can be fit simultaneously using create.model.list() within a function allowing the assumption of mixtures to be compared directly along with other models that vary in structure for p and Phi. However, the results output shows both mixture and no mixture models as pi(~1) models in the model name, requiring some manual parsing within the model results or careful comparison of parameter counts to distinguish between model types. This can be burdensome when many models are being compared.

Any ideas on an an elegant way to fit the no mixture model and have it show up in the model results in a more conspicuous way?

Thanks in advance!

Below is an example that uses model.table with the argument model.name=FALSE which uses the name of the formula object rather than the formula.

Code: Select all

do.dip=function()
{
  dp=process.data(dipper,model="CJSMixture")
  ddl=make.design.data(dp)
  pi.intercept=list(formula=~1)
  pi.fixed0=list(formula=~1,fixed=0)
  p.time=list(formula=~time)
  cml=create.model.list("CJSMixture")
  results=mark.wrapper(cml,data=dp,ddl=ddl)
  results$model.table=model.table(results,model.name=FALSE)
  return(results)
}
data(dipper)
do.dip()


The result will look like:

Code: Select all

               model npar     AICc DeltaAICc    weight Deviance
1    pi.fixed0.p.time    8 680.8255  0.000000 0.7395603 55.80039
2 pi.intercept.p.time    9 682.9129  2.087369 0.2604397 55.80039


Thanks.worked well.