emmGrid objectsThese are the primary methods for obtaining numerical or tabular results
from an emmGrid object.
summary.emmGrid is the general function for summarizing emmGrid objects.
It also serves as the print method for these objects; so for convenience,
summary() arguments may be included in calls to functions such as
emmeans and contrast that construct emmGrid
objects. Note that by default, summaries for Bayesian models are
diverted to hpd.summary.
# S3 method for emmGrid
summary(object, infer, level, adjust, by, type, df, calc,
null, delta, side, frequentist,
bias.adjust = get_emm_option("back.bias.adj"), sigma, ...)# S3 method for emmGrid
confint(object, parm, level = 0.95, ...)
test(object, null, ...)
# S3 method for emmGrid
test(object, null = 0, joint = FALSE, verbose = FALSE,
rows, by, status = FALSE, ...)
# S3 method for emmGrid
predict(object, type, interval = c("none", "confidence",
"prediction"), level = 0.95,
bias.adjust = get_emm_option("back.bias.adj"), sigma, ...)
# S3 method for emmGrid
as.data.frame(x, row.names = NULL, optional = FALSE, ...)
# S3 method for summary_emm
[(x, ..., as.df = TRUE)
An object of class "emmGrid" (see emmGrid-class)
A vector of one or two logical values. The first determines whether confidence intervals are displayed, and the second determines whether t tests and P values are displayed. If only one value is provided, it is used for both.
Numerical value between 0 and 1. Confidence level for confidence
intervals, if infer[1] is TRUE.
Character value naming the method used to adjust \(p\) values
or confidence limits; or to adjust comparison arrows in plot. See
the P-value adjustments section below.
Character name(s) of variables to use for grouping into separate tables. This affects the family of tests considered in adjusted P values.
Character: type of prediction desired. This only has an effect if
there is a known transformation or link function. "response"
specifies that the inverse transformation be applied. "mu" (or
equivalently, "unlink") is usually the same as "response",
but in the case where the model has both a link function and a response
transformation, only the link part is back-transformed. Other valid values
are "link", "lp", and "linear.predictor"; these are
equivalent, and request that results be shown for the linear predictor,
with no back-transformation. The default is "link", unless the
"predict.type" option is in force; see emm_options,
and also the section below on transformations and links.
Numeric. If non-missing, a constant number of degrees of freedom to
use in constructing confidence intervals and P values (NA
specifies asymptotic results).
Named list of character value(s) or formula(s).
The expressions in char are evaluated and appended to the
summary, just after the df column. The expression may include
any names up through df in the summary, any additional names in
object@grid (such as .wgt. or .offset.), or any
earlier elements of calc.
Numeric. Null hypothesis value(s), on the linear-predictor scale,
against which estimates are tested. May be a single value used for all, or
a numeric vector of length equal to the number of tests in each family
(i.e., by group in the displayed table).
Numeric value (on the linear-predictor scale). If zero, ordinary
tests of significance are performed. If positive, this specifies a
threshold for testing equivalence (using the TOST or two-one-sided-test
method), non-inferiority, or non-superiority, depending on side. See
Details for how the test statistics are defined.
Numeric or character value specifying whether the test is
left-tailed (-1, "-", code"<", "left", or
"nonsuperiority"); right-tailed (1, "+", ">",
"right", or "noninferiority"); or two-sided (0,
2, "!=", "two-sided", "both",
"equivalence", or "="). See the special section below for
more details.
Ignored except if a Bayesian model was fitted. If missing
or FALSE, the object is passed to hpd.summary. Otherwise,
a logical value of TRUE will have it return a frequentist summary.
Logical value for whether to adjust for bias in
back-transforming (type = "response"). This requires a value of
sigma to exist in the object or be specified.
Error SD assumed for bias correction (when
type = "response" and a transformation
is in effect), or for constructing prediction intervals. If not specified,
object@misc$sigma is used, and an error is thrown if it is not found.
Note: sigma may be a vector, as long as it conforms to the number of rows
of the reference grid.
Optional arguments such as scheffe.rank
(see “P-value adjustments”).
In as.data.frame.emmGrid, confint.emmGrid,
predict.emmGrid, and
test.emmGrid, these arguments are passed to
summary.emmGrid.
(Required argument for confint methods, but not used)
Logical value. If FALSE, the arguments are passed to
summary.emmGrid with infer=c(FALSE, TRUE). If joint =
TRUE, a joint test of the hypothesis L beta = null is performed, where L
is object@linfct and beta is the vector of fixed effects estimated
by object@betahat. This will be either an F test or a
chi-square (Wald) test depending on whether degrees of freedom are
available. See also joint_tests.
Logical value. If TRUE and joint = TRUE, a table
of the effects being tested is printed.
Integer values. The rows of L to be tested in the joint test. If
missing, all rows of L are used. If not missing, by variables are
ignored.
logical. If TRUE, a note column showing status
flags (for rank deficiencies and estimability issues) is displayed even
when empty. If FALSE, the column is included only if there are
such issues.
Type of interval desired (partial matching is allowed):
"none" for no intervals,
otherwise confidence or prediction intervals with given arguments,
via confint.emmGrid.
object of the given class
passed to as.data.frame
passed to as.data.frame
Logical value. With x[..., as.df = TRUE], the result is
object is coerced to an ordinary data.frame; otherwise, it is left as a
summary_emm object.
summary.emmGrid, confint.emmGrid, and
test.emmGrid return an object of class "summary_emm", which
is an extension of data.frame but with a special print
method that displays it with custom formatting. For models fitted using
MCMC methods, the call is diverted to hpd.summary (with
prob set to level, if specified); one may
alternatively use general MCMC summarization tools with the
results of as.mcmc.
predict returns a vector of predictions for each row of object@grid.
The as.data.frame method returns a plain data frame,
equivalent to as.data.frame(summary(.)).
The misc slot in object may contain default values for
by, calc, infer, level, adjust,
type, null, side, and delta.
These defaults vary depending
on the code that created the object. The update method may be
used to change these defaults. In addition, any options set using
emm_options(summary = ...) will trump those stored in the object's
misc slot.
With type = "response", the transformation assumed can be found in
object@misc$tran, and its label, for the summary is in
object@misc$inv.lbl. Any \(t\) or \(z\) tests are still performed
on the scale of the linear predictor, not the inverse-transformed one.
Similarly, confidence intervals are computed on the linear-predictor scale,
then inverse-transformed.
When bias.adjust is TRUE, then back-transformed estimates
are adjusted by adding
\(0.5 h''(u)\sigma^2\), where \(h\) is the inverse transformation and
\(u\) is the linear predictor. This is based on a second-order Taylor
expansion. There are better or exact adjustments for certain specific
cases, and these may be incorporated in future updates.
The adjust argument specifies a multiplicity adjustment for tests or
confidence intervals. This adjustment always is applied separately
to each table or sub-table that you see in the printed output (see
rbind.emmGrid for how to combine tables).
The valid values of adjust are as follows:
"tukey"Uses the Studentized range distribution with the number of means in the family. (Available for two-sided cases only.)
"scheffe"Computes \(p\) values from the \(F\)
distribution, according to the Scheffe critical value of
\(\sqrt{rF(\alpha; r, d)}\), where \(d\) is
the error degrees of freedom and \(r\) is the rank of the set of linear
functions under consideration. By default, the value of r is
computed from object@linfct for each by group; however, if the
user specifies an argument matching scheffe.rank, its value will
be used instead. Ordinarily, if there are \(k\) means involved, then
\(r = k - 1\) for a full set of contrasts involving all \(k\) means, and
\(r = k\) for the means themselves. (The Scheffe adjustment is available
for two-sided cases only.)
"sidak"Makes adjustments as if the estimates were independent (a conservative adjustment in many cases).
"bonferroni"Multiplies \(p\) values, or divides significance levels by the number of estimates. This is a conservative adjustment.
"dunnettx"Uses our ownad hoc approximation to the
Dunnett distribution for a family of estimates having pairwise
correlations of \(0.5\) (as is true when comparing treatments with a
control with equal sample sizes). The accuracy of the approximation
improves with the number of simultaneous estimates, and is much faster
than "mvt". (Available for two-sided cases only.)
"mvt"Uses the multivariate \(t\) distribution to assess the
probability or critical value for the maximum of \(k\) estimates. This
method produces the same \(p\) values and intervals as the default
summary or confint methods to the results of
as.glht. In the context of pairwise comparisons or comparisons
with a control, this produces “exact” Tukey or Dunnett adjustments,
respectively. However, the algorithm (from the mvtnorm package) uses a
Monte Carlo method, so results are not exactly repeatable unless the same
random-number seed is used (see set.seed). As the family
size increases, the required computation time will become noticeable or even
intolerable, making the "tukey", "dunnettx", or others more
attractive.
"none"Makes no adjustments to the \(p\) values.
For tests, not confidence intervals, the Bonferroni-inequality-based adjustment
methods in p.adjust are also available (currently, these
include "holm", "hochberg", "hommel",
"bonferroni", "BH", "BY", "fdr", and
"none"). If a p.adjust.methods method other than
"bonferroni" or "none" is specified for confidence limits, the
straight Bonferroni adjustment is used instead. Also, if an adjustment method
is not appropriate (e.g., using "tukey" with one-sided tests, or with
results that are not pairwise comparisons), a more appropriate method
(usually "sidak") is substituted.
In some cases, confidence and \(p\)-value adjustments are only approximate
-- especially when the degrees of freedom or standard errors vary greatly
within the family of tests. The "mvt" method is always the correct
one-step adjustment, but it can be very slow. One may use
as.glht with methods in the multcomp package to obtain
non-conservative multi-step adjustments to tests.
When delta = 0, test statistics are the usual tests of significance.
They are of the form
(estimate - null)/SE. Notationally:
\(H_0: \theta = \theta_0\) versus
\(H_1: \theta < \theta_0\) (left-sided), or
\(H_1 \theta > \theta_0\) (right-sided), or
\(H_1: \theta \ne \theta_0\) (two-sided)
The test statistic is
\(t = (Q - \theta_0)/SE\)
where \(Q\) is our estimate of \(\theta\);
then left, right, or two-sided \(p\) values are produced,
depending on side.
When delta is positive, the test statistic depends on side as
follows.
\(H_0: \theta \ge \theta_0 + \delta\) versus \(H_1: \theta < \theta_0 + \delta\) \(t = (Q - \theta_0 - \delta)/SE\) The \(p\) value is the lower-tail probability.
\(H_0: \theta \le \theta_0 - \delta\) versus \(H_1: \theta > \theta_0 - \delta\) \(t = (Q - \theta_0 + \delta)/SE\) The \(p\) value is the upper-tail probability.
\(H_0: |\theta - \theta_0| \ge \delta\) versus \(H_1: |\theta - \theta_0| < \delta\) \(t = (|Q - \theta_0| - \delta)/SE\) The \(p\) value is the lower-tail probability. Note that \(t\) is the maximum of \(t_{nonsup}\) and \(-t_{noninf}\). This is equivalent to choosing the less significant result in the two-one-sided-test (TOST) procedure.
When the model is rank-deficient, each row x of object's
linfct slot is checked for estimability. If sum(x*bhat)
is found to be non-estimable, then the string NonEst is displayed for the
estimate, and associated statistics are set to NA.
The estimability check is performed
using the orthonormal basis N in the nbasis slot for the null
space of the rows of the model matrix. Estimability fails when
\(||Nx||^2 / ||x||^2\) exceeds tol, which by default is
1e-8. You may change it via emm_options by setting
estble.tol to the desired value.
A growing consensus in the statistical and scientific community is that the term “statistical significance” should be completely abandoned, and that criteria such as “p < 0.05” never be used to assess the importance of an effect. These practices are just too misleading and prone to abuse. See the “basics” vignette for more discussion.
confint.emmGrid is equivalent to summary.emmGrid with
infer = c(TRUE, FALSE). When called with joint = FALSE, test.emmGrid
is equivalent to summary.emmGrid with infer = c(FALSE, TRUE).
With joint = TRUE, test.emmGrid calculates the Wald test of the
hypothesis linfct %*% bhat = null, where linfct and
bhat refer to slots in object (possibly subsetted according to
by or rows). An error is thrown if any row of linfct is
non-estimable. It is permissible for the rows of linfct to be linearly
dependent, as long as null == 0, in which case a reduced set of
contrasts is tested. Linear dependence and nonzero null cause an
error.
# NOT RUN {
warp.lm <- lm(breaks ~ wool * tension, data = warpbreaks)
warp.emm <- emmeans(warp.lm, ~ tension | wool)
warp.emm # implicitly runs 'summary'
confint(warp.emm, by = NULL, level = .90)
# --------------------------------------------------------------
pigs.lm <- lm(log(conc) ~ source + factor(percent), data = pigs)
pigs.emm <- emmeans(pigs.lm, "percent", type = "response")
summary(pigs.emm) # (inherits type = "response")
summary(pigs.emm, calc = c(n = ".wgt.")) # Show sample size
# For which percents is EMM non-inferior to 35, based on a 10% threshold?
# Note the test is done on the log scale even though we have type = "response"
test(pigs.emm, null = log(35), delta = log(1.10), side = ">")
con <- contrast(pigs.emm, "consec")
test(con)
test(con, joint = TRUE)
# default Scheffe adjustment - rank = 3
summary(con, infer = c(TRUE, TRUE), adjust = "scheffe")
# Consider as some of many possible contrasts among the six cell means
summary(con, infer = c(TRUE, TRUE), adjust = "scheffe", scheffe.rank = 5)
# Show estimates to more digits
print(test(con), digits = 7)
# }
Run the code above in your browser using DataLab