Standardized effect sizes are typically calculated using pairwise differences of estimates,
divided by the SD of the population providing the context for those effects.
This function calculates effect sizes from an emmGrid
object,
and confidence intervals for them, accounting for uncertainty in both the estimated
effects and the population SD.
Arguments
- object
an
emmGrid
object, typically one defining the EMMs to be contrasted. If instead,class(object) == "emm_list"
, such as is produced byemmeans(model, pairwise ~ treatment)
, a message is displayed; the contrasts already therein are used; andmethod
is replaced by"identity"
.- sigma
numeric scalar, value of the population SD.
- edf
numeric scalar that specifies the equivalent degrees of freedom for the
sigma
. This is a way of specifying the uncertainty insigma
, in that we regard our estimate ofsigma^2
as being proportional to a chi-square random variable withedf
degrees of freedom. (edf
should not be confused with thedf
argument that may be passed via...
to specify the degrees of freedom to use in \(t\) statistics and confidence intervals.)- method
the contrast method to use to define the effects. This is passed to
contrast
after the elements ofobject
are scaled.- ...
Additional arguments passed to
contrast
Value
an emmGrid
object containing the effect sizes
Details
Any by
variables specified in object
will remain in force in the returned
effects, unless overridden in the optional arguments.
For models having a single random effect, such as those fitted using
lm
; in that case, the stats::sigma
and
stats::df.residual
functions may be useful for specifying sigma
and edf
. For models with more than one random effect, sigma
may
be based on some combination of the random-effect variances.
Specifying edf
can be rather unintuitive but is also relatively
uncritical; but the smaller the value, the wider the confidence intervals for
effect size. The value of sqrt(2/edf)
can be interpreted as the
relative accuracy of sigma
; for example, with edf = 50
,
\(\sqrt(2/50) = 0.2\), meaning that sigma
is accurate to plus or
minus 20 percent. Note in an example below, we tried two different edf
values as kind of a bracketing/sensitivity-analysis strategy. A value of
Inf
is allowable, in which case you are assuming that sigma
is
known exactly. Obviously, this narrows the confidence intervals for the
effect sizes – unrealistically if in fact sigma
is unknown.
Note
The effects are always computed on the scale of the linear-predictor;
any response transformation or link function is completely ignored. If you
wish to base the effect sizes on the response scale, it is not enough
to replace object
with regrid(object)
, because this
back-transformation changes the SD required to compute effect sizes.
Paired data: Be careful with paired-data situations, where Cohen's d is typically referenced to
the SD of the paired differences rather than the residual SD.
You may need to enlarge sigma
by a factor of sqrt(2)
to obtain
comparable results with other software.
Disclaimer: There is substantial disagreement among practitioners on
what is the appropriate sigma
to use in computing effect sizes; or,
indeed, whether any effect-size measure is appropriate for some
situations. The user is completely responsible for specifying
appropriate parameters (or for failing to do so).
Cohen effect sizes do not even exist for generalized linear models or other models lacking an additive residual error term.
The examples here illustrate a sobering message that effect sizes are often not nearly as accurate as you may think.
Computation
This function uses calls to regrid
to put the estimated
marginal means (EMMs) on the log scale. Then an extra element is added to
this grid for the log of sigma
and its standard error (where we assume
that sigma
is uncorrelated with the log EMMs). Then a call to
contrast
subtracts log{sigma}
from each of the log EMMs,
yielding values of log(EMM/sigma)
.
Finally, the results are re-gridded back to the original scale and the
desired contrasts are computed using method
. In the log-scaling
part, we actually rescale the absolute values and keep track of the signs.
Examples
fiber.lm <- lm(strength ~ diameter + machine, data = fiber)
emm <- emmeans(fiber.lm, "machine")
eff_size(emm, sigma = sigma(fiber.lm), edf = df.residual(fiber.lm))
#> contrast effect.size SE df lower.CL upper.CL
#> A - B -0.650 0.650 11 -2.081 0.781
#> A - C 0.993 0.726 11 -0.604 2.590
#> B - C 1.643 0.800 11 -0.118 3.405
#>
#> sigma used for effect sizes: 1.595
#> Confidence level used: 0.95
# or equivalently:
eff_size(pairs(emm), sigma(fiber.lm), df.residual(fiber.lm), method = "identity")
#> contrast effect.size SE df lower.CL upper.CL
#> (A - B) -0.650 0.650 11 -2.081 0.781
#> (A - C) 0.993 0.726 11 -0.604 2.590
#> (B - C) 1.643 0.800 11 -0.118 3.405
#>
#> sigma used for effect sizes: 1.595
#> Confidence level used: 0.95
### Mixed model example:
if (require(nlme)) withAutoprint({
Oats.lme <- lme(yield ~ Variety + factor(nitro),
random = ~ 1 | Block / Variety,
data = Oats)
# Combine variance estimates
VarCorr(Oats.lme)
(totSD <- sqrt(214.4724 + 109.6931 + 162.5590))
# I figure edf is somewhere between 5 (Blocks df) and 51 (Resid df)
emmV <- emmeans(Oats.lme, ~ Variety)
eff_size(emmV, sigma = totSD, edf = 5)
eff_size(emmV, sigma = totSD, edf = 51)
}, spaced = TRUE)
#> Loading required package: nlme
#>
#> > Oats.lme <- lme(yield ~ Variety + factor(nitro), random = ~1 | Block/Variety,
#> + data = Oats)
#>
#> > VarCorr(Oats.lme)
#> Variance StdDev
#> Block = pdLogChol(1)
#> (Intercept) 214.4722 14.64487
#> Variety = pdLogChol(1)
#> (Intercept) 109.6928 10.47343
#> Residual 162.5591 12.74987
#>
#> > (totSD <- sqrt(214.4724 + 109.6931 + 162.559))
#> [1] 22.06183
#>
#> > emmV <- emmeans(Oats.lme, ~Variety)
#>
#> > eff_size(emmV, sigma = totSD, edf = 5)
#> contrast effect.size SE df lower.CL upper.CL
#> Golden Rain - Marvellous -0.240 0.330 5 -1.087 0.608
#> Golden Rain - Victory 0.312 0.336 5 -0.551 1.174
#> Marvellous - Victory 0.551 0.365 5 -0.387 1.490
#>
#> Results are averaged over the levels of: nitro
#> sigma used for effect sizes: 22.06
#> Degrees-of-freedom method: inherited from containment when re-gridding
#> Confidence level used: 0.95
#>
#> > eff_size(emmV, sigma = totSD, edf = 51)
#> contrast effect.size SE df lower.CL upper.CL
#> Golden Rain - Marvellous -0.240 0.322 5 -1.067 0.587
#> Golden Rain - Victory 0.312 0.322 5 -0.517 1.140
#> Marvellous - Victory 0.551 0.325 5 -0.285 1.388
#>
#> Results are averaged over the levels of: nitro
#> sigma used for effect sizes: 22.06
#> Degrees-of-freedom method: inherited from containment when re-gridding
#> Confidence level used: 0.95