Model averaging in RMark with multiple grouping variables

[kellyweintraub]

Hello,

I am modeling survival using Known Fate models in RMark and I need to model average my results. I have 4 factor variables and 5 covariates. My factors are site, habitat type, year, and population level. (I did not include site and habitat in the same model.)

I have 11 sites, 2 years, and 4 population levels. Some of my sites were only active in 1 of the 2 years. My covariates are nest height, water depth, nest density, day of season, and season-squared.

I have two questions:
1) What command do I use (and what would the code look like) to model average with more than one grouping variable? I know I can use covariate.predictions to model average with multiple covariates, but can I use this same command to model average with multiple factors? When I tried to use it I did not get the results I expected.

2) Is it possible to obtain one survival estimate for each site in each year rather than getting an estimate for each site at each population level in each year?

Here is an example of some of my code:

Code: Select all

trbl=read.table("TRBL.txt",header=T,colClasses=c("character","factor","factor","numeric","numeric","numeric","numeric","numeric","factor","factor"))

#Create process data
trbl.proc=process.data(trbl,model="Known",groups=c("site","habitat","population","year"))

#Create design data
trbl.ddl=make.design.data(trbl.proc)

do.analysis=function()
{
S.site.season.season2.height.water.nests.pop.year=list(formula=~site+season+season2+height+water+population+nests+year)
S.site.season.season2.height.nests.water.year=list(formula=~site+season+season2+height+nests+water+year)
S.site.season.season2.height.nests.pop.year=list(formula=~site+season+season2+height+nests+population+year)
S.habitat.season.season2.height.nests.water.year=list(formula=~habitat+season+season2+height+nests+water+year)
S.habitat.season.season2.height.nests.pop.year=list(formula=~habitat+season+season2+height+nests+population+year)
cml=create.model.list("Known")

model.list=mark.wrapper(cml,data=trbl.proc,ddl=trbl.ddl)

return(model.list)
}

myresults=do.analysis()

myresults


I ran more models than this, but this should give you an idea of what the analysis looked like.

Thank you very much!

[kellyweintraub]

I believe I have figured out a way to do this. Thank you.

[jlaake]

Was a little busy last week so I wasn't able to answer your question right away. Factor variables differ from other numeric individual covariates in that there is a real parameter for each occasion for each level of the factor variable (eg for 11 sites there will be 11 real parameters for each occasion) and it is the real parameters that are model averaged with model.average. But since you have individual covariates, you want to use covariate.predictions which averages the real parameters at one or more values of the individual covariates. So the answer to question 1 is to use covariate.predictions. For question 2, I'm not sure what you mean by getting just one estimate for a site rather than a site-population level. If a site is always at single population level then just select one of the 4 that are created for each site. If however, sites have different population levels then you have to define what it is you want to estimate (eg mean of values at different population levels at a site) and that will require additional computation.

--jeff

[kellyweintraub]

Jeff,

Thank you for your reply. I thought I had figured out how to do this but I was wrong. I would really appreciate some help in writing the correct code for the covariate.predictions fuction. My first problem is determining which indices to specify. Below I use one of my sites as an example. Each site can, and usually does, have more than one population level as the season progresses, but in my example I chose a site that only had one population level in each year of the study.

This is what my design data looks like for the Burbank site in 2011:

Code: Select all

> trbl.ddl
$S
par.index model.index               group age time Age Time      site  habitat pop year
       1           1   BURBANKCATTAIL12011  1    2   1    0   BURBANK  CATTAIL   1 2011
       2           2   BURBANKCATTAIL12011  2    3   2    1   BURBANK  CATTAIL   1 2011
       3           3   BURBANKCATTAIL12011  3    4   3    2   BURBANK  CATTAIL   1 2011
       4           4   BURBANKCATTAIL12011  4    5   4    3   BURBANK  CATTAIL   1 2011
       5           5   BURBANKCATTAIL12011  5    6   5    4   BURBANK  CATTAIL   1 2011
       6           6   BURBANKCATTAIL12011  6    7   6    5   BURBANK  CATTAIL   1 2011
       7           7   BURBANKCATTAIL12011  7    8   7    6   BURBANK  CATTAIL   1 2011
       8           8   BURBANKCATTAIL12011  8    9   8    7   BURBANK  CATTAIL   1 2011
       9           9   BURBANKCATTAIL12011  9    10  9    8   BURBANK  CATTAIL   1 2011
      10          10   BURBANKCATTAIL12011  10   11  10   9   BURBANK  CATTAIL   1 2011
      11          11   BURBANKCATTAIL12011  11   12  11   10  BURBANK  CATTAIL   1 2011
      12          12   BURBANKCATTAIL12011  12   13  12   11  BURBANK  CATTAIL   1 2011
      13          13   BURBANKCATTAIL12011  13   14  13   12  BURBANK  CATTAIL   1 2011
      14          14   BURBANKCATTAIL12011  14   15  14   13  BURBANK  CATTAIL   1 2011
      15          15   BURBANKCATTAIL12011  15   16  15   14  BURBANK  CATTAIL   1 2011
      16          16   BURBANKCATTAIL12011  16   17  16   15  BURBANK  CATTAIL   1 2011
      17          17   BURBANKCATTAIL12011  17   18  17   16  BURBANK  CATTAIL   1 2011
      18          18   BURBANKCATTAIL12011  18   19  18   17  BURBANK  CATTAIL   1 2011
      19          19   BURBANKCATTAIL12011  19   20  19   18  BURBANK  CATTAIL   1 2011
      20          20   BURBANKCATTAIL12011  20   21  20   19  BURBANK  CATTAIL   1 2011
      21          21   BURBANKCATTAIL12011  21   22  21   20  BURBANK  CATTAIL   1 2011
