RMark MS model: non-sense estimates of Psi

[rkhamede]

Hello,

I have recently moved from Mark to RMark and I am having troubles for estimating Psi in a multi state model.
My data set consist in capture histories of Tasmanian devils from two different sites regarded as healthy sub adults (A), healthy adults (B) and disease adults (C).
So the data file looks like this:
00000000AAB0000C 0 1;
00000000AAB0BB0C 0 1;
00000AA0B000C000 1 0;
00000AA00B000000 1 0;

Where the first two capture histories are from site 1 and the last two from site 2. There are about 250 individuals at each site.
I have specified extra parameters in the design data (i.e. disease prevalence, time from disease emergence, seasons) and added constrains in some transitions (i.e. BtoA and CtoA, transit from adults to sub adults is impossible, and CtoB, the disease is 100% lethal and there is no recovery).
I also fixed recapture rate for stratum C = 1 as they can’t move to any other state, and transitions from C to C = 1 (absorbing state).
Then I separated specific transitions:
• The toB stratum is coded for AtoB only as CtoB has been eliminated and BtoB is estimated by subtraction
• The toC stratum codes for AtoC, BtoC and CtoC and as I need to constrain each rate to different variables I separated them in three different strata

When I run a preliminary models to check parameter estimates I can get outputs for recapture and survival rates but the transition rate parameter (Psi) estimates make no sense (massive standard errors, i.e. 2598.3300 and exact estimate values for all transitions).
I have discussed this with a colleague that have run and published these models for Tasmanian devils with a very similar data set but we haven’t been able to identify the problem. Even when I run the simplest model with no time variation or disease parameters and with only one transition rate from (A to B) the estimates of Psi make no sense. I thought it might be a bug in the program, as I have tripled checked the data set and the code that I have used is the same of my colleague.

I appreciate any help or thoughts you could share. I can send you the TinnR code and the data set if you think is useful.

Thanks in advance
Rodrigo

[jlaake]

I'll contact you offlist and you can send me data and R code. It sounds like you have manipulated the design data quite a bit and that is probably the source of the problem. Have you tried running a simple model prior to removing design data? If there are no transitions then MARK will set estimate them as 0.

--jeff

[rkhamede]

Thanks a lot jeff,
yep, I tried to run a quick model with transition toB only, prior to fixing all the other 'tostratrum' and estimates of Psi seem reliable. So might be something with the design data code, but the weird thing is that the same code worked for other data set.
cheers
Rodrigo

[jlaake]

I believe your Psi estimates were nonsensical due to the contraints and formula you specified. But first upgrade to 1.9.7 of RMark. Not certain if that is an issue in this case but that should be your first reaction if you are having a problem as it may have been something I fixed. Also, you may want to look at the function TransitionMatrix as it will compute the subtracted transition values (and std errors) as well as the others. Below is the script that I used and the results look quite sensical although I know little of your data and problem. Note that even though p for stratum C is high 0.87 it is not 1. Where I changed things in your script I put Jeff: and then my comment about the change. It looks like you are making good use of RMark. From your comments in the script you had partially indentified two of the three problems.

regards --jeff

library(RMark)

#Multi-state analysis for comparing transition rates between two sites
#monthly intervals (May-Aug, Aug/Nov, Nov/May. Feb trips have been deleted from data)

#data is grouped by site only (no sex effects)
MSbysiteFENTWPP<-convert.inp("MS_FENT_WPP.inp", group.df=data.frame(site=c("FENT", "WPP")))

#Process the data - time intervals are in months. So age.unit has to be in months as well =1/12.
MSbysiteFENTWPP.processed<-process.data(MSbysiteFENTWPP, model="Multistrata", groups=("site"),
time.intervals=c(6,3,3,6,3,3,6,3,3,6,6,6,3,3,6), initial.ages=c(1,1), begin.time=1, age.unit=1/12)

MSbysiteFENTWPP.processed

#Create the design data making sure to change the strata that MARK calculates by default (these are A->A, B->B, C->C). We have to change these defaults so that we can apply constraints to the strata we need (see below)
# Jeff: The default would have worked fine
MSbysiteFENTWPP.ddl<-make.design.data(MSbysiteFENTWPP.processed, parameters=list(Psi=list(subtract.stratum=c("A", "B", "C"))))
summary(MSbysiteFENTWPP.ddl$Psi)

