www.phidot.org

Dear RMark User,
I am wondering whether you could give me some suggestion on my models.
For example, I want to build two models using a co-variate ("+Area")

1. Mod1: Movement is differed by starta and time and a function of Area, where survival and p are constant
2. Mod2: movement is differed by starta and function of Area, Survival is also differed by starta and time and

I used following code, I am wondering whether I have been doing Ok or not in terms of building model.

####
ms.pr=process.data(ms, model = "Multistrata")
summary(ms.pr)

ms.ddl=make.design.data(ms.pr)
head(ms.ddl$S)
head(ms.ddl$p)
head(ms.ddl$Psi)
#
ms.ddl$p$Area=NA
#ms.ddl$p$temperature=NA
#ms.ddl$p$flow=NA
ms.ddl$p$Area[ms.ddl$p$stratum=="A"]=467829
ms.ddl$p$Area[ms.ddl$p$stratum=="B"]=299104
ms.ddl$p$Area[ms.ddl$p$stratum=="C"]=236567
ms.ddl$p$Area[ms.ddl$p$stratum=="D"]=450147
ms.ddl$p$Area[ms.ddl$p$stratum=="E"]=104052
ms.ddl$p$Area[ms.ddl$p$stratum=="F"]=10405200
#
ms.ddl$S$Area[ms.ddl$S$stratum=="A"]=467829
ms.ddl$S$Area[ms.ddl$S$stratum=="B"]=299104
ms.ddl$S$Area[ms.ddl$S$stratum=="C"]=236567
ms.ddl$S$Area[ms.ddl$S$stratum=="D"]=450147
ms.ddl$S$Area[ms.ddl$S$stratum=="E"]=104052
ms.ddl$S$Area[ms.ddl$S$stratum=="F"]=10405200
#
ms.ddl$Psi$Area[ms.ddl$Psi$stratum=="A"]=467829
ms.ddl$Psi$Area[ms.ddl$Psi$stratum=="B"]=299104
ms.ddl$Psi$Area[ms.ddl$Psi$stratum=="C"]=236567
ms.ddl$Psi$Area[ms.ddl$Psi$stratum=="D"]=450147
ms.ddl$Psi$Area[ms.ddl$Psi$stratum=="E"]=104052
ms.ddl$Psi$Area[ms.ddl$Psi$stratum=="F"]=10405200
#
Sdot=list(formula=~1)
S.stratum=list(formula=~stratum)
S.stratumxtime.Area=list(formula=~stratum*time+Area)

pdot=list(formula=~1)
p.stratumxtime=list(formula=~stratum*time)
p.stratumxtime.Area=list(formula=~stratum*time+Area)

Psi.s.dot=list(formula=~1)
Psi.s.time.area=list(formula=~-1+stratum:tostratum*time+Area)
Psi.s.area=list(formula=~-1+stratum:tostratum+Area)

# I want to build two models,
#Mod1 is: movement is differed by starta and time and a function of Area, where survival and p are consatnt
#Mod2: movement is differed by starta and function of Area, Survival is also differed by starta and time and function of Area

Mod1<-mark(ms.pr,ms.ddl,model.parameters=list(Psi=Psi.s.time.area,S=Sdot,p=pdot))
Mod2<-mark(ms.pr,ms.ddl,model.parameters=list(Psi=Psi.s.area,S=S.stratumxtime.Area,p=pdot))

Are these procedures correct?

Thanks for your help.

It looks fine to me except for the following suggestions.

1) Change your definition of Area by either subtracting off the mean Area or the minimum Area such that the estimated intercept makes sense. It would be the value for the mean Area or the minimum Area if you did those changes above. Without these changes, the intercept is for a strata with Area=0 which will not make as much sense.

2) Also, MARK does scale covariates but I suggest once you do the above that you divide each value of Area by some constant so the covariate is in smaller units. For example, if it is currently defined as hectares, you could divide by 100 to put in units of 100s of hectares.

3) For formula Psi.s.time.area=list(formula=~-1+stratum:tostratum*time+Area) this one is a bit tricky. As written it will have one too many parameters because stratum:tostratum*time is equivalent to -1+time+ stratum:tostratum + time:stratum:tostratum. Because time is a factor variable -1+time has the same number of parameters as time because it redefines the DM from an intercept + 2 parameters to 3 separate parameters. You can see that by running the following example:

Code: Select all

library(RMark)
data(mstrata)
dp=process.data(mstrata,model="Multistrata")
ddl=make.design.data(dp)
mm=model.matrix(~-1+stratum:tostratum*time,ddl$Psi)
head(mm)
ncol(mm)
mm=model.matrix(~stratum:tostratum*time,ddl$Psi)
head(mm)
ncol(mm)
mm=model.matrix(~-1+stratum:tostratum*time,ddl$Psi)


Both have 30 columns and you only want 29. Actually it will be fewer because RMark removes the subtract.stratum transitions (eg A to A). So when you have an additive factor variable like time, don't use -1 to remove the intercept. Instead use remove.intercept=TRUE as in

Code: Select all

Psi.s.time.area=list(formula=~stratum:tostratum*time+Area,remove.intercept=TRUE)


The argument remove.intercept=TRUE is an RMark argument and will not work with model.matrix but it is equivalent to:

Code: Select all

mm=model.matrix(~stratum:tostratum*time,ddl$Psi)
# remove first column which is the intercept
mm=mm[,-1]
head(mm)
ncol(mm)


4) Unless you have lots of data you may not be able to estimate the time interactions for S and Psi so you may want to consider additive time models if that makes sense for your situation.

regards --jeff