Solved – Confusing LSMEANS results

glmmlsmeansmultiple-comparisonsr

I am having trouble interpreting results from a multiple comparison test. Here is the model I ran:

m1 <-glmer(resp~Grp-1+offset(LArea)+(1|Crk/Meso)+(1|resid),data=df, family=binomial)

And here are the results from the lsmeans function:

> lsmeans(m1, pairwise ~ Grp, adjustSigma = TRUE, adjust = "sidak")
$lsmeans
 Grp      lsmean        SE df asymp.LCL  asymp.UCL
 90U   -0.207353 0.6345219 NA -1.450993  1.0362871
 90I   -2.214079 0.7946987 NA -3.771660 -0.6564978
 96U   -1.605763 0.6551178 NA -2.889770 -0.3217555
 96I   -3.168827 0.7957100 NA -4.728390 -1.6092644
Confidence level used: 0.95

OK, everything is great so far. I can see that the lsmean order is 90U>96U>90I>96I. So now I want to check the observed differences are significant: 

> cld(lsmeans(m1, pairwise ~ Grp, adjustSigma =TRUE, adjust = "sidak"))
 Grp      lsmean        SE df asymp.LCL  asymp.UCL .group
 96I   -3.168827 0.7957100 NA -4.728390 -1.6092644  1    
 90I   -2.214079 0.7946987 NA -3.771660 -0.6564978  12   
 96U   -1.605763 0.6551178 NA -2.889770 -0.3217555  1    
 90U   -0.207353 0.6345219 NA -1.450993  1.0362871   2   

Confidence level used: 0.95 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05

But now I am confused. How can 90I be grouped with 90U while excluding 96U? It is saying -2.2 = -0.2 while -2.2!=-1.6 and -1.6 != -0.2 If it helps, here is an image with the 95%CI:
Plot of estimated coefficients

The explanations I have thought of but find unconvincing are:

  • lsmean ignores the random effects, effectively running a plain-jane pairwise comparison.
  • This is a perfectly reasonable conclusion.

My running hypothesis:

  • I_have_no_idea_what_I_am_doing.jpg

And a brief aside: Why does it perform a Tukey adjustment when I ask for a Sidak?

Best Answer

Short answer -- statistical [non]significance is not transitive.

It will help you if you look at confidence intervals for the differences. The following will show what's going on with the default tukey-adjusted results:

confint(pairs(lsmeans(m1, "Grp")))

You will find that some comparisons have wider CIs than others. The reason Tukey overrode your Sidak is that your lsmeans call asks for both LSmeans and pairwise comparisons thereof -- two different sets of linear functions. If you put the adjust spec in the cld call instead of the lsmeans call, it'll do the right thing.

Related Question