# 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?

Edit 2021-07-26: I think what I probably conceptually wanted to do at the time I wrote this was compare the relative magnitude of the estimates. I'm pretty sure I could have done that by doing t-tests with the standards and their standard errors, doing some hand-waiving due to the repeated measures nature of the data. I certainly did not write this question that way though.

• 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. Jan 21, 2016 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. Jan 22, 2016 at 2:23
• Right, but those lsmeans often comprise a linear combination of several coefficients -- not just one. Jan 22, 2016 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)