I am using a generalized Linear Mixed-Effects model to look at the effects of different treatments on a density of trichomes.

The model is :

fitPoisson = glmer(Count_trichomes ~ Treatment1*Treatment2*Treatment3 + 
                         (1 | Block/Code) + offset(log(Length)), family=poisson(), data=dataset)

Treatment 1 and 2 has 2 levels (0 and 1) and Treatment 3 has 3 levels (0,1,2). Block accounts for the replicates and Code, for each individual. Length is in cm.

An anova(fitPoisson) told me that treatments 1 and 3 are significant and that there is no interactions. What I want now is to know what the density is for each level of treatments.

So I used a lsmeans to look at the differences :

    > lsmeans(fitPoisson, ~ Treatment1)

     Treatment1   lsmean         SE df asymp.LCL asymp.UCL
     0           5.309106 0.06113705 NA  5.189280  5.428933
     1           5.471452 0.06114033 NA  5.351619  5.591285

     Results are averaged over the levels of: Treatment2, Treatment3
     Results are given on the log (not the response) scale. 
     Confidence level used: 0.95

I can see that the density of level 0 is lower than the density of level 1, but I dont understand what are the units used. It doesn't seems like it is for trichomes/cm, since the mean for level 0 is 107 trichomes/cm and the mean for level 1 is 131 trichomes/cm (calculated in excel).

When I transform back from the log scale, it gives me :

    > summary(lsmeans(fitPoisson, ~ Treatment1), type = "response")

     Treatment1   rate       SE df asymp.LCL asymp.UCL
     0           202.1694 12.36004 NA  179.3393  227.9058
     1           237.8053 14.53949 NA  210.9496  268.0799

    Results are averaged over the levels of: Treatment2, Treatment3 
    Confidence level used: 0.95 
    Intervals are back-transformed from the log scale 

Which is still far from the means I found in excel.

Maybe I just don't understand the information lsmeans is giving me, or I am not using the right function.

  • $\begingroup$ With respect, I'm not sure that answers the question, rvl. I believe that the OP may be having trouble with lsmeans recognizing the offsets used in the glmer call. Instead of reporting average rates lsmeans is reporting average counts. I don't have an answer, and I'm also experiencing this problem. $\endgroup$
    – brendan
    Aug 27 '16 at 15:04

If you do

lsmeans(fitPoisson, ~ Treatment1 * Treatment2 * Treatment3)

you will see the predictions (on the log scale) made by the model for each combination of the three factors. The results you show in your question are the averages of these predictions, averaged (with equal weights) over the levels of Treatment2 and Treatment3. The results with type = "response" are the antilogs of these results, and the standard errors are obtained using the delta method. Note that the averaging is still done on the log scale, and the confidence intervals are computed on the log scale, then back-transformed.

If you want the averaging to be done on the raw count scale, that is possible too. Do:

rg = regrid(ref.grid(fitPoisson), transform = TRUE)
lsmeans(rg, ~ Treatment1)

(The regrid function creates a new reference grid for the model based on the back-transformed predictions. You can use summary(rg) to see the individual predictions.)

These results can differ markedly from the raw averages you obtained from the data when there is imbalance in the data, so that the raw averages give far from equal weights to the levels of those two factors.


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.