Selecting appropriate indices for time varying covariates

[rcscott]

Hi all,

I am running RDMultScal models to determine factors impacting prevalence and detection of Bd and ranavirus in larval amphibians. In my top model, there is an effect of larvae species and temperature of the pond on detection of the disease. In this case, temperature is a time varying covariate while species is coded in to the .ddl as 0 or 1 for both theta (prevalence) and p as a parameter index(?).

I have run in to a snag when using covariate.predictions(). I wanted to use covariate.predictions() to estimate p for each species at a range of temperatures. So I orginially used

Code: Select all

pXTemp <- covariate.predictions(psi1eps1gamma1$Psi.null.Epsilon.year.Gamma.null.Theta.SppXCond.p.SppXTemp,
                                data = temp_data,
                                indices = c(534,2093))


Where indices 534 and 2093 should correspond to species = 0 and 1 respectively. However, when I plotted that, one species had a straight line even though the betas indicate it should have a negative slope. The only difference I realized was in the session: model.index 534 had session = 1 while model.index 2093 had session = 5. When I change the indices so that both are from session 1, then the estimates from covariate.predictions() match what I expect.

I'm confused because according to my supervisor as long as we choose any 2 indices where species = 0 and 1, then it should give the correct output, but it doesn't. Are we misunderstanding how indices work in a dynamic model?

Thanks!
Reed

[jlaake]

Depends in how your model for p was specified. If it only included temperature then I would agree with your supervisor. But you need to make sure you are using the correct indices. It sounds like one of the indices may be from a parameter that doesn't include a temperature effect. I can't tell without seeing the ddl data. Please check.

[rcscott]

Hi Jeff,

Thanks for the quick reply! I specified the model for p as:

p.SppXTemp <- list(formula =~ pSpecies*temp)

So I was looking for an interaction between the Species X temperature. I've checked the indices and am positive that I'm calling the right ones, but I'll also paste them here:

Code: Select all

> head(RV.ddl$p)
par.index model.index group time session Time primary pSpecies
1         1         534     1    1       1    0       1        0
2         2         535     1    2       1    1       1        0
3         3         536     1    3       1    2       1        0
4         4         537     1    4       1    3       2        0
5         5         538     1    5       1    4       2        0
6         6         539     1    6       1    5       2        0

> tail(RV.ddl$p)
par.index model.index group time session Time primary pSpecies
1555      1555        2088     1  307       5  306     103        1
1556      1556        2089     1  308       5  307     103        1
1557      1557        2090     1  309       5  308     103        1
1558      1558        2091     1  310       5  309     104        1
1559      1559        2092     1  311       5  310     104        1
1560      1560        2093     1  312       5  311     104        1


So model.index = 534 should have pSpecies = 0 and model.index = 2093 should have pSpecies = 1. Despite this, covariate.predictions() doesn't work with these indices. It does work however if I set indices = c(534,690), where 690 is the 1st case when pSpecies = 1.

[jlaake]

Please send me a .rdata file containing processed data, ddl and model you are working on to jefflaake@gmail.com. That is the only way I can debug the problem. In the meantime it appears you have a solution.

[jlaake]

I worked with Reed on this offlist. I missed that this was a time varying individual covariate which are contained in the data with the capture history and for his situation were named temp1,...temp6 for the 6 sessions in the robust design. As it turns out, the data he supplied to covariate.predictions only contained the variable temp1 and that is sufficient if you only use indices for session1. When he used an index for session 6, it gave constant predictions because it did not find temp6 in the data argument. When covariate.predictions doesn't find a covariate in the data argument that is used in the design matrix, it uses the average covariate value for the prediction and thus the predictions will be constant for the values used in data.

Below I create some dummy time varying individual covariates for Phi in the dipper data to demonstrate.

Code: Select all

library(RMark)
data(dipper)
dipper$temp1=runif(294,0,1)
dipper$temp2=runif(294,0,1)
dipper$temp3=runif(294,0,1)
dipper$temp4=runif(294,0,1)
dipper$temp5=runif(294,0,1)
dipper$temp6=runif(294,0,1)
dp=process.data(dipper,model="CJS")
ddl=make.design.data(dp)
model=mark(dp,ddl,model.parameters = list(Phi=list(formula=~temp)))


If you look at the design matrix you will see that the entries are temp1,...temp6 for the 6 possible Phi values in the simplified model structure.

Code: Select all

