Output Covariate Specific Survival Estimates

[alittle]

I'm working on a known-fate analysis for an avian species and want to evaluate the effect of reproductive activity (initiated a nest vs. did not initiate a nest). I've added the covariate (Y = Yes bird initiated a nest; and N = No bird did not initiate a nest). I ran a set of models (see output below):

model npar AICc DeltaAICc weight Deviance
2 S(~repro.cov) 2 172.1767 0.000000 7.184835e-01 168.16174
1 S(~1) 1 174.4164 2.239699 2.344620e-01 34.31768
4 S(~yr) 3 177.6293 5.452550 4.703386e-02 33.50554
3 S(~time) 14 193.0915 20.914843 2.064507e-05 26.46615

I would like to model average the results and output weekly survival estimates (data was set-up in a weekly fashion). Specifically, I would like to obtain survival estimates for reproductively active vs. inactive. However, I cannot seem to figure out the correct coding. I tried the following code:

S.avg<-model.average(repro.results,vcv=T)$estimates
S.avg<-subset(S.avg,select=c(repro,estimate,lcl,ucl))
write.csv(S.avg,file="estimates.w.repro.csv",row.names=F)

but I receive the following error: Error in eval(expr, envir, enclos) : object 'repro' not found. Maybe it's something very obvious that I'm missing. Any assistance would be much appreciated! Thanks!

[jlaake]

This is a very specific question to RMark and should have been posted to that forum. Please post to the correct forum so Evan doesn't have to move them around.

When you are reporting an error message show the error interwined with the code so it is clear where it happened. I'm guessing that it happened at:

Code: Select all

S.avg<-subset(S.avg,select=c(repro,estimate,lcl,ucl))

Did you check to see whether repro was a field in S.avg? My guess is that it is called repro.cov from the formula in the model.

--jeff

[alittle]

This is a very specific question to RMark and should have been posted to that forum. Please post to the correct forum so Evan doesn't have to move them around.

When you are reporting an error message show the error interwined with the code so it is clear where it happened. I'm guessing that it happened at:

Code: Select all

S.avg<-subset(S.avg,select=c(repro,estimate,lcl,ucl))

Did you check to see whether repro was a field in S.avg? My guess is that it is called repro.cov from the formula in the model.

--jeff

I apologize for not adding this question to the RMark section and will make sure to do so in the future. The following output is what I obtained from the S.avg dataframe, which does not include repro as a field:

par.index estimate se lcl ucl
1 1 0.9778946 0.006371892 0.9612777 0.9874736
2 2 0.9778948 0.006371200 0.9612802 0.9874730
3 3 0.9778951 0.006370934 0.9612813 0.9874730
4 4 0.9778948 0.006371391 0.9612795 0.9874732
5 5 0.9778951 0.006370947 0.9612813 0.9874730
6 6 0.9778948 0.006371391 0.9612795 0.9874732
7 7 0.9778947 0.006371430 0.9612794 0.9874732
8 8 0.9778951 0.006370954 0.9612812 0.9874730
9 9 0.9778951 0.006370970 0.9612812 0.9874730
10 10 0.9778947 0.006371472 0.9612792 0.9874732
11 11 0.9778954 0.006371281 0.9612804 0.9874736
12 12 0.9778954 0.006371281 0.9612804 0.9874736
13 13 0.9778948 0.006371757 0.9612783 0.9874736
14 14 0.9778954 0.006371281 0.9612804 0.9874736
15 15 0.9784176 0.005350470 0.9650226 0.9867533
16 16 0.9784179 0.005349621 0.9650254 0.9867525
17 17 0.9784182 0.005349278 0.9650267 0.9867523
18 18 0.9784178 0.005349855 0.9650246 0.9867527
19 19 0.9784181 0.005349294 0.9650267 0.9867523
20 20 0.9784178 0.005349855 0.9650246 0.9867527
21 21 0.9784178 0.005349903 0.9650245 0.9867528
22 22 0.9784181 0.005349303 0.9650266 0.9867523
23 23 0.9784181 0.005349324 0.9650266 0.9867524
24 24 0.9784178 0.005349953 0.9650243 0.9867528
25 25 0.9784185 0.005349661 0.9650258 0.9867531
26 26 0.9784185 0.005349661 0.9650258 0.9867531
27 27 0.9784179 0.005350284 0.9650234 0.9867533
28 28 0.9784185 0.005349661 0.9650258 0.9867531
29 29 0.9785430 0.005542476 0.9645193 0.9870980
30 30 0.9785432 0.005541651 0.9645221 0.9870972
31 31 0.9785435 0.005541314 0.9645235 0.9870971
32 32 0.9785432 0.005541878 0.9645214 0.9870974
33 33 0.9785435 0.005541329 0.9645234 0.9870971
34 34 0.9785432 0.005541878 0.9645214 0.9870974
35 35 0.9785432 0.005541924 0.9645212 0.9870975
36 36 0.9785435 0.005541338 0.9645234 0.9870971
37 37 0.9785435 0.005541359 0.9645233 0.9870971
38 38 0.9785432 0.005541974 0.9645210 0.9870975
39 39 0.9785438 0.005541676 0.9645226 0.9870978
40 40 0.9785438 0.005541676 0.9645226 0.9870978
41 41 0.9785433 0.005542290 0.9645201 0.9870980
42 42 0.9785438 0.005541676 0.9645226 0.9870978

