The Problem

You have two predictors in your model. One seems to have a stronger coefficient than the other. But is it significant?

Example: when predicting a worker’s salary, is the standardized coefficient of number of extra hours (xtra_hours) really larger than that of number of compliments given the to boss n_comps?

library(parameters)
library(effectsize)

data("hardlyworking", package = "effectsize")

hardlyworkingZ <- standardize(hardlyworking)

m <- lm(salary ~ xtra_hours + n_comps, data = hardlyworkingZ)

model_parameters(m)

#> Parameter   | Coefficient |   SE |        95% CI |   t(497) |      p
#> --------------------------------------------------------------------
#> (Intercept) |    2.76e-16 | 0.02 | [-0.03, 0.03] | 1.59e-14 | > .999
#> xtra hours  |        0.81 | 0.02 | [ 0.78, 0.84] |    46.60 | < .001
#> n comps     |        0.41 | 0.02 | [ 0.37, 0.44] |    23.51 | < .001

Here are 4 6 methods to test coefficient equality in R.

Note

If we were interested in the unstandardized coefficient, we would not need to first standardize the data.
Note that if one parameter was positive, and the other was negative, one of the terms would need to be first reversed (-X) to make this work.

Method 1: As Model Comparisons

Based on this awesome tweet.

Since:

\[ \hat{Y} = a \times \text{xtra_hours} + a \times \text{n_comps} \\ = a \times (\text{xtra_hours} + \text{n_comps}) \]

We can essentially force a constraint on the coefficient to be equal by using a composite variable.

m0 <- lm(salary ~ I(xtra_hours + n_comps), data = hardlyworkingZ)

model_parameters(m0)

#> Parameter            | Coefficient |   SE |        95% CI |   t(498) |      p
#> -----------------------------------------------------------------------------
#> (Intercept)          |    2.80e-16 | 0.02 | [-0.04, 0.04] | 1.31e-14 | > .999
#> xtra hours + n comps |        0.61 | 0.01 | [ 0.58, 0.64] |    41.09 | < .001

We can then compare how this model compares to our first model without this constraint:

anova(m0, m)

#> Analysis of Variance Table
#> 
#> Model 1: salary ~ I(xtra_hours + n_comps)
#> Model 2: salary ~ xtra_hours + n_comps
#>   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
#> 1    498 113.662                                  
#> 2    497  74.942  1     38.72 256.78 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can conclude that the unconstrained model is significantly better than the constrained model - meaning that \(\beta_{\text{xtra_hours}} > \beta_{\text{n_comps}}\).

Method 1b: Composite Variable + Difference

Edit (Feb-17, 2021): Thanks to @joejps84 for pointing this out!

We can achieve the same thing in a single model. If we say that the slope of n_comps as equal to the slope of xtra_hours + some change:

\[ \hat{Y} = a \times \text{xtra_hours} + (a + b) \times \text{n_comps} \\ = a \times \text{xtra_hours} + a \times \text{n_comps} + b \times \text{n_comps} \\ = a \times (\text{xtra_hours} + \text{n_comps}) + b \times \text{n_comps} \]

So the slope of n_comps in such a model is the difference between the slope of n_comps and xtra_hours:

m1 <- lm(salary ~ I(xtra_hours + n_comps) + n_comps, data = hardlyworkingZ)

model_parameters(m1)

#> Parameter            | Coefficient |   SE |         95% CI |   t(497) |      p
#> ------------------------------------------------------------------------------
#> (Intercept)          |    2.56e-16 | 0.02 | [-0.03,  0.03] | 1.47e-14 | > .999
#> xtra hours + n comps |        0.81 | 0.02 | [ 0.78,  0.84] |    46.60 | < .001
#> n comps              |       -0.40 | 0.03 | [-0.45, -0.35] |   -16.02 | < .001

This give the exact same result (\(t^2 = (-16)^2 = 256\) is the same as the F-value from the anova() test)!

Method 2: Paternoster et al (1998)

According to Paternoster et al. (1998), we can compute a t-test to compare the coefficients:

b <- coef(m)
V <- vcov(m)

tibble::tibble(
  diff_estim = b[2] - b[3],
  diff_SE = sqrt(V[2, 2] + V[3, 3] - 2 * V[2, 3]),
  t_stat = diff_estim / diff_SE,
  df = df.residual(m),
  p_value = 2 * pt(abs(t_stat), df = df, lower.tail = FALSE)
)

#> # A tibble: 1 × 5
#>   diff_estim diff_SE t_stat    df  p_value
#>        <dbl>   <dbl>  <dbl> <int>    <dbl>
#> 1      0.402  0.0251   16.0   497 6.83e-47

This gives the exact same results as the first method!

