Predicted Nest Survival at Specified Covariate Values

[mcallen]

Hi,
Can anyone provide an example of how to generate predicted Nest Survival values in RMark for varying levels of a continuous covariate. There is a nice example in the RMark appendix of the MARK book of how to do this with varying levels of "Nest Age" and a continuous covariate, but I've had no luck applying the example to models that just include a continuous variable (not Nest Age). Here's what I've tried, but I'm pretty sure I have the structure of the design matrix all wrong. Any help would be appreciated.

Code: Select all

grsp.veg.height=mark(grsp,nocc=95,model="Nest",model.parameters=list(S=list(formula=~veg.height)))

nocc=95

design=matrix(cbind(rep(1,(nocc-1)),seq(min(grsp$veg.height),max(grsp$veg.height),length=nocc-1)),nocc-1,2)

get.real(grsp$veg.height,parameter="S",beta=grsp.vegheight$results$beta$estimate,se=T,
   design=design)


At this point, I should probably just use regular MARK, but I like working in R!

Thanks,
Mike

[jlaake]

The mallard example was written before I created covariate.predictions which was designed to do what you want. You should be able to pass with arguments data=data.frame(veg.height=...) and indices=1, where ... is the set of values you want. I believe for your example that is correct but I may have overlooked something and if it doesn't work send me the data and R code and I'll look at it.

--jeff

[mcallen]

Jeff,
I'll give covariate.predictions a try with the arguments you suggested. If I'm still stuck I'll send along the data.

Thanks,
Mike

[mcallen]

Hi Jeff,
That worked great - thanks! But one last (related) problem.

I'm now trying to graph both Year (2 dummy variables representing 2009, 2010, and 2012) AND a continuous covariate (% bare ground) simultaneously. It works well predicting DSR at various levels of the continuous variable ("bare"), but is ignoring my arguments for Year2009 (0 or 1) and Year2010 (0 or 1). It just returns the DSR as if both were 0, no matter what I put in. Here is the model and prediction code:

Code: Select all

#model
gp1year.bare=mark(gp1,nocc=103,model="Nest",model.parameters=list(S=list(formula=~year2009+year2010+bare)), groups=c("year2009","year2010"))

#beta estimates
                estimate        se        lcl       ucl
S:(Intercept)  3.0221740 0.5396361  1.9644872 4.0798608
S:year20091    0.3838193 0.6793767 -0.9477590 1.7153977
S:year20101    0.6545762 0.3754082 -0.0812239 1.3903763
S:bare        -1.2772674 0.9206077 -3.0816586 0.5271237

#covariate predictions
covariate.predictions(gp1year.bare,data=data.frame(bare=seq(min(gp1$bare),max(gp1$bare),length=10),year2009=rep(1,10),year2010=rep(0,10)),indices=1)

Again, no matter what value of Year2009 or Year2010 I put into covariate.predictions, it always returns the same set of dsr:

$estimates
vcv.index model.index par.index bare year2009 year2010 estimate se lcl ucl fixed
1 1 1 1 0.0000000 1 0 0.9535659 0.02389400 0.8770178 0.9833714
2 2 1 1 0.1230165 1 0 0.9460914 0.02231640 0.8815761 0.9764007
3 3 1 1 0.2460330 1 0 0.9374926 0.02007459 0.8845753 0.9670532
4 4 1 1 0.3690496 1 0 0.9276272 0.01762251 0.8845571 0.9554379
5 5 1 1 0.4920661 1 0 0.9163436 0.01645984 0.8779144 0.9434555
6 6 1 1 0.6150826 1 0 0.9034839 0.01929397 0.8584993 0.9352462
7 7 1 1 0.7380991 1 0 0.8888872 0.02748818 0.8225804 0.9324483
8 8 1 1 0.8611157 1 0 0.8723946 0.04038604 0.7705184 0.9329778
9 9 1 1 0.9841322 1 0 0.8538565 0.05737670 0.7034877 0.9350138
10 10 1 1 1.1071487 1 0 0.8331403 0.07825519 0.6235711 0.9376939

Any help is appreciated!

Thanks,
Mike

[jlaake]

Factor variables are handled with different real parameters. So you have to use the indices for those years to get different real estimates for them. Currently I suggested indices=1 because I didn't know you had differences for years. You can get those values from the design data or using the PIMS function with the simplified=FALSE argument.

--jeff

[mcallen]

Hi Jeff,
Thanks. I switched datasets on you and am kind of playing around as a learning experience. Below is what I got from the PIMS function. I'm still not clear on how to use that info to define the indices for each year. Is it something like this (since the index values for 2009 seem to be 103-204 based on PIMS)?:

Code: Select all

covariate.predictions(gp1year.bare,data=data.frame(bare=seq(min(gp1$bare),max(gp1$bare),length=10),year2009=rep(1,10),year2010=rep(0,10)),indices=103:204)

Sorry if this is too basic. I still haven't totally wrapped my head around how MARK handles data, design matrices, etc. behind the scenes, and have been heavily reliant on the example code.

Thanks again.
Mike

Code: Select all

> PIMS(gp1year.bare,simplified=F,parameter="S")
group = year20090.year20100;stratum = 1
    1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19
1   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19
   20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38
1  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38
   39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57
1  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57
   58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76
1  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76
   77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95
1  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95
   96  97  98  99 100 101 102
1  96  97  98  99 100 101 102
group = year20091.year20100;stratum = 1
    1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19
1 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
   20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38
1 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
   39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57
1 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
   58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76
1 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178
   77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95
1 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197
   96  97  98  99 100 101 102
1 198 199 200 201 202 203 204
group = year20090.year20101;stratum = 1
    1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19
1 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
   20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38
1 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242
   39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57
1 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261
   58  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76
1 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
   77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95
1 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299
   96  97  98  99 100 101 102
1 300 301 302 303 304 305 306

[jlaake]

No worries it is confusing. If 1 is an index to 2007 and 103 is an index to 2009, then you only need indices=c(1,103) to get the values for 2007 and 2009. You don't need to specify each real parameter unless you had differences based on the values of those real parameters.

--jeff

[mcallen]

Thanks, that worked perfectly.