Manually set PIMs in RMark

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

Manually set PIMs in RMark

Postby drobichaud » Mon Jul 12, 2021 5:09 pm

Hi forum,

How do you manually set the PIMs in a CJS analysis in rmark?? I need one of the Phis to be shared between two occasions.

In MARK, I was able to run a four-occasion CJS, with time intervals set to {1,8,8}, a sin link, and with the second and third Phis being shared. Phi and p PIMs looked like this
1__2__2____3__4__5
___2__2_______4__5
______2__________5

I am hoping the fact that I cannot find documentation for this is evidence simply of it being an odd request (or that I am bad at finding answers within this forum), rather than it being impossible. I am trying to replicate an analysis that I did in MARK, but this time using rmark so that I can bootstrap the input data and see how decreases in sampling could affect the resulting confidence bounds. But before I can code the bootstrapping, I first need to make sure the R code for the base analysis produces that same results as when I ran it manually in MARK.

Here is what I have so far:
Code: Select all
formark = data.frame(DH.2017W$DH)
colnames(formark)=c("ch")
tapply(formark$ch,formark$ch, length)


1000_1001_1100_1101
5083___43____4____1


Code: Select all
time.intervals = c(1, 8, 8)
m.proc = process.data(data = formark, model = "CJS", time.intervals = time.intervals)
Phi.time=list(formula=~time)
p.time=list(formula=~time)
m.ddl = make.design.data(data = m.proc, parameters = list(Phi = list(pim.type = "time"),p =list(pim.type = "time")))
m.ddl$Phi$model.index = c(1,2,2)
m.ddl$p$model.index = c(3,4,5)
m.ddl$Phi


__par.index_model.index_group_time_Time
1______ 1____________1_____1____1____0
2______ 2____________2_____1____1____1
3______ 3____________2_____1___10____9

Code: Select all
m.model = mark(data = m.proc,
                 ddl = m.ddl,
                model.name = "Phi.time.occ2eq3.p.time",
                input.links = rep("sin", 6),
                model.parameters=list(Phi=Phi.time, p = p.time) )


But when I run this, it ignores my attempt to set the PIMS, and I get this

Code: Select all
> m.model$pims

$Phi
$Phi[[1]]
$Phi[[1]]$pim
___[,1]_[,2]_[,3]
[1,]__1___2___3
[2,]__0___2___3
[3,]__0___0___3


$p
$p[[1]]
$p[[1]]$pim
___[,1]_[,2]_[,3]
[1,]__4___5___6
[2,]__0___5___6
[3,]__0___0___6


What am I doing wrong?
Any help would be much apprecaited.
Thanks in advance,
Dave
drobichaud
 
Posts: 2
Joined: Mon Jul 12, 2021 4:19 pm

Re: Manually set PIMs in RMark

Postby jhines » Tue Jul 13, 2021 11:09 am

You don't need to manually change the PIMs. You can let RMark
create them using the design data.

The trick is to create a new design-data variable which is
like the built-in variable, 'time', then change the new
variable to only contain the values 1 and 2, instead of
1, 2 and 3.

Using the new variable in the formula for Phi will cause
the PIM to be numbered just like Mark does.

To get RMark to use the Sin link, like Mark, add the
input.links argument to the Mark function call.

Here is some sample R-code:

Code: Select all
#Sample R/RMark code to analyze expected value data from GENCAPH1

a='0001 2.1773
0010 4.5965
0011 1.4515
0100 11.3165
0101 1.4515
0110 3.0643
0111 0.9677
1000 23.3165
1001 1.4515
1010 3.0643
1011 0.9677
1100 7.5443
1101 0.9677
1110 2.0429
1111 0.6451'

#  read histories and frequencies and split into a matrix
#  split into separate fields
a<-unlist(strsplit(unlist(strsplit(a,'\n')),' '))   
a<-matrix(a,ncol=2,byrow=T)  #  format as a matrix with 2 columns

library(RMark)        #  Load RMark package

#  create data-frame with 2 fields: ch (capture-histories) and freq (frequencies)
data<-data.frame(ch=a[,1],freq=as.numeric(a[,2]))

# create RMark processed data object...
pd<-process.data(data,model='CJS')

#  create design data
dd <- make.design.data(pd)
         # make new design data variable, 'time2' which is the same
         #  as the 'time' variable, for now
dd$Phi$time2 <- as.numeric(dd$Phi$time)
         # change new variable, 'time2' to '2' whenever it's current
         # value is '3'
dd$Phi$time2[dd$Phi$time2==3]=2
         # change time2 to a factor variable
dd$Phi$time2=as.factor(dd$Phi$time2)

#  run mark model, phi(.)p(.) - constant survival and capture probs
m1<-mark(data=pd,ddl=dd,model.parameters=list(
  Phi=list(formula=~time2),p=list(formula=~1)),
  time.intervals=c(1,8,8),
  input.links=rep('Sin',3))
jhines
 
Posts: 599
Joined: Fri May 16, 2003 9:24 am
Location: Laurel, MD, USA

Re: Manually set PIMs in RMark

Postby drobichaud » Tue Jul 13, 2021 11:38 pm

Hi jhines,
Thanks for the reply. I neer would have thought of that as the way to set up PIMs, so I thank you again. I tried your script, and it did not properly implement the time intervals. I moved that piece into the 'process.data' call, and saw that it worked. This is what is working for me now:

Code: Select all
  time.ints = c(1, 8, 8)
  m.proc = process.data(formark, model="CJS", time.intervals=time.ints)
  m.ddl = make.design.data(m.proc)
  m.ddl$Phi$time2 = as.numeric(m.ddl$Phi$time)
  m.ddl$Phi$time2[m.ddl$Phi$time2==3] = 2
  m.ddl$Phi$time2 = as.factor(m.ddl$Phi$time2)
  n.param = length(as.numeric(m.ddl$Phi$time2)) + length(as.numeric(m.ddl$p$time))
  m.model = mark(data = m.proc,
                 ddl = m.ddl,
                 model.parameters=list(Phi=list(formula=~time2), p = list(formula=~time)),
                 #time.intervals=c(1,8,8),
                 model.name = "Phi.time.occ2eq3.p.time",
                 parm.specific = T,
                 input.links = rep("sin",n.param ) )
drobichaud
 
Posts: 2
Joined: Mon Jul 12, 2021 4:19 pm


Return to RMark

Who is online

Users browsing this forum: No registered users and 8 guests