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

[Ben Bolker]

  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!
>>>
>>>
>>
>>
> 


-- 
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / www.zoology.ufl.edu/bolker
GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc



-- 
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / www.zoology.ufl.edu/bolker
GPG key: www.zoology.ufl.edu/bolker/benbolker-publickey.asc