#Some transitions in our model aren't possible (i.e H.Ads to Subads ; D.Ads to Subadults and D.Ads to H.Ads). The easiest way to fix these transitions to zero is to remove the design data. The following code does that
MSbysiteFENTWPP.ddl$Psi=MSbysiteFENTWPP.ddl$Psi[!(MSbysiteFENTWPP.ddl$Psi$stratum=="B"&MSbysiteFENTWPP.ddl$Psi$tostratum=="A"), ]
MSbysiteFENTWPP.ddl$Psi=MSbysiteFENTWPP.ddl$Psi[!(MSbysiteFENTWPP.ddl$Psi$stratum=="C"&MSbysiteFENTWPP.ddl$Psi$tostratum=="A"), ]
MSbysiteFENTWPP.ddl$Psi=MSbysiteFENTWPP.ddl$Psi[!(MSbysiteFENTWPP.ddl$Psi$stratum=="C"&MSbysiteFENTWPP.ddl$Psi$tostratum=="B"), ]
summary(MSbysiteFENTWPP.ddl$Psi)

# Jeff: absorbing state only implies can't leave stratum and not that p(stratum)=1
# Jeff: Also note that these need to be included in the model specification to be used. They were not
# Jeff: in the model specification you were running.
#Also we need to fix the recapture rates of stratum C to equal 1 at both sites. Because diseased individuals can't recover (that is they can't move to any other state) in this model they can't be recaptured in any other stratum and will always be recaptured in stratum C so we need to fix their recapture rate for this stratum to 1.
#pDAds.fixes<-as.numeric(row.names (MSbysiteFENTWPP.ddl$p[MSbysiteFENTWPP.ddl$p$stratum=="C",]))
#pDAds.values=rep(1,length(pDAds.fixes))

# Jeff: no need to do this because if C is done by subtraction and all others set to 0 then it will be 1.
#And we also have to fix transitions C->C to equal 1 - they can't move away from state C - its known as an absorbing state
#PsiCCfixes<-as.numeric(row.names(MSbysiteFENTWPP.ddl$Psi[MSbysiteFENTWPP.ddl$Psi$stratum=="C"
# &MSbysiteFENTWPP.ddl$Psi$tostratum=="C",] ))
#PsiCCvalues<-rep(1,length(PsiCCfixes))

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>CREATE AND ADD EXTRA DESIGN DATA

#For RECAPTURE rates - we will need dummy variables to separate out FENT and WPP within 'site'
MSbysiteFENTWPP.ddl$p$FENT=0
MSbysiteFENTWPP.ddl$p$FENT[MSbysiteFENTWPP.ddl$p$site=="FENT"]=1
MSbysiteFENTWPP.ddl$p$WPP=0
MSbysiteFENTWPP.ddl$p$WPP[MSbysiteFENTWPP.ddl$p$site=="WPP"]=1

#For SURVIVAL - we will need dummy variables to separate out FENT and WPP within 'site'
MSbysiteFENTWPP.ddl$S$FENT=0
MSbysiteFENTWPP.ddl$S$FENT[MSbysiteFENTWPP.ddl$S$site=="FENT"]=1
MSbysiteFENTWPP.ddl$S$WPP=0
MSbysiteFENTWPP.ddl$S$WPP[MSbysiteFENTWPP.ddl$S$site=="WPP"]=1
# and a dummy variable that allows us to model strata AB together versus C
MSbysiteFENTWPP.ddl$S$AB=0
MSbysiteFENTWPP.ddl$S$AB[MSbysiteFENTWPP.ddl$S$stratum%in%c("A","B")]=1


#For TRANSITIONS - we will also need dummy variables to separate out FENT and WPP within 'site'
MSbysiteFENTWPP.ddl$Psi$FENT=0
MSbysiteFENTWPP.ddl$Psi$FENT[MSbysiteFENTWPP.ddl$Psi$site=="FENT"]=1
MSbysiteFENTWPP.ddl$Psi$WPP=0
MSbysiteFENTWPP.ddl$Psi$WPP[MSbysiteFENTWPP.ddl$Psi$site=="WPP"]=1

