problem reproducing RMark multi-strata example

posts related to the RMark library, which may not be of general interest to users of 'classic' MARK

problem reproducing RMark multi-strata example

Postby ctsolomon » Fri Nov 19, 2021 4:31 pm

Hi-

I am having trouble reproducing one of the examples from the RMark appendix of the MARK book, and am hoping someone can suggest what I might be doing wrong.

The example I'm working on is in section C.17 (Multi-strata example); the discussion and (most of) the code begin on page C-81. It's the example of setting PsiAA=PsiBB=PsiCC=0.

In the book, the fitted estimates of Psi show up as 0 for PsiAA, PsiBB, and PsiCC, and as non-zero for the other transitions. These results are on p. C-81 to C-82.

When I try to reproduce the example, the PsiAA, PsiBB, and PsiCC transitions don't show up at all in the summary results; PsiAC and PsiCB show up as 0; and PsiBC shows up as non-zero.

The code that generates my results is pasted below, along with the relevant portion of the output of summary(). The code is copied from Appendix C of the book, and as far as I can tell should produce the same output that the book shows on p. C-81 to C-82. I'd be very grateful for any ideas about what I'm doing wrong - I've searched the forum a few different ways and haven't been able to turn up an answer.

Thanks
Chris

<<Code>>

Code: Select all
library(RMark)
data(mstrata)
mstrata.processed=process.data(mstrata,model="Multistrata")
mstrata.ddl=make.design.data(mstrata.processed,parameters=list(Psi=list(subtract.stratum=c("B","A","A"))))
mstrata.ddl$Psi=mstrata.ddl$Psi[!(mstrata.ddl$Psi$stratum=="A"&mstrata.ddl$Psi$tostratum=="A"),]
mstrata.ddl$Psi=mstrata.ddl$Psi[!(mstrata.ddl$Psi$stratum=="B"&mstrata.ddl$Psi$tostratum=="B"),]
mstrata.ddl$Psi=mstrata.ddl$Psi[!(mstrata.ddl$Psi$stratum=="C"&mstrata.ddl$Psi$tostratum=="C"),]
mymodel=mark(mstrata.processed,mstrata.ddl)
summary(mymodel,show.fixed=T)


<<Abbreviated output>>

Code: Select all
Real Parameter Psi
 Stratum:A To:C
  1 2 3
1 0 0 0
2   0 0
3     0

 Stratum:B To:C
          1         2         3
1 0.2214776 0.2214776 0.2214776
2           0.2214776 0.2214776
3                     0.2214776

 Stratum:C To:B
  1 2 3
1 0 0 0
2   0 0
3     0
ctsolomon
 
Posts: 3
Joined: Fri Nov 19, 2021 2:36 pm

Re: problem reproducing RMark multi-strata example

Postby jlaake » Fri Nov 19, 2021 8:58 pm

Clearly I need to review what is in Appendix C and remove that example. You have stumbled on an approach that I wish I had never added. I added the option to delete design data to set values to their default values as shown in that example. It was originally intended for an example where Steller sea lions were only marked and released every other year so the design data for years in which there were no releases was removed. Unfortunately this has made my code more complicated than it should be. Somewhere along the line it broke for Psi and other mlogit parameters. Had you looked closely at the output of the run you would have seen

Code: Select all
Warning message:
In make.mark.model(data.proc, title = title, parameters = model.parameters,  :
 
Missing rows in design dataframe for parameter Psi
 Deleting rows from design data is still allowed but see warning in help for make.design.data

The help for make.design.data contains the following:

Code: Select all
######WARNING########
Deleting design data for mlogit parameters like Psi in the multistate
model can fail if you do things like delete certain transitions.  It is better
to add the field fix. It should be assigned the value NA for parameters that
are estimated and a fixed real value for those that are fixed. Here is an example
with the mstrata data example:

data(mstrata)
# deleting design data approach to fix Psi A to B to 0 (DON'T use this approach)
dp=process.data(mstrata,model="Multistrata")
ddl=make.design.data(dp)
ddl$Psi=ddl$Psi[!(ddl$Psi$stratum=="A" & ddl$Psi$tostratum=="B"),]
ddl$Psi
summary(mark(dp,ddl,output=FALSE,delete=TRUE),show.fixed=TRUE)
#new approach using fix to set Phi=1 for time 2 (USE this approach)
ddl=make.design.data(dp)
ddl$Psi$fix=NA
ddl$Psi$fix[ddl$Psi$stratum=="A" & ddl$Psi$tostratum=="B"]=0
ddl$Psi
summary(mark(dp,ddl,output=FALSE,delete=TRUE),show.fixed=TRUE)



I may remove design data deletion because the warning doesn't appear to be sufficient.
jlaake
 
Posts: 1258
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: problem reproducing RMark multi-strata example

Postby ctsolomon » Mon Nov 22, 2021 1:06 pm

Thanks Jeff for that quick and helpful reply. I should have looked more carefully at warnings(). I think I understand now how to make that example (and the actual analysis that I'm working on) run. I appreciate the help.

A suggestion FWIW: If you do delete that example from Appendix C, it might be useful to replace it with at least a sentence or two pointing folks towards the right way to do this kind of thing (as indicated in the help for make.design.data). The ability to fix some of the transition probabilities seems pretty useful, and I know I was glad to have the appendix point me towards how to do it.

Chris
ctsolomon
 
Posts: 3
Joined: Fri Nov 19, 2021 2:36 pm

Re: problem reproducing RMark multi-strata example

Postby jlaake » Mon Nov 22, 2021 3:31 pm

Most definitely. You may want to read through the sticky notes at the top of the RMark forum. The one on fixing real parameters is relevant here. Also, the workshop notes (see documentation link when you type library(RMark)) were written more recently and discuss fixing real parameters.

--jeff
jlaake
 
Posts: 1258
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: problem reproducing RMark multi-strata example

Postby ctsolomon » Tue Nov 23, 2021 8:43 am

Great, thanks!!
ctsolomon
 
Posts: 3
Joined: Fri Nov 19, 2021 2:36 pm


Return to RMark

Who is online

Users browsing this forum: No registered users and 1 guest