www.phidot.org

Analysis Forum
http://phidot.org/forum/

question about multi-state model in RMark versus MARK

http://phidot.org/forum/viewtopic.php?f=21&t=2910

Page 1 of 1

question about multi-state model in RMark versus MARK

PostPosted: Fri Jan 23, 2015 12:41 pm
by markmiller
I am having trouble reproducing estimates from MARK when using a multi-state model in RMark and I am hoping someone might spot my mistake and correct it. I do obtain the same correct estimates with MARK and RMark with an initial model, but not with a second, more-complex model.

I feel bad asking such a detailed question, but I am at a loss and appreciate any help you can offer.

I am using the following set of fake data:

Code: Select all
111333     9  ;
111334     3  ;
111344     4  ;
111222    16  ;
112222    32  ;
122222    64  ;


There are four states and six visits. For simplicity I constrain survival and detection probability to one.

All animals begin in State 1. Animals in State 1 can move to States 1, 2 or 3. Animals in State 2 must remain there. Animals in State 3 can move to States 3 or 4. Animals in State 4 must remain there.

The simplest model under this scenario has three constant parameters. Here are the estimates from Program MARK. These estimates are correct, as far as I can tell:

Code: Select all
                Real Function Parameters of {initial DM}
                                                             95% Confidence Interval
Parameter         Estimate       Standard Error      Lower           Upper
----------------  --------------  --------------  --------------  --------------
1:S 1:State 1    1.0000000       0.0000000    1.0000000       1.0000000   Fixed               
2:Psi 1 to 2      0.5000000       0.0334077     0.4348928       0.5651072                           
3:Psi 1 to 3      0.0714286       0.0172076     0.0442159       0.1134020                           
4:Psi 1 to 4      0.0000000       0.0000000     0.0000000       0.0000000  Fixed               
5:Psi 3 to 4      0.2500000       0.0818317     0.1241150       0.4394984   


Model.1 of the RMark code below returns the same estimates as those above.

However, if I try to constrain Psi 1 to 2, Psi 1 to 3, and Psi 3 to 4 to zero between some visits and Psi 1 to 3 is allowed to take on two different estimates (using the following design matrix) I get the correct answer with MARK, but not with RMark.

Here is the design matrix:

Code: Select all
Survival: 0 0 0 0;
Psi 1-2: 1 0 0 0;
Psi 1-2: 1 0 0 0;
Psi 1-2: 1 0 0 0;
Psi 1-2: 0 0 0 0;
Psi 1-2: 0 0 0 0;
Psi 1-3: 0 1 0 0;
Psi 1-3: 0 1 0 0;
Psi 1-3: 0 0 1 0;
Psi 1-3: 0 0 0 0;
Psi 1-3: 0 0 0 0;
Psi 3-4: 0 0 0 0;
Psi 3-4: 0 0 0 0;
Psi 3-4: 0 0 0 0;
Psi 3-4: 0 0 0 1;
Psi 3-4: 0 0 0 1;
Psi 1-4: 0 0 0 0;


Here are the correct estimates from Program MARK using the design matrix immediately above:

Code: Select all
Parameter   Estimate       Standard Error      Lower           Upper
 --------------------  --------------  --------------  --------------  --------------
1:S 1:State 1  1.0000000       0.0000000       1.0000000       1.0000000  Fixed               
2:Psi 1 to 2     0.5000000       0.6705308E-005  0.4999868       0.5000131                           
3:Psi 1 to 2     0.5000000       0.6705308E-005  0.4999868       0.5000131                           
4:Psi 1 to 2     0.5000000       0.0883884       0.3333739       0.6666261                           
5:Psi 1 to 2     0.0000000       0.0000000       0.0000000       0.0000000  Fixed               
6:Psi 1 to 2     0.0000000       0.0000000       0.0000000       0.0000000  Fixed               
7:Psi 1 to 3     0.1067929E-009  0.3003999E-006  -.5886770E-006  0.5888906E-006                     
8:Psi 1 to 3     0.1067929E-009  0.3003999E-006  -.5886770E-006  0.5888906E-006                     
9:Psi 1 to 3     0.5000000       0.0883884       0.3333739       0.6666261                           
10:Psi 1 to 3   0.0000000       0.0000000       0.0000000       0.0000000  Fixed               
11:Psi 1 to 3   0.0000000       0.0000000       0.0000000       0.0000000  Fixed               
12:Psi 3 to 4   0.0000000       0.0000000       0.0000000       0.0000000  Fixed               
13:Psi 3 to 4   0.0000000       0.0000000       0.0000000       0.0000000  Fixed               
14:Psi 3 to 4   0.0000000       0.0000000       0.0000000       0.0000000  Fixed               
15:Psi 3 to 4   0.2499999       0.0818317       0.1241149       0.4394982                           
16:Psi 3 to 4   0.2499999       0.0818317       0.1241149       0.4394982                           
17:Psi 1 to 4   0.0000000       0.0000000       0.0000000       0.0000000  Fixed