Therefore, I changed the coding to repro.cov as stated earlier from the formula of the model.
> S.repro.output<-model.average(repro.results,vcv=T)$estimates
> S.final<-subset(S.repro.output,select=c(repro.cov,estimate,lcl,ucl))
Error in eval(expr, envir, enclos) : object 'repro.cov' not found

I've also included my entire function for running the MARK models:
#Run MARK Models
repro.models=function()
{
#Define range of models for S
S.dot=list(formula=~1)
S.repro=list(formula=~repro.cov)
S.yr<-list(formula=~yr)
S.time<-list(formula=~time)

#Create Model List
model.list=create.model.list("Known")

#Run MARK Wrapper
results=mark.wrapper(model.list,data=turkey.processed,ddl=turkey.ddl,invisible=FALSE,threads=2)

# Return model table and list of models
return(results)
}
repro.results=repro.models()
repro.results

I'm not entirely sure what is going on as I felt like this coding would've worked to allow me to output the survival estimates.

[jlaake]

I'm not sure why you thought that was going to work given that you showed the dataframe and it doesn't contain a variable (column) with the name repro.cov.

You didn't include all of the pertinent code. What is repro.cov? A factor variable used to define groups or an individual covariate. If it is an individual covariate you should be using covariate.predictions but given you defined it as Y/N then I assume it is a factor variable. When you use model.average without specifying a parameter it does not include the design data fields because the design data can vary across parameters so there is no way to create a single dataframe with variable columns. Now, this model is an exception with only one parameter but the code doesn't make an exception. If you want to include the design data use

Code: Select all

S.avg<-model.average(repro.results,parameter="S",vcv=T)$estimates


See ?model.average.marklist where that is explained. You can help yourself work this out by stepping through the code and seeing what is going wrong.

[alittle]

I'm not sure why you thought that was going to work given that you showed the dataframe and it doesn't contain a variable (column) with the name repro.cov.

You didn't include all of the pertinent code. What is repro.cov? A factor variable used to define groups or an individual covariate. If it is an individual covariate you should be using covariate.predictions but given you defined it as Y/N then I assume it is a factor variable. When you use model.average without specifying a parameter it does not include the design data fields because the design data can vary across parameters so there is no way to create a single dataframe with variable columns. Now, this model is an exception with only one parameter but the code doesn't make an exception. If you want to include the design data use

Code: Select all

S.avg<-model.average(repro.results,parameter="S",vcv=T)$estimates


See ?model.average.marklist where that is explained. You can help yourself work this out by stepping through the code and seeing what is going wrong.

I'm learning....the code previously provided was built by another individual for my data as I'm learning about R and RMark. Repro.cov is an individual-specific covariate. In the Excel data set, repro.cov was listed as Y=Yes or N=No but was converted to Y=1 and N=0 for the analysis (see below).

I've included an example of the code for the analysis in case I'm missing something pertinent:

*Note not all code is provided due to the lengthiness of the code.

breeding.3<-subset(breeding,yr==3)

