recreating multistate MARK model in RMark
Posted: Mon Mar 27, 2023 6:50 am
I created fake data (provided below at the very bottom) for a 2-state, multistate, live-recapture CJS model with four capture occasions.
I analyzed those data in MARK and in RMark. I obtained the expected real estimates with both applications with the caveat that the last detection was constrained = 1 in both states.
However, the transition betas from RMark seem off and have huge SE's perhaps because RMark uses three levels of time in the transition design matrix but I only used two levels of time for the transition design matrix in MARK.
I cannot figure out how to remove Psi:time1 from the design matrix in RMark. Rows cannot be deleted from the RMark design matrix. How can I more closely approximately the MARK beta estimates while using RMark? Or should I be trying to more closely approximately the RMark betas from within MARK? Or should I just not worry about it given that the real estimates seem fine? I would feel more comfortable if I could reproduce the betas with both applications.
Here are the true values of all parameters used to create the data:
Here is the design matrix from within MARK: survival, detection and then transition. I tried to create this design matrix so the parameter estimates would appear in the same order as in RMark and match the values returned by RMark:
Here are the betas and real estimates from MARK:
Here is the RMark code:
Here are the beta and real estimates from RMark. Note the large SE's and confidence intervals of the beta transition estimates:
Here are my fake data:
Thank you for any help.
I analyzed those data in MARK and in RMark. I obtained the expected real estimates with both applications with the caveat that the last detection was constrained = 1 in both states.
However, the transition betas from RMark seem off and have huge SE's perhaps because RMark uses three levels of time in the transition design matrix but I only used two levels of time for the transition design matrix in MARK.
I cannot figure out how to remove Psi:time1 from the design matrix in RMark. Rows cannot be deleted from the RMark design matrix. How can I more closely approximately the MARK beta estimates while using RMark? Or should I be trying to more closely approximately the RMark betas from within MARK? Or should I just not worry about it given that the real estimates seem fine? I would feel more comfortable if I could reproduce the betas with both applications.
Here are the true values of all parameters used to create the data:
- Code: Select all
# survival
phiA1 <- 0.80
phiA2 <- 0.70
phiA3 <- 0.60
phiB1 <- 0.40
phiB2 <- 0.30
phiB3 <- 0.20
# detection
pA1 <- 0.75
pA2 <- 0.85
pA3 <- 0.95
pB1 <- 0.65
pB2 <- 0.55
pB3 <- 0.45
# movement
psiAB1 <- 0.22
psiAB2 <- 0.32
psiAB3 <- 0.42
psiBA1 <- 0.37
psiBA2 <- 0.47
psiBA3 <- 0.57
Here is the design matrix from within MARK: survival, detection and then transition. I tried to create this design matrix so the parameter estimates would appear in the same order as in RMark and match the values returned by RMark:
- Code: Select all
INPUT --- design matrix constraints=18 covariates=18;
INPUT --- 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
INPUT --- 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
INPUT --- 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
INPUT --- 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
INPUT --- 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
INPUT --- 1 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
INPUT --- 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
INPUT --- 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0;
INPUT --- 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0;
INPUT --- 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0;
INPUT --- 0 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0;
INPUT --- 0 0 0 0 0 0 1 0 1 1 0 1 0 0 0 0 0 0;
INPUT --- 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0;
INPUT --- 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0;
INPUT --- 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0;
INPUT --- 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0;
INPUT --- 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0;
INPUT --- 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1;
Here are the betas and real estimates from MARK:
- Code: Select all
PARM-SPECIFIC Link Function Parameters of {phi(g,t)p(g,t)psi(g,t) DM5 interactions6 mlogit B}
95% Confidence Interval
Parameter Beta Standard Error Lower Upper
------------------------- -------------- -------------- -------------- --------------
1: 1.3862943 0.1444037 1.1032630 1.6693257
2: -0.5389987 0.2834436 -1.0945481 0.0165507
3: -1.6112381 0.1516505 -1.9084731 -1.3140031
4: -1.7917598 0.1402091 -2.0665696 -1.5169501
5: 0.0971644 0.1887085 -0.2727043 0.4670331
6: 0.2583769 0.1665575 -0.0680758 0.5848295
7: 1.0986121 0.1543161 0.7961526 1.4010716
8: 0.6359843 0.4614205 -0.2683999 1.5403685
9: 0.0000000 0.0000000 0.0000000 0.0000000
10: -0.4795723 0.5244719 -1.5075373 0.5483928
11: -1.0543457 1.2519194 -3.5081077 1.3994164
12: 0.0000000 0.0000000 0.0000000 0.0000000
13: -0.5322165 0.2148391 -0.9533011 -0.1111320
14: -1.2656666 0.1972307 -1.6522388 -0.8790943
15: 0.5118905 0.4458855 -0.3620452 1.3858261
16: 0.1956788 0.2118506 -0.2195483 0.6109060
17: -0.0998140 0.8758751 -1.8165292 1.6169012
18: 1.3656032 0.4160351 0.5501744 2.1810321
Real Function Parameters of {phi(g,t)p(g,t)psi(g,t) DM5 interactions6 mlogit B}
95% Confidence Interval
Parameter Estimate Standard Error Lower Upper
-------------------------- -------------- -------------- -------------- --------------
1:S 1:State 1 0.8000000 0.0231046 0.7508710 0.8414859
2:S 1:State 1 0.6999995 0.0498495 0.5943635 0.7879408
3:S 1:State 1 0.4440000 0.0114345 0.4217170 0.4665106
4:S 2:State 2 0.3999999 0.0314172 0.3402817 0.4628456
5:S 2:State 2 0.2999996 0.0413431 0.2256272 0.3866465
6:S 2:State 2 0.1470000 0.0096619 0.1290535 0.1669638
7:p 1:State 1 0.7500000 0.0289343 0.6891509 0.8023539
8:p 1:State 1 0.8499994 0.0555446 0.7069748 0.9301150
9:p 1:State 1 1.0000000 0.0000000 1.0000000 1.0000000 Fixed
10:p 2:State 2 0.6500001 0.0935853 0.4533279 0.8061700
11:p 2:State 2 0.5500019 0.1793169 0.2280478 0.8348950
12:p 2:State 2 1.0000000 0.0000000 1.0000000 1.0000000 Fixed
13:Psi 1 to 2 0.2200000 0.0338448 0.1608066 0.2933655
14:Psi 1 to 2 0.3199991 0.0865173 0.1775451 0.5063788
15:Psi 1 to 2 0.2554054 0.0147048 0.2276653 0.2852770
16:Psi 2 to 1 0.3700001 0.0500790 0.2782214 0.4722456
17:Psi 2 to 1 0.4700011 0.0995850 0.2882893 0.6600297
18:Psi 2 to 1 0.7367347 0.0301691 0.6735304 0.7914907
Here is the RMark code:
- Code: Select all
### Bring in the Data
library(RMark)
ms <- convert.inp("fake_deterministic_multistate.inp")
head(ms)
# Process data
mstrata.processed = process.data(ms, begin.time = 1, model = "Multistrata")
# Create default design data
mstrata.ddl = make.design.data(mstrata.processed)
# examine the data
head(mstrata.ddl$S)
head(mstrata.ddl$p)
head(mstrata.ddl$Psi)
# mwm: is removing this row wise?
# Deleting rows from design data is no longer allowed
# mstrata.ddl$Psi <- mstrata.ddl$Psi[!(mstrata.ddl$Psi$time == 3),]
# mwm: is changing the value of occ 3 to occ 0 wise? Is it allowed?
# mwm: This does not help.
# mstrata.ddl$Psi$occ[mstrata.ddl$Psi$occ == 3] <- 0
# mwm: can I eliminate time1 this way?
mstrata.ddl$Psi$mytime = 0
mstrata.ddl$Psi$mytime[((mstrata.ddl$Psi$time == 2) | (mstrata.ddl$Psi$time == 3)) & mstrata.ddl$Psi$tostratum == 1] = 1
mstrata.ddl$Psi
### Build Function for Creating Models
# Set up a function that contains the structures for 'S', 'p', and 'Psi'.
run.ms = function() {
# survival probability
S.time = list(formula = ~ time)
S.site = list(formula = ~ stratum)
#S.site.time = list(formula = ~ time : stratum)
S.site.time = list(formula = ~ time * stratum)
# detection probability
p.time = list(formula = ~ time)
p.site = list(formula = ~ stratum)
#p.site.time = list(formula = ~ time : stratum)
p.site.time = list(formula = ~ time * stratum)
# transition probs
Psi.site = list(formula = ~ 1 + stratum:tostratum)
#Psi.site.time = list(formula = ~ -1 + time + stratum:tostratum)
#Psi.site.time = list(formula = ~ -1 + time * stratum:tostratum)
#Psi.site.time = list(formula = ~ -1 + stratum:tostratum * time)
#Psi.site.time = list(formula = ~ -1 + stratum:tostratum * Time)
#Psi.site.time = list(formula = ~ -1 + stratum:tostratum * occ)
#Psi.site.time = list(formula = ~ -1 + stratum:tostratum * (time:mytime))
Psi.site.time = list(formula = ~ -1 + stratum:tostratum * time)
# Fit models
S.dot.p.site.psi.dot = mark(mstrata.processed, mstrata.ddl,
model.parameters = list(S = S.site, p = p.site, Psi = Psi.site),
output=FALSE, silent=TRUE, delete=TRUE)
S.dot.p.time.psi.dot = mark(mstrata.processed, mstrata.ddl,
model.parameters = list(S = S.site, p = p.site.time, Psi = Psi.site),
output=FALSE, silent=TRUE, delete=TRUE)
S.sitetime.p.sitetime.psi.sitetime = mark(mstrata.processed, mstrata.ddl,
model.parameters = list(S = S.site.time, p = p.site.time, Psi = Psi.site.time),
output=FALSE, silent=TRUE, delete=TRUE)
head(mstrata.ddl$p)
mstrata.ddl$p$fix=NA
mstrata.ddl$p$fix[mstrata.ddl$p$time==4] = 1
S.sitetime.p.sitetime.psi.sitetime.lastpis1 = mark(mstrata.processed, mstrata.ddl,
model.parameters = list(S = S.site.time, p = p.site.time, Psi = Psi.site.time),
output=FALSE, silent=TRUE, delete=TRUE)
# Return model table and list of models
return(collect.models())
}
### Run the models and examine the output
ms.results = run.ms()
ms.results
names(ms.results)
# examine the output from the best model
ms.results$S.dot.p.site.psi.dot$results$beta
# examine the output from the most-general model
ms.results$S.sitetime.p.sitetime.psi.sitetime.lastpis1$results$beta
ms.results$S.sitetime.p.sitetime.psi.sitetime.lastpis1$results$real
Here are the beta and real estimates from RMark. Note the large SE's and confidence intervals of the beta transition estimates:
- Code: Select all
# examine the output from the most-general model
ms.results$S.sitetime.p.sitetime.psi.sitetime.lastpis1$results$beta
estimate se lcl ucl
S:(Intercept) 1.3862957 0.1444036 1.1032647 1.6693266
S:time2 -0.5389980 0.2834460 -1.0945522 0.0165561
S:time3 -1.6112395 0.1516503 -1.9084741 -1.3140048
S:stratum2 -1.7917599 0.1402088 -2.0665692 -1.5169506
S:time2:stratum2 0.0971645 0.1887084 -0.2727039 0.4670329
S:time3:stratum2 0.2583765 0.1665573 -0.0680758 0.5848287
p:(Intercept) 1.0986124 0.1543152 0.7961546 1.4010701
p:time3 0.6359883 0.4614243 -0.2684033 1.5403798
p:stratum2 -0.4795764 0.5244659 -1.5075296 0.5483769
p:time3:stratum2 -1.0543539 1.2519246 -3.5081261 1.3994184
Psi:time1 -0.6222362 46.3370340 -91.4428240 90.1983520
Psi:time2 -0.1055139 54.3813650 -106.6929900 106.4819600
Psi:time3 0.1726090 97.8467210 -191.6069700 191.9521900
Psi:stratum2:tostratum1 0.0900185 46.3370960 -90.7306920 90.9107290
Psi:stratum1:tostratum2 -0.6434277 46.3376500 -91.4652240 90.1783690
Psi:stratum2:tostratum1:time2 -0.1046489 53.5548750 -105.0722100 104.8629100
Psi:stratum1:tostratum2:time2 -0.0048308 53.5567860 -104.9761300 104.9664700
Psi:stratum2:tostratum1:time3 0.7664387 96.6687580 -188.7043300 190.2372100
Psi:stratum1:tostratum2:time3 -0.5991695 96.6682550 -190.0689500 188.8706100
ms.results$S.sitetime.p.sitetime.psi.sitetime.lastpis1$results$real
estimate se lcl ucl fixed note
S s1 g1 c1 a0 o1 t1 0.8000002 0.0231046 0.7508713 0.8414860
S s1 g1 c1 a1 o2 t2 0.6999999 0.0498497 0.5943634 0.7879415
S s1 g1 c1 a2 o3 t3 0.4440000 0.0114345 0.4217170 0.4665106
S s2 g1 c1 a0 o1 t1 0.4000002 0.0314170 0.3402823 0.4628456
S s2 g1 c1 a1 o2 t2 0.3000000 0.0413433 0.2256274 0.3866473
S s2 g1 c1 a2 o3 t3 0.1470000 0.0096619 0.1290534 0.1669638
p s1 g1 c1 a1 o1 t2 0.7500000 0.0289341 0.6891513 0.8023536
p s1 g1 c1 a2 o2 t3 0.8499999 0.0555448 0.7069746 0.9301156
p s2 g1 c1 a1 o1 t2 0.6499993 0.0935843 0.4533293 0.8061679
p s2 g1 c1 a2 o2 t3 0.5499999 0.1793163 0.2280474 0.8348931
p s1 g1 c1 a3 o3 t4 1.0000000 0.0000000 1.0000000 1.0000000 Fixed
Psi s1 to2 g1 c1 a0 o1 t1 0.2200004 0.0338446 0.1608073 0.2933656
Psi s1 to2 g1 c1 a1 o2 t2 0.3199999 0.0865176 0.1775455 0.5063801
Psi s1 to2 g1 c1 a2 o3 t3 0.2554053 0.0147048 0.2276653 0.2852770
Psi s1 to2 g1 c2 a0 o2 t2 0.3199999 0.0865176 0.1775455 0.5063801
Psi s1 to2 g1 c2 a1 o3 t3 0.2554053 0.0147048 0.2276653 0.2852770
Psi s1 to2 g1 c3 a0 o3 t3 0.2554053 0.0147048 0.2276653 0.2852770
Psi s2 to1 g1 c1 a0 o1 t1 0.3699998 0.0500786 0.2782218 0.4722445
Psi s2 to1 g1 c1 a1 o2 t2 0.4700000 0.0995851 0.2882883 0.6600290
Psi s2 to1 g1 c1 a2 o3 t3 0.7367348 0.0301691 0.6735306 0.7914908
Psi s2 to1 g1 c2 a0 o2 t2 0.4700000 0.0995851 0.2882883 0.6600290
Psi s2 to1 g1 c2 a1 o3 t3 0.7367348 0.0301691 0.6735306 0.7914908
Psi s2 to1 g1 c3 a0 o3 t3 0.7367348 0.0301691 0.6735306 0.7914908
Here are my fake data:
- Code: Select all
1000 314.01016428 ;
1001 6.293399868 ;
1002 2.189955852 ;
1010 39.19820016 ;
1011 23.307419016 ;
1012 7.994740824 ;
1020 20.98902036 ;
1021 2.664842796 ;
1022 0.952256844 ;
1100 199.2186144 ;
1101 16.15605264 ;
1102 5.61493296 ;
1110 105.2801568 ;
1111 62.60003568 ;
1112 21.47260752 ;
1120 49.1819328 ;
1121 6.24431808 ;
1122 2.23134912 ;
1200 88.40735332 ;
1201 1.686376692 ;
1202 0.591149988 ;
1210 7.62322704 ;
1211 4.532803704 ;
1212 1.554809256 ;
1220 8.53365084 ;
1221 1.083463524 ;
1222 0.387165636 ;
2000 683.91040431 ;
2001 2.577455811 ;
2002 0.899679879 ;
2010 14.20077932 ;
2011 8.443844682 ;
2012 2.896345998 ;
2020 10.46758097 ;
2021 1.329002367 ;
2022 0.474906663 ;
2100 47.2505688 ;
2101 3.83188428 ;
2102 1.33174692 ;
2110 24.9702936 ;
2111 14.84744436 ;
2112 5.09286204 ;
2120 11.6649456 ;
2121 1.48102416 ;
2122 0.52923024 ;
2200 126.58325589 ;
2201 2.414584809 ;
2202 0.846419301 ;
2210 10.91507508 ;
2211 6.490150758 ;
2212 2.226204162 ;
2220 12.21863643 ;
2221 1.551322773 ;
2222 0.554350797 ;
0100 425.6808 ;
0101 34.52148 ;
0102 11.99772 ;
0110 224.9576 ;
0111 133.76076 ;
0112 45.88164 ;
0120 105.0896 ;
0121 13.34256 ;
0122 4.76784 ;
0200 772.79155 ;
0201 14.741055 ;
0202 5.167395 ;
0210 66.6366 ;
0211 39.62241 ;
0212 13.59099 ;
0220 74.59485 ;
0221 9.470835 ;
0222 3.384315 ;
0010 556 ;
0011 330.6 ;
0012 113.4 ;
0020 853 ;
0021 108.3 ;
0022 38.7 ;
Thank you for any help.