Replace compIntercepts() with emmeans()
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.
Introduction
The compIntercepts()
function in FSA
(prior to v0.9.0) was used to statistically compare intercepts for all pairs of groups with the same slope in an indicator/dummy variable regression (I/DVR). However, the excellent emmeans()
function in the emmmeans
package is a more general approach that follows principals similar to those of emtrends()
, which I demonstrated in a post yesterday. As such, compIntercepts()
will be removed from the next version (0.9.0) of FSA
.
In this post I demonstrate how to use emmeans()
for the same purpose as compIntercepts()
was used (prior to FSA v0.9.0). Note, however, that the results will not be identical because compSlopes()
and emtrends()
use different methods to correct for multiple comparisons when comparing pairs of slopes.
The Mirex
data described here, but filtered to include just three years, will be used in this post.^{1}
The lm()
below fits the I/DVR to determine if the relationship between mirex concentration and weight of the salmon differs by year.
The weight:year
interaction term pvalue suggests that the slopes (i.e., relationship between mirex concentration and salmon weight) do NOT differ among the three years. However, the year
term pvalue suggests that the intercepts of at least one pair of these parallel lines DO differ.^{2}
A next step is to determine which pair(s) of intercepts differ significantly. In preparation for this next step, I fit a model that does not include the insignificant interaction term.^{3}
compIntercepts()
from FSA
The procedure coded in compItercepts()
is described in more detail in this supplement to the Introductory Fisheries Analyses with R book. The results of compIntercepts()
applied to the saved lm()
object are assigned to an object below.^{4}
The $comparisons
component in this saved object contains the results from comparing all pairs of intercepts. Each paired comparison is a row in these results with the groups being compared under comparison
, the differences in sample intercepts under diff
, 95% confidence intervals for the difference in intercepts under 95% LCI
and 95% UCI
, and adjusted (for multiple comparisons) pvalues for the hypothesis test comparing the intercepts under p.adj
.
For example, these results suggest that the intercepts for 1982 and 1977 ARE statistically different (first row), but the intercepts for 1986 and 1982 are NOT statistically different (last row).
The $smeans
component in this saved object contains the mean value of the response variable predicted at the mean value of the covariate. For example, the results below show the predicted mean mirex concentration at the overall mean salmon weight (i.e., 3.782083 kg).
Because the lines are known to be parallel, differences in intercepts also represent differences in predicted means of the response at all other values of the covariate. compIntercepts()
defaulted to show these means at the mean (i.e., center) of the covariate. This could be adjusted with common.cov=
in compIntercepts()
. For example, the actual intercepts are shown below.
emmeans()
from emmeans
Similar results can be obtained with emmeans()
from emmeans
using the fitted lm()
object (without the interaction term) as the first argument and a specs=
argument with pairwise~
followed by the name of the factor variable from the lm()
model (year
in this case).
The object saved from emmeans()
is then given as the first argument to summary()
, which also requires infer=TRUE
if you would like pvalues to be calculated.^{5}
The $contrasts
component in this saved object contains the results for comparing all pairs of predicted means at the overall mean of the covariate. Each paired comparison is a row with the groups compared under contrast
, the difference in predicted means under estimate
, the standard error of the difference in predicted means under SE
, the degreesoffreedom under df
, a 95% confidence interval for the difference in predicted means under lower.CL
and upper.CL
, and the t test statistic and pvalue adjusted for multiple comparisons for testing a difference in predicted means under t.ratio
and p.value
, respectively.
Comparing these results to the $comparison
component from compIntercepts()
shows that the difference in sample intercepts or predicted means are the same, though the signs differ because the subtraction was reversed. The confidence interval values and pvalues are slightly different. Again, this is due to emmeans()
and compIntercepts()
using different methods of adjusting for multiple comparisons. These differences did not result in different conclusions in this case, but they could, especially if the pvalues are near the rejection criterion.
The $emmeans
component contains results for predicted means for each group with the groups under the name of the factor variable (year
in this example), the predicted means under emmean
, standard errors of the predicted means under SE
, degreesoffreedom under df
, 95% confidence intervals for the predicted mean under lower.CL
and upper.CL
, and t test statistics and pvalues adjusted for multiple comparisons for testing that the predicted mean is not equal to zero under t.ratio
and p.adj
, respectively. While it is not obvious here, these predict means of the response variable are at the mean of the covariate, as they were for compIntercepts()
.
Here the predicted means match exactly (within rounding) with those in the $means
component of compIntercepts()
.
The means can be predicted at any other “summary” value of the covariate using cov.reduce=
in emmeans()
. For example, the predicted values at the minimum value of the covariate are obtained below.
The following will compute predicted means that represent the actual intercepts.
Conclusion
The emmeans()
function in the emmeans
package provides a more general solution to comparing multiple intercepts (or predicted means on parallel lines) than what was used in compIntercepts()
in the FSA
package (prior to v0.9.0). Thus, compIntercepts()
will be removed from FSA with v0.9.0. You should use emmeans()
instead.
The emmeans
package has extensive vignettes that further explain its use. Their “Basics” vignette is very useful.
Note that this change to FSA
does not affect anything in the published Introductory Fisheries Analyses with R book. However, the specific analysis in this supplement to the book will no longer work as described. The use of compIntercepts()
there will need to be replaced with emmeans()
.
In a previous post I demonstrated how to use emtrends()
from the emmeans
package to replace compSlopes()
, which will also be removed from the next version of FSA
.
Footnotes

I filtered the data to only three years to reduce the amount of output below to make it easier to follow the main concepts. ↩

The
weight
term pvalue suggests that there is a significant relationsip between mirex concentration and salmon weight, regardless of which year is considered. ↩ 
Note that use of the additive ‘+’ in this model formula rather than the multiplicative ‘*’. ↩

compIntercepts()
had aprint()
function for nicely printing the results. However, here we will look at each component separately to ease comparison with theemtrends()
results. ↩ 
emmeans
does not compute pvalues by default. ↩
Rbloggers.com offers daily email updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/datascience job.
Want to share your content on Rbloggers? click here if you have a blog, or here if you don't.