Model averaging on nest survival data

[fltt]

Hello forum,

I am running my nest survival data on RMark and have encountered some troubles with the model.average function. My data consists of continuous and categorical (2 levels) variables, where categorical variables are "grouped" in the model. I followed the Mallard example by Rotella to do this. I am able to run the models, and obtained beta estimates for each variable, as well as DSRs without troubles. The problem arises when I try to conduct model averaging, as I get the following error message:

"Error in model.average.marklist(YHPA.results.Final, parameter = "S", vcv = TRUE) :
Cannot model average models with different structures".

Any idea what this means? My code looks as follows:

Code: Select all

run.YHPA.Poaching.Final <- function()
{ 
  ##modelo con S constante- null model
  Null <- mark(data = YHPA_Poaching, nocc = 125,model = "Nest", model.parameters = list(S=list(formula=~1)))
 
  Temporal <- mark(YHPA_Poaching,nocc=125,model="Nest",model.parameters=list(S=list(formula=~Time + I(Time^2) + Year)), groups = c("Year"))
 
  Anthropogenic <- mark(YHPA_Poaching, nocc = 125, model = "Nest", model.parameters = list(S=list(formula=~SteOpnnss + Site + Settlement + Road + PrevPchng)), groups = c("SteOpnnss","Site","PrevPchng"))
 
  Global <- mark(YHPA_Poaching, nocc = 125, model = "Nest", model.parameters = list(S=list(formula=~Time + I(Time^2) + Year + PrevPchng + SteOpnnss + Site + Settlement + Road)), groups = c("SteOpnnss","Site","PrevPchng","Year"))
 
  return(collect.models())
}
YHPA.results.Final <- run.YHPA.Poaching.Final()
YHPA.results.Final


Avg.Estimates <- model.average(YHPA.results.Final, parameter = "S", vcv = TRUE)

# Model.average produces the error below
# Error in model.average.marklist(YHPA.results.Final, parameter = "S", vcv = TRUE) :
# Cannot model average models with different structures

After reading this post (viewtopic.php?f=21&t=3258&hilit=Cannot+model+average+models+with+different+structures), I thought the issue came from the structure of my code, so I changed it as to use process data and make design data.Now I encounter the following error message:

"Using default formula for S
Error in create.model.list("Nest") :
No model specifications found. Use case sensitive parameter.description notation (e.g., Phi.time)"

In this case, the function with the DSR models did not run. My code, after modifying it to use process.data and make.design.data looks as follows:

Code: Select all

# what are the parameters of the model
# There is only one parameter, the daily survival probality (S)
setup.parameters("Nest", check=TRUE)


# 1. Process the data.
# The nocc variable is the data at which hatching occurs
YHPAData.proc <- process.data(YHPA_Poaching, model="Nest", nocc=125,
            groups = c("SteOpnnss","Site","PrevPchng","Year"))
YHPAData.proc

# 2. Examine and/or modify the ddl. Here you could standardize covariates
# Modify the ddl for all TIME VARYING covariates!!!
YHPAData.ddl <- make.design.data(YHPAData.proc)
YHPAData.ddl$S$Time2 <- (YHPAData.ddl$S$Time)^2 #quadratic time variable
YHPAData.ddl

# 3. Develop nest survival function
run.YHPA.Poaching.Final <- function()
{
 
  Null <- list(formula=~1)
 
  Temporal <- list(formula=~Time + I(Time^2) + Year)
 
  Anthropogenic <-list(formula=~SteOpnnss + Site + Settlement + Road + PrevPchng)
 
  Global <- list(formula=~Time + I(Time^2) + Year + PrevPchng + SteOpnnss + Site + Settlement + Road)
 
  cml <- create.model.list("Nest")
 
  return(mark.wrapper(cml,data=YHPAData.proc,ddl=YHPAData.ddl)) 
 
}
YHPA_Poaching_Results <- run.YHPA.Poaching.Final()

# The following error message appears
#Using default formula for S
#Error in create.model.list("Nest") :
#No model specifications found. Use case sensitive parameter.description notation (e.g., Phi.time)



Any comments/suggestions on how to obtain model averaged DSR estimates for the first sample code OR a potential solution to the second sample code would be of great help. Thanks in advance!

[jlaake]

The mallard example isn't a particularly good example for how I'd suggest using RMark. It is inefficient in that each call to mark function is processing data and making design data. It is better to call process.data and then make.design data and then using those two results in the call to mark or mark.wrapper. in your calls to mark you have used different groups arguments which changes the underlying structure of the pims. That prevents you from using model.average because parameter x is not the same across all models Use a single call to process.data with your largest set of variables defining groups so the pim structure will be the same across all models.

[Rotella]

I agree that the mallard example in RMark should be updated. Although it certainly shows how to run nest survival models in RMark, it is not doing so as efficiently as it could. It will be awhile before I get to that, but it's on my to-do list.

[jlaake]

A lot of my documentation and examples need updating as well. Didn't mean my comment as a criticism of your example. You were one of the first to provide examples for RMark which has been very useful for a lot of folks. It does have this one side effect discovered here.

[fltt]

Thanks for your response!

It is better to call process.data and then make.design data and then using those two results in the call to mark or mark.wrapper.

Yes, so this is what I was attempting to do with the second line of code I provided. But I encountered a problem with the create.model.list function, as it is not recognizing the "Nest" model. Have you encountered this problem before?

[jlaake]

Sorry, I was on my phone and I stopped reading once I saw what the problem was. Have you read the workshop notes or the help for create.model.list? Problem wasn't Nest model, the problem was that you didn't specify the parameter models with the correct naming convention. It is how formulas are recognized. With the Nest model you only have survival but for others you have multiple parameters so a naming convention was needed to know how to put the models together from the various parameter specifications. You have to use parameter.something where parameter is the name of the parameter (S in this case) followed by a dot and then a numbering scheme or words. I mean this in a helpful way - looking at help and reading the documentation will typically be a faster way to solve a problem. I have written 100s of pages of documentation and while it may not be perfect and needs updating it will generally solve most problems.