However, model.2 in RMark, does not constrain Psi 1 to 2 to a constant value for the first three time periods. Rather, RMark allows Psi 1 to 2 to differ between the second and third time period even though I think I specifically constrain Psi 1 to 2 to have the same value for each of the first three time periods.

Why does RMark allow the estimate for Psi 1 to 2 to differ among the first three visits versus between the third and fourth visit even though I, apparently, have constrained Psi 1 to 2 to be constant among the first four visits?

Here is the R code.

library(RMark)

setwd('C:/Users/markm/simple R programs/')
my.data = convert.inp("fakeJan23.inp")
n.visits <- dim(my.data)[1]

mstrat.processed = process.data(my.data, model="Multistrata", time.intervals = rep(1, n.visits-1) )
mstrat.ddl = make.design.data(mstrat.processed)

#
# Additions to the design matrix
#
mstrat.ddl$Psi$mystrat12 = 0
mstrat.ddl$Psi$mystrat12[ mstrat.ddl$Psi$stratum == 1 & mstrat.ddl$Psi$to2 == 1 & mstrat.ddl$Psi$time %in% c(1,2,3)] = 1

mstrat.ddl$Psi$mystrat13a = 0
mstrat.ddl$Psi$mystrat13a[mstrat.ddl$Psi$stratum == 1 & mstrat.ddl$Psi$to3 == 1 & mstrat.ddl$Psi$time %in% c(1:2)] = 1

mstrat.ddl$Psi$mystrat13b = 0
mstrat.ddl$Psi$mystrat13b[mstrat.ddl$Psi$stratum == 1 & mstrat.ddl$Psi$to3 == 1 & mstrat.ddl$Psi$time %in% c(3) ] = 1

mstrat.ddl$Psi$mystrat34 = 0
mstrat.ddl$Psi$mystrat34[ mstrat.ddl$Psi$stratum == 3 & mstrat.ddl$Psi$to4 == 1 & mstrat.ddl$Psi$time %in% c(4:5)] = 1

#############################################################
#
# Constrain survival and p = 1 and constrain some transitions = 0
#

pre.S = as.numeric(row.names(mstrat.ddl$S[ mstrat.ddl$S$time %in% c(1:(n.visits-1)), ]))
pre.p = as.numeric(row.names(mstrat.ddl$p[ mstrat.ddl$p$time %in% c(2: n.visits ), ]))

pre.P12 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 1 & mstrat.ddl$Psi$to2 == 1 & mstrat.ddl$Psi$time %in% c(4:5), ]))
pre.P13 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 1 & mstrat.ddl$Psi$to3 == 1 & mstrat.ddl$Psi$time %in% c(4:5), ]))
pre.P14 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 1 & mstrat.ddl$Psi$to4 == 1, ]))

pre.P21 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 2 & mstrat.ddl$Psi$to1 == 1, ]))
pre.P23 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 2 & mstrat.ddl$Psi$to3 == 1, ]))
pre.P24 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 2 & mstrat.ddl$Psi$to4 == 1, ]))

pre.P31 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 3 & mstrat.ddl$Psi$to1 == 1, ]))
pre.P32 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 3 & mstrat.ddl$Psi$to2 == 1, ]))
pre.P34 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 3 & mstrat.ddl$Psi$to4 == 1 & mstrat.ddl$Psi$time %in% c(1:3), ]))

pre.P41 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 4 & mstrat.ddl$Psi$to1 == 1, ]))
pre.P42 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 4 & mstrat.ddl$Psi$to2 == 1, ]))
pre.P43 = as.numeric(row.names(mstrat.ddl$Psi[ mstrat.ddl$Psi$stratum == 4 & mstrat.ddl$Psi$to3 == 1, ]))

#####################################################################
#
# This works
#

