RMark: Plotting the effects of time-varying covariates

[bootzies]

Hello list.

I am running a multi-state analysis, in which I have multiple 'design' (time-varying, population-level) enviromental covariates in my model for S (apparent survival). I am trying to work out the best way to present my results, but struggling a bit.

I have found quite a bit of info online about plotting the effect of individual covariates, but not much on how to do the same for time-varying/design covariates. I would like to show plots of predicted survival against a range of values of my environmental covariates. But, I can't work out a way to do this when I have multiple covariates in my models.

For example, if my model for S has the covariates winter temperature and winter rainfall: I have two covariates affecting survival at each time step, so I can't just plot the real estimate of survival against the value of one covariate (eg. rain) at the same time step, as the S value will also be representative of the effects of the other covariate (temperature) at that same time step (so won't be representative of only the rain effect). Besides I think I have read that the real parameter estimates in the model output represent estimates at the mean values for the covariates.....

So, how can I create plots of predicted S where one covariate (e.g. temp) is held at the mean and the other (e.g. rain) is varied within the range of sampled values? This is a pretty straight forward thing to do in regression type models, but I can't seem to work it out for RMark!

I have read that you can use the function covariate.predictions (documented for use with individual covariates) to produce these types of estimates; but I don't understand how (and haven't been able to find a worked example of this). Would I need to go back to the encounter history data frame and add in the values for each of my covariates alongside each individual's ch data in columns for each year (with the same value for each year among individuals)? This seems like an extremely long-winded work around (especiallywhen I am using multiple covariates and 30 years worth of data on >5000 individuals)!

Hopefully someone can help me out with this.

Thanks very much,

Alice

[birdman]

The RMark documentation for the nest survival models in Appendix C helped me immensely in figuring out just this scenario. You have to extract the covariate values from your design matrix using find.covariates and then set up a loop which builds sequential matrices for generating your output at each time step and covariate value of interest. The loop uses fill.covariates to create the necessary matrices at each time step, and compute.real then calculates your output values and appends those to an output dataframe you've established. It then loops back to rebuild the loop for the next set of covariate values.

It is a bit cumbersome, but work through the nest survival examples and you should be able to figure it out if you are familiar with other regression approaches.

[bootzies]

Thank you for your quick response!!

However, I am still having some problems (stumbling at the first hurdle). I am trying to use the 'fill.covariates' function on my model, but get an error message:
test <- find.covariates(mod)
"Error in find.covariates(mod) :
This is an old model object. Re-run model or use the data argument must be specified and be a dataframe. Use data and not processed data list."

I have no idea what is meant by "this is an old model object"! And, it isn't clear to me from the function documentation what the 'data' argument should be (except that it shouldn't be a processed data object). I have tried using the dataframe containing my time-varying covariate data in the data argument. When I do this I don't get the error message anymore, but the command only produces an empty object (called 'test' in my case).

In the mallard example, the covariates are included in the initial ch data dataframe. Because I am using time varying (design) covariates in my models (as opposed to individual or grouping covariates), the covariate data were not included in the original ch dataframe, but instead I merged them into the design data using 'merge.design.covariates'. So, my ch data data frame doesn't contain any covariate data and is not of the same format as that used in the mallard example.

Any further tips would be great!

[jlaake]

I've been away and now in the field so I have limited time to answer this question. If the covariates are in your design data then you can use those values to plot against the reals directly but at present there is no good way to plot a different set of covariate values without using the coefficients and computing the values yourself. That is because the covariate values are in the design data and in the design matrix(dm) which is fixed for the model fit. I've not yet written code that lets you change those values.

If the covariates are in the data, I call these individual covariates because they can differ for each individual, then you should use covariate.predictions which replaces all of the need to use find.covariates and fill.covariates. Those functions find names of individual covariates in the dm and then replaces them with values. With design covariates there will not be any names in the dm.

There is discussion of this previously on the forum.

--jeff

[bacollier]

I've been away and now in the field so I have limited time to answer this question. If the covariates are in your design data then you can use those values to plot against the reals directly but at present there is no good way to plot a different set of covariate values without using the coefficients and computing the values yourself. That is because the covariate values are in the design data and in the design matrix(dm) which is fixed for the model fit. I've not yet written code that lets you change those values.

--jeff

Along the lines of what jeff said for using the coefficients, here is a simple example building a 3d plot in R where you have 2 continuous covariates you want to build predictions for the range of both. It would be pretty simple to do the same approach using the mean value for one covariate is that is your choice.

Code: Select all

