MSORD and parameter identifiability

Forum for discussion of general questions related to study design and/or analysis of existing data - software neutral.

MSORD and parameter identifiability

Postby violetblue » Wed Oct 02, 2019 5:43 am

Hi, thanks for your help. I have some turtle recapture data and am hoping to get some advice/comments on our model and whether we should resign ourselves to something simpler.

Ideally the model will estimate Survival, Psi and the derived parameters Residency Time and Abundance all time-varying because we would like to monitor these parameters over time. I've tried several models and this one appears the best supported by the data, but as it has many poorly estimated parameters,  I think the AIC may be incorrect because the inestimable parameters aren't counted. Can you suggest how to find AIC in this situation?

Once I have found the 'best' (lowest AIC/as close to the ideal as we can get/no estimation problems) I will use MCMC to get estimates robust to sampling variation.

I think that since the poorly estimated parameters are close to 0 or 1 this show problems with estimation near the boundary, and not confounding, which I think the constant Psi in 'nest to skip' has taken care of. Then finding the profile likelihood CI of the cloned data would be the next step.

This is the model specification:  
    a. MSORD: multi-state open robust design. Two states, nest and skip.
    a. S: 2 ageclass survival (transients & residents), time-varying.
    b. Psi: time-varying 'nest to skip' transition probabilities but constant 'skip to nest' to avoid parameter confounding.
    c. p, pent & phi probability of zero in the skipped state.
    d. all other parameters time-varying.

In RMark:
Code: Select all
tagging_test<-convert.inp("tagging.inp") # 15 years, 5 secondary occasions each
t.int<-c(rep(c(0,0,0,0,1),14),c(0,0,0,0))
ordms.process<-process.data(tagging_test,model="ORDMS",time.interval=t.int,begin.time=2005,           strata.labels=c("1","2")) # 2 states; 1 nest; 2 skip
ordms.ddl=make.design.data(ordms.process)
# transient/resident effect on Survival
ordms.ddl=add.design.data(ordms.process,ordms.ddl,parameter="S",type="age",right=FALSE,bins=c(0,0.5,13),name="ageclass",replace=TRUE)

# parameter indexes for fixing values to 0.
up=as.numeric(row.names(ordms.ddl$p[ordms.ddl$p$stratum=="2",]))
upent=as.numeric(row.names(ordms.ddl$pent[ordms.ddl$pent$stratum=="2",]))
uphi=as.numeric(row.names(ordms.ddl$Phi[ordms.ddl$Phi$stratum=="2",]))

# Survival                                   
S.a2.time=list(formula=~-1+ageclass:time)  # time-varying 2 ageclass survival

# Transition probabilities
Psi.dot.time=list(formula=~-1+to1+to2:time)  

p.time.time=list(formula=~-1+time:session:stratum,fixed=list(index=up,value=0))  
pent.time.time=list(formula=~-1+time:session:stratum,fixed=list(index=upent,value=0))    
Phi.time.time=list(formula=~-1+time:session:stratum,fixed=list(index=uphi,value=0))   
 

These are the poorly estimated parameters:
95% Confidence Interval
[list=]  Parameter                  Estimate       Standard Error      Lower           Upper
 --------------------------  --------------  --------------  --------------  --------------
 11:S s1 g1 c2005 a10 t2   0.9996639       0.0053729       0.7291067E-10   1.0000000                    
 13:S s1 g1 c2005 a12 t2   1.0000000       0.2552875E-06   0.9999995       1.0000005                          
 27:Phi s1 g1 c1 a0 s200   1.0000000       0.1681379E-06   0.9999997       1.0000003                          
 28:Phi s1 g1 c1 a0 s200   1.0000000       0.1443151E-04  0.1239199E-289   1.0000000                         
 30:Phi s1 g1 c1 a0 s200   1.0000000       0.0000000       1.0000000       1.0000000     
 28:Phi s1 g1 c1 a0 s200   1.0000000       0.1443151E-04  0.1239199E-289   1.0000000                       
 30:Phi s1 g1 c1 a0 s200   1.0000000       0.0000000       1.0000000       1.0000000                          
 38:Phi s1 g1 c1 a0 s200   1.0000000       0.2656015E-05   0.9999948       1.0000052    
 47:Phi s1 g1 c1 a0 s201   1.0000000       0.0000000       1.0000000       1.0000000                        
 54:Phi s1 g1 c1 a0 s201   1.0000000       0.4262330E-05   0.9999916       1.0000083                          
 55:Phi s1 g1 c1 a0 s201   1.0000000       0.0000000       1.0000000       1.0000000                          
 56:Phi s1 g1 c1 a0 s201   1.0000000       0.0000000       1.0000000       1.0000000                          