#and we need dummy variables to separate out each of the specific transitions.
#NOTE: the "toB" column in the design data already codes for subadult to adult transitions only A->B so this is good. This is because we eliminated all the other toB data (ie C->B) and B->B is estimated by subtraction.
#HOWEVER the "toC" column codes for A->C ; B->C and C->C transitions. Not good when we want to constrain each rate to different variables. So let's separate them.
MSbysiteFENTWPP.ddl$Psi$AtoC=0
MSbysiteFENTWPP.ddl$Psi$AtoC[MSbysiteFENTWPP.ddl$Psi$stratum=="A"&MSbysiteFENTWPP.ddl$Psi$tostratum=="C"]=1

MSbysiteFENTWPP.ddl$Psi$BtoC=0
MSbysiteFENTWPP.ddl$Psi$BtoC[MSbysiteFENTWPP.ddl$Psi$stratum=="B"&MSbysiteFENTWPP.ddl$Psi$tostratum=="C"]=1


#And now we need to add Disease "Prevalence" fields to the Transition rate design data
##Prevalence is the 2-year rolling mean of #diseased/#totalcaught per trip (prev t1+ prev t2)/2) to smooth.

df.prev<-data.frame(group=c(rep("FENT",15), rep("WPP", 15)),
time=rep(c(1, 7, 10, 13, 19, 22, 25, 31, 34, 37, 43, 49, 55, 58, 61),2),
prev=c(0.014, 0.028, 0.105, 0.222, 0.109, 0.214, 0.184, 0.129, 0.368, 0, 0, 0, 0, 0, 0, # for FENT
0, 0, 0, 0, 0.014, 0.035, 0.027, 0.025, 0.035, 0.054, 0.057, 0.064, 0.104, 0.136, 0.121)) # for WPP

#becuase prevalence is different at the two sites we need to merge the prevalence data by group into our design table.
MSbysiteFENTWPP.ddl<-merge.occasion.data(MSbysiteFENTWPP.processed, MSbysiteFENTWPP.ddl,"Psi", df.prev, bygroup=T)

# I also use a "time since disease emergence" variable to look at variation in transition rates as an increasing trend from disease outbreak #the following codes for this field...
df.tsince<-data.frame(group=c(rep("FENT",15), rep("WPP", 15)),
time=rep(c(1, 7, 10, 13, 19, 22, 25, 31, 34, 37, 43, 49, 55, 58, 61),2), #SHELLY: this code include the first trip but not the last one (15 intervals). That was the foramt you had for FNP-FOR
tsince=c(1,2,3,4,5,6,7,8,9,0,0,0,0,0,0,0,0,0,0,1,2,3,4,5,6,7,8,9,10,11)) #SHELLY:last 6 trips not done in FENT - first 4 trips not done in WPP, so the first 15 numbers correspond to FENT, the next 15 to WPP.

MSbysiteFENTWPP.ddl<-merge.occasion.data(MSbysiteFENTWPP.processed, MSbysiteFENTWPP.ddl,"Psi", df.tsince, bygroup=T)

#We know from the single site analysis that repcature rates at least at FNP might be best explained by a seasonal effect between different trapping Intervals. So we need to Add a "season" field to recapture rates
#"Season" allows rates to vary between intervals that are entirely within the mating season (April-June, "a"), entirely outside mating season (June to April, "b"); (June to Nov, "d") and intervals spanning an entire year ("c").
#note that the time here has to include trips 1 through 67 (that is all 16 occasions) - the 1st trip isn't actually measurable in terms of recapture rate though
# SHELLY: Seasons differ in my FENT-WPP data so here "a" is Nov to May (6 months), "b" is May to Aug (3 months), and "c" Aug to Nov (3 months)
df.season<-data.frame(time=c(1, 7, 10, 13, 19, 22, 25, 31, 34, 37, 43, 49, 55, 58, 61, 67),
season=c("c", "a", "b", "c", "a", "b", "c", "a", "b", "c", "a", "c", "a", "b", "c", "a"))

