I have bird count data and use classical poisson loglinear model, i.e. we have counts obs(i,j)
- observed count for site i
and year j
, and the model is:
ln(model(i,j)) = site_effect(i) + year_effect(j)
where model(i,j)
is expected count for site i
and year j
, and it is assumed that counts are poisson distributed, i.e. obs(i,j) ~ Poiss(lambda = model(i,j))
. Serial correlation and overdispersion was taken into account.
Now I need to analyse residuals of such model in another regression with other explaining variables analyse residuals of this model as a dependent variable in another regression with some explaining variables (Please don't tell me I should put these variables in the above model - to be brief, there are some technical issues behind that.)
The problem now is:
- how to define the residuals
- what is their distribution
- which regression use to explain them.
Possible solutions of residual definition and analysis:
A) obs/model with another poisson loglinear regression
B) log((obs+1)/model), normal distribution -> normal linear regression
C) log((obs+1)/(model+1)), normal distribution -> normal linear regression
D) ... any other??
Ad Solution A: use obs/model
,i.e. divide observed value by the expected value from model (for each i,j). But what is their distribution? I would say Poisson, but the numbers are not integerwhole numbers! So I multiplied by 20 and rounded to whole numbers to be able to fit poisson. Is this a good idea, or is it nonsense? Or how should I fit poisson to real number continuous random variable? Here is the distribution:
Note that the Poisson distribution doesn't fit (used goodfit() test in R).
So there is a problem that obs/model are real numbers - not whole natural, as the poisson distribution. Moreover, I'd like to avoid another poisson loglinear regression if possible, as I'm quite infamiliar with how to analyze the explaining variables, the explained/unexplained variability etc. as in the normal regression.
Ad Solution B/C: use something like log((obs+1)/model)
or log((obs+1)/(model+1))
. There is a problem of zero observed counts (there are no zero model counts), so I can't just use log(obs/model)
as I would like to. How to solve this problem?? It is better to use solution B or C?
Here are the distributions.
Note that also in this case, the normal distributions don't fit. Nothing really fits, but this is quite normal with real count data. :-)
So, summary:
- which solution you think fits best?
- how to overcome the problems each solution has?
- is it a big mistake to use B/C and go for normal linear regression, if I'd like to avoid the poisson? Should I use B or C?
Thanks a lot!