59:Phi s1 g1 c1 a0 s201   1.0000000       0.0000000       1.0000000       1.0000000                          
 65:Phi s1 g1 c1 a0 s201   1.0000000       0.0000000       1.0000000       1.0000000                          
 71:Phi s1 g1 c1 a0 s201   1.0000000       0.0000000       1.0000000       1.0000000                        
 73:Phi s1 g1 c1 a0 s201   1.0000000       0.4715585E-06   0.9999991       1.0000009                         
 78:Phi s1 g1 c1 a0 s201   0.9999999       0.4220431E-04  0.7254796E-301   1.0000000                          
 79:Phi s1 g1 c1 a0 s201   0.9999709       0.0030372       0.5979220E-84   1.0000000                          
89:p s1 g1 s2006 t3       0.5214649       0.0000000       0.5214649       0.5214649                          
 91:p s1 g1 s2006 t5       1.0000000       0.0000000       1.0000000       1.0000000                          
 147:p s1 g1 s2018 t1       1.0000000       0.7081960E-05  0.3379634E-300   1.0000000 
336:pent s1 g1 a0 s2005    0.5686495E-56   0.0000000       0.5686495E-56   0.5686495E-56                    
 338:pent s1 g1 a0 s2005    0.4545653E-63   0.7364965E-60  -0.1443079E-59   0.1443988E-59                      
 342:pent s1 g1 a0 s2006    0.4259980E-11   0.1122624E-08  -0.2196084E-08   0.2204604E-08                      
 346:pent s1 g1 a0 s2007    0.1141971E-20   0.0000000       0.1141971E-20   0.1141971E-20                      
 354:pent s1 g1 a0 s2009    0.1142761E-03   0.0000000       0.1142761E-03   0.1142761E-03                      
 356:pent s1 g1 a0 s2010    0.3884196E-05   0.5456457E-03  0.1026143E-124   1.0000000                          
 358:pent s1 g1 a0 s2010    0.1767426E-07   0.4032951E-05  0.1034506E-201   1.0000000                          
 360:pent s1 g1 a0 s2011    0.7211397E-40   0.1853445E-36  -0.3632031E-36   0.3633474E-36                      
368:pent s1 g1 a0 s2013    0.6005180E-10   0.0000000       0.6005180E-10   0.6005180E-10                         372:pent s1 g1 a0 s2014    0.2564095E-07   0.8951283E-05  0.1770696E-304   1.0000000                          
 373:pent s1 g1 a0 s2014    0.5904186E-12   0.0000000       0.5904186E-12   0.5904186E-12                      
 376:pent s1 g1 a0 s2015    0.3842264E-06   0.1344566E-03  0.5111796E-304   1.0000000                          
 378:pent s1 g1 a0 s2015    0.5558775E-10   0.0000000       0.5558775E-10   0.5558775E-10                      
 380:pent s1 g1 a0 s2016    0.1954978E-07   0.4965941E-05  0.1172356E-223   1.0000000                          
 384:pent s1 g1 a0 s2017    0.1519375E-20   0.6520015E-18  -0.1276403E-17   0.1279442E-17                      
 390:pent s1 g1 a0 s2018    0.6985258E-64   0.1156272E-60  -0.2265594E-60   0.2266991E-60                      

                                                    95% Confidence Interval
 Grp. Str. Ses.    N-hat          Standard Error     Lower           Upper
 ---- ---- ----   --------------  --------------  --------------  --------------
   1     1    2    1468.0061       0.0000000       1468.0061       1468.0061  [/list]

Many thanks!
violetblue
 
Posts: 6
Joined: Tue Nov 06, 2018 11:16 pm

Re: MSORD and parameter identifiability

Postby Bill Kendall » Sun Oct 06, 2019 11:29 pm

