multi-season multi-state occupancy model

[markmiller]

I am attempting to use the multi-season multi-state occupancy model (MacKenzie et al. 2009, Ecology) in MARK. Below is a small data set of 9 capture histories that I have been using to understand parameter definitions in MARK.

Code: Select all

0000 100 ;
00AA 200 ;
00BB 700 ;
AA00 300 ;
AAAA 300 ;
AABB 400 ;
BB00 600 ;
BBAA 300 ;
BBBB 100 ;

The data set consists of two primary sampling periods and two secondary sampling visits in each primary period.

If I use the 'Robust Design Multiple-State Occupancy Estimation Conditional Binomial' model I obtain the estimates I expect, without using any constraints or fixed values.

However, if I use the 'Robust Design Multiple-State Occupancy Estimation General' model the parameter estimates are confusing. Perhaps I do not understand the model correctly, but I am expecting the following estimates of initial occupancy:

PsiA = 0.33
PsiB = 0.33

For the transition matrix (top of right-hand column of page 826 of MacKenzie et al. 2009) I am expecting:

0.1, 0.2, 0.7
0.3, 0.3, 0.4
0.6, 0.3, 0.1

or using parameter names from MARK,

Psi 0 to A = 0.2
Psi 0 to B = 0.7
Psi A to 0 = 0.3
Psi A to B = 0.4
Psi B to 0 = 0.6
Psi B to A = 0.3

For detection probability I am expecting:

p[A,A] = 1.0
p[A,B] = 0.0
p[B,B] = 1.0

I actually get these estimates if I fix all of the detection probabilities to the above values. Although, that model has a delta AIC of > 9400.

If I use a model with only three detection probability parameters (p[A,A], p[A,B], p[B,B]), none of which are fixed, I obtain the following estimates:

Code: Select all

                         Real Function Parameters
                                                               95% Confidence Interval
  Parameter                  Estimate       Standard Error      Lower           Upper
 --------------------------  --------------  --------------  --------------  --------------
     1:Phi0                  0.0000000       0.0000000       0.0000000       0.0000000
     2:Phi0                  0.6666667       0.0086066       0.6495898       0.6833170
     3:Psi 0 to A            1.0000000       0.9424320E-008  1.0000000       1.0000000
     4:Psi 0 to B            1.0000000       0.0000000       1.0000000       1.0000000
     5:Psi A to 0            0.3419188       657.95658       0.5122786E-304  1.0000000
     6:Psi A to B            0.7695708       63.043177       0.8078489E-302  1.0000000
     7:Psi B to 0            0.4002814       84.596508       0.7128231E-300  1.0000000
     8:Psi B to A            0.0497186       84.596510       0.5158571E-305  1.0000000
     9:p[A,A] A:Strata 1 Se  0.3885781E-015  0.1395841E-008  -.2735849E-008  0.2735850E-008
    10:p[A,B] B:Strata 2 Se  1.0000000       0.0000000       1.0000000       1.0000000 
    11:p[B,B] B:Strata 2 Se  1.0000000       0.0000000       1.0000000       1.0000000

This model has a delta AIC of approximately 3800.

1. Are my expected values correct? In other words, do I understand the model correctly?

2. Should I be concerned that the estimates MARK returns immediately above are so different from what I expect unless I fix every detection probability?

3. Should I be concerned that the model that returns the values I expect has a delta AIC of > 9400?

I am particularly confused as to why MARK estimates p[A,A] = 0 in the last model and why MARK estimates any of the initial occupancies to be 0 in that model.

I tried the above example data set after a more realistic simulated data set yielded parameter estimates that were more-or-less equally far removed from what I expected them to be. That concerns me more than do the model results above because with the more realistic simulated data set I could not fix any of the detection probabilities to a known value.

Thank you for any advice. I will try to provide additional information if that will help.

[markmiller]

I have written my own likelihood function in R for the general multi-state multi-season occupancy model (MacKenzie et al. 2009) as I understand it. When I analyze the above data set, reproduced here:

Code: Select all

0000 100 ;
00AA 200 ;
00BB 700 ;
AA00 300 ;
AAAA 300 ;
AABB 400 ;
BB00 600 ;
BBAA 300 ;
BBBB 100 ;

that likelihood function returns the following point estimates, which match my expected values to several decimal places:

Code: Select all

psi.A = 0.333333812     # initial occupancy in state A
psi.B = 0.333333499     # initial occupancy in state B

phi.00_12 = 0.099999621  # transition prob from state 0 in year 1 to state 0 in year 2
phi.0A_12 = 0.199997508  # transition prob from state 0 in year 1 to state A in year 2

phi.A0_12 = 0.299999159  # transition prob from state A in year 1 to state 0 in year 2
phi.AA_12 = 0.299999255  # transition prob from state A in year 1 to state A in year 2

phi.B0_12 = 0.599999448  # transition prob from state B in year 1 to state 0 in year 2
phi.BA_12 = 0.299999830  # transition prob from state B in year 1 to state A in year 2

p.yr1_visit1_1of1 = 0.999999999
p.yr1_visit1_1of2 = 0.000000001
p.yr1_visit1_2of2 = 0.999999999

p.yr1_visit2_1of1 = 0.999999999
p.yr1_visit2_1of2 = 0.000000001
p.yr1_visit2_2of2 = 0.999999999

p.yr2_visit1_1of1 = 0.999999999
p.yr2_visit1_1of2 = 0.000000001
p.yr2_visit1_2of2 = 0.999999999

p.yr2_visit2_1of1 = 0.999999999
p.yr2_visit2_1of2 = 0.000000001
p.yr2_visit2_2of2 = 0.999999999

These estimates were obtained without using any constrained or fixed values. Standard errors were not estimated presumably because some point estimates were effectively 0 or 1. Standard errors were obtained with a different data set that did not return point estimates of 0 or 1.

