1
$\begingroup$

I'm having trouble with lsmeans outputs from my glmmadmb & GLMER models for poisson (count data) with LOTS of zeros. Fixed effect is period (late or early), response is no. of birds per km, and random term is degree square sampled (block) Here is an example of my code:

ztpoiss <- glmmadmb(YBK ~ (PERIOD) + offset(logdist) + (1 | BLOCK), 
  data = counts, zeroInflation = TRUE, family = "poisson")

summary(ztpoiss)
lsmeans (ztpoiss, pairwise ~ PERIOD)

Back-transformed means do not match up to the raw data (way out). Using estimates generated by the main models to manually calculate back-transformed means do work. Is it something to do with the offset?

Here is a snippet of my data, just using one particular species (yellow-billed kite = YBK).

PERIOD BLOCK      DATE DAY MONTH YEAR SEASON    BIOME PA HA TIME DIST YBK
     L  2221 29-Apr-15  29   Apr 2015    Dry Kalahari UN NO   59 63.0   0      
     L  2226 23-Apr-15  23   Apr 2015    Dry Kalahari UN NO  162 93.5   0      
     L  2226 30-Jul-15  30   Jul 2015    Dry Kalahari UN NO  100 73.5   0      
     L  2226 31-Jul-15  31   Jul 2015    Dry Kalahari UN NO  123 71.0   0      
     L  2226 26-Jan-16  26   Jan 2016    Wet Kalahari UN NO  102 73.5   4      
     L  2226 27-Jan-16  27   Jan 2016    Wet Kalahari UN NO   88 77.5   7  
     E  1724  8-Aug-92   8   Aug 1992    Dry SubTrop  PR NO   UNK 152   0          
     E  1724 14-Jan-93  14   Jan 1993    Wet SubTrop  UN NO   UNK  44   0      
     E  1724 15-Jul-93  15   Jul 1993    Dry SubTrop  UN NO   UNK  51   0      
     E  1724 24-Feb-94  24   Feb 1994    Wet SubTrop  UN NO    51  42   2      
     E  1724 28-Feb-94  28   Feb 1994    Wet SubTrop  UN NO    36  32   0      
     E  1724 28-Feb-94  28   Feb 1994    Wet SubTrop  UN NO    18  14   0         
$\endgroup$
  • 1
    $\begingroup$ It's going to be hard to answer this without a small reproducible example ... $\endgroup$ – Ben Bolker Sep 20 '16 at 14:53
  • $\begingroup$ The offset should be handled correctly, but note that with the code shown, the results will be on the log scale. Try adding type = "response" to the lsmeans call. $\endgroup$ – rvl Sep 21 '16 at 2:53
  • $\begingroup$ @BenBolker thanks for your comment, apologies for the late response, it didn't flag that I had a response to my question. I have edited post to include a snapshot of the data I'm using. $\endgroup$ – Bop Nov 16 '16 at 15:02
  • $\begingroup$ @rvl thanks for your input, the command you suggested doesn't solve the problem. This just back-transforms the lsmeans outputs which are not correct. $\endgroup$ – Bop Nov 16 '16 at 15:04
  • 1
    $\begingroup$ Define "correct". The model has an offset, which amounts to having a covariate with a coefficient forced to be 1. Thus, the lsmeans are in essence like adjusted means in an ancova setting, and there is no reason to believe that they should reproduce the raw averages of the counts. $\endgroup$ – rvl Nov 16 '16 at 15:35
1
$\begingroup$

To illustrate what is going on, here is an illustration using a simple dataset. I'm not showing one with zero inflation or random effects, but that doesn't matter because that isn't the issue behind the question.

Here is a fake dataset:

> fake
   treat   x count
1      1 3.7     8
2      1 2.6    13
3      1 1.8     7
4      1 1.1    11
5      2 1.7     9
6      2 3.5    17
7      2 3.3    15
8      2 1.3    15
9      3 0.7     3
10     3 4.0     4
11     3 2.6     5
12     3 2.7     6

I'll fit a Poisson regression model:

> fake.glm = glm(count ~ treat + offset(log(x)), data = fake, family = "poisson")

... and the LS means for each treat level, on the response scale:

> lsmeans(fake.glm, "treat", type = "response")
 treat     rate       SE df asymp.LCL asymp.UCL
 1     10.24457 1.640405 NA  7.485062 14.021409
 2     13.80952 1.845373 NA 10.627527 17.944244
 3      4.35000 1.025298 NA  2.740694  6.904274

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

You can verify that the raw means of count are 9.75, 14, and 4.5, respectively. These definitely do differ from the LS means. Why? Because the mean values of x are different for each treatment too. To get a fair comparison of treatments, we need to compare them with the same value of x, not different ones. The lsmeans function uses by default the mean values of covariates (including offset variables):

> mean(fake$x)
[1] 2.416667

Consider making a prediction for each treatment where x is set to this value...

> newfake
  treat        x
1     1 2.416667
2     2 2.416667
3     3 2.416667

> predict(fake.glm, newdata = newfake, type = "response")
       1        2        3 
10.24457 13.80952  4.35000

Voila! These predictions are the same as you get from lsmeans, and they are not the same as the raw mean counts. Again, what they do do is equalize for the same value of x.

Footnote

In some instances, a covariate (or offset) may be a mediating variable (also affected by some or all of the predictors) -- in which case you do want to make predictions at different values of that covariate. In this example, if you think that x is affected by treat, then you do:

> lsmeans(fake.glm, "treat", type = "response", cov.reduce = x ~ treat)
 treat  rate       SE df asymp.LCL asymp.UCL
 1      9.75 1.561213 NA  7.123714 13.344513
 2     14.00 1.870827 NA 10.774114 18.191751
 3      4.50 1.060653 NA  2.835200  7.142353

Confidence level used: 0.95 
Intervals are back-transformed from the log scale

... and in fact you now do get the raw mean counts.

$\endgroup$
  • $\begingroup$ Hi @rvl, I have veered away from trying to solve this for a while and have continued with other thesis work but alas to get back to it! Using the predict function for each treatment (in my case, each transect) doesn't produce what I am looking for. Essentially I want to get just 2 estimates and lsmeans from my data: one for the E period and one for the L period to tell me mean no. of birds per km (hence logdist as an offset, because all transect lengths were different) in each period. $\endgroup$ – Bop Feb 8 '17 at 8:31
  • $\begingroup$ Your footnote seems to head more in the right direction but from reading up, using an offset and alternatively including it as a covariate produces quite different results (and different interpretations). Using the cov.reduce argument still produces lsmeans which I know will not be the same as my raw data, but are around 60 times higher than. After reading your response to @Pharcyde I used lsm@grid and it looked as though lsm is using the mean of the offset (logged distance) in its calculations, but this value was same for both periods. I just can't seem to figure out what is going wrong! $\endgroup$ – Bop Feb 8 '17 at 8:36
  • $\begingroup$ Can you find a statistical consulting service on your campus? This isn't something that's going to get reliably solved via social networks like this one. $\endgroup$ – rvl Feb 8 '17 at 23:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.