# Summaries, predictions, intervals, and tests for `emmGrid`

objects

Source: `R/summary.R`

, `R/test.R`

`summary.emmGrid.Rd`

These 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`

.

## Usage

```
# S3 method for class 'emmGrid'
summary(object, infer, level, adjust, by,
cross.adjust = "none", type, df, calc, null, delta, side, frequentist,
bias.adjust = get_emm_option("back.bias.adj"), sigma, ...)
# S3 method for class 'emmGrid'
confint(object, parm, level = 0.95, ...)
test(object, null, ...)
# S3 method for class 'emmGrid'
test(object, null = 0, joint = FALSE, verbose = FALSE,
rows, by, status = FALSE, ...)
# S3 method for class 'emmGrid'
predict(object, type, interval = c("none", "confidence",
"prediction"), level = 0.95,
bias.adjust = get_emm_option("back.bias.adj"), sigma, ...)
# S3 method for class 'emmGrid'
as.data.frame(x, row.names = NULL, optional,
check.names = TRUE, destroy.annotations = FALSE, ...)
# S3 method for class 'summary_emm'
x[..., as.df = FALSE]
```

## Arguments

- object
An object of class

`"emmGrid"`

(see emmGrid-class)- infer
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.- level
Numerical value between 0 and 1. Confidence level for confidence intervals, if

`infer[1]`

is`TRUE`

.- adjust
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.- by
Character name(s) of variables to use for grouping into separate tables. This affects the family of tests considered in adjusted

*P*values.- cross.adjust
Character: \(p\)-value adjustment method to additionally apply

*across*the`by`

groups. See the section on P-value adjustments for details.- type
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.- df
Numeric. If non-missing, a constant number of degrees of freedom to use in constructing confidence intervals and

*P*values (`NA`

specifies asymptotic results).- calc
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`