Nevertheless, this done not mean my likelihood function accurately represents the actual MacKenzie et al. 2009 model. It only suggests that my likelihood function is doing what I intended it to do. My understanding of the MacKenzie et al. 2009 model might still be incorrect.

However, if my understanding of the MacKenzie et al. 2009 model is correct, that I obtained estimates matching my expected values increases my concern regarding why I have not yet been able to obtain similar estimates with Program MARK (unless I fix the value of every detection probability to a known true value).

[markmiller]

My first question in my first post in this thread was whether my expected values were correct.

I have now used MacKenzie et al's 2009 Bayesian code for the general model to analyze several fake data sets created according to my understanding of how the model works. Two of these data set consisted of expected values such that the data matched the model without error (if my understanding of the model was correct).

In those two cases all parameter point estimates matched expected values to within 0.02. This suggests to me that my understanding of the general model and my expected values are correct.

My likelihood function in R appears to be less stable (sensitive to initial values) and slower than the Bayesian version run with JAGS. So, I might try using JAGS with my real data.

I am not sure whether I am using the model correctly in Program MARK. If I obtain parameter estimates with MARK that more-or-less match those from JAGS I will post a follow-up.

[markmiller]

In my last post I wrote that I would post an update if I was able to obtain parameter estimates with MARK that match those obtained with JAGS (or match my expected values). I have now been able to do so with a model that has only three detection probability parameters (p[A,A], p[A,B], and p[B,B]), none of which are constrained to a fixed value.

To do this using the data set in the original post:

Code: Select all

0000 100 ;
00AA 200 ;
00BB 700 ;
AA00 300 ;
AAAA 300 ;
AABB 400 ;
BB00 600 ;
BBAA 300 ;
BBBB 100 ;

Step 1. I first created the most general model using a PIM Chart (there are 20 parameters).

Step 2. Then I set all p[A,A] equal, all p[A,B] equal and all p[B,B] equal (there are now 3 detection probability parameters and 11 parameters total).

Step 3. Then I ran the model in Step 2 but using the 'Parm-Specific' link and specifying five different MLogit links as follows:

1. Phi0: MLogit(1)
2. Phi0: MLogit(1)
3. Psi 0 to A: MLogit(2)
4. Psi 0 to B: MLogit(2)
5. Psi A to 0: MLogit(3)
6. Psi A to B: MLogit(3)
7. Psi B to 0: MLogit(4)
8. Psi B to A: MLogit(4)
9. p[A,A] A Strata 1 Session 1: Logit
10. p[A,B] B Strata 2 Session 1: MLogit(5)
11. p[B,B] B Strata 2 Session 1: MLogit(5)

The parameter estimates are as follows:

Code: Select all

            PARM-SPECIFIC Link Function Parameters of {general PIM Chart with 3 parameters}
                                                              95% Confidence Interval
 Parameter                    Beta         Standard Error      Lower           Upper
 -------------------------  --------------  --------------  --------------  --------------
    1:Phi0                  0.7331316E-007  0.0447214       -0.0876538      0.0876540     
    2:Phi0                  0.3865503E-007  0.0447214       -0.0876538      0.0876539     
    3:Psi 0 to A            0.6931461       0.1224746       0.4530960       0.9331963     
    4:Psi 0 to B            1.9459094       0.1069045       1.7363765       2.1554423     
    5:Psi A to 0            0.4904879E-006  0.0816497       -0.1600330      0.1600340     
    6:Psi A to B            0.2876823       0.0763763       0.1379847       0.4373799     
    7:Psi B to 0            1.7917600       0.1080124       1.5800558       2.0034643     
    8:Psi B to A            1.0986129       0.1154701       0.8722915       1.3249342     
    9:p[A,A] A:Strata 1 Se  27.562751       0.0000000       27.562751       27.562751     
   10:p[A,B] B:Strata 2 Se  0.8522012       0.0000000       0.8522012       0.8522012     
   11:p[B,B] B:Strata 2 Se  34.321182       0.0000000       34.321182       34.321182     


            Real Function Parameters of {general PIM Chart with 3 parameters}
                                                               95% Confidence Interval
  Parameter                  Estimate       Standard Error      Lower           Upper
 --------------------------  --------------  --------------  --------------  --------------
     1:Phi0                  0.3333333       0.0086066       0.3166830       0.3504102                           
     2:Phi0                  0.3333333       0.0086066       0.3166830       0.3504102                           
     3:Psi 0 to A            0.1999999       0.0126491       0.1763542       0.2259464                           
     4:Psi 0 to B            0.7000000       0.0144914       0.6708529       0.7276140                           
     5:Psi A to 0            0.3000001       0.0144914       0.2723860       0.3291472                           
     6:Psi A to B            0.4000000       0.0154919       0.3700546       0.4307118                           
     7:Psi B to 0            0.6000000       0.0154919       0.5692882       0.6299454                           
     8:Psi B to A            0.3000000       0.0144914       0.2723860       0.3291471                           
     9:p[A,A] A:Strata 1 Se  1.0000000       0.0000000       1.0000000       1.0000000                           
    10:p[A,B] B:Strata 2 Se  0.2914783E-014  0.0000000       0.2914783E-014  0.2914783E-014                     
    11:p[B,B] B:Strata 2 Se  1.0000000       0.0000000       1.0000000       1.0000000                           

 Attempted ordering of parameters by estimatibility:
   2  1  7  6  4  5  8  3 10  9 11
 Beta number 11 is a singular value.

Step 4. Run the model in Step 3 using a reduced design matrix and the same MLogit links.

If I find any errors in the above steps I will post a correction. I am not sure why this model has a delta AICc of > 5000.