Predict Time-Varying Covariate

[ctlamb]

Hello,

I am looking to predict the effect of annual hunter harvest on the survival and recruitment of a population. I am using a Pradel Robust Design Model.

I have 8 years of data (occasions), with 3 secondary sessions within each year.

Originally I constructed the model with hunter harvest rate (Kill_F) in the design data, with a unique harvest rate for each year.

However, when I went to predict this effect using "covariate.prediction" I realized the hunter harvest rate must be included as in individual covariate.

I really just want to make a graph that shows the effect of harvest on survival. I currently have the following information from the model where the kill data was in the design matrix:

Code: Select all

                    estimate           se          lcl        ucl
Phi:(Intercept)    1.9951422 4.812250e-01  1.051941300  2.9383431
Phi:Kill_F        -0.2170472 1.308897e-01 -0.473591000  0.0394966
f:(Intercept)      1.3437004 1.258996e+00 -1.123931200  3.8113321
f:Kill_F          -1.5799839 7.908900e-01 -3.130128300 -0.0298394
p:(Intercept)     -2.1753491 1.920780e-01 -2.551822000 -1.7988763
p:SexM             0.2173607 1.394550e-01 -0.055971100  0.4906925
p:TrapTypeRS      -1.5453538 3.501957e-01 -2.231737400 -0.8589701
p:TrapNights       0.0011118 8.217706e-05  0.000950719  0.0012729
p:time2            0.0437163 1.002392e-01 -0.152752500  0.2401851
p:time3           -0.2339337 2.941970e-01 -0.810559900  0.3426925
p:time4           -0.2015211 3.028178e-01 -0.795043900  0.3920018
p:time5           -1.1301934 4.259638e-01 -1.965082400 -0.2953044
p:session2007     -0.1969021 1.828471e-01 -0.555282500  0.1614784
p:session2008      0.0408710 1.852748e-01 -0.322267500  0.4040096
p:session2009      0.0181858 1.963324e-01 -0.366625700  0.4029973
p:session2010     -0.6265472 2.243603e-01 -1.066293400 -0.1868009
p:session2011     -0.6089051 2.292836e-01 -1.058301000 -0.1595093
p:session2012      0.7997242 3.234028e-01  0.165854600  1.4335938
p:session2013     -0.5813095 2.762739e-01 -1.122806400 -0.0398127
p:SexM:TrapTypeRS  1.0341984 1.836248e-01  0.674293800  1.3941030


Can I simply back transform the Phi:Kill_F = -0.2170472 to get a slope?


I attempted to include the hunter kill as an individual covariate following the instructions in the MARK book, but I was stumped when it came to including continuous time-varying individual covariates. I kept getting an error when I tried to run my models

Code: Select all


##Import captures, specify column classes
SR<- read.csv("SR_GBPU_DATA_Short.txt" ,
                 colClasses=c( Individual="factor",ch="character", Sex="factor", freq="numeric"),
                 strip.white=FALSE)

##Add individual covariates
SR$KillFem2006<-2
SR$KillFem2007<-3
SR$KillFem2008<-6
SR$KillFem2009<-1
SR$KillFem2010<-5
SR$KillFem2011<-8
SR$KillFem2012<-1

##################
##Process Data
SR.proc=process.data(SR, model="RDPdfHuggins", groups="Sex", begin.time=2006,time.intervals=c(0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,1,0,0))

####################################### 

###SR MODEL
##Sex + TrapType + TrapNights + time + session + (Sex * TrapType)
SR1<-mark(
               data=SR.proc,
               ddl=SR.ddl,
               model.parameters =
               list(
               p=list(formula=~Sex + TrapType + TrapNights + time + session+ (Sex * TrapType),share=TRUE),
               Phi=list(formula=~KillFem),
               f=list(formula=~1)
               ),
               output=FALSE,model="RDPdfHuggins",
               delete=TRUE)




