Stage specific nest survival in RMark

[bacollier]

OK, as an update, I just ran 2 models, one under 2.1.0 and most updated MARK with a long covariate name (canhtmean100) and the same model under 2.1.1. and most updated MARK with a short covariate name (chm). Same model in both cases (simple single season occupancy model)

occurrence=BaseCC + canhtmean100
detection=SurveyDayn

BaseCC and canhtmean100/chm are both continuous covariates, SurveyDayn is a time-specific covariate for day of sampling from a start date. Other than my changing the column name in .csv from canhtmean100 to chm, everything else was the same. As far as I can tell, all values for beta, reals, -2LL, AIC, etc. match up the same in both cases, output files look identical other than covariate names.

So, if it does not change the parameter estimates (granted this is only one test), what is going on? One thing it will change is that I am using short names in the future unless Gary decides to lengthen it

Bret

[jlaake]

From a test Murray suggested I believe it uses the first variable with that name in each part of the DM. You can test by running the model with

~canhtmean###

instead of ~canhtmean

where ### is the first time prefix. When we did that with the nest stage example, the outputs were identical. Note that this primarily a problem with time-varying individual covariates and it will not affect individual covariates that are not time dependent unless you use two covariates that are the same for the first 10 characters (eg myweightofpup and myweightofadult).

--jeff

[gwhite]

MARK is only using the first 10 characters, so you can have longer names, but the extra characters are ignored. Where you get in trouble is if the distinction between 2 variables is only after 10 characters, the value used in the design matrix will only be the value from the first of these covariates.

Gary

[bacollier]

MARK is only using the first 10 characters, so you can have longer names, but the extra characters are ignored. Where you get in trouble is if the distinction between 2 variables is only after 10 characters, the value used in the design matrix will only be the value from the first of these covariates.

Gary

Thanks Gary, now this all makes sense. The way (I think) Jeff has RMark is that the distinction for the time-specific variables is the number at the end of the covariate name. So, if the time-specific covariate name for 3 time-specific variables is: Iammorethanten1, Iammorethanten2, and Iammorethanten3, then it won't see anything after Iammoretha and the distinction between times is lost and hence bad results occur.

It seems that if one does not have 'standard' covariates (e.g., non-time dependent) with the same name, then any length is fine as long as they are not of equivalent text.

Jeff, is there a way to change the blanket '>10" error in RMark to only be on those covariates that have duplicate names? That way the 'check' occurs, but not on variables that only occur once?

Bret

Bret

[jlaake]

Jeff, is there a way to change the blanket '>10" error in RMark to only be on those covariates that have duplicate names? That way the 'check' occurs, but not on variables that only occur once?

Done and posted to github.

[B.K. Sandercock]

I'm interested in calculating stage-specific nest survival with RMark. The previous series of posts showed how to fit the model but I'm wondering about the best way to extract the real parameter estimates for the different stages from the model. The attached code with predict_real seems to work based on back-calculation from the beta coefficients, but I'm just wondering if this is the best way. I tried the functions for compute.real and covariate.predictions but was not able to sort out the correct arguments. Thanks in advance if anybody can help with some sample code to look at. Thanks, Brett.

Code: Select all

library(RMark)

# stage is 1 if less than cutoff value, 0 if greater than cutoff value
create.stage.var=function(data,agevar.name,stagevar.name,time.intervals,cutoff)
{
nocc=length(time.intervals)
age.mat=matrix(data[,agevar.name],nrow=dim(data)[1],ncol=nocc-1)
age.mat=t(t(age.mat)+cumsum(c(0,time.intervals[1:(nocc-2)])))
stage.mat=age.mat<cutoff
stage.mat=t(apply(stage.mat,1,function(x) as.numeric(x)))
stage.mat=data.frame(stage.mat)
names(stage.mat)=paste(stagevar.name,1:(nocc-1),sep="")
return(stage.mat)
}

# input data
data(mallard)
cutoff = 10
x=create.stage.var(mallard,"AgeDay1","Inc",rep(1,90),cutoff)
mallard = cbind(mallard,x)

# process data and run model
dp = process.data(mallard, model='Nest', nocc=90)
ddl = make.design.data(dp)
model = mark(dp,ddl,model.parameters=list(S=list(formula=~Inc)),nocc=90)

# extract estimates for Inc stage equals 1 (<cutoff) or 0 (>cutoff)
pred = predict_real(model, ddl$S[1,,drop=FALSE], "S", replicate=TRUE, data=data.frame(Inc=c(1,0)))

# make a little table
stage= c(paste("<", cutoff, sep=""), paste(">", cutoff, sep=""))
data.frame(stage, pred$real, pred$se, pred$lcl, pred$ucl)

With the following result where the cutoff value was 10 days:

Code: Select all

  stage pred.real     pred.se  pred.lcl  pred.ucl
1   <10 0.9325094 0.009598268 0.9110864 0.9490593
2   >10 0.9557878 0.002761649 0.9500477 0.9608954

[jlaake]

What you did was fine. The following also works. Obviously it would be more complicated if you had a more complex model that involved nest age.

Code: Select all

 covariate.predictions(model,data=data.frame(Inc1=c(0,1),index=c(1,1)))

$estimates
  vcv.index model.index par.index Inc1 index  estimate          se       lcl       ucl fixed
1         1           1         1    0     1 0.9557878 0.002761649 0.9500477 0.9608954     
2         1           1         1    1     1 0.9325094 0.009598267 0.9110864 0.9490593     

$vcv
              [,1]          [,2]
[1,]  7.626703e-06 -2.638670e-06
[2,] -2.638670e-06  9.212674e-05



[JClarke]

I am currently working through a similar issue as the original post except that I have multiple stages. I know it was mentioned that the function could be modified to incorporate multiple stages (for instance a nest that was in the incubation stage for ages 1-10 and nestling stage for ages 10-22). However, I've tried editing the code myself to get the multiple stages to work but it doesn't seem to be functioning correctly. Could someone point me in the correct direction?

Thanks

[jlaake]

Been awhile since I wrote that function but i believe it creates a time varying individual covariate that has 0/1 values. 0 when not in the stage and 1 when it is in the stage. So a multi-stage variable would be a factor variable and individual covariates (time varying in this case) have to be numeric. If you have k stages then you need to create k-1 numeric 0/1 variables. The stage you leave out will be the intercept value. The function i wrote only uses a single cutoff value. You need to modify it to have an upper and lower cutoff and then run it k-1 times once for each stage using a different stage name. Be aware of the limitation of 10 characters on individual covariate names which has to include the time suffix part of the name. You will need to use each of the k-1 variable names for a stage specific model.

Good on you for scouring phidot posts for partial answers to your questions. You used phidot as a resource as it was intended. Kudos. Let me know if you need more help writing the code.

Jeff

[JClarke]

Thanks for the quick reply. I believe that I misunderstood the original post and it seems designed for two stages already. Your explanation of using k-1 numeric 0/1 variables for k stages cleared things up. It sounds like, if I'm only working with two stages (incubation and nestling), I can actually run the code as is.

test <- create.stage.var(CCSP.surv,
"AgeDay1",
"Incub",
rep(1,max(CCSP.surv$LastChecked)),
12)

If I understand this correctly, any nest with an AgeDay1 value of >12 (Incubation) will be assigned a value of 1 and anything <12 will be assigned a value of 0 (Nestling).

My last question is, that a nest with an AgeDay1 value of -2 will have 14 columns in the incubation stage. When I run this through the mark model, will it function similar to incorporating NestAge where nests will be excluded until "the add(AgeDay1,x) value is zero"?