> model$design.matrix
                Phi:(Intercept) Phi:temp p:(Intercept)
Phi g1 c1 a0 t1 "1"             "temp1"  "0"         
Phi g1 c1 a1 t2 "1"             "temp2"  "0"         
Phi g1 c1 a2 t3 "1"             "temp3"  "0"         
Phi g1 c1 a3 t4 "1"             "temp4"  "0"         
Phi g1 c1 a4 t5 "1"             "temp5"  "0"         
Phi g1 c1 a5 t6 "1"             "temp6"  "0"         
p g1 c1 a1 t2   "0"             "0"      "1"         


If you pick a model.index from ddl.Phi that corresponds to time 1 it will use the first row in the simplified structure which contains temp1 and if you use a model.index for time 6 it will the 6th row containing temp6. If you were to choose those, you would need to include values for temp1 and temp6. Below I show what happens when you do it incorrectly and then correctly.

Code: Select all

> covariate.predictions(model,data=data.frame(temp1=c(0,.5,1)),indices=c(1,6))
$estimates
  vcv.index model.index par.index covdata  estimate         se       lcl       ucl fixed
1         1           1         1     0.0 0.5082236 0.05316857 0.4051506 0.6106021     
2         2           6         6     0.0 0.5585437 0.02523609 0.5086633 0.6072700     
3         3           1         1     0.5 0.5599992 0.02519812 0.5101794 0.6086390     
4         4           6         6     0.5 0.5585437 0.02523609 0.5086633 0.6072700     
5         5           1         1     1.0 0.6105015 0.05103428 0.5071930 0.7047597     
6         6           6         6     1.0 0.5585437 0.02523609 0.5086633 0.6072700     

$vcv
              [,1]         [,2]         [,3]         [,4]          [,5]         [,6]
[1,]  0.0028268967 0.0006932227 0.0006316270 0.0006932227 -0.0014704436 0.0006932227
[2,]  0.0006932227 0.0006368601 0.0006350426 0.0006368601  0.0005661585 0.0006368601
[3,]  0.0006316270 0.0006350426 0.0006349454 0.0006350426  0.0006245740 0.0006350426
[4,]  0.0006932227 0.0006368601 0.0006350426 0.0006368601  0.0005661585 0.0006368601
[5,] -0.0014704436 0.0005661585 0.0006245740 0.0005661585  0.0026044975 0.0005661585
[6,]  0.0006932227 0.0006368601 0.0006350426 0.0006368601  0.0005661585 0.0006368601

> covariate.predictions(model,data=data.frame(temp1=c(0,.5,1),temp6=c(0,0.5,1)),indices=c(1,6))
$estimates
  vcv.index model.index par.index temp1 temp6  estimate         se       lcl       ucl fixed
1         1           1         1   0.0   0.0 0.5082236 0.05316857 0.4051506 0.6106021     
2         2           6         6   0.0   0.0 0.5082236 0.05316857 0.4051506 0.6106021     
3         3           1         1   0.5   0.5 0.5599992 0.02519812 0.5101794 0.6086390     
4         4           6         6   0.5   0.5 0.5599992 0.02519812 0.5101794 0.6086390     
5         5           1         1   1.0   1.0 0.6105015 0.05103428 0.5071930 0.7047597     
6         6           6         6   1.0   1.0 0.6105015 0.05103428 0.5071930 0.7047597     

$vcv
             [,1]         [,2]         [,3]         [,4]         [,5]         [,6]
[1,]  0.002826897  0.002826897 0.0006316270 0.0006316270 -0.001470444 -0.001470444
[2,]  0.002826897  0.002826897 0.0006316270 0.0006316270 -0.001470444 -0.001470444
[3,]  0.000631627  0.000631627 0.0006349454 0.0006349454  0.000624574  0.000624574
[4,]  0.000631627  0.000631627 0.0006349454 0.0006349454  0.000624574  0.000624574
[5,] -0.001470444 -0.001470444 0.0006245740 0.0006245740  0.002604498  0.002604498
[6,] -0.001470444 -0.001470444 0.0006245740 0.0006245740  0.002604498  0.002604498


Now in this dummy example they are all the same and you wouldn't choose a model.index for both 1 and 6 but in more complex examples you might have a separate intercept per time or a separate slope per time, in which case you would need to specify all of the indices and a range for all of the tempx values for x=1,...,6.

[rcscott]

Thanks again Jeff!