#Define 7-day periods
first.date<-breeding.3$date.code[3]
first.day<-first.date$yday
day<-breeding.3$date.code$yday-first.day+1
breeding.3$per<-as.integer(day/7)+1

convert.animal.id<-function(x,data)
{
levels<-levels(factor(data$animal_id))
index<-match(x,levels)
return(index)
}

breeding.3$id.no<-convert.animal.id(breeding.3$animal_id,data=breeding.3)

n.animals<-max(breeding.3$id.no)

#n.occasions<-max(breeding.3$per)
n.occasions<-14

#Create empty capture histories
y<-array(0,dim=c(n.animals,2*n.occasions))

data.set<-subset(breeding.3,select=c(id.no,per,Status))
for(i in 1:nrow(data.set))
{

status<-data.set[i,3]
id<-data.set[i,1]
per<-data.set[i,2]
ind1<-(per-1)*2+1
ind2<-per*2
status<-substr(status,1,5) #because some of these had spaces at the end
if (status=="Alive") y[id,ind1]=1
if (status=="Dead") {y[id,ind1]=y[id,ind2]=1 }
}

pasty<-function(x)
{
k<-ncol(x)
n<-nrow(x)
out<-array(dim=n)
for (i in 1:n)
{
out[i]<-paste(x[i,],collapse="")
}
return(out)
}

capt.hist.3<-pasty(y)

#Covariate
breeding.3$repro.cov[substr(breeding.3$repro,1,1)=="N"]=0
breeding.3$repro.cov[substr(breeding.3$repro,1,1)=="Y"]=1

cov<-with(breeding.3,aggregate(repro.cov~id.no,FUN=mean,na.action=NULL))

mark.data.3<-data.frame(ch=capt.hist.3,yr="3",repro.cov=cov$repro.cov)

#Eliminate capture histories with missing values
mark.data.3<-subset(mark.data.3,!is.na(repro.cov))
mark.data.3

##Then combine the data and run through MARK with some models including year effect

mark.data.final<-rbind(mark.data.1,mark.data.2,mark.data.3)

#Process the data
turkey.processed = process.data(data=mark.data.final,model="Known",groups=c("yr"))

#Create default design data
turkey.ddl=make.design.data(turkey.processed)

#Run MARK Models
repro.models=function()
{
#Define range of models for S
S.dot=list(formula=~1)
S.repro=list(formula=~repro.cov)
S.yr<-list(formula=~yr)
S.time<-list(formula=~time)

#Create Model List
model.list=create.model.list("Known")

#Run MARK Wrapper
results=mark.wrapper(model.list,data=turkey.processed,ddl=turkey.ddl,invisible=FALSE,threads=2)

# Return model table and list of models
return(results)
}
repro.results=repro.models()
repro.results

#Model Averaging Estimates
S.avg<-model.average(repro.results,parameter="S",vcv=T)$estimates

S.avg