...


The Burbank site was active for 2 years and had 1 population level in each year. I had 21 intervals in my study. When I look up how many parameters correspond to the Burbank site I find 42 indices for the site, one for each interval in each year:

Code: Select all

> as.numeric(row.names(trbl.ddl$S[trbl.ddl$S$site=="BURBANK",]))
 [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292
[41] 293 294


But when I ask R for the real parameter estimates from one of my models I get only 2 survival estimates for Burbank (one for each year because this site had only one population level in each year):

Code: Select all

 
myresults[[133]]$results$real
                                  estimate        se       lcl       ucl fixed    note
S gBURBANKCATTAIL12011 a1 t2      0.9420701 0.0062921 0.9284332 0.9532393             
S gECLACATTAIL12011 a1 t2         0.9754643 0.0040203 0.9662209 0.9822248             
S gUNIT1CATTAIL12011 a1 t2        0.7658776 0.0638431 0.6194746 0.8679596             
S gTOLEDOCATTAIL22011 a1 t2       0.9181596 0.0112330 0.8932714 0.9376492             
S gUNIT1CATTAIL22011 a1 t2        0.6988002 0.0921607 0.4958451 0.8455088             
S gBAKERSFIELDMESQUITE22011 a1 t2 0.9821960 0.0035103 0.9738376 0.9879172             
S gHACIENDASALT_CEDAR22011 a1 t2  0.9254179 0.0164141 0.8861720 0.9518674             
S gCostaA5SILAGE22011 a1 t2       0.9864633 0.0030187 0.9790747 0.9912663             
S gCostaA7SILAGE22011 a1 t2       0.9805017 0.0038133 0.9714406 0.9867272             
S gECLACATTAIL32011 a1 t2         0.9341202 0.0091114 0.9138504 0.9498824             
S gBAKERSFIELDMESQUITE32011 a1 t2 0.9652076 0.0083744 0.9444948 0.9783681             
S gCostaA5SILAGE32011 a1 t2       0.9734366 0.0071141 0.9553000 0.9843351             
S gCostaA5SILAGE42011 a1 t2       0.9361968 0.0168624 0.8940426 0.9622878             
S gBURBANKCATTAIL12012 a1 t2      0.9580376 0.0084968 0.9378393 0.9718696             
S gECLACATTAIL12012 a1 t2         0.9823992 0.0039238 0.9728082 0.9886467             
S gSEMITROPICCATTAIL12012 a1 t2   0.8776591 0.0212416 0.8295868 0.9135832             
S gTOLEDOCATTAIL12012 a1 t2       0.9569115 0.0053300 0.9451679 0.9662298             
S gSEMITROPICCATTAIL22012 a1 t2   0.8357390 0.0314714 0.7645029 0.8885675             
S gTOLEDOCATTAIL22012 a1 t2       0.9403001 0.0066672 0.9258128 0.9521046             
S gECLACATTAIL32012 a1 t2         0.9521677 0.0047402 0.9419719 0.9606469             
S gTOLEDOCATTAIL32012 a1 t2       0.8878978 0.0156983 0.8532500 0.9151783             
S gCORNERSTONESILAGE32012 a1 t2   0.9679120 0.0052766 0.9557934 0.9767892             
S gRIVERVIEWSILAGE32012 a1 t2     0.9957073 0.0011157 0.9928604 0.9974220             
S gECLACATTAIL42012 a1 t2         0.8885252 0.0165193 0.8518108 0.9170294             
S gSEMITROPICCATTAIL42012 a1 t2   0.5060397 0.0724799 0.3672191 0.6439353             
S gTOLEDOCATTAIL42012 a1 t2       0.7602715 0.0280162 0.7011717 0.8108363             
S gRIVERVIEWSILAGE42012 a1 t2     0.9893477 0.0024471 0.9833107 0.9932160     


If I understand correctly, the "real" estimates above are the probability of surviving from one interval to the next (the daily survival rate).

Which indices would I specify to use covariate.predictions to get the daily survival rate for the Burbank site?

Thank you very much.