# Comparing lsmeans of levels from different factors

I have a glmer model that I am performing contrasts on. My model has four factors (emotion, gender, filter, ageGroup). No interactions. I have been using lsmeans(x, spec, contr='pairwise')) (and) lsm within glht) to get lsmeans and to do pairwise contrasts within-factors. I am also interested in comparing the means from one factor against a level from a different factor.

> summary(m7.only_within)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial  ( logit )

Formula: FaceStimulus.RESP.binary ~ emotion + emotion.Sample + gender +
gender.Sample + ageGroup + ageGroup.Sample + filter + filter.Sample + (1 | Subject)

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)                 -0.823451   0.264442  -3.114  0.00185 **
emotiondisgust               1.504896   0.092532  16.264  < 2e-16 ***
emotionfear                  1.555492   0.091832  16.938  < 2e-16 ***
emotion.Sample              -0.106948   0.044364  -2.411  0.01592 *
gendermale                   0.012718   0.074198   0.171  0.86391
gender.Sample                0.007021   0.010046   0.699  0.48462
ageGroupold                  0.364239   0.076379   4.769 1.85e-06 ***
ageGroup.Sample             -0.013519   0.009577  -1.412  0.15808
filtercropAdaptiveThreshold  0.419095   0.075906   5.521 3.37e-08 ***
filter.Sample                0.070147   0.029658   2.365  0.01802 *
---


I understand that lsmeans calculates standard errors using variable variances and covariances from vcov() to do pairwise comparisons. Can I just compare one model variable against another? Like so:

summary(glht(m3.only_within, linfct = c("filtercropAdaptiveThreshold - ageGroupold = 0")))

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
filtercropAdaptiveThreshold - ageGroupold == 0  0.05486    0.10798   0.508        1


That contrast is the test of the difference between the betas in the model, not between the lsmeans. Naturally, all my other contrasts compare the lsmeans since I'm using lsm(). Can I just manually take the difference between the lsmeans, and then use the variances and covariance SE's from the model vcov for the two model effects to get the SE for the lsmeans difference? Naturally, those SEs are much lower than the SEs reported from the lsmeans output.

$lsmeans ageGroup lsmean SE df asymp.LCL asymp.UCL young 0.01043529 0.2466025 NA -0.4728967 0.4937673 old 0.37467465 0.2474967 NA -0.1104101 0.8597594$lsmeans
filter                     lsmean        SE df   asymp.LCL asymp.UCL
cropDesat37           -0.01699274 0.2462214 NA -0.49957782 0.4655923
cropAdaptiveThreshold  0.40210268 0.2478032 NA -0.08358268 0.8877880


So it would be: $$\displaystyle \begin{array}{l} \bar{X_1}-\bar{X_2}= lsmean\_of\_cropAdaptiveThreshold - lsmean\_of\_old \\ = 0.40210268 - 0.37467465 = 0.0274\\ \\ SE_{\bar{X_1}-\bar{X_2}}= \sqrt{SE_1^2 + SE_2^2 + 2*cov}\\ = \sqrt{0.075906^2 + 0.076379^2 + 2*-3.225794e-05}\\ = 0.10798 \end{array}$$

...which is the same SE as the SE in the beta-contrast.

Is it legal to compare cross-factor level lsmeans like this?

• I think this is like comparing apples and oranges. The lsmeans for each factor are averaged over the levels of the other. And the math you have worked out is not correct, since the lsmeans are not just single coefficients. – rvl Jan 21 '16 at 22:50
• @rvl Yes, the averaging over the other factors is what makes me uncomfortable doing this. But then again, seems like there should be some way to compare apples and oranges if they're measured on the same scale. Is there anything I can do to compare the effects? Also, can you elaborate on why the math is wrong? I fixed my notation, I realized I was missing bars over the x's. lsmeans' pairwise contrasts seem to simply take the difference in estimated lsmeans and then compute SEs from the vcov matrix to get the z-score. – deargle Jan 22 '16 at 2:23
• Right, but those lsmeans often comprise a linear combination of several coefficients -- not just one. – rvl Jan 22 '16 at 2:57

Let $\hat y_{ij}$ denote the prediction from your model at the $i$th level of ageGroup $(i=1,2)$ and the $j$th level of filter $(j=1,2)$. Then the LS means you seem to want to compare are $\hat y_{2\cdot} = (\hat y_{21}+\hat y_{22})/2$ and $\hat y_{\cdot 2} = (\hat y_{12}+\hat y_{22})/2$. Consequently, $$\hat y_{\cdot 2} - \hat y_{2\cdot} = [(\hat y_{12}+\hat y_{22})/2] - [(\hat y_{21}+\hat y_{22})/2] = (\hat y_{12} - \hat y_{21})/2$$ So what you think is a comparison of two marginal LS means is really one half of a diagonal comparison of two cell means. You can get it from lsmeans using code like this:
library(lsmeans)