www.phidot.org

Hi,
I am fitting multistate models to estimate sex-specific survival. The three defined states (stratum): M, F, and U.
I have successfully fitted a series of models to estimate sex-specific survival with and without a time trend. I would like though to fit a model with a time trend on female survival (F), a constant male survival (M) and either constant or a time trend for the unknown state animals (U).

Process data and design data:

Code: Select all

#Process data
sexMS.processed=process.data(markdataMS,model="Multistrata",begin.time=1989)

#Create the design data
sexMS.ddl=make.design.data(sexMS.processed,parameters=list(Psi=list(pim.type="time")))


#add design data to allow survival probabilities to be M=F<>U, M=U<>F,F=U<>M

Code: Select all

sexMS.ddl$S$SFM="FM"
sexMS.ddl$S$SFM[sexMS.ddl$S$stratum=="U"]="SU"
sexMS.ddl$S$SFU="FU"
sexMS.ddl$S$SFU[sexMS.ddl$S$stratum=="M"]="SM"
sexMS.ddl$S$SMU="MU"
sexMS.ddl$S$SMU[sexMS.ddl$S$stratum=="F"]="SF"

Define survival parameters to run models

Code: Select all

#Time trend for each state (sex)
S.sexTime=list(formula=~Time:stratum)
#Time trend for F and M+U (together as one parameter)
S.MUTimeb=list(formula=~Time:SMU)

I wonder whether I need to modify the design matrix for S before I can define the model parameters? Any tips on how to specify the survival parameter would be much appreciated. I just haven't been able to see any examples.
Many thanks,
Monica

Actually, I figured it out by looking at another post with a similar issue.
In case it is of anyone's interest the trick is to add a dummy variable for each of the groups, so that you can call them separately when defining survival. In my case:

Code: Select all

#define dummy variable for each sex
sexMS.ddl$S$Female=0
sexMS.ddl$S$Female[sexMS.ddl$S$stratum=="F"]=1
sexMS.ddl$S$Male=0
sexMS.ddl$S$Male[sexMS.ddl$S$stratum=="M"]=1
sexMS.ddl$S$Unknown=0
sexMS.ddl$S$Unknown[sexMS.ddl$S$stratum=="U"]=1

and then the models can be run, for example for a Time trend in Female survival but keeping male and unknown animals constant this would be:

Code: Select all

S.TF.Mdot.Udot=list(formula=~Time:Female+Male+Unknown)
p.time=list(formula=~time)
#fix some transitions in years I know there were no unknown animals being sexed as males or females.
Psi.time=list(formula=~-1+time,fixed=list(index=c(110,112,113,114,115,116,117,131,132,133,134,135,136,137,138,139,140,141,142,143),value=0))

model=mark(sexMS.processed,sexMS.ddl,model="Multistrata",model.parameters = list(S=S.TF.Mdot.Udot,p=p.time,Psi=Psi.time))


Hopefully I did not make anyone lose time with this.
Have fun Rmarkcoding.