Providing partial results is helpful, but it's difficult to tell what's causing the problem without seeing how sparse the data are. What species of sea turtle is this analysis from. You state on one hand that you had time varying nest to skip transition but constant skip to nest, but then in another place you said you had constant nest to skip (which makes more sense given the low level at which sea turtles nest consecutively). One thing that helps with pent, and which is a reasonable hypothesis for sea turtles, is to assume arrivals are distributed as a quadratic function across the season. Also, phi as a linear function of time since arrival (age) is also a reasonable starting point.
Bill Kendall
 
Posts: 76
Joined: Wed Jun 04, 2003 8:58 am

Re: MSORD and parameter identifiability

Postby violetblue » Tue Oct 08, 2019 3:55 am

Thank you so much for your insight on this.

I've tried what you suggested but am still getting many poorly estimated parameters.

With the model fitting so poorly, I'm unsure if I should conclude that it is simply not possible to estimate both time-varying S and time-varying Psi (skip to nest).

I don't think the data is sparse. There are 4315 different encounter histories with recaptures at every occasion (flatback species,15 years of sampling and 5 fortnights per season), apart from one occasion; the 3rd fortnight in the second season had missing data (there were no recaptures at all) so I've now fixed p recapture probability to 0 for this occasion.
I'm able to send data and results privately if that helps!

I had made Psi skip to nest constant, but have changed it now to the more realistic constant nest to skip.

I've also tried what you suggested for pent and phi by;

(i) assuming arrivals are distributed as a quadratic function across the season, which I coded as
pent.session.sq=list(formula = ~ session + (session)^2,fixed=list(index=upent,value=0)), and

(ii) phi as a linear function of time since arrival (age), which I coded as a 4 ageclass TSM (not sure if this is the usual way but it seemed to give me the results I wanted in the ddl).
ordms.ddl$Phi$TimeSinceArrival=ordms.ddl$Phi$Time-ordms.ddl$Phi$Cohort
ordms.ddl$Phi$four_a=as.factor(ifelse(ordms.ddl$Phi$ TimeSinceArrival >3,4,ordms.ddl$Phi$ TimeSinceArrival +1))
Phi.4a.time=list(formula=~-1+four_a:session:stratum,fixed=list(index=uphi,value=0))

Modelling Phi as a function of TimeSinceArrival
Phi.tsa=list(formula=~-1+TimeSinceArrival,fixed=list(index=uphi,value=0))
produced an error message:
' One or more formulae are invalid because the design matrix has all zero rows for the following non-fixed parameters
Phi s1 g1 c1 a0 s2005 t1'
violetblue
 
Posts: 6
Joined: Tue Nov 06, 2018 11:16 pm

Re: MSORD and parameter identifiability

Postby jlaake » Tue Oct 08, 2019 12:03 pm

This post falls into a crack between RMark and this forum. Many of your problems are due to misspecification of the formula. I'm confused why you think poorly estimated parameters are not being counted. Unless you set adjust=FAALSE, RMark assumes the design matrix is full rank and uses the number of dm columns as the parameter count.

Here are some problems.

1) using ~session + session^2 is probably not working. Create a variable sqsession=session^2 in the design data for pent and then use ~session + sqsession in the formula.

2) ~-1 + TimeSinceArrival for Phi gave you an error because you created a regression with no intercept. -1 should only be used with factor variable combinations.

3) See post at top of RMark forum or the workshop notes for a discussion about fixing real parameter values. What you are doing may be fine but it was the original way of fixing parameters and it is kludgy in comparison to using the fix variable in the design data.

Maybe address these issues and re-post if you still have questions that Bill can address with regard to the model.

--jeff
jlaake
 
Posts: 1120
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: MSORD and parameter identifiability

Postby violetblue » Wed Oct 09, 2019 7:25 am

Thank you very much, that has made a huge difference.
I have only a handful of inestimable parameters now so I'll continue working towards the best solution.
violetblue
 
Posts: 6
Joined: Tue Nov 06, 2018 11:16 pm

Re: MSORD and parameter identifiability

Postby jlaake » Wed Oct 09, 2019 8:25 am

If session is a factor variable then you should change to numeric. Send a snippet of your design data for phi and results of str function on that design data to my email account.
jlaake
 
Posts: 1120
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA


Return to analysis & design questions

Who is online

Users browsing this forum: No registered users and 1 guest