Not a complete answer...
The difference between glht(myfit, mcp(myfactor="Tukey")) and the two other methods is that this way uses a "z" statistic (normal distribution), whereas the other ones use a "t" statistic (Student distribution). The "z" statistic it the same as a "t" statistic with an infinite degree of freedom. This method is an asymptotic one and ...
The short answer is that LS means (or more modernly, estimated marginal means) are incredibly useful with experimental data. With observational data, not so much. A long-winded explanation follows.
The underlying ideas are very old (and predate SAS by at least 50 years). Look at a standard experimental design textbook -- pretty much any of them. Such a book ...
The output represents predictions from your model for each image. With the poison family, the default link function is the natural log - so those values are on the log scale. If you do lsmeans(..., type = "response"), it will back-transform the predictions to the original response scale.
The history of the least squares mean, its appearance in SAS, and its interpretation is discussed in Searle, Milliken, and Speed (1979).
Some discussion of the concepts around least squares means (population marginal means) is found in Searle, Speed, and Milliken (1980). Earliest mention of the concept that they note is Damon et al (1959). They provide ...
Provision in upcoming version of lsmeans
The next update of lsmeans (2.20 or later) will include an rbind method for ref.grid and lsmobj objects. It makes it easy to combine two or ore objects into one family, and defaults to the "mvt" adjustment method. Here is the present example:
> w.t <- pairs(lsmeans(warp.lm, ~ wool | tension))
> t.w <- ...
Both of the lsmeans statements you show generate lists of lsmobjs, and each element of those lists is handled separately. If you want to incorporate an overall adjustment for two or more lists combined, it is technical and it takes a bit of work.
First, save the list:
lsmlist = lsmeans(warp.lm, list(pairwise ~ wool|tension,
pairwise ~ ...
The following code illustrates how this computation is done:
grid <- with(data, expand.grid(f1 = levels(f1), f2 = levels(f2)))
X <- model.matrix(~ f1 * f2, data = grid)
V <- vcov(m)
betas <- fixef(m)
grid$emmean <- c(X %*% betas)
grid$SE <- sqrt(diag(X %*% V %*% t(X)))
It is doable! You have to do it in stages, though. Using the warpbreaks data to illustrate, I'll do such comparisons of wool at each tension:
> warp.lm = lm(breaks ~ wool * tension, data = warpbreaks)
> (warp.emm = emmeans(warp.lm, ~ wool | tension))
tension = L:
wool emmean SE df lower.CL upper.CL
A 44.55556 3....
I disagree strongly with the "only situation" in the OP. EMMs (estimated marginal means, more restrictively known as least-squares means) are very useful for heading off a Simpson's paradox situation in evaluating the effects of a factor. In your example, consider a scenario where these three things are true:
When $x_2$ is held at any fixed level, the ...
It uses $t$ tests, the observed contrast divided by the estimated standard error. It gets this information from the fitted model. Thus, the validity of the result depends on the validity of the model.
If, for example, you fitted a model using lm(), and that the errors are actually normally distributed with common variance (as assumed in the model), and the ...
You can exert more control over what you get if you do the contrasts in a separate step:
lsm = lsmeans(model, ~ factor1 | factor2*factor3)
comps = pairs(lsm)
summary(comps, by = NULL, adjust = "holm")
The by = NULL part of the last statement tells it to ignore the by grouping and thus treat it as one family of tests.
To answer part of the first question, ...
Least-squares means are not means of the data values. They are marginal averages of predictions from the model, taken over a grid consisting of all factor combinations.
If your model had included all the two- and three-way interactions, then those predictions would be the same as the cell means, making the LS means the same as the raw means, provided the ...
I'll try to explain it briefly.
First, it is important to know that one should never use side-by-side confidence intervals to do comparisons. The reason is because the standard errors of means are not the same as the standard errors of differences of means. In fact, they can be hugely different, as when comparing levels of a within-subjects treatment, where ...
It makes perfect sense if your data are balanced, i.e., you have the same number of observations in each Treatment * SubjectType * Year combination. (That's how ANOVA was done in the 1940s, and balance like that was the reason ANOVA was feasible in 1940s ☺ ).
There are a few simple reasons to report LS means or EM means and their related statistics.
In cases of unbalanced designs, the EM means themselves may not equal the arithmetic means, and EM means may be more meaningful. With your toy data, you can compare,
dat <- data.frame(
trt = factor(c(1, 1, 1, 1, 2,2,2,2,2,2)),
eff = c(10, 12, 11, 15, ...
Your t.test() results are each based on only selected portions of the data, whereas the emmeans() results are based on a model that is fitted to all of the data. Consequently, the denominators for the $t$ statistics are based on pooled information in the model. More information, more df.
Let me create your data in a reproducible way, by ...
There is enough information in the question to clear this up. lsmeans is simply using the coefficients to obtain the group predicted probabilities.
So for the GLM, OP's implied model is:
where $d$ is an indicator for membership in group 2. So the predicted probabilities are: ...
contrast(emmeans(m, ~f1*f2), interaction = “consec”)
But @whuber’s comment does apply in a case as simple as this one.
Or you may prefer
pairs(pairs(emmeans(m, ~ f2|f1)), by = NULL)
That is, compute the results you already have, remove the by grouping, and compare the two pairwise differences.
I will explain using a somewhat simpler model, but with the same kind of discrepancy. Consider the pigs dataset in the emmeans package.
pigs$pct = factor(pigs$percent)
I'll fit an ordinary regression model, making the contrast coding explicit so there is no question about how it is parameterized:
mod = lm(conc ...
I'm glad you love the package; but it reminds me of a letter I read once in the newspaper column Hints for Heloise, in which the writer said that she had a new dishwasher and loves it, but it doesn't get her dishes clean. It made me wonder what exactly she loves about it, and whether she feels the same way about her husband.
As the developer of lsmeans, I ...
This seems to be described in the vignette you quote (p. 10, italics added by myself):
Another way to display pairwise comparisons is via the comparisons
argument of plot. When this is set to TRUE, arrows are added to the
plot, with lengths set so that the amount by which they overlap (or
don't overlap) matches as closely as possible to the amounts ...
My view is that the $F$ test of statistical significance of the interaction effect is less important than the subjective nature of the interaction, as evidenced by the plot. The plot tells me that it is reasonably sensible to compare the overall averages of Depression and Top, but it'd be silly to compare those averages with the overall average of Slope -- ...
I'm not sure you understand the implications of the model you fitted. Since X is a quantitative predictor and Fact is a factor, your model has fitted 3 straight lines having different slopes and intercepts. It's almost the same as separately fitting the 3 lines with the data split according to Fact, except that they share common variance estimates.
I think what's most useful to you is that the  operator is defined for lsmeans() results (also emmeans::emmeans()). So you can do:
TL_lsm = lsmeans(TL_interact, ~ Year*Sex|River)
TL_lsm_pairs = pairs(TL_lsm)
wanted = c(1,2,3,5,8,13,21) ## whichever ones you want
summary(TL_lsm_pairs[wanted], adjust = "fdr")
There are also other contrast families ...
There are two extra points:
Because you’re using a mixed effects Poisson regression with a log (nonlinear) link function, the coefficients you obtain have an intepretation conditional on the random effect. I.e., an interpretation conditional on site. This seems to be different than the plot you produce where you average over sites. Hence, it is a bit as ...
Look at the model summary:
Linear mixed model fit by REML ['lmerMod']
Formula: dep ~ f1 * f2 + (1 | sub)
REML criterion at convergence: 371.7578
Groups Name Std.Dev.
sub (Intercept) 1.222
Number of obs: 64, groups: sub, 8
(Intercept) f1Male ...
Your model includes an interaction term, so it invites the following question:
When you designed your study, did you envision conducting a set of planned comparisons on the basis of this model?
A planned comparison is something you choose carefully when designing your study (hence prior to seeing your data). For example, comparing the mean body mass ...