Movatterモバイル変換
[0]ホーム
Comparisons and contrasts in emmeans
emmeans package, Version 2.0.1
Pairwise comparisons
The most common follow-up analysis for models having factors aspredictors is to compare the EMMs with one another. This may be donesimply via thepairs() method foremmGridobjects. In the code below, we obtain the EMMs forsourcefor thepigs data, and then compare the sourcespairwise.
pigs.lm <- lm(log(conc) ~ source + factor(percent), data = pigs)pigs.emm.s <- emmeans(pigs.lm, "source")pairs(pigs.emm.s)
## contrast estimate SE df t.ratio p.value## fish - soy -0.273 0.0529 23 -5.153 <0.0001## fish - skim -0.402 0.0542 23 -7.428 <0.0001## soy - skim -0.130 0.0530 23 -2.442 0.0570## ## Results are averaged over the levels of: percent ## Results are given on the log (not the response) scale. ## P value adjustment: tukey method for comparing a family of 3 estimates
In its out-of-the-box configuration,pairs() sets twodefaults forsummary():adjust = "tukey" (multiplicity adjustment), andinfer = c(FALSE, TRUE) (test statistics, not confidenceintervals). You may override these, of course, by callingsummary() on the result with different values forthese.
In the example above, EMMs for later factor levels are subtractedfrom those for earlier levels; if you want the comparisons to go in theother direction, usepairs(pigs.emm.s, reverse = TRUE).Also, in multi-factor situations, you may specifybyfactor(s) to perform the comparisons separately at the levels of thosefactors.
Matrix displays
The numerical main results associated with pairwise comparisons canbe presented compactly in matrix form via thepwpm()function. We simply hand it theemmGrid object to use inmaking the comparisons:
pwpm(pigs.emm.s)
## fish soy skim## fish [3.39] <.0001 <.0001## soy -0.273 [3.67] 0.0570## skim -0.402 -0.130 [3.80]## ## Row and column labels: source## Upper triangle: P values adjust = "tukey"## Diagonal: [Estimates] (emmean) ## Lower triangle: Comparisons (estimate) earlier vs. later
This matrix shows the EMMs along the diagonal,\(P\) values in the upper triangle, and thedifferences in the lower triangle. Options exist to switch off any oneof these and to switch which triangle is used for the latter two. Also,optional arguments are passed. For instance, we can reverse thedirection of the comparisons, suppress the display of EMMs, swap wherethe\(P\) values go, and performnon-inferiority tests with a threshold of 0.05 as follows:
pwpm(pigs.emm.s, means = FALSE, flip = TRUE, # args for pwpm() reverse = TRUE, # args for pairs() side = ">", delta = 0.05, adjust = "none") # args for test()
## fish soy skim## fish 0.273 0.402## soy <.0001 0.130## skim <.0001 0.0013 ## ## Row and column labels: source## Lower triangle: P values side = ">" delta = 0.05## Upper triangle: Comparisons (estimate) later vs. earlier
With all threeP values so small, we have fish, soy, andskim in increasing order of non-inferiority based on the giventhreshold.
When more than one factor is present, an existing or newly specifiedby variables() can split the results into l list ofmatrices.
Effect size
Some users desire standardized effect-size measures. Most popular isprobably Cohen’sd, which is defined as the observeddifference, divided by the population SD; and obviously Cohen effectsizes are close cousins of pairwise differences. They are available viatheeff_size() function, where the user must specify theemmGrid object with the means to be compared, the estimatedpopulation SDsigma, and its degrees of freedomedf. This is illustrated with the current example:
eff_size(pigs.emm.s, sigma = sigma(pigs.lm), edf = 23)
## contrast effect.size SE df lower.CL upper.CL## fish - soy -2.37 0.577 23 -3.56 -1.175## fish - skim -3.49 0.698 23 -4.94 -2.051## soy - skim -1.12 0.490 23 -2.14 -0.112## ## Results are averaged over the levels of: percent ## sigma used for effect sizes: 0.1151 ## Confidence level used: 0.95
The confidence intervals shown take into account the error inestimatingsigma as well as the error in the differences.Note that the intervals are narrower if we claim that we knowsigma perfectly (i.e., infinite degrees of freedom):
eff_size(pigs.emm.s, sigma = sigma(pigs.lm), edf = Inf)
## contrast effect.size SE df lower.CL upper.CL## fish - soy -2.37 0.460 23 -3.32 -1.418## fish - skim -3.49 0.470 23 -4.47 -2.521## soy - skim -1.12 0.461 23 -2.08 -0.172## ## Results are averaged over the levels of: percent ## sigma used for effect sizes: 0.1151 ## Confidence level used: 0.95
Note thateff_size() expects the object with the means,not the differences. If you want to use the differences, use themethod argument to specify that you don’t want to computepairwise differences again; e.g.,
eff_size(pairs(pigs.emm.s), sigma = sigma(pigs.lm), edf = 23, method = "identity")
(results are identical to the first effect sizes shown).
Graphical comparisons
Comparisons may be summarized graphically via thecomparisons argument inplot.emm():
plot(pigs.emm.s, comparisons = TRUE)
The blue bars are confidence intervals for the EMMs, and the redarrows are for the comparisons among them. If an arrow from one meanoverlaps an arrow from another group, the difference is not“significant,” based on theadjust setting (which defaultsto"tukey") and the value ofalpha (whichdefaults to 0.05). See the“xplanations” supplement for detailson how these are derived.
Note: Don’tever use confidence intervals for EMMsto perform comparisons; they can be very misleading. Use the comparisonarrows instead; or better yet, usepwpp().
A caution: it really is not good practice to draw a brightdistinction based on whether or not aP value exceeds somecutoff. This display does dim such distinctions somewhat by allowing theviewer to judge whether aP value is close toalpha one way or the other; but a better strategy is tosimply obtain all theP values usingpairs(), andlook at them individually.
PairwiseP-value plots
In trying to develop an alternative to compact letter displays (seenext subsection), we devised the “pairwiseP-value plot”displaying all theP values in pairwise comparisons:
pwpp(pigs.emm.s)
Each comparison is associated with a vertical line segment that joinsthe scale positions of the two EMMs being compared, and whose horizontalposition is determined by theP value of that comparison.
This kind of plot can get quite “busy” as the number of means beingcompared goes up. For example, suppose we include the interactions inthe model for the pigs data, and compare all 12 cell means:
pigs.lmint <- lm(log(conc) ~ source * factor(percent), data = pigs)pigs.cells <- emmeans(pigs.lmint, ~ source * percent)pwpp(pigs.cells, type = "response")
While this plot has a lot of stuff going on, consider looking at itrow-by-row. Next to each EMM, we can visualize theP values ofall 11 comparisons with each other EMM (along with their color codes).Also, note that we can include arguments that are passed tosummary(); in this case, to display the back-transformedmeans.
If we are willing to forgo the diagonal comparisons (where neitherfactor has a common level), we can make this a lot less cluttered via aby specification:
pwpp(pigs.cells, by = "source", type = "response")
In this latter plot we can see that the comparisons withskim as the source tend to be statistically stronger. Thisis also an opportunity to remind the user that multiplicity adjustmentsare made relative to eachby group. For example, comparingskim:9 versusskim:15 has a Tukey-adjustedP value somewhat greater than 0.1 when all are in one family of12 means, but about 0.02 relative to a smaller family of 4 means asdepicted in the three-paneled plot.
Compact letter displays (CLDs)
Another way to depict comparisons is bycompact letterdisplays, whereby two EMMs sharing one or more grouping symbols arenot “significantly” different. These may be generated by themultcomp::cld() function. I really recommend against thiskind of display, though, and decline to illustrate it. These displayspromote visually the idea that two means that are “not significantlydifferent” are to be judged as being equal; and that is a very wronginterpretation. In addition, they draw an artificial “bright line”betweenP values on either side ofalpha, evenones that are very close.
Back to Contents
Other contrasts
Pairwise comparisons are an example of linear functions of EMMs. Youmay usecoef() to see the coefficients of these linearfunctions:
coef(pairs(pigs.emm.s))
## source c.1 c.2 c.3## 1 fish 1 1 0## 2 soy -1 0 1## 3 skim 0 -1 -1
The pairwise comparisons correspond to columns of the above results.For example, the first pairwise comparison,fish - soy,gives coefficients of 1, -1, and 0 to fish, soy, and skim, respectively.In cases, such as this one, where each column of coefficients sums tozero, the linear functions are termedcontrasts
Thecontrast() function provides for general contrasts(and linear functions, as well) of factor levels. Its second argument,method, is used to specify what method is to be used. Inthis section we describe the built-in ones, where we simply provide thename of the built-in method. Consider, for example, the factorpercent in the modelpigs.lm . It is treatedas a factor in the model, but it corresponds to equally-spaced values ofa numeric variable. In such cases, users often want to computeorthogonal polynomial contrasts:
pigs.emm.p <- emmeans(pigs.lm, "percent")ply <- contrast(pigs.emm.p, "poly")ply
## contrast estimate SE df t.ratio p.value## linear 0.9374 0.2110 23 4.452 0.0002## quadratic -0.0971 0.0883 23 -1.099 0.2830## cubic 0.1863 0.1880 23 0.992 0.3313## ## Results are averaged over the levels of: source ## Results are given on the log (not the response) scale.
coef(ply)
## percent c.1 c.2 c.3## 1 9 -3 1 -1## 2 12 -1 -1 3## 3 15 1 -1 -3## 4 18 3 1 1
We obtain tests for the linear, quadratic, and cubic trends. Thecoefficients are those that can be found in tables in manyexperimental-design texts. It is important to understand that theestimated linear contrast isnot the slope of a line fitted tothe data. It is simply a contrast having coefficients that increaselinearly. Itdoes test the linear trend, however.
There are a number of other named contrast methods, for example"trt.vs.ctrl","eff", and"consec". The"pairwise" and"revpairwise" methods incontrast() are thesame asPairs() andpairs(..., reverse = TRUE). Seehelp("contrast-methods") for details.
Normalized contrasts
Sometimes users like contrast coefficients to be normalized – i.e.,the sum of the squares of each set of contrast coefficients is equal to1. This is possible with any contrast family by specifying"nrmlz" as the method and the desired contrast family asthefamily argument. For example,
nrmlz.emmc(1:4, family = "helmert")
## 2 vs earlier 3 vs earlier 4 vs earlier## 1 -0.7071068 -0.4082483 -0.2886751## 2 0.7071068 -0.4082483 -0.2886751## 3 0.0000000 0.8164966 -0.2886751## 4 0.0000000 0.0000000 0.8660254
contrast(pigs.emm.s, "nrmlz", family = "pairwise")
## contrast estimate SE df t.ratio p.value## fish - soy -0.1929 0.0374 23 -5.153 <0.0001## fish - skim -0.2845 0.0383 23 -7.428 <0.0001## soy - skim -0.0916 0.0375 23 -2.442 0.0570## ## Results are averaged over the levels of: percent ## Results are given on the log (not the response) scale. ## P value adjustment: tukey method for comparing a family of 3 estimates
Note that the above estimates and SEs differ from those in theexample atthe beginning of this vignette, butthe thet ratios andP values are the same, becauseall we have done is rescale those comparisons.
Back to Contents
Formula interface
If you already know what contrasts you will want before callingemmeans(), a quick way to get them is to specify the methodas the left-hand side of the formula in its second argument. Forexample, with theoranges dataset provided in thepackage,
org.aov <- aov(sales1 ~ day + Error(store), data = oranges, contrasts = list(day = "contr.sum"))org.emml <- emmeans(org.aov, consec ~ day)org.emml
## $emmeans## day emmean SE df lower.CL upper.CL## 1 7.87 2.77 29.2 2.21 13.5## 2 7.10 2.77 29.2 1.43 12.8## 3 13.76 2.77 29.2 8.09 19.4## 4 8.04 2.77 29.2 2.37 13.7## 5 12.92 2.77 29.2 7.26 18.6## 6 11.60 2.77 29.2 5.94 17.3## ## Warning: EMMs are biased unless design is perfectly balanced ## Confidence level used: 0.95 ## ## $contrasts## contrast estimate SE df t.ratio p.value## day2 - day1 -0.772 3.78 25 -0.204 0.9997## day3 - day2 6.658 3.78 25 1.763 0.3244## day4 - day3 -5.716 3.78 25 -1.513 0.4679## day5 - day4 4.882 3.78 25 1.293 0.6128## day6 - day5 -1.321 3.78 25 -0.350 0.9965## ## P value adjustment: mvt method for 5 tests
The contrasts shown are the day-to-day changes.
This two-sided formula technique is quite convenient, but it can alsocreate confusion. For one thing, the result is not anemmGrid object anymore; it is alist ofemmGrid objects, called anemm_list. You mayneed to be cognizant of that if you are to do further contrasts or otheranalyzes. For example if you want"eff" contrasts as well,you need to docontrast(org.emml[[1]], "eff") orcontrast(org.emml, "eff", which = 1).
Another issue is that it may be unclear which part of the results isaffected by certain options. For example, if you were to addadjust = "bonf" to theorg.emm call above,would the Bonferroni adjustment be applied to the EMMs, or to thecontrasts? (See the documentation if interested; but the best practiceis to avoid such dilemmas.)
Back to Contents
Custom contrasts and linear functions
The user may write a custom contrast function for use incontrast(). What’s needed is a function having the desiredname with".emmc" appended, that generates the neededcoefficients as a list or data frame. The function should take a vectorof levels as its first argument, and any optional parameters asadditional arguments. It should also always have a...argument to allow for unspecified arguments that may occur in thecall.
As an example, suppose we want to compare every third level of atreatment. The following function provides for this:
skip_comp.emmc <- function(levels, skip = 1, reverse = FALSE, ...) { if((k <- length(levels)) < skip + 1) stop("Need at least ", skip + 1, " levels") coef <- data.frame() coef <- as.data.frame(lapply(seq_len(k - skip - 1), function(i) { sgn <- ifelse(reverse, -1, 1) sgn * c(rep(0, i - 1), 1, rep(0, skip), -1, rep(0, k - i - skip - 1)) })) names(coef) <- sapply(coef, function(x) paste(which(x == 1), "-", which(x == -1))) attr(coef, "adjust") = "fdr" # default adjustment method coef}
To test it, try 5 levels:
skip_comp.emmc(1:5)
## 1 - 3 2 - 4 3 - 5## 1 1 0 0## 2 0 1 0## 3 -1 0 1## 4 0 -1 0## 5 0 0 -1
skip_comp.emmc(1:5, skip = 0, reverse = TRUE)
## 2 - 1 3 - 2 4 - 3 5 - 4## 1 -1 0 0 0## 2 1 -1 0 0## 3 0 1 -1 0## 4 0 0 1 -1## 5 0 0 0 1
(The latter is the same as"consec" contrasts.) Now tryit with theoranges example we had previously:
contrast(org.emml[[1]], "skip_comp", skip = 2, reverse = TRUE)
## contrast estimate SE df t.ratio p.value## 4 - 1 0.17 3.78 25 0.045 0.9645## 5 - 2 5.82 3.78 25 1.542 0.4069## 6 - 3 -2.15 3.78 25 -0.571 0.8601## ## P value adjustment: fdr method for 3 tests
Thecontrast() function may in fact be used to computearbitrary linear functions of EMMs. Suppose for some reason we want toestimate the quantities\(\lambda_1 =\mu_1+2\mu_2-7\) and\(\lambda_2 =3\mu_2-2\mu_3+1\), where the\(\mu_j\) are the population values of thesource EMMs in thepigs example. This may bedone by providing the coefficients in a list, and the added constants intheoffset argument:
LF <- contrast(pigs.emm.s, list(lambda1 = c(1, 2, 0), lambda2 = c(0, 3, -2)), offset = c(-7, 1))confint(LF, adjust = "bonferroni")
## contrast estimate SE df lower.CL upper.CL## lambda1 3.73 0.0827 23 3.53 3.93## lambda2 4.41 0.1340 23 4.09 4.73## ## Results are averaged over the levels of: percent ## Note: contrasts are still on the log scale. Consider using## regrid() if you want contrasts of back-transformed estimates. ## Confidence level used: 0.95 ## Conf-level adjustment: bonferroni method for 2 estimates
Back to Contents
Special properties of log (and logit) transformations
Suppose we obtain EMMs for a model having a response transformationor link function. In most cases, when we compute contrasts of thoseEMMs, there is no natural way to express those contrasts on anythingother than the transformed scale. For example, in a model fitted usingglm() with thegamma() family, the defaultlink function is the inverse. Predictions on such a model are estimatesof\(1/\mu_j\) for various\(j\). Comparisons of predictions will beestimates of\(1/\mu_j - 1/\mu_{k}\)for\(j \ne k\). There is no naturalway to back-transform these differences to some other interpretablescale.
However, logs are an exception, in that\(\log\mu_j - \log\mu_k =\log(\mu_j/\mu_k)\). Accordingly, whencontrast()(orpairs()) notices that the response is on the log scale,it back-transforms contrasts to ratios when results are to be ofresponse type. For example:
pairs(pigs.emm.s, type = "lp")
## contrast estimate SE df t.ratio p.value## fish - soy -0.273 0.0529 23 -5.153 <0.0001## fish - skim -0.402 0.0542 23 -7.428 <0.0001## soy - skim -0.130 0.0530 23 -2.442 0.0570## ## Results are averaged over the levels of: percent ## Results are given on the log (not the response) scale. ## P value adjustment: tukey method for comparing a family of 3 estimates
pairs(pigs.emm.s, type = "response")
## contrast ratio SE df null t.ratio p.value## fish / soy 0.761 0.0403 23 1 -5.153 <0.0001## fish / skim 0.669 0.0362 23 1 -7.428 <0.0001## soy / skim 0.879 0.0466 23 1 -2.442 0.0570## ## Results are averaged over the levels of: percent ## P value adjustment: tukey method for comparing a family of 3 estimates ## Tests are performed on the log scale
As is true of EMM summaries withtype = "response", thetests and confidence intervals are done before back-transforming. Theratios estimated here are actually ratios ofgeometric means.In general, a model with a log response is in fact a model forrelative effects of any of its linear predictors, and thisback-transformation to ratios goes hand-in-hand with that.
In generalized linear models, this behaviors will occur in two commoncases: Poisson or count regression, for which the usual link is the log;and logistic regression, because logits are logs of odds ratios.
Back to Contents
Index of all vignette topics
[8]ページ先頭