www.phidot.org

I'm super new to RMark. I have a CJS model with p(.) and phi that depends on age at release (decimal years), weight at release, and their interaction. I want to predict the survival at a range of weights at 2-3 ages (that differ from the ages used for the model). The manuals seem to provide examples with a single continuous variable, where PIMS are used for intercept/categorical variable states.

My release ages are 0.83 and 1.08, but say I want to calculate survival at age 0.5, 0.7, and 1, across weights 10:300. How do I get that using covariate.predictions?

Thank you in advance for your patience....


try <- structure(list(ch = c("10000000000", "00001100000", "00010010010",
"00010000010", "00010100000", "01000000101", "10000000000", "00000101100",
"10000000000", "10000000000", "10000000000", "00001100000", "10000000000",
"10000000000", "10000000000", "10000000000", "10000000000", "00010100000",
"10000000000", "10000000000", "10000000000", "10000000000", "10000000000",
"10000000000", "00010100000", "10000000000", "10000000000", "10000000000",
"00010000010", "10000000000", "00010000010", "10000000000", "10000000000",
"00000011010", "10000000000", "10000000000", "10000000000", "00010000010",
"00010100000", "10000000000", "00000101100", "10000000000", "00010011010",
"10000000000", "10000000000", "00001100000", "00101000100", "00001100000",
"00010000010", "00001100000"), Release.Weight = c(123.1, 49.2,
37.1, 36.9, 121, 57.1, 95.7, 113.4, 96, 91.1, 34.3, 74.8, 36.7,
97.4, 48.9, 122.9, 69, 67.2, 137, 52.2, 43.5, 111, 120.4, 137.9,
63.3, 70.7, 133.5, 32.4, 86.2, 113, 19.7, 163.6, 50.4, 47.9,
32.8, 38.8, 35.8, 84, 86.6, 57.1, 45.8, 103.9, 41.2, 96.5, 53.7,
64, 43, 78, 55, 34), Release.Age = c(1.08, 0.83, 0.83, 0.83,
1.08, 0.83, 1.08, 1.08, 1.08, 0.83, 0.83, 0.83, 0.83, 1.08, 0.83,
1.08, 0.83, 0.83, 1.08, 0.83, 0.83, 1.08, 1.08, 1.08, 0.83, 1.08,
1.08, 0.83, 0.83, 1.08, 0.83, 1.08, 0.83, 0.83, 0.83, 0.83, 0.83,
0.83, 0.83, 0.83, 0.83, 1.08, 0.83, 1.08, 0.83, 0.83, 0.83, 0.83,
0.83, 0.83)), .Names = c("ch", "Release.Weight", "Release.Age"
), row.names = c(64L, 81L, 82L, 102L, 122L, 160L, 179L, 194L,
195L, 205L, 211L, 222L, 229L, 251L, 279L, 283L, 320L, 355L, 371L,
410L, 427L, 542L, 548L, 645L, 649L, 686L, 727L, 760L, 782L, 799L,
809L, 828L, 831L, 859L, 882L, 918L, 950L, 954L, 973L, 988L, 1030L,
1102L, 1105L, 1176L, 1187L, 1225L, 1242L, 1417L, 1419L, 1740L
), class = "data.frame")

try.proc <- process.data(try, model="CJS", begin.time=1990)
try.ddl <- make.design.data(try.proc)

try.analysis=function(){
Phi.weight.x.age.release=list(formula = ~Release.Weight*Release.Age)
p.1=list(formula=~1)
cml=create.model.list("CJS")
mark.wrapper(cml,data=try.proc,ddl=try.ddl,output=FALSE)
}

try.results <- try.analysis()

covariate.predictions only works on individual covariates. If you are using Age in the design data, it will not work. I have not gotten to coding in design covariates like Age into covariate.predictions. To do what you want, you will have to create a series of variables for age at each occasion. Look at documentation for time varying individual covariates. Then you can use covariate.predictions with both. With the way you have it setup you can only get predictions at the ages you observed but you can use it with a range of release.weights because that was entered as an individual covariate. I guess I don't see the problem with using the observed ages because it will provide a nice range if you have sufficient number of occasions.

--jeff

