Skip to content

'adjusting for bias' in a Poisson model

3 messages · Henrik Pärn, Ben Bolker, Douglas Bates

#
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!
#
"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:

  
    
#
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: