covariate prediction confusion

[ked]

Hi, I am having some trouble using the covariate.prediction function. I've run CJS models with several weather covariates for survival of different groups, and I'd like to predict survival rates under certain weather conditions. The problem I'm having is that the output using covariate.prediction is exactly the same as the original survival estimates over a number of years (through the whole range of weather conditions), rather than a single estimate for survival in the one particular set of weather conditions. I'm sure it's just a coding error somewhere, but I haven't figured it out yet.

Here is some (abbreviated) code, where A1, A, H, and L are the different groups, and djf1, avg, rain, and rain1 are different weather covariates:

model=function() {
Phi.L.djf1_A1_A1.avg_A_A.rain_H_H.rain1=list(formula=~L:djf1+A1+A1:avg+A+A:rain+H+H:rain1)
p.time=list(formula=~time)
cml=create.model.list("CJS")
results=mark.wrapper(cml,data=sosp.process,ddl=sosp.ddl,output=FALSE)
return(results)
}
results=model()

A1pims=c(1,30,58,85,111,136,160,183,205,226,246,265,283,300,316,331,345,358,370,381,391,400,408,415,421,426,430,433,435)

#prediction for A1 group, for avg=1.59
predict.A1=covariate.predictions(results,data=data.frame(avg=1.59),indices=A1pims)

The covariate is labelled the same as in the original model, and it doesn't seem to matter whether I use only one model or a marklist of models, or if I try the prediction for multiple groups/covariates at the same time.

I would appreciate any insights you might have. Thanks!

[jlaake]

covariate.predictions is for individual covariates and what you are using is an occasion covariate. An individual covariate can be different for each animal which is not the case for your covariates. --jeff

[ked]

Ah, I see. Thanks for the clarification. I didn't realize it was limited to individual covariates. So is there not a function for an occasion covariate? I suppose I can try writing one!

[jlaake]

I can understand your confusion as I don't use the word individual in describing the covariates but those are the only ones that can have a value plugged in as the other covariates have their values specified in the design data which is then translated to the design matrix where there values are fixed. Individual covariates on the other hand are represented by a covariate name rather than a value and this function plugs in specific values.

For design covariates like those for occasion or group you can specify new values in the design data, apply the model formula, multiply by the covariates and then transform. It might go something like this:

Phidf=sosp.ddl$Phi
Phidf$avg=1.59
beta=summary(results[[1]])$beta$estimate
dm=model.matrix(~L:djf1+A1+A1:avg+A+A:rain+H+H:rain1,Phidf)
plogis(dm%*%beta)

I'm doing this from memory and may have some synatx wrong but that should get you started. --jeff

[jlaake]

WIth regard to the code I sent, you'll need to restrict beta to the coefficients specific to Phi because what I'd would get all of the betas for Phi and p.

--jeff

[ked]

First, let me send a (very belated) thank you for this response. I'm revisiting this analysis, and I'm now wondering if there's an easy way to generate a prediction interval around these estimates?

[jlaake]

What I sent previously was an explanation of how it is computed. What you have is occasion/interval covariates and summary computes the estimates and confidence intervals for each time interval and each time interval is associated with a particular value of the covariate, so you can bind the design data with the dataframe from summary which you get by setting se=T.

Here is a simple example with dipper data.

data(dipper)
model=mark(dipper)
summary(model,se=T)$reals$Phi
cbind(summary(model,se=T)$reals$Phi,model$design.data$Phi)

In this case there are no new fields being added. Also you can use the R subset function to extract the fields you want from the design data to attach.

--jeff

[ked]

Ok, thanks, but I think this is just binding observed interval covariate values to survival estimates, right? What I was trying to do is predict survival under a specific value of the interval covariate (or combination of covariates), along with an associated prediction interval, similar to the way you can predict survival for a specific value of an individual covariate.

[jlaake]

Look at documentation for compute.real. You can provide a design matrix with the values you want substituted in and it will compute the reals with std errors and conf intervals. The function fill.covariates will not be useful because it is intended for individual covariates.

Here is a simple (but nonsensical) example with the dipper data, which fits a model with the 0/1 flood covariate and then predicts its values at 0.5 and 0.75 for flood:

Code: Select all

data(dipper)
dp=process.data(dipper,groups="sex")
ddl=make.design.data(dp)
ddl$Phi$Flood=0
ddl$Phi$Flood[ddl$Phi$time==2 | ddl$Phi$time==3]=1
mod=mark(dp,ddl,model.parameters=list(Phi=list(formula=~Flood)))
dm=mod$design.matrix
dm[1:2,2]=c(0.5,.75)
compute.real(mod,design=dm)


This uses the simplified design matrix (dm) which is model specific and could become a little confusing if you try to model average values across models because the rows in the dm aren't consistent across models. You'll just have to keep track of them.

A couple of alternatives. The simplest is to turn your occasion covariates into time-varying individual covariates. The values will be the same for each animal but will vary across time. That will let you use covariate.predictions which does the model averaging for you. For my example, above the covariates would be named Flood1,...,Flood6 for survival because I let begin.time default to 1.

The second is to use model.matrix with each parameter in the ddl (eg Phi and p), construct the complete dm and then use that with compute.real. That will use the all different PIM so each real parameter is the same in all models and model averaging is straightforward with model.average. Below is an example with a single model.

Code: Select all

phidm=model.matrix(~Flood,ddl$Phi)
pdm=model.matrix(~1,ddl$p)
dm=matrix(0,nrow=nrow(phidm)+nrow(pdm),ncol=ncol(phidm)+ncol(pdm))
dm[1:nrow(phidm),1:ncol(phidm)]=phidm
dm[(nrow(phidm)+1):nrow(dm),(ncol(phidm)+1):ncol(dm)]=pdm
compute.real(mod,design=dm)


That will compute the real values at the observed values of Flood at each time and you would have to change them as in:

Code: Select all

dm[1:42,2]=0.5
compute.real(mod,design=dm)


[gstauffer]

Jeff,
This make sense, I think. But, I'm not sure how to take the next step and get model-averaged estimates for specified occasion-specific covariate values. Is that possible?

Glenn