coding a time effect on Survival for just one age class

posts related to the RMark library, which may not be of general interest to users of 'classic' MARK

coding a time effect on Survival for just one age class

Postby splash1199 » Mon Feb 04, 2019 4:49 pm

I am using RMark to run live-dead models. I have set up my data with the following code:

Code: Select all
LD.proc = process.data(data = df,
   model       = "Burnham",
   groups      = c("age", "sex"),
   age.var     = 1,
   initial.age = c(0, 1, 2))
#for age, HY = 0, SY = 1, ASY = 2
#for sex, male = 0, female = 1

LD.ddl=make.design.data(LD.proc)

LD.ddl=add.design.data(data = LD.proc,
   ddl = LD.ddl,
   parameter="S",
   type = "age",
   bins = c(0.5, 1.5, 2.5, 20),
   right = FALSE,
   name = "age",
   replace = TRUE)


My data has two covariates: sex and age (young-of-the-year, 1 year olds, and +1 birds).
So far, I have run the following models (using sin link function) for S:

Code: Select all
run.final = function()
{
    r.dot = list(formula = ~ 1,                    link="sin")
    p.dot = list(formula = ~ 1,                    link="sin")
    F.dot = list(formula = ~ 1,                    link="sin")

    S.a3xsxt = list(formula = ~ -1 + age:sex:time, link="sin")
    S.a3xs   = list(formula = ~ -1 + age:sex,      link="sin")
    S.a3xt   = list(formula = ~ -1 + age:time,     link="sin")
    S.sxt    = list(formula = ~ -1 + sex:time,     link="sin")
    S.a3     = list(formula = ~ -1 + age,          link="sin")
    S.s      = list(formula = ~ -1 + sex,          link="sin")
    S.t      = list(formula = ~ -1 + time,         link="sin")
    S.dot    = list(formula = ~ 1,                 link="sin")

    final.model.list = create.model.list("Burnham")
    final.results = mark.wrapper(final.model.list,
              data = LD.proc, ddl = LD.ddl)
    #
    # Return model table and list of models
    #
    return(final.results)
}

final.results=run.final()


However, I would like to add another model where survival is constant for young-of-the-year (HY) birds, but can but can vary over time for older age classes. In other words, if age == 0, the model for S is ~1, if age >0, the model is age:time:sex
Can anyone tell me how I would code this model?
splash1199
 
Posts: 5
Joined: Mon Feb 04, 2019 2:22 pm

Re: coding a time effect on Survival for just one age class

Postby jlaake » Mon Feb 04, 2019 7:03 pm

I'm thinking you have not read the workshop notes in the documentation zip file or Appendix C of Cooch and White or skipped over the parts dealing with age.

First of all I think you have something incorrect on the process.data unless you aren't giving the full story about how you handled your data. Your code was as follows:

Code: Select all
LD.proc = process.data(data = df,
   model       = "Burnham",
   groups      = c("age", "sex"),
   age.var     = 1,
   initial.age = c(0, 1, 2))
#for age, HY = 0, SY = 1, ASY = 2


Unless you specified the factor levels for age to be HY,SY,ASY in that order then you have the initial.age vector wrong because the default ordering of a factor variable is alphabetical so the initial.age argument should be c(2,0,1). You can tell if you have it correct by looking at the age values by group in the design data for one of the parameters. This issue is covered in several places in various documents. For example, see C15 in Appendix C of Cooch and White electronic book.

Second, how to constrain a temporal effect to a particular age group is also covered in C15 and the workshop notes. Essentially you create a 0/1 variable and then interact the 0/1 variable with time so it is only used where the variable is 1.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: coding a time effect on Survival for just one age class

Postby splash1199 » Wed Feb 27, 2019 3:32 pm

Thanks very much for pointing me towards those resources.

I have a similar, follow-up question as I expand my analysis. How would I code an a constant survival for HY, but an additive time effect for SY and ASY? I have tried the following code for a Brownie recovery model, but the model.matrix() still shows an additive time effect being applied to the HYs.

Code: Select all
pre.df <- as.data.frame(read_csv("./output/LD_recovery_data_1998-2007.csv"))
pre.df <- pre.df[,-1]
#coerce age and sex into factors
pre.df$age <- as.factor(pre.df$age)
pre.df$sex <- as.factor(pre.df$sex)

pre.proc = process.data(data = pre.df,
   model       = "Brownie",
   groups      = c("age", "sex"),
   age.var     = 1,
   initial.age = c(2, 0, 1))
#ASY (2), HY (0), SY (1)

# make the design data from the process data above
pre.ddl = make.design.data(pre.proc)

pre.ddl=add.design.data(data = pre.proc,
   ddl = pre.ddl,
   parameter="S",
   type = "age",
   bins = c(0, 0.5, 1.5, 30),
   right = FALSE,
   name = "age",
   replace = TRUE)
