5
$\begingroup$

this is my first post, so I hope everything is in the right format. I have some problems with glmer and don't know how to fix it, so I hope somebody can help me out with this. I could not find an answer to this anywhere.

My experiment: The experimental setup has eight sites, each site has a central area in which organisms get marked. From the central area, organisms have the choice to go into three different areas which are equipped with traps. The three areas are crossed with sites. The recapture rates in these three areas are very low and the sampling effort (trap days) is very high. Due to external influences, traps got destroyed to different degrees in all the sampling areas. In addition to the marked organisms, unmarked organisms of the same species get caught, too.

The dataset looks something like this:

library(lme4)
data <- structure(list(site = c("A", "A", "A", "B", "B", "B", "C", "C", "C", "D", "D", "D", "E", "E", "E", "F", "F", "F", "G", "G", "G", "H", "H", "H"), area = c("I", "II", "III", "I", "II", "III", "I", "II", "III", "I", "II", "III", "I", "II", "III", "I", "II", "III", "I", "II", "III", "I", "II", "III"), marked = c(2, 6, 3, 5, 3, 9, 0, 8, 1, 1, 1, 18, 3, 0, 0, 1, 5, 6, 3, 0, 2, 2, 4, 5), unmarked = c(38, 78, 104, 1, 6, 10, 1, 13, 0, 13,7, 85, 7, 1, 0, 9, 4, 36, 3, 4, 3, 10, 20, 29), sampl_effort = c(9300, 9100, 8700, 9900, 9600, 8600, 9800, 9400, 10800, 11600, 11000, 13950, 10300, 9700, 9800, 10450, 10100, 10800, 9600, 9900, 9300, 11800, 11250, 9450)), .Names = c("site", "area", "marked", "unmarked", "sampl_effort"), row.names = c(NA,-24), class = "data.frame")

Now I wanted to fit a glmer. Because the amount of marked organisms may be related to the total amount of organisms, I chose to take a binomial approach, with cbind(marked, unmarked). I use area with the three treatments as explanatory variable, and site as random factor. Because the sampling effort differs between the different areas, I want to include it as an offset. The code looks like this:

mod.glmer1= glmer(cbind(marked,unmarked) ~ area + (1 | site) + offset(sampl_effort), family=binomial, data=data)

Then, I get the error:

Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

If I try the glmer without the offset, everything works out fine:

mod.glmer2= glmer(cbind(marked,unmarked) ~ area + (1 | site), family=binomial, data=data)

Out of interest, I tried a glm without the random factor and with offset:

mod.glm1= glm(cbind(marked,unmarked) ~ area + offset(sampl_effort), family=binomial, data=data)

and I get the following Warnings:

Warning messages: 1: glm.fit: fitted probabilities numerically 0 or 1 occurred 2: glm.fit: fitted probabilities numerically 0 or 1 occurred

By looking at the fitted values

round(cbind(data[,3:5],fits=fitted(mod.glm1)),8)

I can see that all the fitted values are 0. I now thought that the offset is just too large to get reasonable fitted value. As the sampling effort is arbitrary (instead of days, I could have taken sampling hours, weeks, etc.), I decided to divide the sampling effort by 1000. By doing this, the glm now works:

mod.glm2= glm(cbind(marked,unmarked) ~ area + offset(sampl_effort/1000), family=binomial, data=data)

However, for the glmer

mod.glmer3= glmer(cbind(marked,unmarked) ~ area + (1 | site) + offset(sampl_effort/1000), family=binomial, data=data)

I still get

Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

So my questions are basically:

  1. Am I allowed to divide the sampling effort by 1000? In my mind it should lead to the same results, as the relative differences in sampling effort stay the same (and sampling days are a arbitrary measurement). However, I of course tried it with e.g. dividing by 10000 and I get different results.

  2. How can I include the sampling effort in the glmer? The sampling effort is just too important to keep it out, however, I also need the sites as random factor. Is the offset the right approach and if yes, why doesn't it work.

P.S.: My session information:

R version 3.0.2 (2013-09-25) Platform: i386-w64-mingw32/i386 (32-bit)

locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages: [1] stats graphics grDevices utils datasets methods base

other attached packages: [1] car_2.0-19 lme4_1.0-6 Matrix_1.1-2 lattice_0.20-23

loaded via a namespace (and not attached): [1] grid_3.0.2 MASS_7.3-29 minqa_1.2.3 nlme_3.1-111 nnet_7.3-7
[6] Rcpp_0.11.0 splines_3.0.2 tools_3.0.2

$\endgroup$
  • $\begingroup$ this is a little bit tricky. Don't have time for an answer right now, but the basic issue is that the offset gets added to the model on the linear predictor scale. Adding offsets to a binomial model generally makes the most sense if you are using a complementary log-log link (family=binomial(link="cloglog")) (google "cloglog offset" or see rpubs.com/bbolker/4082 for an example). This approach may or may not make sense for your case -- I haven't had time to think about it. $\endgroup$ – Ben Bolker Mar 5 '14 at 23:14
  • $\begingroup$ @ Ben Bolker: Thank you so much for your fast answer. I do not know if it makes sense or not (still trying to think this through), the only thing I can say so far is that mod.glmer4= glmer(cbind(marked,unmarked) ~ area + (1 | site) + offset(sampl_effort), family=binomial(link="cloglog"), data=data) gives me the Error Error in (function (fr, X, reTrms, family, nAGQ = 1L, verbose = 0L, control = glmerControl(), : Downdated VtV is not positive definite $\endgroup$ – John Mar 6 '14 at 15:28
  • $\begingroup$ If you did want to use an offset with a cloglog link in the standard way I believe you would need offset(log(sample_effort)) rather than offset(sample_effort) -- see the link I posted. $\endgroup$ – Ben Bolker Mar 6 '14 at 18:02
  • $\begingroup$ Oh yes, my bet. However, still get the same error. $\endgroup$ – John Mar 6 '14 at 18:36
2
$\begingroup$

I don't think you do want to add sampling_effort as a offset at all.

You are estimating the logit of marked over unmarked catch probabilities - call that quantity $\eta$ - and you think, quite reasonably, that the quantity $N=$ marked+unmarked will be a function of how long the traps remain undamaged, i.e. a function of sampling_effort. However, it's not obvious why that relationship is relevant; your binomial model already treats $N$ as known and conditions on it.

The circumstance I can imagine where it would be relevant is when $\eta$ really depends on sampling_effort. To examples would be a) if the rate of unmarked captures were constant but the rate of marked captures is increasing because of some feature of the experiment such as marked animals being expected to take a random walk from a common origin, or b) if capturing one kind of animal makes it more or less likely to capture the other type. But even in these cases sampling_effort would be a regular covariate and not have its coefficient fixed at 1 as an offset would have.

Here's some possibly relevant discussion of this question elsewhere on the site.

$\endgroup$
  • $\begingroup$ @concugateprior That is an interesting thought. So basically what you say is that if sampling effort decreases, it should (under normal circumstances) decrease the likelihood to catch marked as well as unmarked individuals and thereby not affect the marked over unmarked catch probabilities? And by adding sampling_effort as covariate I would actually see if it affects the probabilities or not, right? $\endgroup$ – John Mar 6 '14 at 22:58
  • $\begingroup$ Exactly. Either sampling_effort affects the precision of the estimate of the marked/unmarked balance by increasing marked+unmarked, or sampling_effort affects the balance itself. If the former you don't need to add it to the model explicitly because its already conditioned on, and if the latter then it ought to be a regular covariate. $\endgroup$ – conjugateprior Mar 6 '14 at 23:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.