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

[Douglas Bates]

Thanks for pointing this out, Ben.

At the risk of muddying the waters, I will point out that there are
two ways of specifying an offset and the way that Ben suggests is the
preferred one (the other is with an offset() term in the model
formula).  Using the offset argument is simpler and easier to modify.

Also, you may find it handy to assign, say

myoff <- (x==1)*log(0.8)

and use offset = myoff.  Creation of a model frame and the
corresponding model matrices is more subtle than it may seem at first,
and sometimes evaluation of the offset as an expression is tricky.

On Sat, Jan 3, 2009 at 9:23 AM, Ben Bolker <bolker at ufl.edu> 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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>