.- null
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).- delta
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.- side
Numeric or character value specifying whether the test is left-tailed (

`-1`

,`"-"`

,`"<"`

,`"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.- frequentist
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.- bias.adjust
Logical value for whether to adjust for bias in back-transforming (

`type = "response"`

). This requires a valid value of`sigma`

to exist in the object or be specified.- sigma
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 a warning is issued if it is not found or not valid.*Note:*`sigma`

may be a vector, but be careful that it correctly corresponds (perhaps after recycling) to the order of the reference grid.- ...
Optional arguments such as

`scheffe.rank`

(see “P-value adjustments”). In`confint.emmGrid`

,`predict.emmGrid`

, and`test.emmGrid`

, these arguments are passed to`summary.emmGrid`

.- parm
(Required argument for

`confint`

methods, but not used)- joint
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`

.- verbose
Logical value. If

`TRUE`

and`joint = TRUE`

, a table of the effects being tested is printed.- rows
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.- status
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.- interval
Type of interval desired (partial matching is allowed):

`"none"`

for no intervals, otherwise confidence or prediction intervals with given arguments, via`confint.emmGrid`

. Note: prediction intervals are not available unless the model family is`"gaussian"`

.- x
object of the given class

- row.names
passed to

`as.data.frame`

- optional
required argument, but ignored in

`as.data.frame.emmGrid`

- check.names
passed to

`data.frame`

- destroy.annotations
Logical value. If

`FALSE`

, an object of class`summary_emm`

is returned (which inherits from`data.frame`

), but if displayed, details like confidence levels, P-value adjustments, transformations, etc. are also shown. But unlike the result of`summary`

, the number of digits displayed is obtained from`getOption("digits")`

rather than using the optimal digits algorithm we usually use. Thus, it is formatted more like a regular data frame, but with any annotations and groupings still intact. If`TRUE`

(not recommended), a “plain vanilla” data frame is returned, based on`row.names`

and`check.names`

.- as.df
Logical value. With

`x[..., as.df = TRUE]`

, the result is object is coerced to a`data.frame`

before the subscripting is applied. With`as.df = FALSE`

, the result is returned as a`summary_emm`

object when possible.

## Value

`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 an object that inherits
from `"data.frame"`

.

## Details

`confint.emmGrid`

is equivalent to ```
summary.emmGrid with
infer = c(TRUE, FALSE)
```

. The function `test.emmGrid`

, when called with
`joint = FALSE`

, 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. The returned object has an additional `"est.fcns"`

attribute, which
is a list of the linear functions associated with the joint test.

## Note

In doing testing and a transformation and/or link is in force, any
`null`

and/or `delta`

values specified must always be on the
scale of the linear predictor, regardless of the setting for `type`. If
`type = "response"`

, the null value displayed in the summary table
will be back-transformed from the value supplied by the user. But the
displayed `delta`

will not be changed, because there (often) is
not a natural way to back-transform it.

When we have `type = "response"`

, and `bias.adj = TRUE`

,
the `null`

value displayed in the output is both back-transformed
and bias-adjusted, leading to a rather non-intuitive-looking null value.
However, since the tests themselves are performed on the link scale,
this is the response value at which a *P* value of 1 would be obtained.

The default `show`

method for `emmGrid`

objects (with the
exception of newly created reference grids) is `print(summary())`

.
Thus, with ordinary usage of `emmeans`

and such, it is
unnecessary to call `summary`

unless there is a need to
specify other than its default options.

If a data frame is needed, `summary`

, `confint`

,
and `test`

serve this need. `as.data.frame`

routes to
`summary`

by default; calling it with `destroy.annotations = TRUE`

is not recommended for exactly that reason.
If you want to see more digits in the output, use
`print(summary(object), digits = ...)`

; and if you *always* want
to see more digits, use `emm_options(opt.digits = FALSE)`

.

## Defaults

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.

## Transformations and links

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.

Be aware that only univariate transformations and links are
supported in this way. Some multivariate transformations are supported by
`mvregrid`

.

## Bias adjustment when back-transforming

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.

Note: In certain models, e.g., those with non-gaussian families,
`sigma`

is initialized as `NA`

, and so by default, bias adjustment
is skipped and a warning is issued. You may override this by specifying a
value for `sigma`

. However, *with ordinary generalized linear models,
bias adjustment is inappropriate* and you should not try to do it. With GEEs and GLMMs,
you probably should *not* use `sigma(model)`

, and instead you should create an
appropriate value using the estimated random effects, e.g., from `VarCorr(model)`

.
An example is provided in the “transformations” vignette.

## P-value adjustments

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). If there are non-estimable
cases in a `by`

group, those cases are *excluded* before determining
the adjustment; that means there could be different adjustments in different groups.

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 own

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

*Warning:* Non-estimable cases are *included* in the family to which adjustments
are applied. You may wish to subset the object using the `[]`

operator
to work around this problem.

The `cross.adjust`

argument is a way of specifying a multiplicity
adjustment across the `by`

groups (otherwise by default, each group is
treated as a separate family in regards to multiplicity adjustments). It
applies only to \(p\) values. Valid options are one of the
`p.adjust.methods`

or `"sidak"`

. This argument is ignored unless
it is other than `"none"`

, there is more than one `by`

group, and
they are all the same size. Under those conditions, we first use
`adjust`

to determine the within-group adjusted \(p\) values.
Imagine each group's adjusted \(p\) values arranged in side-by-side
columns, thus forming a matrix with the number of columns equal to the
number of `by`

groups. Then we use the `cross.adjust`

method to
further adjust the adjusted \(p\) values in each row of this matrix. Note
that an *overall* Bonferroni (or Sidak) adjustment is obtainable by
specifying *both* `adjust`

and `cross.adjust`

as
`"bonferroni"`

(or `"sidak"`

). However, less conservative (but
yet conservative) overall adjustments are available when it is possible to
use an “exact” within-group method (e.g., `adjust = "tukey"`

for pairwise comparisons) and `cross.adjust`

as a conservative
adjustment. [`cross.adjust`

methods other than `"none"`

,
`"bonferroni"`

, or `"sidak"`

do not seem advisable, but other
`p.adjust`

methods are available if you can make sense of them.]

## Tests of significance, nonsuperiority, noninferiority, or equivalence

When `delta = 0`

, test statistics are the usual tests of significance.
They are of the form
`(estimate - null)/SE`. Notationally:

- Significance
\(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.

- Left-sided (nonsuperiority)
\(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.- Right-sided (noninferiority)
\(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.- Two-sided (equivalence)
\(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.

## Non-estimable cases

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.

See the warning above that non-estimable cases are still included when
determining the family size for *P*-value adjustments.

## Warning about potential misuse of P values

Some in the statistical and scientific community argue 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 can be too misleading and are prone to abuse. See the “basics” vignette for more discussion.

## Examples

```
warp.lm <- lm(breaks ~ wool * tension, data = warpbreaks)
warp.emm <- emmeans(warp.lm, ~ tension | wool)
warp.emm # implicitly runs 'summary'
#> wool = A:
#> tension emmean SE df lower.CL upper.CL
#> L 44.6 3.65 48 37.2 51.9
#> M 24.0 3.65 48 16.7 31.3
#> H 24.6 3.65 48 17.2 31.9
#>
#> wool = B:
#> tension emmean SE df lower.CL upper.CL
#> L 28.2 3.65 48 20.9 35.6
#> M 28.8 3.65 48 21.4 36.1
#> H 18.8 3.65 48 11.4 26.1
#>
#> Confidence level used: 0.95
confint(warp.emm, by = NULL, level = .90)
#> tension wool emmean SE df lower.CL upper.CL
#> L A 44.6 3.65 48 38.4 50.7
#> M A 24.0 3.65 48 17.9 30.1
#> H A 24.6 3.65 48 18.4 30.7
#> L B 28.2 3.65 48 22.1 34.3
#> M B 28.8 3.65 48 22.7 34.9
#> H B 18.8 3.65 48 12.7 24.9
#>
#> Confidence level used: 0.9
# --------------------------------------------------------------
pigs.lm <- lm(log(conc) ~ source + factor(percent), data = pigs)
pigs.emm <- emmeans(pigs.lm, "percent", type = "response")
summary(pigs.emm) # (inherits type = "response")
#> percent response SE df lower.CL upper.CL
#> 9 31.4 1.28 23 28.8 34.1
#> 12 37.5 1.44 23 34.7 40.6
#> 15 39.0 1.70 23 35.6 42.7
#> 18 42.3 2.24 23 37.9 47.2
#>
#> Results are averaged over the levels of: source
#> Confidence level used: 0.95
#> Intervals are back-transformed from the log scale
summary(pigs.emm, calc = c(n = ".wgt.")) # Show sample size
#> percent response SE df n lower.CL upper.CL
#> 9 31.4 1.28 23 8 28.8 34.1
#> 12 37.5 1.44 23 9 34.7 40.6
#> 15 39.0 1.70 23 7 35.6 42.7
#> 18 42.3 2.24 23 5 37.9 47.2
#>
#> Results are averaged over the levels of: source
#> Confidence level used: 0.95
#> Intervals are back-transformed from the log scale
# 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 = ">")
#> percent response SE df null t.ratio p.value
#> 9 31.4 1.28 23 35 -0.360 0.6390
#> 12 37.5 1.44 23 35 4.295 0.0001
#> 15 39.0 1.70 23 35 4.635 0.0001
#> 18 42.3 2.24 23 35 5.384 <.0001
#>
#> Results are averaged over the levels of: source
#> Statistics are tests of noninferiority with a threshold of 0.09531
#> P values are right-tailed
#> Tests are performed on the log scale
con <- contrast(pigs.emm, "consec")
test(con)
#> contrast ratio SE df null t.ratio p.value
#> percent12 / percent9 1.20 0.0671 23 1 3.202 0.0109
#> percent15 / percent12 1.04 0.0604 23 1 0.650 0.8613
#> percent18 / percent15 1.09 0.0750 23 1 1.194 0.5200
#>
#> Results are averaged over the levels of: source
#> P value adjustment: mvt method for 3 tests
#> Tests are performed on the log scale
test(con, joint = TRUE)
#> df1 df2 F.ratio p.value
#> 3 23 7.981 0.0008
#>
# default Scheffe adjustment - rank = 3
summary(con, infer = c(TRUE, TRUE), adjust = "scheffe")
#> contrast ratio SE df lower.CL upper.CL null t.ratio p.value
#> percent12 / percent9 1.20 0.0671 23 1.011 1.42 1 3.202 0.0343
#> percent15 / percent12 1.04 0.0604 23 0.872 1.24 1 0.650 0.9344
#> percent18 / percent15 1.09 0.0750 23 0.882 1.34 1 1.194 0.7027
#>
#> Results are averaged over the levels of: source
#> Confidence level used: 0.95
#> Conf-level adjustment: scheffe method with rank 3
#> Intervals are back-transformed from the log scale
#> P value adjustment: scheffe method with rank 3
#> Tests are performed on the log scale
# Consider as some of many possible contrasts among the six cell means
summary(con, infer = c(TRUE, TRUE), adjust = "scheffe", scheffe.rank = 5)
#> contrast ratio SE df lower.CL upper.CL null t.ratio p.value
#> percent12 / percent9 1.20 0.0671 23 0.976 1.47 1 3.202 0.1090
#> percent15 / percent12 1.04 0.0604 23 0.841 1.28 1 0.650 0.9940
#> percent18 / percent15 1.09 0.0750 23 0.845 1.40 1 1.194 0.9165
#>
#> Results are averaged over the levels of: source
#> Confidence level used: 0.95
#> Conf-level adjustment: scheffe method with rank 5
#> Intervals are back-transformed from the log scale
#> P value adjustment: scheffe method with rank 5
#> Tests are performed on the log scale
# Show estimates to more digits
print(test(con), digits = 7)
#> contrast ratio SE df null t.ratio p.value
#> percent12 / percent9 1.196684 0.06710564 23 1 3.202 0.0110
#> percent15 / percent12 1.038570 0.06042501 23 1 0.650 0.8613
#> percent18 / percent15 1.085945 0.07499759 23 1 1.194 0.5201
#>
#> Results are averaged over the levels of: source
#> P value adjustment: mvt method for 3 tests
#> Tests are performed on the log scale
# --------------------------------------------------------------
# Cross-adjusting P values
prs <- pairs(warp.emm) # pairwise comparisons of tension, by wool
test(prs, adjust = "tukey", cross.adjust = "bonferroni")
#> wool = A:
#> contrast estimate SE df t.ratio p.value
#> L - M 20.556 5.16 48 3.986 0.0013
#> L - H 20.000 5.16 48 3.878 0.0018
#> M - H -0.556 5.16 48 -0.108 1.0000
#>
#> wool = B:
#> contrast estimate SE df t.ratio p.value
#> L - M -0.556 5.16 48 -0.108 1.0000
#> L - H 9.444 5.16 48 1.831 0.3407
#> M - H 10.000 5.16 48 1.939 0.2777
#>
#> P value adjustment: tukey method for comparing a family of 3 estimates
#> Cross-group P-value adjustment: bonferroni
# Same comparisons taken as one big family (more conservative)
test(prs, adjust = "bonferroni", by = NULL)
#> contrast wool estimate SE df t.ratio p.value
#> L - M A 20.556 5.16 48 3.986 0.0014
#> L - H A 20.000 5.16 48 3.878 0.0019
#> M - H A -0.556 5.16 48 -0.108 1.0000
#> L - M B -0.556 5.16 48 -0.108 1.0000
#> L - H B 9.444 5.16 48 1.831 0.4396
#> M - H B 10.000 5.16 48 1.939 0.3504
#>
#> P value adjustment: bonferroni method for 6 tests
```