S.stratum = list(formula =~ 1, link="sin",
fixed = list(index = c(pre.S), value = 1))

p.stratum = list(formula =~ 1, link="sin",
fixed = list(index = c(pre.p), value = 1))

Psi.s = list(formula =~ -1 + stratum:tostratum,
fixed = list(index = c( pre.P14,
pre.P21, pre.P23, pre.P24,
pre.P31, pre.P32,
pre.P41, pre.P42, pre.P43), value = 0))

model.1 <- mark(mstrat.processed, mstrat.ddl,
model.parameters=list(S = S.stratum, p = p.stratum, Psi = Psi.s))

########################################################################
#
# This does not works
#

S.stratum = list(formula =~ 1, link="sin",
fixed = list(index = c(pre.S), value = 1))

p.stratum = list(formula =~ 1, link="sin",
fixed = list(index = c(pre.p), value = 1))

Psi.s = list(formula =~ -1 + mystrat12 + mystrat13a + mystrat13b + mystrat34,
fixed = list(index = c(pre.P12, pre.P13, pre.P14,
pre.P21, pre.P23, pre.P24,
pre.P31, pre.P32, pre.P34,
pre.P41, pre.P42, pre.P43), value = 0))

model.2 <- mark(mstrat.processed, mstrat.ddl,
model.parameters=list(S = S.stratum, p = p.stratum, Psi = Psi.s))

########################################################################

Re: question about multi-state model in RMark versus MARK

PostPosted: Fri Jan 23, 2015 12:51 pm
by jlaake
Mark-

You are using the old approach to fixing parameters and it is messy and too difficult to debug. That is one of the reasons I came up with the new approach. See http://www.phidot.org/forum/viewtopic.php?f=21&t=2887 which is a sticky at the top of the RMark forum. My guess is that you have a bug somewhere in your code. If you use the newer approach, it will be easier to specify and debug. If you find it still doesn't work that way let me know. Either way post back here to close the loop.

Re: question about multi-state model in RMark versus MARK

PostPosted: Fri Jan 23, 2015 1:04 pm
by jlaake
Ps

Did you look at the design matrix that RMark created to see if it is the same as what you produced in MARK?

Re: question about multi-state model in RMark versus MARK

PostPosted: Fri Jan 23, 2015 1:25 pm
by markmiller
Thank you, Jeff, for the replies. I will look at the newer approach for specifying constraints. Here is the design matrix that RMark returns for model.2. I think it is the same as what I created in MARK, except that there are more lines to account for a potentially more general structure to transition probabilities.

Code: Select all
INPUT ---    design matrix constraints=181 covariates=4;
  INPUT ---        0 0 0 0;
  INPUT ---        1 0 0 0;
  INPUT ---        1 0 0 0;
  INPUT ---        1 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        1 0 0 0;
  INPUT ---        1 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        1 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 1 0 0;
  INPUT ---        0 1 0 0;
  INPUT ---        0 0 1 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 1 0 0;
  INPUT ---        0 0 1 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 1 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 1;
  INPUT ---        0 0 0 1;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 1;
  INPUT ---        0 0 0 1;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 1;
  INPUT ---        0 0 0 1;
  INPUT ---        0 0 0 1;
  INPUT ---        0 0 0 1;
  INPUT ---        0 0 0 1;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;
  INPUT ---        0 0 0 0;

Re: question about multi-state model in RMark versus MARK

PostPosted: Fri Jan 23, 2015 1:36 pm
by jlaake
If you are using the Mlogit link in MARK you have to be careful about how you construct the PIMS. See pages 10-24 - 10-27 of the MARK book for an explanation. That restriction is the reason I can't simplify the parameter structure for Mlogit parameters in RMark. If you are not using the MLogit link in MARK that could explain the difference in the results. See the other sticky message about Reasons for differences bt MARK and RMark. Bottom line is that if the PIM,DM,links and data are all the same, then the results will be the same as well.

Re: question about multi-state model in RMark versus MARK

PostPosted: Fri Jan 23, 2015 1:43 pm
by markmiller
Thank you, Jeff. I will look at those sections of the MARK book as well. Yes, I am using the mlogit link within MARK. I might not post again until tomorrow though. I appreciate your feedback.
All times are UTC - 5 hours [ DST ]
Page 1 of 1
Powered by phpBB® Forum Software © phpBB Group
http://www.phpbb.com/