Thank you for the super-quick reply. Age at release is currently set up as an individual covariate - it does not update with time. Does that change your answer? As for the ages - I just wanted more points within the (min, max) range of ages. And it would be weird to report survival at 0.83 years ("well, that certainly seems oddly specific...")

Thank you (again)

Your formula shows Age and not Release.Age. If it as you say then you just need to create a dataframe with the combination of release weights and ages that you want to get computations for and specify the parameter index for which they should be computed. Note that the values can vary depending on what you have in the model (ie time effect). You may want to look at expand.grid in R to create combinations of 2 sets of values.

--jeff

Yes, mistake in the model copying, it should be Release.Age. I edited the original post, thank you for catching that. There are no time effects, model specification is just like in the post...

What you're suggesting is exactly what I'm after, I just fail on figuring out the parameter index. I can understand how the indices work for non-continuous params, but I just don't get how to figure out the index in this (I'm sure, simple) case.

Say, weight/age data for prediction are:

pred.data <- expand.grid(Release.Age = c(0.5, 0.75, 1), Release.Weight = seq(10, 300, 10))
preds <- covariate.predictions(try.results,
data = data.frame(Release.Weight = pred.data$Release.Weight, Release.Age = pred.data$Release.Age), indices = ?????)

When I experimented and entered different indices (within the range of my PIMS), it got me the graphs I wanted, but if I understand correctly, these are for modelled ages, not the ages I want to predict at. Also, at any given index value, I got identical graphs across the three ages I used.

I'm just missing something super-basic here.

Here is an example with the dipper data where I create two dummy covariates wt and length.
Here the indices could be 1 to 21 because there are 21 possible Phi estimates. In this case the estimates are all the same because there is no time, age, cohort or other variation, so I could have just chosen 1. If I passed covariate.predictions a marklist of models it would compute the predictions for each model and then model average the predictions. The se and v-c matrix include model uncertainty in that case. You have to use all different PIMS because the smplified PIMS differ between models. --jeff

Code: Select all

data(dipper)
dipper$wt=rnorm(294)
dipper$length=rnorm(294)
mod=mark(dipper,model.parameters=list(Phi=list(formula=~wt+length)))
PIMS(mod,"Phi",simplified=F)
df=expand.grid(wt=-3:3,length=-3:3)
pred=covariate.predictions(mod,df,indices=1:5)
pred$estimates


Jeff (or anyone else) - I have a follow-up question regarding this thread that I would love to get feedback on.

As Jeff's said in multiple posts, covariate.predictions only works with individual covariates and not with design covariates. I'm trying to predict capture probability (in RDHuggins) for an interaction between a factor with 3 levels (currently a group in the design data) and an individual covariate (age).

Can I use covariate.predictions if I re-code the 3 level factor as an individual covariate (2 columns of dummy variables with a reference level)? I would then need to specify values for all 3 individual covariates (2 fields representing the factor and the 1 age field) in the "data" argument of covariate.predictions?

Thanks!
Joe

Joe-

You are missing the point with regard to individual and design covariates. The folks that asked that question wanted to predict at values of the design covariates that were not used in fitting the model. In the case you are describing of using a factor variable for grouping, there is a parameter for each group and you can make predictions of values for the individual covariate for each parameter - thus each level of the factor variable.

--jeff

Thanks Jeff.
I think I follow...cov.preds. can be used with a mix of design and individual covariates but the "predictions" for the design covariate are restricted to observed values for that covariate.

Any insights on this error? Do I need to specify the values for the covariates differently?

Code: Select all

 # make dataframe for cov.pred
df = expand.grid(ad_ave = c(0,1), Tpref = c("rand","hav", "sherm"))

covariate.predictions(model, df, indices = 484:1173) # indices for "p"

Warning message:
In fixedvalues[fixedparms] = model$fixed$value[match(indices[fixedparms],  :
  number of items to replace is not a multiple of replacement length

The function still runs and returns predicted values, but the values are not what I would expect (they are the same for every level of the factor).

Thanks again.
Joe

Hard to tell. Why are you putting factor covariates in the data frame. If that is the variable used for the groups then it is not needed because some of the indices (parameters) are for each group. Did the model you run have different values of p for each group?

--jeff