www.phidot.org

Hi,

I'm currently using RMark to evaluate seasonal survival for a species of interest. My data is currently organized in a monthly fashion and I need to collapse the months into the following seasons: 1 April–30 June (spring), 1 July–30 September (summer), 1 October–31 December (fall), 1 January–31 March (winter). However, I cannot find an easy way to collapse the monthly data into this structure. Any suggestions would be much appreciated.

Thanks!

You need to explain in more detail what you want. Are you asking how to collapse your capture history into seasonal occasions from monthly occasions? I would advise against doing that but that is something you can do with R code or with your database. Or are you asking how to make survival estimates constant across season or combine monthly survival estimates to seasonal estimates?

--jeff

Sorry for the confusion. I have 3 objectives for a project I'm working on as they each build on each other.

1) Estimate monthly and annual survival for a species of interest
-My data is currently structured with monthly occasions (n=32 occasions). I evaluated the following 2 models: constant survival and time-varying survival

#Define range of models for S
S.dot=list(formula=~1) #Constant Survival
S.time=list(formula=~time) #Time-varying Survival
#
#Create model list and run assortment of models
model.list=create.model.list("Known")
results=mark.wrapper(model.list,data=animal.processed,ddl=animal.ddl,invisible=FALSE,threads=2)

The model output was as follows:
model npar AICc DeltaAICc weight Deviance
1 S(~1) 1 240.3178 0.00000 9.999763e-01 4.404026e+01
2 S(~time) 32 261.6186 21.30077 2.369121e-05 2.271242e-06

Survival was constant; therefore, I output the monthly survival estimate: 0.96 (SE = 0.01). To estimate annual survival, I used the monthly survival (0.96 ^ 12) = 0.59 (annual survival).

2) Estimate seasonal survival during the following 4 seasons (1 April–30 June (spring), 1 July–30 September (summer), 1 October–31 December (fall), 1 January–31 March (winter)).

- Since survival was constant, I guess I could just calculate seasonal survival based on the product of the monthly survival estimate above. Would that be correct? If so, is there an easy way to obtain standard errors, lower and upper confidence intervals for the derived estimates (e.g., seasonal and annual)?

3) Evaluate the effects of reproductive status (reproductively active vs. inactive) on seasonal survival
-This requires the addition of a covariate for reproductive status. I've read about adding a covariate into the R code but I'm struggling with how to deal with animals monitored across multiple seasons/years. For example, 1 year an animal was reproductively active and the following year it was not. Therefore, it's not as simple as just adding a covariate that is not time-varying.

Thanks for the assistance!

alittle wrote:2) Estimate seasonal survival during the following 4 seasons (1 April–30 June (spring), 1 July–30 September (summer), 1 October–31 December (fall), 1 January–31 March (winter)).

- Since survival was constant, I guess I could just calculate seasonal survival based on the product of the monthly survival estimate above. Would that be correct? If so, is there an easy way to obtain standard errors, lower and upper confidence intervals for the derived estimates (e.g., seasonal and annual)?

I already answered this for you earlier - Delta method - http://www.phidot.org/software/mark/doc ... /app_2.pdf

Simply read the first paragraph, which describes virtually the exact same thing you're after.

cooch wrote:

alittle wrote:2) Estimate seasonal survival during the following 4 seasons (1 April–30 June (spring), 1 July–30 September (summer), 1 October–31 December (fall), 1 January–31 March (winter)).

- Since survival was constant, I guess I could just calculate seasonal survival based on the product of the monthly survival estimate above. Would that be correct? If so, is there an easy way to obtain standard errors, lower and upper confidence intervals for the derived estimates (e.g., seasonal and annual)?

I already answered this for you earlier - Delta method - http://www.phidot.org/software/mark/doc ... /app_2.pdf

Simply read the first paragraph, which describes virtually the exact same thing you're after.

Thanks for the assistance! I forgot to add that I was wondering if you could easily do the calculation in RMark coding--seamless process.

alittle wrote:

cooch wrote:

alittle wrote:2) Estimate seasonal survival during the following 4 seasons (1 April–30 June (spring), 1 July–30 September (summer), 1 October–31 December (fall), 1 January–31 March (winter)).

- Since survival was constant, I guess I could just calculate seasonal survival based on the product of the monthly survival estimate above. Would that be correct? If so, is there an easy way to obtain standard errors, lower and upper confidence intervals for the derived estimates (e.g., seasonal and annual)?

I already answered this for you earlier - Delta method - http://www.phidot.org/software/mark/doc ... /app_2.pdf

Simply read the first paragraph, which describes virtually the exact same thing you're after.

Thanks for the assistance! I forgot to add that I was wondering if you could easily do the calculation in RMark coding--seamless process.

I just found what I was looking for in regards to the Delta method....

