0
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – rvl Jan 21 '16 at 22:50
  • $\begingroup$ @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. $\endgroup$ – deargle Jan 22 '16 at 2:23
  • $\begingroup$ Right, but those lsmeans often comprise a linear combination of several coefficients -- not just one. $\endgroup$ – rvl Jan 22 '16 at 2:57
1
$\begingroup$

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)
m3.lsm = lsmeans(m3.only_within, c("ageGroup", "filter"))
contrast(m3.lsm, list(diag.comp = c(0, -1/2, 1/2, 0)))

I still, however, think you are going to confuse other people, if not yourself, by calling this a comparison of marginal LS means of different factors.

$\endgroup$
  • $\begingroup$ Beautiful, thank you for spelling out how the contrast works. Enlightening to me but probably not surprising to you was when I saw that the mean difference is the same but the SE from the contrast is 1/2 the one I got by hand. I'll stop calling it a comparison of cross-factor marginal LS means. $\endgroup$ – deargle Jan 22 '16 at 19:42

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.