row.names par.index estimate se lcl ucl fixed note group age time Age Time yr
1 S g1 a1 t2 1 0.9778946 0.006371892 0.9612777 0.9874736 1 1 2 1 0 1
2 S g1 a2 t3 2 0.9778948 0.006371200 0.9612802 0.9874730 1 2 3 2 1 1
3 S g1 a3 t4 3 0.9778951 0.006370934 0.9612813 0.9874730 1 3 4 3 2 1
4 S g1 a4 t5 4 0.9778948 0.006371391 0.9612795 0.9874732 1 4 5 4 3 1
5 S g1 a5 t6 5 0.9778951 0.006370947 0.9612813 0.9874730 1 5 6 5 4 1
6 S g1 a6 t7 6 0.9778948 0.006371391 0.9612795 0.9874732 1 6 7 6 5 1
7 S g1 a7 t8 7 0.9778947 0.006371430 0.9612794 0.9874732 1 7 8 7 6 1
8 S g1 a8 t9 8 0.9778951 0.006370954 0.9612812 0.9874730 1 8 9 8 7 1
9 S g1 a9 t10 9 0.9778951 0.006370970 0.9612812 0.9874730 1 9 10 9 8 1
10 S g1 a10 t11 10 0.9778947 0.006371472 0.9612792 0.9874732 1 10 11 10 9 1
11 S g1 a11 t12 11 0.9778954 0.006371281 0.9612804 0.9874736 1 11 12 11 10 1
12 S g1 a12 t13 12 0.9778954 0.006371281 0.9612804 0.9874736 1 12 13 12 11 1
13 S g1 a13 t14 13 0.9778948 0.006371757 0.9612783 0.9874736 1 13 14 13 12 1
14 S g1 a14 t15 14 0.9778954 0.006371281 0.9612804 0.9874736 1 14 15 14 13 1
15 S g2 a1 t2 15 0.9784176 0.005350470 0.9650226 0.9867533 2 1 2 1 0 2
16 S g2 a2 t3 16 0.9784179 0.005349621 0.9650254 0.9867525 2 2 3 2 1 2
17 S g2 a3 t4 17 0.9784182 0.005349278 0.9650267 0.9867523 2 3 4 3 2 2
18 S g2 a4 t5 18 0.9784178 0.005349855 0.9650246 0.9867527 2 4 5 4 3 2
19 S g2 a5 t6 19 0.9784181 0.005349294 0.9650267 0.9867523 2 5 6 5 4 2
20 S g2 a6 t7 20 0.9784178 0.005349855 0.9650246 0.9867527 2 6 7 6 5 2
21 S g2 a7 t8 21 0.9784178 0.005349903 0.9650245 0.9867528 2 7 8 7 6 2
22 S g2 a8 t9 22 0.9784181 0.005349303 0.9650266 0.9867523 2 8 9 8 7 2
23 S g2 a9 t10 23 0.9784181 0.005349324 0.9650266 0.9867524 2 9 10 9 8 2
24 S g2 a10 t11 24 0.9784178 0.005349953 0.9650243 0.9867528 2 10 11 10 9 2
25 S g2 a11 t12 25 0.9784185 0.005349661 0.9650258 0.9867531 2 11 12 11 10 2
26 S g2 a12 t13 26 0.9784185 0.005349661 0.9650258 0.9867531 2 12 13 12 11 2
27 S g2 a13 t14 27 0.9784179 0.005350284 0.9650234 0.9867533 2 13 14 13 12 2
28 S g2 a14 t15 28 0.9784185 0.005349661 0.9650258 0.9867531 2 14 15 14 13 2
29 S g3 a1 t2 29 0.9785430 0.005542476 0.9645193 0.9870980 3 1 2 1 0 3
30 S g3 a2 t3 30 0.9785432 0.005541651 0.9645221 0.9870972 3 2 3 2 1 3
31 S g3 a3 t4 31 0.9785435 0.005541314 0.9645235 0.9870971 3 3 4 3 2 3
32 S g3 a4 t5 32 0.9785432 0.005541878 0.9645214 0.9870974 3 4 5 4 3 3
33 S g3 a5 t6 33 0.9785435 0.005541329 0.9645234 0.9870971 3 5 6 5 4 3
34 S g3 a6 t7 34 0.9785432 0.005541878 0.9645214 0.9870974 3 6 7 6 5 3
35 S g3 a7 t8 35 0.9785432 0.005541924 0.9645212 0.9870975 3 7 8 7 6 3
36 S g3 a8 t9 36 0.9785435 0.005541338 0.9645234 0.9870971 3 8 9 8 7 3
37 S g3 a9 t10 37 0.9785435 0.005541359 0.9645233 0.9870971 3 9 10 9 8 3
38 S g3 a10 t11 38 0.9785432 0.005541974 0.9645210 0.9870975 3 10 11 10 9 3
39 S g3 a11 t12 39 0.9785438 0.005541676 0.9645226 0.9870978 3 11 12 11 10 3
40 S g3 a12 t13 40 0.9785438 0.005541676 0.9645226 0.9870978 3 12 13 12 11 3
41 S g3 a13 t14 41 0.9785433 0.005542290 0.9645201 0.9870980 3 13 14 13 12 3
42 S g3 a14 t15 42 0.9785438 0.005541676 0.9645226 0.9870978 3 14 15 14 13 3


S.final<-subset(S.avg,select=c(repro.cov,estimate,lcl,ucl))
write.csv(S.final,file="estimates.w.repro.csv",row.names=F)

Using the code that you provided, it still doesn't contain a variable (column) with the name repro.cov. Am I missing something obvious?

[jlaake]