This example is excluded from testing to reduce package check time
#
# The following are examples only to demonstrate selecting different
# model sets for adjusting chat and showing model selection table.
# It is not a realistic analysis.
#
data(dipper)
mod=mark(dipper,model.parameters=list(Phi=list(formula=~time)))
rr=get.real(mod,"Phi",se=TRUE,vcv=TRUE)
deltamethod.special("prod",rr$estimates$estimate[1:6],rr$vcv.real)
deltamethod.special("cumprod",rr$estimates$estimate[1:6],rr$vcv.real,ses=FALSE)
deltamethod.special("sum",rr$estimates$estimate[1:6],rr$vcv.real)
deltamethod.special("cumsum",rr$estimates$estimate[1:6],rr$vcv.real,ses=FALSE)

Thanks!

Thanks for the assistance! I forgot to add that I was wondering if you could easily do the calculation in RMark coding--seamless process.

Sure, see deltamethod.special:

Code: Select all

data(dipper)
  mod=mark(dipper,model.parameters=list(Phi=list(formula=~time)))
  rr=get.real(mod,"Phi",se=TRUE,vcv=TRUE)
  deltamethod.special("prod",rr$estimates$estimate[1:6],rr$vcv.real)
  deltamethod.special("cumprod",rr$estimates$estimate[1:6],rr$vcv.real,ses=FALSE)
  deltamethod.special("sum",rr$estimates$estimate[1:6],rr$vcv.real)
  deltamethod.special("cumsum",rr$estimates$estimate[1:6],rr$vcv.real,ses=FALSE)


But, I'd suggest doing one by hand first, so you understand what the code is doing.

deltamethod.special is overkill for your problem but it should work fine. I do agree with Evan that you should try this simple problem by hand so you understand how it is computed.

You should probably consider fitting a seasonal model. Right now you have compared models with 1 and 32 parameters. Something intermediate like a seasonal model might be supported.

With regard to your question about repro covariate, see the section on time varying individual covariates in the workshop notes and the RMark Appendix in Cooch and White book.

--jeff

Thanks for the assistance!

I have 2 issues associated to the suggestions below:

1) If I'm wanting to obtain the SE estimate for the annual survival for my species of interest, would you use the brackets to define [1:12} for the first 12 time occasions (i.e., months)? My interest is defining an annual survival estimate during the 2.5 year study.

#Delta Method to Derive SE
#Listed below models
deltamethod.special("prod",animal.results$estimates$estimate[1:12])


#Main Code

#Create Data Frame for Known Fate Analysis
mark.data<-data.frame(ch=capt.hist,Id=Id)

run.mark.data=function()
{
#Process the data
animal.proc = process.data(mark.data,model="Known")

# Create default design data
animal.ddl=make.design.data(animal.proc)
#
#Define range of models for S
S.dot=list(formula=~1) #Constant Survival
S.time=list(formula=~time) #Time-varying Survival
#
#Create model list and run assortment of models
model.list=create.model.list("Known")
animal.results=mark.wrapper(model.list,data=animal.proc,ddl=animal.ddl,invisible=FALSE,threads=2)

# Return model table and list of models
}
animal.results=run.mark.data()
animal.results

#Derived Parameters
animal.results$S.dot$results$real$estimate^12 #Annual Survival
animal.results$S.dot$results$derived #Survival across 32 months of the study
animal.results$S.dot$results #Monthly Survival

#Delta Method to Derive SE
deltamethod.special("prod",animal.results$estimates$estimate[1:12])


2) How would you about grouping the months into seasons to create a seasonal model? I recognize that you need to change the PIMS but I don't know how to code this in R. I've spent most of today trying to solve this problem with little luck in figuring out how to change the PIMs to create a seasonal model.

Thanks for the information regarding time varying individual covariates as well as the other information!

The answer depends on the model. With a constant survival, the monthly survivals are all the same, so it is S^12. Using deltamethod.special with "prod" and estimates 1 to 12 (all being the same) is equivalent to S^12. But you can't just give it the estimates to get the se. You need to specify the var-cov values as well. Please try the delta method by hand so you can see what it is doing.

To create seasonal values you need to add a variable to the design data that associates each interval with a season. Examples like this are in the Appendix and Workshop notes. For example see the creation of the flood covariate in ?dipper. With known fate data, the design data are simple. Let's say that each season is 3 months and your first interval is for Jan. Then something like

ddl$S$Season=factor(c(rep(c(rep("Winter",3),rep("Spring",3),rep("Summer",3),rep("Fall",3)),2),c(rep("Winter",3),rep("Spring",3),rep("Summer",1))))

where ddl is the name of your design data list. You'll have to play around with it to get the order correct for what you want.

Please spend some more time reading the documentation and examples. You'll find the answers to most of your questions.