[R-sig-ME] 'adjusting for bias' in a Poisson model

[Henrik Parn]

Thanks for the comment!

I would also be concerned. However, although the test.data and the real 
data are similar in the sense that both have two random grouping 
variables, the real data has, among other things, more levels in both 
groups. I just wanted to keep the test.data small, but apparently I lost 
some (too much?) realism on the way (I wish I had taken the short course 
"How to efficiently produce frequently used minimal, self contained, 
reproducible code and example data in R").

Henrik

Ben Bolker wrote:
>   PS I would be somewhat concerned with your number of
> random-effects levels (4 for group 1, 2 for group 2).  It's
> not at all surprising that you got an estimated value
> of zero variance for group 2 ... there are various schools
> of thought on this, but I would suggest that you might
> alternatively try (1) fitting the grouping variables as
> fixed effects and (??) (2) Bayesian analysis with appropriate
> non-informative priors on the variances [this is a big
> can of worms, though].
> 
>   [Forwarding back to the list]
> 
>   Ben Bolker
> 
> 
> Henrik Parn wrote:
>> Thanks a lot Ben for your very rapid answer! I will try it out!
>>
>> Cheers,
>>
>> Henrik
>>
>> Ben Bolker wrote:
>>>   "offset" is you want to look for.
>>>   offsets apply on the scale of the linear predictor
>>> (log scale in this case), so I think adding something
>>> like offset=((x==1)*log(0.8)) to the model will do
>>> the trick.
>>>
>>> Henrik Parn wrote:
>>>> Dear all,
>>>>
>>>> Basically, I have fitted a mixed model with poisson error to analyse how
>>>> number of offspring (y) depend on a fixed factor (x). The data is
>>>> grouped by two random factors (gr1, gr2).
>>>>
>>>> # Some test data with at least a similar structure:
>>>>
>>>> set.seed(100)
>>>> test.data <- data.frame(
>>>> y = c(rpois(40, 5.5), rpois(40, 3.5)),
>>>> x = factor(rep(0:1, each = 40)),
>>>> gr1 = factor(rep(1:4, each = 10)),
>>>> gr2 = factor(rep(1:2, each = 5)))
>>>>
>>>> # The model
>>>> model <- lmer(y ~ x + (1|gr1) + (1|gr2),
>>>> data = test.data, family = poisson)
>>>>
>>>> summary(model)
>>>> Generalized linear mixed model fit by the Laplace approximation
>>>> Formula: y ~ x + (1 | gr1) + (1 | gr2)
>>>>    Data: test.data
>>>>    AIC   BIC logLik deviance
>>>>  61.01 70.54 -26.51    53.01
>>>> Random effects:
>>>>  Groups Name        Variance   Std.Dev.
>>>>  gr1    (Intercept) 0.00045437 0.021316
>>>>  gr2    (Intercept) 0.00000000 0.000000
>>>> Number of obs: 80, groups: gr1, 4; gr2, 2
>>>>
>>>> Fixed effects:
>>>>             Estimate Std. Error z value Pr(>|z|)
>>>> (Intercept)  1.75331    0.06666  26.302  < 2e-16 ***
>>>> x1          -0.48660    0.10665  -4.563 5.05e-06 ***
>>>>
>>>>
>>>> Thus, y is significantly higher for x = 0 than for x = 1. However, it
>>>> has been suggested that the estimate of x may be biased (for biological
>>>> reasons) and actually be less negative than what is found above.
>>>> Specifically, even in absence of an effect of x, the bias could cause y
>>>> for x = 1 to be 20% lower than y for x = 0.
>>>> (y for x=0 - y for x=1)/ y for x=0 = 20%; y for x=1 could be 0.8*y for
>>>> x=0 due to bias.
>>>>
>>>> Thus, the estimate of x would be (1.75331 - 0.2 * 1.75331) - 1.75331 =
>>>> -1.75331 * 0.2 = -0.350662 just due to the bias.
>>>>
>>>> I wish to test if there is an effect of x on y over and above the
>>>> potential 20% bias.
>>>>
>>>> My naive starting point: Instead of testing if the estimate of x
>>>> (-0.48660) differs from zero, I would need to test if -0.48660 -
>>>> (-0.350662) = -0.135938 differs fro zero.
>>>>
>>>> Or could I somehow make the adjustment already in the data set, e.g.
>>>> adjusting the y's for x = 1? But I assume that I cannot 'just add ?? %
>>>> to y in group x', because the response has to be integers?
>>>>
>>>> Does anyone have a suggestion of a convenient way of performe test for
>>>> significance of x on y while taking the bias into account?
>>>>
>>>>
>>>> Thanks a lot in advance!
>>>>
>>>>
>>>
> 
> 

-- 
Henrik Pärn
Centre for Conservation Biology
Department of Biology
Norwegian University of Science and Technology
NO-7491 Trondheim
Norway

Office: +47 73596285
Fax: +47 73596100
Mobile: +47 90989255

E-mail: henrik.parn at bio.ntnu.no