MSbysiteFENTWPP.ddl<-merge.occasion.data(MSbysiteFENTWPP.processed, MSbysiteFENTWPP.ddl,"p", df.season)

#Add season to Survival rates also
MSbysiteFENTWPP.ddl<-merge.occasion.data(MSbysiteFENTWPP.processed, MSbysiteFENTWPP.ddl,"S", df.season)



# Note: this model has no time effect, takes only a minute to run
p.prelim<-list(formula=~site*stratum)
S.prelim<-list(formula=~site*stratum)
# Jeff: Note that I changed Psi by getting rid of CtoC and removed the intercept because you only
# Jeff: have 3 transitions and you specified 3 effects, so by having an intercept you had one too many parameters. That was probably part of the reason for the nonsense.
Psi.prelim<-list(formula=~site:toB +site:AtoC +site:BtoC,remove.intercept=TRUE)#SHELLY: this code spits an error, if I remove the last bit +site:CtoC it works, otherwise it tells me that CtoC variable is not specified in the data( cause we fixed that transition, I guess)
prelim.model<-mark(MSbysiteFENTWPP.processed,MSbysiteFENTWPP.ddl,model.parameters=list(S=S.prelim,p=p.prelim,Psi=Psi.prelim)
,invisible=F)

[rkhamede]

Hi Jeff,
Thanks a lot for your help.
I did try your code but when I run the model it spits an error saying the formulae has zero rows for a whole bunch of non fixed parameters.
Not sure why...the design matrix looks fine to me,
this is the message that R spits:

"Error in make.mark.model(data.proc, title = title, covariates = covariates, :
One or more formulae are invalid because the design matrix has all zero rows for the following non-fixed parameters
Psi sA toC gFENT c43 a2.25 t58,Psi sA toC gFENT c43 a2.5 t61,Psi sA toC gFENT c49 a1 t49,Psi sA toC gFENT c49 a1.5 t55,Psi sA toC gFENT c49 a1.75 t58,Psi sA toC gFENT c49 a2 t61,Psi sA toC gFENT c55 a1 t55,Psi sA toC gFENT c55 a1.25 t58,Psi sA toC gFENT c55 a1.5 t61,Psi sA toC gFENT c58 a1 t58,Psi sA toC gFENT c58 a1.25 t61,Psi sA toC gFENT c61 a1 t61,Psi sB toA gFENT c1 a1 t1,Psi sB toA gFENT c1 a1.5 t7,Psi sB toA gFENT c1 a1.75 t10,Psi sB toA gFENT c1 a2 t13,Psi sB toA gFENT c1 a2.5 t19,Psi sB toA gFENT c1 a2.75 t22,Psi sB toA gFENT c1 a3 t25,Psi sB toA gFENT c1 a3.5 t31,Psi sB toA gFENT c1 a3.75 t34,Psi sB toA gFENT c1 a4 t37,Psi sB toA gFENT c1 a4.5 t43,Psi sB toA gFENT c1 a5 t49,Psi sB toA gFENT c1 a5.5 t55,Psi sB toA gFENT c1 a5.75 t58,Psi sB toA gFENT c1 a6 t61,Psi sB toA gFENT c7 a1 t7,Psi sB toA gFENT c7 a1.25 t10,Psi sB toA gFENT c7 a1.5 t13,Psi sB toA"

thanks again for your thoughts
cheers
Rodrigo

[rkhamede]

Hey Jeff,
I did upgrade RMark to 1.9.7 and it solved the problem.
I run it in the old (2.9.0) version of R and my colleague also run it in R version 2.11.1 and it worked as well
best
Rodrigo

[rkhamede]

Hi Jeff,
For some reason that I haven't been able to identify I keep having non sense estimates of Psi. I have adapted the code with your corrections and run the model for S and p successfully but after running the longer global model with time variation for Psi I get non sense values.
I have run the code with and without the link=“logit” function but both yield non sense outputs. I have tried a dozen times, I have upgraded RMark and R but still can’t get Psi values.

I can send you the code I am running (with the modifications you suggested) if you want.
I appreciate your time and help
Best
Rodrigo