# Tag Info

18

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 ...

9

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.

8

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 <- ...

7

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 ~ ...

7

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 ...

6

It is a lot easier to do using the lsmeans package library(lsmeans) lsmeans(mod, pairwise ~ tension | wool) lsmeans(mod, pairwise ~ wool | tension)

6

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))) grid

5

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 ...

4

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, ...

4

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 ...

4

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 ...

4

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 ☺ ).

4

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) > library(emmeans) > (warp.emm = emmeans(warp.lm, ~ wool | tension)) tension = L: wool emmean SE df lower.CL upper.CL A 44.55556 3....

4

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 ...

4

NA is just a code that no degrees of freedom are needed. The tests and CIs are based on the standard normal rather than the t distribution. Note that the headings say z rather than t

4

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, ...

4

OK, let us dissect this model. First, the model itself: > getOption("contrasts") unordered ordered "contr.treatment" "contr.poly" > warp.lm <- lm(breaks ~ wool*tension, data = warpbreaks) > summary(warp.lm) ... (some output omitted) ... Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) ...

4

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. More details Let me create your data in a reproducible way, by ...

4

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: $$\hat\pi=\frac{1}{1+e^{-(-1.9684+0.2139\times d)}}$$ where $d$ is an indicator for membership in group 2. So the predicted probabilities are: ...

3

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 ...

3

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 ...

3

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. Given ...

3

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 ...

3

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 ...

3

Look at the model summary: > m Linear mixed model fit by REML ['lmerMod'] Formula: dep ~ f1 * f2 + (1 | sub) Data: data REML criterion at convergence: 371.7578 Random effects: Groups Name Std.Dev. sub (Intercept) 1.222 Residual 4.768 Number of obs: 64, groups: sub, 8 Fixed Effects: (Intercept) f1Male ...

3

Yes. 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.

3

People seem to have an assortment of pretty rigid practices. I'll just comment that I think ANOVA F tests are kind of overrated. They can be useful for deciding what model is suitable for a dataset; but they do not help much, in my opinion, for doing any meaningful inference. In particular, if you have fitted a two-way factorial model including interaction,...

3

Note that in both sets of comparisons, the estimates and standard errors for the comparisons within the same exposure match up. It is only the $P$ values that differ. The reason is that the $P$ values in the first set are adjusted for multiplicity, and the $P$ values in the second set are not. The first set of EMM comparisons consists of ${8\choose2}=28$ ...

3

You are indeed right that this difference can be explained from the missing data you have. In particular, when you have missing data that are of the missing at random type, then the observed data are not a representative sample of your target population. In this case, the simple sample means will be biased and should not be trusted. The mixed model, on the ...

3

EMMs are a linear function of the regression coefficients. You can always tell what that linear function is by looking at the linfct slot of the object that emmeans() creates: > mod = lm(formula = y ~ x1 + x2, data = d) > (emm = emmeans(mod, "x1")) x1 emmean SE df lower.CL upper.CL T1 10.30 0.820 66 8.66 11.94 T2 13.03 0.542 66 11.95 ...

Only top voted, non community-wiki answers of a minimum length are eligible