I realize you are learning which is why I suggested that you step through and examine the results at each step. If something goes wrong, then look at the help for that function which was the base R function subset in that case. That is how you learn but you have to be patient.

If your animals are in either the repro.cov=0 or 1 for the whole study then you should it as a grouping variable (groups=c("yr","repro.cov")). Since it is numeric you'll get a warning that it is being converted to a factor variable or if you simply left it as Y/N then it would be a factor variable. If you do that, estimates will be constructed for each value of repro.cov and it will be in the design data and it will appear in the model.average dataframe. If you leave it as an individual covariate then you need to use covariate.predictions instead of model.average and specify that you want to compute the value of S at both repro.cov=0 and repro.cov=1. It computes those values and model.averages over the set of models. Using it as a grouping variable is easier. But if repro.cov is changing for each animal then you have an entirely different issue because you are venturing into time-varying individual covariates.

--jeff

[alittle]

I realize you are learning which is why I suggested that you step through and examine the results at each step. If something goes wrong, then look at the help for that function which was the base R function subset in that case. That is how you learn but you have to be patient.

If your animals are in either the repro.cov=0 or 1 for the whole study then you should it as a grouping variable (groups=c("yr","repro.cov")). Since it is numeric you'll get a warning that it is being converted to a factor variable or if you simply left it as Y/N then it would be a factor variable. If you do that, estimates will be constructed for each value of repro.cov and it will be in the design data and it will appear in the model.average dataframe. If you leave it as an individual covariate then you need to use covariate.predictions instead of model.average and specify that you want to compute the value of S at both repro.cov=0 and repro.cov=1. It computes those values and model.averages over the set of models. Using it as a grouping variable is easier. But if repro.cov is changing for each animal then you have an entirely different issue because you are venturing into time-varying individual covariates.

--jeff

The data was collected over 2.5 years with some individuals being reproductively active in one year and not the next; therefore, I would suspect that this is a time-varying covariate case. For example, individual 1 was caught in 2011 and monitored until 2013. She was reproductively active the first 2 years but not the last year.

I attempted the following code based on the RMark manual but obtain survival estimates for every individual in the dataset. Is there a way to just extract the survival estimate for reproductively active (value=1) and inactive (value=0)?
Sbyrepro<-covariate.predictions(repro.results,data=mark.data.final,indices=c(1))

[jlaake]

I can't tell for sure but it looks like you have split the capture history into different years and then used yr as a grouping variable. So the animal observed over 3 years has 3 capture histories and for the first 2 years it has the value repro.cov=1 and in the last year it is 0. If that is the case then for any one year the individual is either repro.cov=0 or 1 and doesn't vary within the year so you could use it as grouping variable.

[bacollier]

I'm working on a known-fate analysis for an avian species and want to evaluate the effect of reproductive activity (initiated a nest vs. did not initiate a nest). I've added the covariate (Y = Yes bird initiated a nest; and N = No bird did not initiate a nest). I ran a set of models (see output below):

model npar AICc DeltaAICc weight Deviance
2 S(~repro.cov) 2 172.1767 0.000000 7.184835e-01 168.16174
1 S(~1) 1 174.4164 2.239699 2.344620e-01 34.31768
4 S(~yr) 3 177.6293 5.452550 4.703386e-02 33.50554
3 S(~time) 14 193.0915 20.914843 2.064507e-05 26.46615

I would like to model average the results and output weekly survival estimates (data was set-up in a weekly fashion). Specifically, I would like to obtain survival estimates for reproductively active vs. inactive. Thanks!

Andy,
Isn't this the same problem I walked you through a while back on how to set up the models in MARK? Are you just trying to get year specific estimates of survival and a function of whether they nested or not each year?

Are you just trying to rebuild those analyses in RMark or is there something different you are trying to do as based on the results of this thread so far I cannot seem to figure out what you are trying to figure out?

\bret

[alittle]

Bret,

Thanks for the message! The problem that I was dealing with has been addressed and I was able to output the weekly survival estimates per reproductively active vs. inactive covariate. I ended up taking an RMark and MARK class this spring to help me build a solid foundation for survival and population estimation analyses. In order to beef up my R skills, I forced myself to conduct the same analysis that you assisted me with in RMark. This was a great experience and I now am more R savvy because of it. Thanks again for all your help!