Different DMs, same model, different results

[TimG]

I know this topic has been discussed elsewhere on this forum, and I have checked that I am not using different link functions, different data/models, or a confounded parameter as a reference level.

I am using the Barker robust design manatee example project that comes with the MARK installation (MEE-13-04-179 MARK EXAMPLE FILE_KENDALL ET AL BARKER-RD 20 YEAR MANATEE). The deviance for the top model in this project = 26248.4710. Except for Survival, the design matrix is set up as an Identity matrix. I attempted to recreate this model but using Intercepts in the design matrix for all parameters, with the 1st year as the reference level (I am not including the DM in this post for the sake of brevity, but others can attempt it or I can provide mine). This resulted in a deviance of 26235.6490. Also, the real parameter estimates were generally close but not identical to those in the Identity model, with greater differences near the end of the time series.

Interestingly, I have tried many different approaches to this same model (reference level in the middle of the time series, some parameters as Identity others as Intercept, simulated annealing) and each give a different model deviance. Is this possibly related to the large number of parameters in the Barker RD model, especially here when nearly everything is fully time-dependent? (I don’t encounter the same issues with a simpler model/dataset e.g. phi(t)p(t) using the dipper data) Are there possibly optimization issues? (although simulated annealing didn't shed much light, and I am still investigating MCMC results)

[cooch]

I know this topic has been discussed elsewhere on this forum, and I have checked that I am not using different link functions, different data/models, or a confounded parameter as a reference level.

I am using the Barker robust design manatee example project that comes with the MARK installation (MEE-13-04-179 MARK EXAMPLE FILE_KENDALL ET AL BARKER-RD 20 YEAR MANATEE). The deviance for the top model in this project = 26248.4710. Except for Survival, the design matrix is set up as an Identity matrix. I attempted to recreate this model but using Intercepts in the design matrix for all parameters, with the 1st year as the reference level (I am not including the DM in this post for the sake of brevity, but others can attempt it or I can provide mine). This resulted in a deviance of 26235.6490. Also, the real parameter estimates were generally close but not identical to those in the Identity model, with greater differences near the end of the time series.

Interestingly, I have tried many different approaches to this same model (reference level in the middle of the time series, some parameters as Identity others as Intercept, simulated annealing) and each give a different model deviance. Is this possibly related to the large number of parameters in the Barker RD model, especially here when nearly everything is fully time-dependent? (I don’t encounter the same issues with a simpler model/dataset e.g. phi(t)p(t) using the dipper data) Are there possibly optimization issues? (although simulated annealing didn't shed much light, and I am still investigating MCMC results)

Quick observations --

1\ Barked RD models are 'twitchy', especially with near/full time-dependence., especially if you have lots of parms estimated near [0,1] (as is the case with the example data set).

2\ the structure for S is an additive age+time TSM structure. Unless you're *really* careful, its easy to get structures with different reference levels out of synch for such models. Estimates from said models are often *close*, but not exactly the same. I'd suggest that your results indicate a problem with the DM construction. Given the size of that example data set, it is very easy to gt things off by even a little a bit;. Perfect candidate for RMark sorts of inefficiencies....

[TimG]

Thanks for the quick response.
1) Yes, I’m getting the impression that this model is twitchy/unstable and thinking that the different DMs are possibly affecting the optimization results? The R’ parameter seems especially troublesome, and given the high survival and moderate recovery rates in this example, I suspect there’s little data informing this parameter.
My concern is that the difference in deviance/AIC between what I believe to be the same models is greater than that between 2 truly different models. Also, even the small differences in S I’m observing can be significant, e.g. in a PVA for a long-lived species like this. So what would be the recommendation for dealing with model selection?

2) I think I have the DM correct, but wouldn’t mind having someone else review it again. What’s the easiest way to share one of this size, sending an Excel file offline?
Yes, right now I’m building up my confidence in RMark because at the very least I was curious to see how it would create the DM for this particular model. For the sake of moving forward, I was considering using the default DMs built by RMark for use in model selection.

[TimG]

Just following up on this in case others encounter the same issue. I built this model in RMark, which uses the 'Intercept' approach for the design matrix as default. Similar to what I had been observing, this gave slightly different deviance and real parameter estimates than the model in the example MARK project, which primarily used the 'Identity' approach.

Maybe not a totally satisfying answer, but I'm going to attribute these differences to the 'twichiness' of this model, as many parameters are near 0,1 and may have sparse data for full time-dependence. I still suspect the different design matrices are resulting in convergence on different maxima. For the purpose of model selection, hopefully AIC differences between models that are actually different will be greater than those between the 'same' model.

I haven't seen an RMark example for the Barker RD model, so some code is below if others are interested.

Code: Select all

#Closed Barker robust design for manatees
library(RMark)
#Read in data. Note, this function only allows 1 set of /*comments*/ in the .inp file
manatee=convert.inp("C:/Program Files (x86)/MARK/Examples/MEE-13-04-179 MARK example file_Kendall et al Barker-RD 20 year manatee.inp", use.comments=TRUE, covariates='days')
dim(manatee)
head(manatee)

#20 primary periods each with 2 secondary occasions
t.int<-c(rep(c(0,1),20)) #Barker RD a bit different since you have data for last interval
t.int

#process data for RMark
manatee.data<-process.data(manatee,model="RDBarkHug",time.interval=t.int)
manatee.ddl<-make.design.data(manatee.data)
manatee.ddl$S #design data for S

#create grouping index to fix all F's=1
F.fix=as.numeric(row.names(manatee.ddl$F))


##################################
# Fully time-dependent model
#Parameter specifications
S.time=list(formula=~time)  #fully time-dependent
r.time=list(formula=~time)
R.time=list(formula=~time)
RPrime.time=list(formula=~time)
aDoublePrime.time=list(formula=~time)
aPrime.time=list(formula=~time)
F.dot=list(formula=~1, fixed=list(index=F.fix,value=1)) #constant; fixed=1
p.session.times.time=list(formula=~session*time, share=TRUE)  #p=c; different for each session(year) and time (secondary period) - note different meaning for 'time' for this parameter compared to others
#c.session=list(formula=~session, share=TRUE)  #p=c; varies with session; shouldn't need if p=c

#run this model
m1=mark(manatee.data,manatee.ddl,model="RDBarkHug", model.parameters=list(S=S.time,r=r.time,R=R.time,RPrime=RPrime.time,aDoublePrime=aDoublePrime.time,
   aPrime=aPrime.time,F=F.dot,p=p.session.times.time))
summary(m1) #review model
m1$results$real
#export Design Matrix
library(MASS)
write.matrix(m1$design.matrix,"C:\\m1dm.txt")