Error in make.mark.model(data.proc, title = title, parameters = model.parameters, :

Following covariates are duplicates of another covariate within the first 10 characters
KillFem2007, KillFem2008, KillFem2009, KillFem2011, KillFem2012
Error in mark(data = SR.proc, ddl = SR.ddl, model.parameters = list(p = list(formula = ~Sex + :
Misspecification of model or internal error in code

Any help would be much appreciated!!!

[ctlamb]

I was able to fix the error above by shortening the column names of individual covariates, such that the names were <10 characters, as required my MARK

Code: Select all


##Import captures, specify column classes
SR<- read.csv("SR_GBPU_DATA_Short.txt" ,
                 colClasses=c( Individual="factor",ch="character", Sex="factor", freq="numeric"),
                 strip.white=FALSE)

##Add individual covariates, ****less than 10 characters long****
SR$Kill2006<-2
SR$Kill2007<-3
SR$Kill2008<-6
SR$Kill2009<-1
SR$Kill2010<-5
SR$Kill2011<-8
SR$Kill2012<-1[/b]

##################
##Process Data
SR.proc=process.data(SR, model="RDPdfHuggins", groups="Sex", begin.time=2006,time.intervals=c(0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,1,0,0))

####################################### 

###SR MODEL
##Sex + TrapType + TrapNights + time + session + (Sex * TrapType)
SR1<-mark(
               data=SR.proc,
               ddl=SR.ddl,
               model.parameters =
               list(
               p=list(formula=~Sex + TrapType + TrapNights + time + session+ (Sex * TrapType),share=TRUE),
               Phi=list(formula=~KillFem),
               f=list(formula=~1)
               ),
               output=FALSE,model="RDPdfHuggins",
               delete=TRUE)

Now my issue is back to making the predictions:

My design data looks for Phi like this:

Code: Select all

summary(SR1,se=TRUE)$real$Phi


all.diff.index par.index estimate se lcl ucl fixed note group age time Age Time Sex
Phi gF a0 t2006 1 1 0.8423074 0.0392169 0.7496610 0.9050113 F 0 2006 0 0 F
Phi gF a1 t2007 2 2 0.7310661 0.0319770 0.6640192 0.7889851 F 1 2007 1 1 F
Phi gF a2 t2008 3 3 0.8404231 0.0385227 0.7499571 0.9024171 F 2 2008 2 2 F
Phi gF a3 t2009 4 4 0.8261047 0.0330206 0.7517084 0.8817171 F 3 2009 3 3 F
Phi gF a4 t2010 5 5 0.7434315 0.0255843 0.6901849 0.7903071 F 4 2010 4 4 F
Phi gF a5 t2011 6 6 0.7059933 0.0470437 0.6062991 0.7892191 F 5 2011 5 5 F
Phi gF a6 t2012 7 7 0.8458767 0.0405073 0.7490638 0.9098345 F 6 2012 6 6 F
Phi gM a0 t2006 8 1 0.8423074 0.0392169 0.7496610 0.9050113 M 0 2006 0 0 M
Phi gM a1 t2007 9 2 0.7310661 0.0319770 0.6640192 0.7889851 M 1 2007 1 1 M
Phi gM a2 t2008 10 3 0.8404231 0.0385227 0.7499571 0.9024171 M 2 2008 2 2 M
Phi gM a3 t2009 11 4 0.8261047 0.0330206 0.7517084 0.8817171 M 3 2009 3 3 M
Phi gM a4 t2010 12 5 0.7434315 0.0255843 0.6901849 0.7903071 M 4 2010 4 4 M
Phi gM a5 t2011 13 6 0.7059933 0.0470437 0.6062991 0.7892191 M 5 2011 5 5 M
Phi gM a6 t2012 14 7 0.8458767 0.0405073 0.7490638 0.9098345 M 6 2012 6 6 M

All I really want to do is get a general feel for how kill rate affects survival, so can I just provide a data frame with varying kill rates and predict the effect on a single parameter, i.e., parameter 1, Female survival in 2006.

Like this, but not sure how to specify the time covariate, as per the MARK the help " If a time-varying covariate is used then you need to specify the covariate name with the time index included as it is specified in the data":

Code: Select all

data<-data.frame(Kill=runif(10,min=min(SR.ddl$Phi$Kill_F), max=max(SR.ddl$Phi$Kill_F)),
           Time=1,
           time=2006,
           Sex="F"
           )
           
           
 Phibykill<-covariate.predictions(
                                 SR1,
                                 data=data,
                                 indices=c(1)
                                 )

> Phibykill$estimates$estimate
 [1] 0.8423073 0.8423073 0.8423073 0.8423073 0.8423073 0.8423073 0.8423073 0.8423073 0.8423073 0.8423073


The estimates don't vary with varying kill rates, even if I make the time variables span the range of dates sampled, such as runif(10,2006,2013).

Any thoughts?

[jlaake]

I have answered this question several times on this forum. When I did a search on time varying covariates it returned 13 pages of relevant posts. In particular see http://www.phidot.org/forum/viewtopic.php?f=21&t=3044&p=9773&hilit=time+varying+covariates#p9773. Also this issue is discussed in sections 6.7 and 12 of the RMark workshop notes at http://www.phidot.org/software/mark/rmark/RMarkDocumentation.zip.

If your kill covariate is a design covariate you can create the graph yourself using the estimates, link function and the beta variance-covariance matrix. See ch 6 of Evan and Gary's book and sec 6.7 of the RMark workshop notes and the first url above. If you use a time-varying individual covariate, the you have to specify the exact name of the covariate (Kill2008) and the parameter index where Kill2008 is used. Using Kill will not work because it doesn't match any of the covariate names(e.g. KillYYYY covariates). Using Kill in the formula works because it creates KillYYYY in the design matrix for each year. It doesn't work in covariate.predicitions because you have to specify which covariate you are assigning values. If there is no other time component to the model using the covariate (eg. using annual weather as well) then you can choose any one of the values.

[ctlamb]

Thanks, Jeff.

To conclude this thread, as stated by Jeff above, I needed to specify the exact name of the covariate I wanted to plot. Initially I had just coded the prediction of "Kill" and then added in a second column depicting what year, but what I needed to do was code the prediction of "KillYYYY", YYYY being some year, and also include the year index (time) in a subsequent column.

I was able to produce the desired graph with the following code using the individual kill covariates for each year:

Code: Select all

 # Compute and plot survival and recruitment values

data<-data.frame(Kill2007=runif(100,min=min(SR.ddl$Phi$Kill_F), max=max(SR.ddl$Phi$Kill_F)),
           Time=2,
           time=2007
           )
           
           
 Phibykill<-covariate.predictions(
   TopSR,
                                 data=data,
                                 indices=c(2)
                                 )


 # Plot predicted model estimates
# with pointwise confidence interval
Phi<-Phibykill$estimates

 plot(  Phi$Kill2007,  Phi$estimate,
             type="l",lwd=2,xlab="Female Kill (% Fem Pop)",ylab="Apparent Survival",
             ylim=c(0.5,1),las=1)
lines(Phi$Kill2007,Phi$lcl,lty=2)
lines(Phi$Kill2007,Phi$ucl,lty=2)