Occasion and vegetation specific estimates of g0 and sigma

[mjgould]

Hey All,

I'm comparing density estimates between two different areas - NSC and SSC. I conducted 4 sampling occasions and used a categorical vegetation type covariate with 9 factors. My top supported model for the NSC = D~1, g0~t:veg, sigma~t:veg, while the top model for the SSC = D~1, g0~t+veg, sigma~t+veg. I would like to see what the occasion-specific estimates of g0 and sigma are for each vegetation type. When the two covariates are in independent models I simply use the following:
predict(NSC.V, newdata = data.frame(veg=factor(1:9)))
predict(NSC.Ti, newdata = data.frame(t=factor(1:4))).

Can I call up the g0 estimate for a specific vegetation type in each sampling occasion? Any help would be great!

[murray.efford]

I think you just provide more complex newdata with both predictors:

Code: Select all

predict(NSC.V, newdata = expand.grid(veg = factor(1:9), t = 1:4))

(I find it a bit scary that you have so much going on in the detection model, but presumably the data are strong enough to support that.)
Murray

[mjgould]

Thanks for the help, Murray. I really appreciate it. The code just needed a simple tweak [predict(NSC.V, newdata = expand.grid(veg = factor(1:9), t = factor(1:4)))], but otherwise it worked great.

Also, thanks for your input on detection model. Here's what the two are were comprised of:

NSC
N animals : 379
N detections : 455

SSC
N animals : 130
N detections : 154

I used an a priori stepwise modeling approach for each parameter. I held g0 constant while modeling sigma and vice versa. I then combined the top model from g0 and sigma into one model and fitted.

[murray.efford]

The code just needed a simple tweak

You're so discrete! I stuffed up (~t implies occasion is categorical, so newdata needs t=factor(1:4)), but I'm glad you could recover the essence.

On the question of model complexity - I don't have the answer on how best to navigate the model maze, but I am still a bit suspicious of where you ended up. The good news is that slight variations in the detection model often have little effect on density estimates if that's what you care about.
Murray

[mjgould]

That's great to hear. The defensibility of the density estimates is definitely the most important aspect of the project. Just to add to the discussion, here a few number from the AIC results:

NSC
g0
Model AICwt
t:veg 1 0.89
t + veg 0.11

sigma
Model AICwt
t:veg 1.00
t + veg 0.00

SSC
g0
Model AICwt
t + veg 0.99
t 0.01

sigma
Model AICwt
t + veg 0.94
t 0.06

[murray.efford]

Unfortunately, the magnitude of AIC differences bears little relation to differences in density estimates. I predict that the variation in density estimates among these models is much less than the sampling error of density, but I may be wrong.
Murray

[mjgould]

Once again, thanks for the insight, Murray. A quick follow up question. After I've extracted the veg and time specific parameter estimates, I'm presented with this type of heading: "$`veg = 1, t = 1". I'm not surprised that "veg" was coded/labeled as a number, however, I'm having trouble matching that number to the factor level name. I've tried the usual investigatory arguments, e.g. levels(), but I can't seem to get anything to work. Additionally, I've also tried using these arguments on just the capture history file. I feel like there's a simple solution that I'm not thinking of?

Best,

Matt

[murray.efford]

Hi Matt
I inferred from your original post that veg was coded 1:9. If the levels of veg were something else then you will want to provide these as the factor levels in newdata (veg = factor('A','B',,)). It's probably important to match the order of levels. If 'trps' is a traps object with a 'veg' covariate, this should be safe

Code: Select all

..., newdata = data.frame(expand.grid(veg = factor(levels(covariates(trps)$veg)), t=factor(1:5))), ...

Or maybe I'm not understanding the problem...
Murray

[mjgould]

My fault, Murray. I was thinking I had to specify the number of levels in the factor instead of their names. I may have wrongly assumed that looking at a previous post.

I tried listing the factor levels, however, R spat back and error listing the last three cov names as unused arguments. Therefore, I'm going to try the second route you suggest, but I'm a little confused by where "trps" is coming from? Would that be the name of the model, capt.hist, or other?

[murray.efford]

My 'trps' was just my shorthand for the name of an arbitrary 'traps' object. In the case of the dataset captdata, trps <- traps(captdata). I'm assuming your 'veg' covariate was recorded at detector sites, hence exists as a covariate at the level of the traps object.
Murray