levels(pre.ddl$S$age) = c("HY", "SY", "ASY")
#add a dummy variables; all three age classes
pre.ddl$S$hy     = ifelse(pre.ddl$S$age == "HY", 1, 0)
pre.ddl$S$sy     = ifelse(pre.ddl$S$age == "SY", 1, 0)
pre.ddl$S$asy    = ifelse(pre.ddl$S$age == "ASY", 1, 0)

#repeat for F
pre.ddl=add.design.data(data = pre.proc,
   ddl = pre.ddl,
   parameter="f",
   type = "age",
   bins = c(0, 0.5, 1.5, 30),
   right = FALSE,
   name = "age",
   replace = TRUE)
levels(pre.ddl$f$age) = c("HY", "SY", "ASY")
#dummy variables for all three age classes
pre.ddl$f$hy     = ifelse(pre.ddl$f$age == "HY", 1, 0)
pre.ddl$f$sy     = ifelse(pre.ddl$f$age == "SY", 1, 0)
pre.ddl$f$asy    = ifelse(pre.ddl$f$age == "ASY", 1, 0)

model1 = mark(data = pre.proc,
   ddl = pre.ddl,
  model.parameters = list(
    S = list(formula = ~ -1 + hy + (sy + asy + time)),
    f = list(formula = ~ -1 + age)),
   invisible = FALSE,
   model = "Brownie")

dm <- model.matrix(~ -1 + hy + (sy + asy + time), pre.ddl$S)
splash1199
 
Posts: 5
Joined: Mon Feb 04, 2019 2:22 pm

Re: coding a time effect on Survival for just one age class

Postby jlaake » Wed Feb 27, 2019 6:59 pm

Create another variable that is sy+asy in your S design data. I'll call it ny for not young then use

~ -1 + hy + sy + ny:time

Putting things in parenthesis like you did, didn't do anything as you found. When you want to limit an effect to a particular group, you create a 0/1 variable and use the colon operator which acts like multiplication. May want to read up on R formulas. There are examples of this in the RMark documentation.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: coding a time effect on Survival for just one age class

Postby PPalencia » Tue Jun 21, 2022 10:20 am

Dear,

I'm newbie on survival analysis.

I'm using Rmark to apply CSJ analysis. After an exploratory analysis, I found that the best model is: phi(~ageclass + time)p(~1). My variable "ageclass" has 3 levels: cubs (< 1year), juvenilles (1-2years), adults (>=2years).
Exploring the results, it seems that the variation in phi is more strong in some age groups (in this case cubs) than for others (juvenilles & adults). So, now I’m interested in exploring if phi is constant for some age classes but varies for others. Specifically, I’m thinking in a model in which phi~1 for juveniles and adults, BUT phi~time for cubs.

I found in phidot that this could be analysed by creating a new variable 0/1 that then I interact with time.

So, at this is the first time that I coded and run these type of analysis, I would like to have your confirmation that I did it properly:

Code: Select all
# In the design matrix, I created a new variable
df_lynxOK_ddl$Phi$Cubs <- 0
df_lynxOK_ddl$Phi$Cubs[df_lynxOK_ddl$Phi$age == 0] <- 1
head(df_lynxOK_ddl$Phi)

phiCubs2t <- list(formula=~Cubs:time) #

# run the model
phiCubs2t.p <- mark(df_lynxOK_pro,
                   df_lynxOK_ddl,
                   model.parameters = list(Phi = phiCubs2t, p = p)) # time-age-varying surv & constant p

On the other hand, I have other quesstion:
I run a new model in which I included rabbit relative abundance as covariate of phi. We have a unique rabbit abundance estimate for each year for the entire study area (e.g. year1 - rabbitAbundance = 5; year2 - rabbitAbundance = 3; etc). Considering that, what I expect is very similar values after run a model with phi(~time) or with phi(~rabbitabundance), essentially because each year is associated with a constant rabbit abundance. However, the results of both models are not equivalent… please, can you help to interpret these results? Maybe I coded something wrong…

Finally, in the previous comments of the forum you typed that workshop notes in the documentation zip there are useful materials about that issues. Where can I download these files?

Thanks a lot for your help
All the best
Pablo
PPalencia
 
Posts: 4
Joined: Tue Jun 21, 2022 4:13 am

Re: coding a time effect on Survival for just one age class

Postby jlaake » Tue Jun 21, 2022 6:34 pm

I understand your thought to post as reply to this thread but this really should have been a new thread.

First, your cub:time model looks correct. The rabbit abundance model is not the same because you entered it as numeric so it is fitting a linear model for the survival as a function of rabbit abundance which has just 2 parameters- intercept and slope. Whereas, time model has a separate parameter for each time because it is a factor variable.

The documentation link is given each time you type library (RMark) in a new R session.
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: coding a time effect on Survival for just one age class

Postby PPalencia » Wed Jun 22, 2022 3:17 am

Thanks a lot for your response jlaake! It has been very useful for me
PPalencia
 
Posts: 4
Joined: Tue Jun 21, 2022 4:13 am


Return to RMark

Who is online

Users browsing this forum: No registered users and 12 guests