recreating multistate MARK model in RMark

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

recreating multistate MARK model in RMark

Postby markmiller » 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:

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.
markmiller
 
Posts: 49
Joined: Fri Nov 08, 2013 6:23 pm

Re: recreating multistate MARK model in RMark

Postby jhines » Mon Mar 27, 2023 8:21 am

Hi Mark,

With only 2 stratum, you should only use time*stratum for Psi: since there is only 1 possible transition from each stratum. (eg.,

Psi.site.time=list(formula=~time*stratum)

Also, an easier way to fix the last p is to do it when defining the formula:

p.site.time = list(formula=~time*stratum, fixed=list(time=4, value=1))

Then, you don't need to create another column in the ddl.

Jim
jhines
 
Posts: 599
Joined: Fri May 16, 2003 9:24 am
Location: Laurel, MD, USA

Re: recreating multistate MARK model in RMark

Postby jlaake » Wed Mar 29, 2023 8:04 pm

Jim is correct on both counts. In this case setting p to 1 for the last occasion is easily handled with the original way I handled fixed parameters. However, it doesn't work in more general cases which is why I recommend adding fix to the ddl for the parameter. It also simplifies the formula specification which can be useful if you have a lot of formulas. It is Important to understand that using fix in the ddl will fix the parameter for all of the models and the older format that Jim showed is the only way to do it for just a subset of models.

Now back to the original question about Psi. It looks like you were wildly throwing formula at RMark in hopes that one of them would do what you wanted. Not the best approach. It would be best to use the R function model.matrix with the Psi dataframe (ddl$Psi) and the formula to examine the design matrix it creates. As Jim said, with just 2 strata there is only one estimated transition for each stratum because the other is computed by subtraction. Jim suggested time*stratum which will create 6 parameters and it sounds like you only want 4. I presume you are trying to fit a model that is same for 2 of the years. If that is correct and it is the variable in your ddl$Psi then mytime*stratum should do what you want. If it is correctly defined as a factor variable with only 2 values.

If there are more than 2 states, then you would want to use ~-1 + mytime:stratum:tostratum because the betas vary between stratum. Meaning A moves to B and C but B only moves to A and C etc. That model would create 2×3 parameters with 3 states.

Try using model.matrix to learn more and look at the ?mstrata example in RMark.

Jeff
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: recreating multistate MARK model in RMark

Postby jlaake » Wed Mar 29, 2023 8:33 pm

For my 3 state example the 2x3 was for each time because each of the 3 states would each have 2 transition parameters. Then if mytime had 2 values, there would be 12 parameters in total.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: recreating multistate MARK model in RMark

Postby markmiller » Thu Mar 30, 2023 3:30 am

Thank you, Jim and Jeff. Both of you have been extremely helpful. Jeff is correct that I was pretty much wildly guessing trying to figure out the correct formula for Psi. I have always had a bit of a problem getting the formulae correct in RMark when nested effects are involved. I was not aware of the model.matrix command. It is very powerful and will be very helpful moving forward. Jim sent me some example R code showing how to use it with a variety of models, which is really helpful.

In trying to obtain the design matrix for Psi with model.matrix I was confused initially with what seem to be repeat rows. Then I assumed I could extract just the rows corresponding to Cohort 1 and use those with model.matrix. The formula for Psi in the 4-state model appears to correspond to the identity matrix.

Below is revised RMark code used to estimate the most-general 2-state model with the data set I posted previously as well as the most-general model using a 4-state deterministic data set. The betas and reals look good with RMark and match those with MARK at least to maybe six decimal places. The reals also match the true values used to create the data.

Below are the true values I used to create the 4-state data set, the RMark code and the 4-state data set itself.

Here are the true values of all parameters I used to create the 4-state data set:

Code: Select all
n.released.A1 <- 1000
n.released.B1 <- 1000
n.released.C1 <- 1000
n.released.D1 <- 1000

n.released.A2 <- 1000
n.released.B2 <- 1000
n.released.C2 <- 1000
n.released.D2 <- 1000

n.released.A3 <- 1000
n.released.B3 <- 1000
n.released.C3 <- 1000
n.released.D3 <- 1000

# survival
phiA1 <- 0.80
phiA2 <- 0.70
phiA3 <- 0.60
phiB1 <- 0.40
phiB2 <- 0.30
phiB3 <- 0.20
phiC1 <- 0.23
phiC2 <- 0.33
phiC3 <- 0.43
phiD1 <- 0.53
phiD2 <- 0.63
phiD3 <- 0.73

# detection
pA1 <- 0.75
pA2 <- 0.85
pA3 <- 0.95
pB1 <- 0.65
pB2 <- 0.55
pB3 <- 0.45
pC1 <- 0.44
pC2 <- 0.54
pC3 <- 0.64
pD1 <- 0.74
pD2 <- 0.84
pD3 <- 0.94

# movement
psiAB1 <- 0.22
psiAB2 <- 0.25
psiAB3 <- 0.28
psiAC1 <- 0.17
psiAC2 <- 0.16
psiAC3 <- 0.15
psiAD1 <- 0.10
psiAD2 <- 0.12
psiAD3 <- 0.14
psiAA1 <- 1 - (psiAB1 + psiAC1 + psiAD1)
psiAA2 <- 1 - (psiAB2 + psiAC2 + psiAD2)
psiAA3 <- 1 - (psiAB3 + psiAC3 + psiAD3)

psiBA1 <- 0.11
psiBA2 <- 0.13
psiBA3 <- 0.15
psiBC1 <- 0.17
psiBC2 <- 0.19
psiBC3 <- 0.21
psiBD1 <- 0.23
psiBD2 <- 0.25
psiBD3 <- 0.27
psiBB1 <- 1 - (psiBA1 + psiBC1 + psiBD1)
psiBB2 <- 1 - (psiBA2 + psiBC2 + psiBD2)
psiBB3 <- 1 - (psiBA3 + psiBC3 + psiBD3)

psiCA1 <- 0.26
psiCA2 <- 0.24
psiCA3 <- 0.22
psiCB1 <- 0.20
psiCB2 <- 0.18
psiCB3 <- 0.16
psiCD1 <- 0.14
psiCD2 <- 0.12
psiCD3 <- 0.10
psiCC1 <- 1 - (psiCA1 + psiCB1 + psiCD1)
psiCC2 <- 1 - (psiCA2 + psiCB2 + psiCD2)
psiCC3 <- 1 - (psiCA3 + psiCB3 + psiCD3)

psiDA1 <- 0.10
psiDA2 <- 0.12
psiDA3 <- 0.14
psiDB1 <- 0.16
psiDB2 <- 0.21
psiDB3 <- 0.19
psiDC1 <- 0.17
psiDC2 <- 0.15
psiDC3 <- 0.13
psiDD1 <- 1 - (psiDA1 + psiDB1 + psiDC1)
psiDD2 <- 1 - (psiDA2 + psiDB2 + psiDC2)
psiDD3 <- 1 - (psiDA3 + psiDB3 + psiDC3)


Here is the RMark code with selected output shown as comments:

Code: Select all
library(RMark)

# Read 2-state data
ms <- convert.inp("fake_deterministic_multistate_Mar30_2023.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)

# constrain last detection probability = 1 in both states
mstrata.ddl$p$fix=NA
mstrata.ddl$p$fix[mstrata.ddl$p$time==4] = 1

# examine data
head(mstrata.ddl$S)
head(mstrata.ddl$p)
head(mstrata.ddl$Psi)

# survival probability
S.site.time = list(formula = ~ time * stratum)

# detection probability
p.site.time = list(formula = ~ time * stratum)

# transition probs
Psi.site.time = list(formula = ~ time * stratum)

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)

# examine the design matrix for transition probabilites in the above model
# I am a little confused as to why some rows appear to be repeated.
# The MARK design matrix is Rows 1, 2, 3, 7, 8 and 9.
# The repeated rows might be related to treatment contrasts for time & stratum.
# I need to study this more.
# The pattern of row repeats seems to follow the cohort column in mstrata.ddl$Psi.
desmat = model.matrix(object = ~ time * stratum, data=mstrata.ddl$Psi)
desmat
#    (Intercept) time2 time3 stratum2 time2:stratum2 time3:stratum2
# 1            1     0     0        0              0              0
# 2            1     1     0        0              0              0
# 3            1     0     1        0              0              0
# 4            1     1     0        0              0              0
# 5            1     0     1        0              0              0
# 6            1     0     1        0              0              0
# 7            1     0     0        1              0              0
# 8            1     1     0        1              1              0
# 9            1     0     1        1              0              1
# 10           1     1     0        1              1              0
# 11           1     0     1        1              0              1
# 12           1     0     1        1              0              1
# attr(,"assign")
# [1] 0 1 1 2 3 3
# attr(,"contrasts")
# attr(,"contrasts")$time
# [1] "contr.treatment"
#
# attr(,"contrasts")$stratum
# [1] "contr.treatment"

# mem: Try just using Cohort 1
my.mstrata.ddl <- mstrata.ddl$Psi[mstrata.ddl$Psi$cohort == 1,]
my.mstrata.ddl

my.desmat = model.matrix(object = ~ time * stratum, data=my.mstrata.ddl)
my.desmat
#   (Intercept) time2 time3 stratum2 time2:stratum2 time3:stratum2
# 1           1     0     0        0              0              0
# 2           1     1     0        0              0              0
# 3           1     0     1        0              0              0
# 7           1     0     0        1              0              0
# 8           1     1     0        1              1              0
# 9           1     0     1        1              0              1
# attr(,"assign")
# [1] 0 1 1 2 3 3
# attr(,"contrasts")
# attr(,"contrasts")$time
# [1] "contr.treatment"
#
# attr(,"contrasts")$stratum
# [1] "contr.treatment"


# examine output betas
S.sitetime.p.sitetime.psi.sitetime$results$beta
#                      estimate        se        lcl        ucl
# S:(Intercept)       1.3862933 0.1444050  1.1032594  1.6693271
# S:time2            -0.5389957 0.2834660 -1.0945890  0.0165976
# S:time3            -1.6112369 0.1516518 -1.9084745 -1.3139994
# S:stratum2         -1.7917597 0.1402096 -2.0665705 -1.5169489
# S:time2:stratum2    0.0971641 0.1887100 -0.2727075  0.4670357
# S:time3:stratum2    0.2583766 0.1665581 -0.0680773  0.5848304
# p:(Intercept)       1.0986109 0.1543189  0.7961460  1.4010759
# p:time3             0.6359897 0.4614653 -0.2684823  1.5404617
# p:stratum2         -0.4795662 0.5244893 -1.5075652  0.5484329
# p:time3:stratum2   -1.0543625 1.2520545 -3.5083894  1.3996643
# Psi:(Intercept)    -1.2656690 0.1972362 -1.6522520 -0.8790861
# Psi:time2           0.5118968 0.4459330 -0.3621320  1.3859256
# Psi:time3           0.1956814 0.2118557 -0.2195557  0.6109185
# Psi:stratum2        0.7334545 0.3798621 -0.0110752  1.4779842
# Psi:time2:stratum2 -0.0998261 0.8759708 -1.8167289  1.6170767
# Psi:time3:stratum2  1.3655987 0.4160464  0.5501478  2.1810496

# examine output reals
S.sitetime.p.sitetime.psi.sitetime$results$real
#                            estimate        se       lcl       ucl fixed note
# S s1 g1 c1 a0 o1 t1       0.7999998 0.0231048 0.7508703 0.8414861           
# S s1 g1 c1 a1 o2 t2       0.6999999 0.0498535 0.5943549 0.7879474           
# S s1 g1 c1 a2 o3 t3       0.4440000 0.0114345 0.4217170 0.4665106           
# S s2 g1 c1 a0 o1 t1       0.3999997 0.0314178 0.3402804 0.4628466           
# S s2 g1 c1 a1 o2 t2       0.3000000 0.0413464 0.2256222 0.3866541           
# S s2 g1 c1 a2 o3 t3       0.1470000 0.0096619 0.1290535 0.1669638           
# p s1 g1 c1 a1 o1 t2       0.7499997 0.0289348 0.6891495 0.8023546           
# p s1 g1 c1 a2 o2 t3       0.8499999 0.0555493 0.7069604 0.9301200           
# p s2 g1 c1 a1 o1 t2       0.6500013 0.0935882 0.4533225 0.8061749           
# p s2 g1 c1 a2 o2 t3       0.5500003 0.1793318 0.2280259 0.8349103           
# 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.2199995 0.0338457 0.1608048 0.2933672           
# Psi s1 to2 g1 c1 a1 o2 t2 0.3199999 0.0865251 0.1775356 0.5063971           
# Psi s1 to2 g1 c1 a2 o3 t3 0.2554054 0.0147048 0.2276654 0.2852771           
# Psi s1 to2 g1 c2 a0 o2 t2 0.3199999 0.0865251 0.1775356 0.5063971           
# Psi s1 to2 g1 c2 a1 o3 t3 0.2554054 0.0147048 0.2276654 0.2852771           
# Psi s1 to2 g1 c3 a0 o3 t3 0.2554054 0.0147048 0.2276654 0.2852771           
# Psi s2 to1 g1 c1 a0 o1 t1 0.3700005 0.0500804 0.2782195 0.4722490           
# Psi s2 to1 g1 c1 a1 o2 t2 0.4700001 0.0995933 0.2882751 0.6600436           
# Psi s2 to1 g1 c1 a2 o3 t3 0.7367347 0.0301691 0.6735305 0.7914907           
# Psi s2 to1 g1 c2 a0 o2 t2 0.4700001 0.0995933 0.2882751 0.6600436           
# Psi s2 to1 g1 c2 a1 o3 t3 0.7367347 0.0301691 0.6735305 0.7914907           
# Psi s2 to1 g1 c3 a0 o3 t3 0.7367347 0.0301691 0.6735305 0.7914907           



# Read 4-state data
ms.4state <- convert.inp("fake_deterministic_multistate_4states.inp")
head(ms.4state)

# Process data
mstrata.4state.processed = process.data(ms.4state, begin.time = 1,
     model = "Multistrata")

# Create default design data
mstrata.4state.ddl = make.design.data(mstrata.4state.processed)

# constrain last detection probability = 1 in all four states
mstrata.4state.ddl$p$fix=NA
mstrata.4state.ddl$p$fix[mstrata.4state.ddl$p$time==4] = 1

# examine data
head(mstrata.4state.ddl$S)
head(mstrata.4state.ddl$p)
head(mstrata.4state.ddl$Psi)

# survival probability
S.site.time.4state = list(formula = ~ time * stratum)

# detection probability
p.site.time.4state = list(formula = ~ time * stratum)

# transition probs
Psi.site.time.4state = list(formula = ~ -1 + stratum:tostratum:time)


S.sitetime.p.sitetime.psi.sitetime.4state = mark(mstrata.4state.processed,
     mstrata.4state.ddl,
     model.parameters = list(S = S.site.time.4state,
     p = p.site.time.4state, Psi = Psi.site.time.4state),
     output=FALSE, silent=TRUE, delete=TRUE)

# examine the design matrix for transition probabilites in the above model
# This design matrix is too big to copy and paste here as a comment.
# I only copy and paste the bottom text.
desmat.4state = model.matrix(object = ~ -1 + stratum:tostratum:time,
                data=mstrata.4state.ddl$Psi)
dim(desmat.4state)
#[1] 72 48
desmat.4state
# attr(,"assign")
#  [1] 1 1 1 1 1 1 1 1 1 1
#      1 1 1 1 1 1 1 1 1 1
#      1 1 1 1 1 1 1 1 1 1
#      1 1 1 1 1 1 1 1 1 1
#      1 1 1 1 1 1 1 1
# attr(,"contrasts")
# attr(,"contrasts")$stratum
# [1] "contr.treatment"
#
# attr(,"contrasts")$tostratum
# [1] "contr.treatment"
#
# attr(,"contrasts")$time
# [1] "contr.treatment"


# mem: Try just using Cohort 1
my.mstrata.4state.ddl <- mstrata.4state.ddl$Psi[mstrata.4state.ddl$Psi$cohort == 1,]
my.mstrata.4state.ddl

my.desmat.4state = model.matrix(object = ~ -1 + stratum:tostratum:time,
                         data=my.mstrata.4state.ddl)
my.desmat.4state

colnames(my.desmat.4state)
# [1] "stratum1:tostratum1:time1" "stratum2:tostratum1:time1"
#     "stratum3:tostratum1:time1" "stratum4:tostratum1:time1"
#     "stratum1:tostratum2:time1" "stratum2:tostratum2:time1"
#     "stratum3:tostratum2:time1" "stratum4:tostratum2:time1"
#     "stratum1:tostratum3:time1" "stratum2:tostratum3:time1"
#     "stratum3:tostratum3:time1" "stratum4:tostratum3:time1"
#     "stratum1:tostratum4:time1" "stratum2:tostratum4:time1"
#     "stratum3:tostratum4:time1" "stratum4:tostratum4:time1"
#     "stratum1:tostratum1:time2" "stratum2:tostratum1:time2"
#     "stratum3:tostratum1:time2" "stratum4:tostratum1:time2"
#     "stratum1:tostratum2:time2" "stratum2:tostratum2:time2"
#     "stratum3:tostratum2:time2" "stratum4:tostratum2:time2"
#     "stratum1:tostratum3:time2" "stratum2:tostratum3:time2"
#     "stratum3:tostratum3:time2" "stratum4:tostratum3:time2"
#     "stratum1:tostratum4:time2" "stratum2:tostratum4:time2"
#     "stratum3:tostratum4:time2" "stratum4:tostratum4:time2"
#     "stratum1:tostratum1:time3" "stratum2:tostratum1:time3"
#     "stratum3:tostratum1:time3" "stratum4:tostratum1:time3"
#     "stratum1:tostratum2:time3" "stratum2:tostratum2:time3"
#     "stratum3:tostratum2:time3" "stratum4:tostratum2:time3"
#     "stratum1:tostratum3:time3" "stratum2:tostratum3:time3"
#     "stratum3:tostratum3:time3" "stratum4:tostratum3:time3"
#     "stratum1:tostratum4:time3" "stratum2:tostratum4:time3"
#     "stratum3:tostratum4:time3" "stratum4:tostratum4:time3"

# mwm: the 12 columns that sum to 0 are for
# movement from a given state in time t to itself in time t+1
# That all other columns sum to 1 indicates this is the identity matrix for Psi.
table(colSums(my.desmat.4state))
#
#  0  1
# 12 36

# examine output betas
S.sitetime.p.sitetime.psi.sitetime.4state$results$beta
#                                 estimate        se        lcl        ucl
# S:(Intercept)                  1.3862933 0.1823386  1.0289097  1.7436770
# S:time2                       -0.5389928 0.2647895 -1.0579803 -0.0200054
# S:time3                       -1.5576717 0.1888112 -1.9277416 -1.1876018
# S:stratum2                    -1.7917562 0.1634296 -2.1120783 -1.4714341
# S:stratum3                    -2.5946056 0.1799613 -2.9473297 -2.2418815
# S:stratum4                    -1.2661483 0.1770994 -1.6132632 -0.9190334
# S:time2:stratum2               0.0971590 0.2047437 -0.3041387  0.4984567
# S:time3:stratum2               0.1431859 0.1871186 -0.2235666  0.5099384
# S:time2:stratum3               1.0391202 0.2470936  0.5548167  1.5234237
# S:time3:stratum3               1.9393348 0.1957953  1.5555760  2.3230936
# S:time2:stratum4               0.9510654 0.2223413  0.5152764  1.3868544
# S:time3:stratum4               1.8047547 0.1903935  1.4315836  2.1779259
# p:(Intercept)                  1.0986165 0.2266325  0.6544168  1.5428163
# p:time3                        0.6359712 0.4078675 -0.1634492  1.4353916
# p:stratum2                    -0.4795894 0.6085485 -1.6723445  0.7131658
# p:stratum3                    -1.3397703 0.5426158 -2.4032972 -0.2762434
# p:stratum4                    -0.0526512 0.3053609 -0.6511587  0.5458563
# p:time3:stratum2              -1.0543339 1.0318180 -3.0766973  0.9680294
# p:time3:stratum3              -0.2344663 0.7101180 -1.6262977  1.1573651
# p:time3:stratum4              -0.0237016 0.4910494 -0.9861584  0.9387553
# Psi:stratum2:tostratum1:time1 -1.4939320 0.2802415 -2.0432054 -0.9446586
# Psi:stratum3:tostratum1:time1 -0.4307764 0.3177958 -1.0536561  0.1921034
# Psi:stratum4:tostratum1:time1 -1.7404674 0.1852245 -2.1035074 -1.3774274
# Psi:stratum1:tostratum2:time1 -0.8407776 0.2241687 -1.2801482 -0.4014070
# Psi:stratum3:tostratum2:time1 -0.6931379 0.3893794 -1.4563216  0.0700457
# Psi:stratum4:tostratum2:time1 -1.2704577 0.2557514 -1.7717305 -0.7691849
# Psi:stratum1:tostratum3:time1 -1.0986164 0.2859874 -1.6591518 -0.5380810
# Psi:stratum2:tostratum3:time1 -1.0586162 0.3866898 -1.8165282 -0.3007041
# Psi:stratum4:tostratum3:time1 -1.2098423 0.2858949 -1.7701962 -0.6494883
# Psi:stratum1:tostratum4:time1 -1.6292391 0.1581419 -1.9391972 -1.3192811
# Psi:stratum2:tostratum4:time1 -0.7563309 0.2619453 -1.2697436 -0.2429181
# Psi:stratum3:tostratum4:time1 -1.0498148 0.3356086 -1.7076076 -0.3920220
# Psi:stratum2:tostratum1:time2 -1.1962504 0.3534912 -1.8890932 -0.5034076
# Psi:stratum3:tostratum1:time2 -0.6505809 0.1863614 -1.0158493 -0.2853124
# Psi:stratum4:tostratum1:time2 -1.4663336 0.1287831 -1.7187486 -1.2139187
# Psi:stratum1:tostratum2:time2 -0.6312701 0.3203327 -1.2591222 -0.0034179
# Psi:stratum3:tostratum2:time2 -0.9382609 0.3768503 -1.6768874 -0.1996343
# Psi:stratum4:tostratum2:time2 -0.9067176 0.3262446 -1.5461570 -0.2672783
# Psi:stratum1:tostratum3:time2 -1.0775652 0.1745329 -1.4196496 -0.7354808
# Psi:stratum2:tostratum3:time2 -0.8167671 0.3887033 -1.5786256 -0.0549087
# Psi:stratum4:tostratum3:time2 -1.2431962 0.1643356 -1.5652940 -0.9210984
# Psi:stratum1:tostratum4:time2 -1.3652438 0.1237982 -1.6078884 -1.1225992
# Psi:stratum2:tostratum4:time2 -0.5423277 0.3415565 -1.2117784  0.1271231
# Psi:stratum3:tostratum4:time2 -1.3437311 0.2085390 -1.7524677 -0.9349946
# Psi:stratum2:tostratum1:time3 -0.1556513 0.2115730 -0.5703344  0.2590319
# Psi:stratum3:tostratum1:time3 -0.4652077 0.1134613 -0.6875918 -0.2428235
# Psi:stratum4:tostratum1:time3 -1.3393444 0.0875736 -1.5109887 -1.1677002
# Psi:stratum1:tostratum2:time3 -1.1762097 0.1007885 -1.3737552 -0.9786642
# Psi:stratum3:tostratum2:time3 -1.5308767 0.1681893 -1.8605276 -1.2012257
# Psi:stratum4:tostratum2:time3 -1.7811774 0.1050282 -1.9870326 -1.5753221
# Psi:stratum1:tostratum3:time3 -1.4481434 0.1126660 -1.6689686 -1.2273181
# Psi:stratum2:tostratum3:time3 -0.2141745 0.2150470 -0.6356666  0.2073176
# Psi:stratum4:tostratum3:time3 -1.8084459 0.1067827 -2.0177401 -1.5991517
# Psi:stratum1:tostratum4:time3 -1.1327243 0.0994449 -1.3276363 -0.9378122
# Psi:stratum2:tostratum4:time3  0.4215526 0.1832271  0.0624274  0.7806778
# Psi:stratum3:tostratum4:time3 -1.2642468 0.1522390 -1.5626352 -0.9658584

# examine output reals
S.sitetime.p.sitetime.psi.sitetime.4state$results$real
#                            estimate        se       lcl       ucl fixed note
# S s1 g1 c1 a0 o1 t1       0.7999998 0.0291742 0.7367045 0.8511535           
# S s1 g1 c1 a1 o2 t2       0.7000006 0.0388705 0.6188121 0.7703152           
# S s1 g1 c1 a2 o3 t3       0.4572600 0.0121631 0.4335359 0.4811795           
# S s2 g1 c1 a0 o1 t1       0.4000005 0.0297900 0.3432717 0.4595442           
# S s2 g1 c1 a1 o2 t2       0.3000002 0.0318789 0.2414342 0.3659199           
# S s2 g1 c1 a2 o3 t3       0.1394400 0.0092183 0.1223333 0.1585070           
# S s3 g1 c1 a0 o1 t1       0.2299998 0.0205110 0.1922724 0.2726314           
# S s3 g1 c1 a1 o2 t2       0.3300000 0.0212448 0.2897671 0.3728863           
# S s3 g1 c1 a2 o3 t3       0.3043540 0.0126105 0.2802162 0.3296190           
# S s4 g1 c1 a0 o1 t1       0.5300002 0.0234332 0.4839467 0.5755484           
# S s4 g1 c1 a1 o2 t2       0.6300002 0.0312027 0.5670567 0.6888145           
# S s4 g1 c1 a2 o3 t3       0.5907890 0.0120473 0.5669841 0.6141758           
# p s1 g1 c1 a1 o1 t2       0.7500008 0.0424935 0.6580051 0.8238738           
# p s1 g1 c1 a2 o2 t3       0.8499983 0.0434827 0.7438648 0.9170574           
# p s2 g1 c1 a1 o1 t2       0.6499973 0.1172649 0.4034161 0.8360745           
# p s2 g1 c1 a2 o2 t3       0.5499984 0.1554850 0.2629581 0.8072095           
# p s3 g1 c1 a1 o1 t2       0.4400020 0.0987336 0.2637566 0.6327951           
# p s3 g1 c1 a2 o2 t3       0.5400021 0.0592668 0.4237754 0.6520350           
# p s4 g1 c1 a1 o1 t2       0.7399994 0.0446849 0.6435387 0.8177495           
# p s4 g1 c1 a2 o2 t3       0.8400009 0.0336177 0.7627786 0.8955281           
# 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.2200011 0.0400026 0.1515411 0.3081558           
# Psi s1 to2 g1 c1 a1 o2 t2 0.2500007 0.0602949 0.1507296 0.3850122           
# Psi s1 to2 g1 c1 a2 o3 t3 0.1653327 0.0133047 0.1408749 0.1930823           
# Psi s1 to2 g1 c2 a0 o2 t2 0.2500007 0.0602949 0.1507296 0.3850122           
# Psi s1 to2 g1 c2 a1 o3 t3 0.1653327 0.0133047 0.1408749 0.1930823           
# Psi s1 to2 g1 c3 a0 o3 t3 0.1653327 0.0133047 0.1408749 0.1930823           
# Psi s1 to3 g1 c1 a0 o1 t1 0.1699992 0.0407616 0.1041586 0.2651409           
# Psi s1 to3 g1 c1 a1 o2 t2 0.1599991 0.0259538 0.1154002 0.2175946           
# Psi s1 to3 g1 c1 a2 o3 t3 0.1259677 0.0119430 0.1043592 0.1512947           
# Psi s1 to3 g1 c2 a0 o2 t2 0.1599991 0.0259538 0.1154002 0.2175946           
# Psi s1 to3 g1 c2 a1 o3 t3 0.1259677 0.0119430 0.1043592 0.1512947           
# Psi s1 to3 g1 c3 a0 o3 t3 0.1259677 0.0119430 0.1043592 0.1512947           
# Psi s1 to4 g1 c1 a0 o1 t1 0.1000001 0.0146689 0.0746971 0.1326455           
# Psi s1 to4 g1 c1 a1 o2 t2 0.1199998 0.0152318 0.0932026 0.1531999           
# Psi s1 to4 g1 c1 a2 o3 t3 0.1726808 0.0135801 0.1476618 0.2009395           
# Psi s1 to4 g1 c2 a0 o2 t2 0.1199998 0.0152318 0.0932026 0.1531999           
# Psi s1 to4 g1 c2 a1 o3 t3 0.1726808 0.0135801 0.1476618 0.2009395           
# Psi s1 to4 g1 c3 a0 o3 t3 0.1726808 0.0135801 0.1476618 0.2009395           
# Psi s2 to1 g1 c1 a0 o1 t1 0.1099996 0.0217404 0.0740556 0.1603680           
# Psi s2 to1 g1 c1 a1 o2 t2 0.1300003 0.0252870 0.0879301 0.1880489           
# Psi s2 to1 g1 c1 a2 o3 t3 0.2043892 0.0284443 0.1542141 0.2657597           
# Psi s2 to1 g1 c2 a0 o2 t2 0.1300003 0.0252870 0.0879301 0.1880489           
# Psi s2 to1 g1 c2 a1 o3 t3 0.2043892 0.0284443 0.1542141 0.2657597           
# Psi s2 to1 g1 c3 a0 o3 t3 0.2043892 0.0284443 0.1542141 0.2657597           
# Psi s2 to3 g1 c1 a0 o1 t1 0.1699990 0.0474624 0.0957876 0.2836683           
# Psi s2 to3 g1 c1 a1 o2 t2 0.1899992 0.0413896 0.1216245 0.2843685           
# Psi s2 to3 g1 c1 a2 o3 t3 0.1927710 0.0278642 0.1439244 0.2532905           
# Psi s2 to3 g1 c2 a0 o2 t2 0.1899992 0.0413896 0.1216245 0.2843685           
# Psi s2 to3 g1 c2 a1 o3 t3 0.1927710 0.0278642 0.1439244 0.2532905           
# Psi s2 to3 g1 c3 a0 o3 t3 0.1927710 0.0278642 0.1439244 0.2532905           
# Psi s2 to4 g1 c1 a0 o1 t1 0.2299997 0.0357179 0.1674776 0.3072482           
# Psi s2 to4 g1 c1 a1 o2 t2 0.2499996 0.0409687 0.1784516 0.3384174           
# Psi s2 to4 g1 c1 a2 o3 t3 0.3640277 0.0337348 0.3007924 0.4323369           
# Psi s2 to4 g1 c2 a0 o2 t2 0.2499996 0.0409687 0.1784516 0.3384174           
# Psi s2 to4 g1 c2 a1 o3 t3 0.3640277 0.0337348 0.3007924 0.4323369           
# Psi s2 to4 g1 c3 a0 o3 t3 0.3640277 0.0337348 0.3007924 0.4323369           
# Psi s3 to1 g1 c1 a0 o1 t1 0.2600005 0.0443388 0.1827764 0.3556529           
# Psi s3 to1 g1 c1 a1 o2 t2 0.2400007 0.0309502 0.1846353 0.3057429           
# Psi s3 to1 g1 c1 a2 o3 t3 0.2952811 0.0220844 0.2539067 0.3403222           
# Psi s3 to1 g1 c2 a0 o2 t2 0.2400007 0.0309502 0.1846353 0.3057429           
# Psi s3 to1 g1 c2 a1 o3 t3 0.2952811 0.0220844 0.2539067 0.3403222           
# Psi s3 to1 g1 c3 a0 o3 t3 0.2952811 0.0220844 0.2539067 0.3403222           
# Psi s3 to2 g1 c1 a0 o1 t1 0.2000009 0.0498672 0.1195015 0.3153084           
# Psi s3 to2 g1 c1 a1 o2 t2 0.1800009 0.0529874 0.0979719 0.3073118           
# Psi s3 to2 g1 c1 a2 o3 t3 0.1017236 0.0147045 0.0763061 0.1343759           
# Psi s3 to2 g1 c2 a0 o2 t2 0.1800009 0.0529874 0.0979719 0.3073118           
# Psi s3 to2 g1 c2 a1 o3 t3 0.1017236 0.0147045 0.0763061 0.1343759           
# Psi s3 to2 g1 c3 a0 o3 t3 0.1017236 0.0147045 0.0763061 0.1343759           
# Psi s3 to4 g1 c1 a0 o1 t1 0.1400004 0.0303694 0.0903251 0.2106683           
# Psi s3 to4 g1 c1 a1 o2 t2 0.1200000 0.0201047 0.0858348 0.1653050           
# Psi s3 to4 g1 c1 a2 o3 t3 0.1328059 0.0166316 0.1034542 0.1689162           
# Psi s3 to4 g1 c2 a0 o2 t2 0.1200000 0.0201047 0.0858348 0.1653050           
# Psi s3 to4 g1 c2 a1 o3 t3 0.1328059 0.0166316 0.1034542 0.1689162           
# Psi s3 to4 g1 c3 a0 o3 t3 0.1328059 0.0166316 0.1034542 0.1689162           
# Psi s4 to1 g1 c1 a0 o1 t1 0.0999999 0.0168922 0.0714185 0.1383156           
# Psi s4 to1 g1 c1 a1 o2 t2 0.1200003 0.0153504 0.0930172 0.1534862           
# Psi s4 to1 g1 c1 a2 o3 t3 0.1643396 0.0117134 0.1426490 0.1886028           
# Psi s4 to1 g1 c2 a0 o2 t2 0.1200003 0.0153504 0.0930172 0.1534862           
# Psi s4 to1 g1 c2 a1 o3 t3 0.1643396 0.0117134 0.1426490 0.1886028           
# Psi s4 to1 g1 c3 a0 o3 t3 0.1643396 0.0117134 0.1426490 0.1886028           
# Psi s4 to2 g1 c1 a0 o1 t1 0.1600008 0.0346870 0.1030236 0.2400560           
# Psi s4 to2 g1 c1 a1 o2 t2 0.2100006 0.0542959 0.1227771 0.3354911           
# Psi s4 to2 g1 c1 a2 o3 t3 0.1056468 0.0097103 0.0880702 0.1262458           
# Psi s4 to2 g1 c2 a0 o2 t2 0.2100006 0.0542959 0.1227771 0.3354911           
# Psi s4 to2 g1 c2 a1 o3 t3 0.1056468 0.0097103 0.0880702 0.1262458           
# Psi s4 to2 g1 c3 a0 o3 t3 0.1056468 0.0097103 0.0880702 0.1262458           
# Psi s4 to3 g1 c1 a0 o1 t1 0.1699993 0.0408285 0.1040719 0.2653222           
# Psi s4 to3 g1 c1 a1 o2 t2 0.1499995 0.0234562 0.1095654 0.2019708           
# Psi s4 to3 g1 c1 a2 o3 t3 0.1028049 0.0096381 0.0853916 0.1232906           
# Psi s4 to3 g1 c2 a0 o2 t2 0.1499995 0.0234562 0.1095654 0.2019708           
# Psi s4 to3 g1 c2 a1 o3 t3 0.1028049 0.0096381 0.0853916 0.1232906           
# Psi s4 to3 g1 c3 a0 o3 t3 0.1028049 0.0096381 0.0853916 0.1232906

# Examine the estimated transition matrix
Psilist = get.real(S.sitetime.p.sitetime.psi.sitetime.4state,"Psi",vcv=TRUE)
Psivalues = Psilist$estimates
TransitionMatrix(Psivalues[Psivalues$time==1,])
#            1         2         3         4
# 1 0.50999966 0.2200011 0.1699992 0.1000001
# 2 0.10999961 0.4900017 0.1699990 0.2299997
# 3 0.26000050 0.2000009 0.3999982 0.1400004
# 4 0.09999989 0.1600008 0.1699993 0.5700001


Here is the 4-state deterministic data set:

Code: Select all
1100 133.38319580856  ;
1101 6.204165912  ;
1102 2.6888684508  ;
1103 4.02340119552  ;
1104 4.57672863312  ;
1110 46.443835746  ;
1111 20.97391779  ;
1112 6.46931124  ;
1113 4.92899904  ;
1114 6.756836184  ;
1120 25.3456434  ;
1121 0.83939625  ;
1122 0.98076825  ;
1123 0.7916832  ;
1124 1.4950089  ;
1130 12.87423704448  ;
1131 1.6632133056  ;
1132 0.5729730048  ;
1133 2.64840855552  ;
1134 0.7480480896  ;
1140 8.83542201696  ;
1141 2.0963051424  ;
1142 1.3476247344  ;
1143 1.31137284096  ;
1144 8.00063526528  ;
1200 88.806535774528  ;
1201 0.75615327216  ;
1202 0.45028711728  ;
1203 0.729684713472  ;
1204 1.01986712256  ;
1210 2.0582654664  ;
1211 0.929507436  ;
1212 0.286702416  ;
1213 0.218439936  ;
1214 0.2994447456  ;
1220 6.9848901408  ;
1221 0.23132538  ;
1222 0.270285444  ;
1223 0.2181763584  ;
1224 0.4120026768  ;
1230 2.449530955872  ;
1231 0.31645311984  ;
1232 0.10901734272  ;
1233 0.503902384128  ;
1234 0.14232819744  ;
1240 2.9492655192  ;
1241 0.699747048  ;
1242 0.449837388  ;
1243 0.4377364992  ;
1244 2.6706135456  ;
1300 44.9170319710746  ;
1301 0.632161584384  ;
1302 0.2600391074688  ;
1303 0.70493194371072  ;
1304 0.44671155336192  ;
1310 2.186389446912  ;
1311 0.98736789888  ;
1312 0.30454921728  ;
1313 0.23203749888  ;
1314 0.318084738048  ;
1320 1.682371392768  ;
1321 0.0557167248  ;
1322 0.06510059424  ;
1323 0.052549668864  ;
1324 0.099234419328  ;
1330 3.41228587569408  ;
1331 0.4408307266176  ;
1332 0.1518651307008  ;
1333 0.70195438190592  ;
1334 0.1982683650816  ;
1340 0.81454176308736  ;
1341 0.1932593693184  ;
1342 0.1242381659904  ;
1343 0.12089608667136  ;
1344 0.73758237493248  ;
1400 28.3613611584832  ;
1401 0.797536470528  ;
1402 0.441465964128  ;
1403 0.6901390282752  ;
1404 1.4857485785856  ;
1410 2.06468716608  ;
1411 0.9324074592  ;
1412 0.2875969152  ;
1413 0.2191214592  ;
1414 0.30037900032  ;
1420 3.70702398528  ;
1421 0.122769108  ;
1422 0.1434460104  ;
1423 0.11579065344  ;
1424 0.21865824288  ;
1430 2.101529870496  ;
1431 0.27149511312  ;
1432 0.09352941696  ;
1433 0.432313749504  ;
1434 0.12210784992  ;
1440 6.6664125335808  ;
1441 1.581682781952  ;
1442 1.016796074112  ;
1443 0.9894436651008  ;
1444 6.0365577452544  ;
1000 359.411903316858  ;
1001 3.559999564128  ;
1002 1.6248203046432  ;
1003 2.87370827563008  ;
1004 3.16529717874048  ;
1010 20.097683475408  ;
1011 9.07606260792  ;
1012 2.79947096352  ;
1013 2.13293025792  ;
1014 2.923891895232  ;
1020 15.653310332352  ;
1021 0.5184058572  ;
1022 0.60571631736  ;
1023 0.488938576896  ;
1024 0.923308116192  ;
1030 10.6916752114099  ;
1031 1.3812497322624  ;
1032 0.4758372283392  ;
1033 2.19942541099008  ;
1034 0.6212319369984  ;
1040 7.91214929011584  ;
1041 1.8772481056896  ;
1042 1.2068023536576  ;
1043 1.17433866461184  ;
1044 7.16459502592512  ;
2100 14.38446229308  ;
2101 0.669076716  ;
2102 0.2899760094  ;
2103 0.43389620736  ;
2104 0.49356877416  ;
2110 5.008648953  ;
2111 2.261893095  ;
2112 0.69767082  ;
2113 0.53155872  ;
2114 0.728678412  ;
2120 2.7333537  ;
2121 0.090523125  ;
2122 0.105769125  ;
2123 0.0853776  ;
2124 0.16122645  ;
2130 1.38839811264  ;
2131 0.1793661408  ;
2132 0.0617912064  ;
2133 0.28561268736  ;
2134 0.0806718528  ;
2140 0.95283962928  ;
2141 0.2260721232  ;
2142 0.1453320792  ;
2143 0.14142256128  ;
2144 0.86281360704  ;
2200 98.898187567088  ;
2201 0.84207978036  ;
2202 0.50145610788  ;
2203 0.812603430912  ;
2204 1.13576111376  ;
2210 2.2921592694  ;
2211 1.035133281  ;
2212 0.319282236  ;
2213 0.243262656  ;
2214 0.3334725576  ;
2220 7.7786276568  ;
2221 0.257612355  ;
2222 0.300999699  ;
2223 0.2429691264  ;
2224 0.4588211628  ;
2230 2.727886746312  ;
2231 0.35241370164  ;
2232 0.12140567712  ;
2233 0.561164018688  ;
2234 0.15850185624  ;
2240 3.2844093282  ;
2241 0.779263758  ;
2242 0.500955273  ;
2243 0.4874792832  ;
2244 2.9740923576  ;
2300 22.4585159855373  ;
2301 0.316080792192  ;
2302 0.1300195537344  ;
2303 0.35246597185536  ;
2304 0.22335577668096  ;
2310 1.093194723456  ;
2311 0.49368394944  ;
2312 0.15227460864  ;
2313 0.11601874944  ;
2314 0.159042369024  ;
2320 0.841185696384  ;
2321 0.0278583624  ;
2322 0.03255029712  ;
2323 0.026274834432  ;
2324 0.049617209664  ;
2330 1.70614293784704  ;
2331 0.2204153633088  ;
2332 0.0759325653504  ;
2333 0.35097719095296  ;
2334 0.0991341825408  ;
2340 0.40727088154368  ;
2341 0.0966296846592  ;
2342 0.0621190829952  ;
2343 0.06044804333568  ;
2344 0.36879118746624  ;
2400 32.6155653322557  ;
2401 0.9171669411072  ;
2402 0.5076858587472  ;
2403 0.79365988251648  ;
2404 1.70861086537344  ;
2410 2.374390240992  ;
2411 1.07226857808  ;
2412 0.33073645248  ;
2413 0.25198967808  ;
2414 0.345435850368  ;
2420 4.263077583072  ;
2421 0.1411844742  ;
2422 0.16496291196  ;
2423 0.133159251456  ;
2424 0.251456979312  ;
2430 2.4167593510704  ;
2431 0.312219380088  ;
2432 0.107558829504  ;
2433 0.4971608119296  ;
2434 0.140424027408  ;
2440 7.66637441361792  ;
2441 1.8189351992448  ;
2442 1.1693154852288  ;
2443 1.13786021486592  ;
2444 6.94204140704256  ;
2000 698.090779747143  ;
2001 1.4009856348208  ;
2002 0.7105290416784  ;
2003 1.30963429001216  ;
2004 1.6606803263456  ;
2010 5.129373216352  ;
2011 2.31641186448  ;
2012 0.71448689088  ;
2013 0.54437096448  ;
2014 0.746241863808  ;
2020 7.668047762144  ;
2021 0.2539501734  ;
2022 0.29672072892  ;
2023 0.239515110912  ;
2024 0.452298624624  ;
2030 4.95224777710656  ;
2031 0.6397772828832  ;
2032 0.2204017433856  ;
2033 1.01874583609344  ;
2034 0.2877467205312  ;
2040 5.2980770780784  ;
2041 1.257029511696  ;
2042 0.808090400376  ;
2043 0.7863522960384  ;
2044 4.7975051138112  ;
3100 19.549791934686  ;
3101 0.9093360822  ;
3102 0.39410375823  ;
3103 0.589704390912  ;
3104 0.670804833972  ;
3110 6.80720925885  ;
3111 3.07411834275  ;
3112 0.948198069  ;
3113 0.722436624  ;
3114 0.9903402054  ;
3120 3.714876165  ;
3121 0.12302915625  ;
3122 0.14374985625  ;
3123 0.11603592  ;
3124 0.2191214025  ;
3130 1.886959253088  ;
3131 0.24377489136  ;
3132 0.08397986688  ;
3133 0.388173606912  ;
3134 0.10964038176  ;
3140 1.294995677976  ;
3141 0.30725256744  ;
3142 0.19751950764  ;
3143 0.192206117376  ;
3144 1.172642129568  ;
3200 23.210799122888  ;
3201 0.19763096886  ;
3202 0.11768867838  ;
3203 0.190713050112  ;
3204 0.26655617976  ;
3210 0.5379557469  ;
3211 0.2429394435  ;
3212 0.074933586  ;
3213 0.057092256  ;
3214 0.0782639676  ;
3220 1.8255962868  ;
3221 0.0604600425  ;
3222 0.0706427865  ;
3223 0.0570233664  ;
3224 0.1076825178  ;
3230 0.640218318012  ;
3231 0.08270933814  ;
3232 0.02849316912  ;
3233 0.131701759488  ;
3234 0.03719941524  ;
3240 0.7708307607  ;
3241 0.182888433  ;
3242 0.1175711355  ;
3243 0.1144084032  ;
3244 0.6980012676  ;
3300 30.3850510392563  ;
3301 0.427638718848  ;
3302 0.1759088079936  ;
3303 0.47686572662784  ;
3304 0.30218722727424  ;
3310 1.479028155264  ;
3311 0.66792534336  ;
3312 0.20601858816  ;
3313 0.15696654336  ;
3314 0.215174969856  ;
3320 1.138074765696  ;
3321 0.0376907256  ;
3322 0.04403863728  ;
3323 0.035548305408  ;
3324 0.067129166016  ;
3330 2.30831103355776  ;
3331 0.2982090209472  ;
3332 0.1027322942976  ;
3333 0.47485149364224  ;
3334 0.1341227175552  ;
3340 0.55101354561792  ;
3341 0.1307342792448  ;
3342 0.0840434652288  ;
3343 0.08178264686592  ;
3344 0.49895278304256  ;
3400 11.4154478662895  ;
3401 0.32100842938752  ;
3402 0.17769005056152  ;
3403 0.277780958880768  ;
3404 0.598013802880704  ;
3410 0.8310365843472  ;
3411 0.375294002328  ;
3412 0.115757758368  ;
3413 0.088196387328  ;
3414 0.1209025476288  ;
3420 1.4920771540752  ;
3421 0.04941456597  ;
3422 0.057737019186  ;
3423 0.0466057380096  ;
3424 0.0880099427592  ;
3430 0.84586577287464  ;
3431 0.1092767830308  ;
3432 0.0376455903264  ;
3433 0.17400628417536  ;
3434 0.0491484095928  ;
3440 2.68323104476627  ;
3441 0.63662731973568  ;
3442 0.40926041983008  ;
3443 0.398251075203072  ;
3444 2.4297144924649  ;
3000 831.697436091167  ;
3001 1.06658290875248  ;
3002 0.48105432306288  ;
3003 1.00377850366259  ;
3004 0.961846745354656  ;
3010 4.7331235293188  ;
3011 2.137466516262  ;
3012 0.659291997672  ;
3013 0.502317712512  ;
3014 0.6885938642352  ;
3020 4.1940075573088  ;
3021 0.13889701518  ;
3022 0.162290196684  ;
3023 0.1310018164224  ;
3024 0.2473828944048  ;
3030 4.2087658468716  ;
3031 0.543727393902  ;
3032 0.187312786416  ;
3033 0.8658013238784  ;
3034 0.244547248932  ;
3040 2.49077476400861  ;
3041 0.59096486125152  ;
3042 0.37990598223312  ;
3043 0.369686289143808  ;
3044 2.25544183136294  ;
4100 17.32673867121  ;
4101 0.805933317  ;
4102 0.34928928405  ;
4103 0.52264770432  ;
4104 0.59452602342  ;
4110 6.03314532975  ;
4111 2.72455304625  ;
4112 0.840376215  ;
4113 0.64028664  ;
4114 0.877726269  ;
4120 3.292448775  ;
4121 0.10903921875  ;
4122 0.12740371875  ;
4123 0.1028412  ;
4124 0.1942045875  ;
4130 1.67238863568  ;
4131 0.2160546696  ;
4132 0.0744303168  ;
4133 0.34403346432  ;
4134 0.0971729136  ;
4140 1.14773864436  ;
4141 0.2723141484  ;
4142 0.1750590954  ;
4143 0.17034990336  ;
4144 1.03929820848  ;
4200 42.7886036004544  ;
4201 0.364328394768  ;
4202 0.216956520144  ;
4203 0.3515753619456  ;
4204 0.491390522688  ;
4210 0.99170972472  ;
4211 0.4478535828  ;
4212 0.1381384368  ;
4213 0.1052483328  ;
4214 0.14427792288  ;
4220 3.36544706784  ;
4221 0.111456774  ;
4222 0.1302284412  ;
4223 0.10512133632  ;
4224 0.19851038064  ;
4230 1.1802285514656  ;
4231 0.152472866832  ;
4232 0.052526537856  ;
4233 0.2427893305344  ;
4234 0.068576313312  ;
4240 1.42100975016  ;
4241 0.3371508504  ;
4242 0.2167398324  ;
4243 0.21090940416  ;
4244 1.28675016288  ;
4300 29.7575336808369  ;
4301 0.4188070496544  ;
4302 0.17227590869808  ;
4303 0.467017412708352  ;
4304 0.295946404102272  ;
4310 1.4484830085792  ;
4311 0.654131233008  ;
4312 0.201763856448  ;
4313 0.153724843008  ;
4314 0.2107311389568  ;
4320 1.1145710477088  ;
4321 0.03691233018  ;
4322 0.043129143684  ;
4323 0.0348141556224  ;
4324 0.0657428028048  ;
4330 2.26063939264733  ;
4331 0.29205035638416  ;
4332 0.10061064908928  ;
4333 0.465044778012672  ;
4334 0.13135279186656  ;
4340 0.539633918045376  ;
4341 0.12803433217344  ;
4342 0.08230778496864  ;
4343 0.080093657419776  ;
4344 0.488648323392768  ;
4400 107.099590074722  ;
4401 3.01169709683136  ;
4402 1.66708584703836  ;
4403 2.60613750552422  ;
4404 5.61055806988387  ;
4410 7.7967749109096  ;
4411 3.521003667804  ;
4412 1.086037851024  ;
4413 0.827457410304  ;
4414 1.1343061999584  ;
4420 13.9986493244136  ;
4421 0.463606844085  ;
4422 0.541687996773  ;
4423 0.4372544550528  ;
4424 0.8257081896756  ;
4430 7.93590217346052  ;
4431 1.0252334209194  ;
4432 0.3531904607952  ;
4433 1.63252479656448  ;
4434 0.4611097682604  ;
4440 25.1740403299345  ;
4441 5.97282960534624  ;
4442 3.83967617486544  ;
4443 3.7363866403369  ;
4444 22.7955511855169  ;
4000 574.318407267022  ;
4001 2.05601226589024  ;
4002 1.03824561557556  ;
4003 1.87358144408422  ;
4004 2.81070645742746  ;
4010 7.1279772175612  ;
4011 3.218976334938  ;
4012 0.992878869528  ;
4013 0.756479138688  ;
4014 1.0370068192848  ;
4020 9.2466360577976  ;
4021 0.306229812735  ;
4022 0.357805360143  ;
4023 0.2888230654848  ;
4024 0.5454114138396  ;
4030 6.85843792177975  ;
4031 0.88603659911844  ;
4032 0.30523748869152  ;
4033 1.41087550328525  ;
4034 0.39850449912504  ;
4040 10.6794785380409  ;
4041 2.53382868803232  ;
4042 1.62888987087792  ;
4043 1.58507178078413  ;
4044 9.6704619702647  ;
0100 435.89279676  ;
0101 20.275052  ;
0102 8.7871518  ;
0103 13.14836992  ;
0104 14.95662952  ;
0110 151.777241  ;
0111 68.542215  ;
0112 21.14154  ;
0113 16.10784  ;
0114 22.081164  ;
0120 82.8289  ;
0121 2.743125  ;
0122 3.205125  ;
0123 2.5872  ;
0124 4.88565  ;
0130 42.07267008  ;
0131 5.4353376  ;
0132 1.8724608  ;
0133 8.65492992  ;
0134 2.4446016  ;
0140 28.87392816  ;
0141 6.8506704  ;
0142 4.4040024  ;
0143 4.28553216  ;
0144 26.14586688  ;
0200 776.28090712  ;
0201 6.6097314  ;
0202 3.9360762  ;
0203 6.37836288  ;
0204 8.9149224  ;
0210 17.991831  ;
0211 8.125065  ;
0212 2.50614  ;
0213 1.90944  ;
0214 2.617524  ;
0220 61.056732  ;
0221 2.022075  ;
0222 2.362635  ;
0223 1.907136  ;
0224 3.601422  ;
0230 21.41198388  ;
0231 2.7661986  ;
0232 0.9529488  ;
0233 4.40474112  ;
0234 1.2441276  ;
0240 25.780293  ;
0241 6.11667  ;
0242 3.932145  ;
0243 3.826368  ;
0244 23.344524  ;
0300 750.618849784  ;
0301 10.5641976  ;
0302 4.34557332  ;
0303 11.780279808  ;
0304 7.465099488  ;
0310 36.5372568  ;
0311 16.500132  ;
0312 5.089392  ;
0313 3.877632  ;
0314 5.3155872  ;
0320 28.1144952  ;
0321 0.931095  ;
0322 1.087911  ;
0323 0.8781696  ;
0324 1.6583292  ;
0330 57.023493912  ;
0331 7.36682364  ;
0332 2.53785312  ;
0333 11.730521088  ;
0334 3.31330824  ;
0340 13.611994704  ;
0341 3.22960176  ;
0342 2.07617256  ;
0343 2.020322304  ;
0344 12.325908672  ;
0400 479.077046596  ;
0401 13.47189984  ;
0402 7.45719534  ;
0403 11.657753856  ;
0404 25.097104368  ;
0410 34.8764724  ;
0411 15.750126  ;
0412 4.858056  ;
0413 3.701376  ;
0414 5.0739696  ;
0420 62.6186484  ;
0421 2.0738025  ;
0422 2.4230745  ;
0423 1.9559232  ;
0424 3.6935514  ;
0430 35.49881538  ;
0431 4.5860661  ;
0432 1.5798888  ;
0433 7.30259712  ;
0434 2.0626326  ;
0440 112.608319824  ;
0441 26.71761456  ;
0442 17.17560936  ;
0443 16.713575424  ;
0444 101.968880832  ;
0010 542.74  ;
0011 245.1  ;
0012 75.6  ;
0013 57.6  ;
0014 78.96  ;
0020 860.56  ;
0021 28.5  ;
0022 33.3  ;
0023 26.88  ;
0024 50.76  ;
0030 695.646  ;
0031 89.87  ;
0032 30.96  ;
0033 143.104  ;
0034 40.42  ;
0040 409.211  ;
0041 97.09  ;
0042 62.415  ;
0043 60.736  ;
0044 370.548  ;
markmiller
 
Posts: 49
Joined: Fri Nov 08, 2013 6:23 pm


Return to RMark

Who is online

Users browsing this forum: No registered users and 18 guests