library(rgl)
for(Days_Nesting in 1:55){
   Breeding_Period<-1:118
   DSR_RGWT= exp(6.6205449-0.0589883*Days_Nesting)/(1+
            exp(6.6205449-0.0589883*Days_Nesting))
   DSR_BS<-DSR_RGWT^Breeding_Period
   DSR_Final<-data.frame(cbind(Days_Nesting, Breeding_Period, DSR_BS))
   write(c(t(DSR_Final)), "location.txt", append=T, ncolumns=3)
}
dsr_data<-read.delim("location.txt", header=F, sep="")


Also, you might look at ?expand.grid in Rwhich I did not know about when I wrote the above but is much better and creating dataframes for all possible combinations.

\bret

[bootzies]

Hi Bret,

Thanks so much for the code. I am working with it now, but just have a couple of further questions. I have created a dataframe containing a range of values for one of my covariates (temperature) and a constant value for the other (rainfall). But I am struggling to interpret your code (completely my fault as I am a relative novice with R coding). I am unsure about what the names of your 2 covariates are in the example.

It looks like 6.6205449 is your intercept value and then 0.0589883 is the coefficient for one of your covariates (but is it one of the named variables in the loop)? I am not sure what your 2nd covariate of interest is in this code. Could you just describe what Days_Nesting and Breeding_Period are? And what this line does: DSR_BS<-DSR_RGWT^Breeding_Period?

Thank you very much - I really aprreciate your help with this.

[bootzies]

Jeff - thanks very much for your response to my message.

I am still struggling somewhat with this issue. I have 3 specific queries:

1) If I have more than one covariate in the model for survival, is it still acceptable to apply the coeffcient for one of my model design covariates to the intercept to calculate the effect of that covariate effect on survival? What covariate values does the intercept of S represent? Are all covariates ignored/set to 0/set to their mean value at the overall S intercept?

2) Following on from above, If I have used model averaging, is it ok to use the weights of the models to calculate a weighted average of the covariate coefficients and then apply this average to the model averaged intercept for S, to get an estimate of the effect of the covariate on survival?

3)

jlaake wrote: If the covariates are in your design data then you can use those values to plot against the reals directly

If I wanted to do as you suggest above; how can I extract the survival estimate for each time period if I haven't explicitly modelled S as a function of time?
My model: S ~ state : age group + winter temp + rain.
I am sure that there must be estimates of survival for each time period somewhere in the model output, but I am unsure how to access these (they are not given as coefs/real parameter estimates in the model output, because I did not have time as a parameter in the model formula).

Thanks very much. I know you are busy on fieldwork, so very much appreciate you taking the time to respond.

Alice

[jlaake]

1) If I have more than one covariate in the model for survival, is it still acceptable to apply the coeffcient for one of my model design covariates to the intercept to calculate the effect of that covariate effect on survival?

If I was doing this I would use the mean of all the other covariate values not being plotted or some other value that makes sense for your situation or you can show a set of lines at the quartiles of the other covariate value. Often 0 doesn't make sense (eg temperature data in the tropics) and because the values are transformed typically to 0-1 it will affect the real parameter value. Alternatively you can do as Bret suggested and display 2 variables with a 3-d plot.

1) What covariate values does the intercept of S represent? Are all covariates ignored/set to 0/set to their mean value at the overall S intercept?

The answer depends on how you treated your covariate data. The intercept is the value of 0 for any numeric covariate. If you centered the covariate by subtracting off the mean then 0 represents the mean covariate value. You can chose whatever values you want to use and make sense for your application.

2) Following on from above, If I have used model averaging, is it ok to use the weights of the models to calculate a weighted average of the covariate coefficients and then apply this average to the model averaged intercept for S, to get an estimate of the effect of the covariate on survival?

This question comes up frequently in regard to model averaging beta values (coefficients). It can work if you are very careful with the design matrices to avoid changing the meaning of coefficients but typically it is safer and more direct to compute the real parameter values and model average those.

If I wanted to do as you suggest above; how can I extract the survival estimate for each time period if I haven't explicitly modelled S as a function of time?
My model: S ~ state : age group + winter temp + rain.
I am sure that there must be estimates of survival for each time period somewhere in the model output, but I am unsure how to access these (they are not given as coefs/real parameter estimates in the model output, because I did not have time as a parameter in the model formula).

The functions that compute real parameter values like summary or get.real provide an estimate for each row in the design data with a 1-1 correspondence. So you can simply cbind the real parameter estimates to the design data which contain the design covariate values and you can use those to plot. For any one value of say winter temp, you'll may get a scatter of values because you can have different values of rain for the same winter temp.

I hope this helps. Eventually me or someone else we'll get this coded generally so it is easier to work with design covariates for prediction at some specified set of values. Your question has been asked before so you may want to read what is in the forum about this issue. It is possible that I have made other suggestions that I'm not remembering now.

regards --jeff