Edit (Feb-17, 2021): As @bmwiernik pointed out, this can also be done with some fancy matrix multiplication:

contr <- c(0, 1, -1)

tibble::tibble(
  diff_estim = b %*% contr,
  diff_SE = sqrt(contr %*% V %*% contr),
  t_stat = diff_estim / diff_SE,
  df = df.residual(m),
  p_value = 2 * pt(abs(t_stat), df = df, lower.tail = FALSE)
)

#> # A tibble: 1 × 5
#>   diff_estim[,1] diff_SE[,1] t_stat[,1]    df p_value[,1]
#>            <dbl>       <dbl>      <dbl> <int>       <dbl>
#> 1          0.402      0.0251       16.0   497    6.83e-47

All of the following solutions are essentially this method, wrapped in some nice functions.

Method 3: `emmeans` <3

We can estimate the slopes from the model using emmeans, and then, with some trickery, compare them!

library(emmeans)

trends <- rbind(
  emtrends(m, ~1, "xtra_hours"),
  emtrends(m, ~1, "n_comps")
)

# clean up so it does not error later
trends@grid$`1` <- c("xtra_hours", "n_comps")
trends@misc$estName <- "trend"

trends

#>  1          trend     SE  df lower.CL upper.CL
#>  xtra_hours 0.811 0.0174 497    0.772    0.850
#>  n_comps    0.409 0.0174 497    0.370    0.448
#> 
#> Confidence level used: 0.95 
#> Conf-level adjustment: bonferroni method for 2 estimates

Compare them:

pairs(trends)

#>  contrast             estimate     SE  df t.ratio p.value
#>  xtra_hours - n_comps    0.402 0.0251 497  16.024  <.0001

Once again - exact same results!

Method 4: `lavaan`

The lavaan package for latent variable analysis and structural equation modeling allows for parameter constraining. So let’s do just that:

library(lavaan)

L <- sem("salary ~ xtra_hours + n_comps", data = hardlyworkingZ)
L0 <- sem("salary ~ a * xtra_hours + a * n_comps", data = hardlyworkingZ)

anova(L0, L)

#> 
#> Chi-Squared Difference Test
#> 
#>    Df    AIC    BIC  Chisq Chisq diff   RMSEA Df diff Pr(>Chisq)    
#> L   0 475.99 488.63   0.00                                          
#> L0  1 682.25 690.68 208.26     208.26 0.64383       1  < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This too yields similar results! (Only slightly different due to different estimation methods.)

We can also directly estimate the difference using lavaan’s “defined parameters” - defines with the := operator:

L <- sem("salary ~ b1 * xtra_hours + b2 * n_comps
         diff := b1 - b2", 
         data = hardlyworkingZ)

parameterEstimates(L, output = "text")[7,]

#> 
#> Defined Parameters:
#>                    Estimate  Std.Err  z-value  P(>|z|) ci.lower ci.upper
#>     diff              0.402    0.025   16.073    0.000    0.353    0.451

Which, again, gives the same results.

Method 5: `car`

Edit (Feb-17, 2021): Thanks to @DouglasKGAraujo for his suggestion!

Even more methods!

library(car)

linearHypothesis(m, c("xtra_hours - n_comps"))

#> Linear hypothesis test
#> 
#> Hypothesis:
#> xtra_hours - n_comps = 0
#> 
#> Model 1: restricted model
#> Model 2: salary ~ xtra_hours + n_comps
#> 
#>   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
#> 1    498 113.662                                  
#> 2    497  74.942  1     38.72 256.78 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Which once again gives the same result!

Method 6: `multcomp`

Edit (Feb-17, 2021): Thanks to Stefan Th. Gries for his suggestion!

Even more more methods!

library(multcomp)

summary(glht(m, matrix(c(0, 1, -1), nrow = 1)))

#> 
#>   Simultaneous Tests for General Linear Hypotheses
#> 
#> Fit: lm(formula = salary ~ xtra_hours + n_comps, data = hardlyworkingZ)
#> 
#> Linear Hypotheses:
#>        Estimate Std. Error t value Pr(>|t|)    
#> 1 == 0  0.40171    0.02507   16.02   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> (Adjusted p values reported -- single-step method)

Which once again again gives the same result!

Summary

That’s it really - these 4 6 simple methods have wide applications to GL(M)M’s, SEM, and more.

Enjoy!

The Problem

Method 1: As Model Comparisons

Method 1b: Composite Variable + Difference

Method 2: Paternoster et al (1998)

Method 3: emmeans <3

Method 4: lavaan

Method 5: car

Method 6: multcomp

Summary

Method 3: `emmeans` <3

Method 4: `lavaan`

Method 5: `car`

Method 6: `multcomp`