Comparison of β of two linear regressions

Question

In my empirical research I use data from two waves of a survey: wave $i$ and wave $j$. My aim is to find out whether there is a difference in the magnitude of the influence of certain explanatory variables on the dependent variable. To test this question I run a separate regression for each wave ($i$ and $j$):
\begin{align} Y_i &= β_{1i}X_{1i} + β_{2i}X_{2i} \tag{1} \\ Y_j &= β_{1j}X_{1j} + β_{2j}X_{2j} \tag{2} \end{align} Next I test the significance of ($β_{1i} - β_{1j}$) and ($β_{2i} – β_{2j}$). I would like to make sure that this is the right way to test the significance of the changes in the influence.

Answer 1

This is straightforward using SEM software, like lavaan, to run two separate simultaneous regression models with and without constraints. You can use a Wald test to test whether individual parameters (e.g., the difference in coefficients) is different from 0, or you can do a likelihood ratio test to compare constrained and unconstrained models.

The unconstrained model:

model1 <- '
    Yi ~ ai*X1i + bi*X2i
    Yj ~ aj*X1j + bj*X2j
    Yi~~Yj
    adiff := aj - ai
    bdiff := bj - bi'
fit1 <- sem(model1, data = data)
summary(fit1)

The tests for adiff and bdiff are tests of whether the individual parameter estimates differ between the two models. These tests allow the residual variance of each regression to differ. Note that you need to allow the residual variances to correlate, since they probably do if they are of the same outcome at different times.

Next, the constrained model, which holds the coefficients equal between the two models:

model2 <- '
    Yi ~ a*X1i + b*X2i
    Yj ~ a*X1j + b*X2j
    Yi~~Yj'
fit2 <- sem(model2, data = data)
summary(fit2)

Finally, compare the fit of the two models using a likelihood ratio test.

lavTestLRT(fit2, fit1)

If significant, there is evidence that at least one of the parameters differs between the models. Otherwise, there is not enough evidence that the parameters differ; unconstraining the coefficients does not add much to the fit of the model.

Answer 2

Having seen the OP's comment on a response stating a reviewer's preference for regression with dummies, I demonstrate that below. Essentially, I repeat what others users have said about using dummies and interactions to control predictor entry in the model. However, we can use generalized least squares or any multilevel software that will allow a model for the variances to account for potential heteroskedastic nature of the responses from different waves. This approach is equivalent to using SEM software.

library(lavaan) # For comparison to regression only
library(nlme) # For GLS

# Data prep
dat <- clubSandwich::MortalityRates[, c(
  "year", "state", "mrate", "cause", "totpercap", "spiritpercap")]
dat <- dat[dat$cause == "Motor Vehicle" & dat$year %in% c(1970, 1980), ]
dat <- dat[, c("year", "state", "mrate", "totpercap", "spiritpercap")]
head(dat)

    year state    mrate totpercap spiritpercap
2   1970     1 62.18053      1.38         0.70
42  1980     1 42.21172      1.86         0.74
110 1970     2 61.80470      3.62         2.07
150 1980     2 30.05453      3.76         1.70
218 1970     4 84.74664      2.73         1.05
258 1980     4 68.25127      3.03         0.98

We now have data from each state in 1970 and 1980. mrate will serve as the response variable. totpercap and spiritpercap as predictors.

dat$year.f <- factor(dat$year) # We make the year a dummy variable

Now we can run our regression models. First, a model where all predictors make year specific predictions. We use a general correlation structure (corSymm) to permit the response to correlate since two responses from the same state by different years are not independent. We permit the response to have different variances by year using the varIdent argument:

(m0 <- gls(
  mrate ~ 0 + year.f + year.f:totpercap + year.f:spiritpercap, dat,
  corSymm(form = ~ 1 | state), varIdent(form = ~ 1 | year.f),
  method = "ML"))

Generalized least squares fit by maximum likelihood
  Model: mrate ~ 0 + year.f + year.f:totpercap + year.f:spiritpercap 
  Data: dat 
  Log-likelihood: -435.8944

Coefficients:
             year.f1970              year.f1980    year.f1970:totpercap    year.f1980:totpercap 
             67.5682105              47.8940439              -0.9308490              12.9756910 
year.f1970:spiritpercap year.f1980:spiritpercap 
              0.9838328             -21.1029073 

Correlation Structure: General
 Formula: ~1 | state 
 Parameter estimate(s):
 Correlation: 
  1    
2 0.822
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | year.f 
 Parameter estimates:
     1970      1980 
1.0000000 0.8592329 
Degrees of freedom: 102 total; 96 residual
Residual standard error: 24.84121 

The predictors enter the model as interactions with the year dummies. The residual standard error is the value for 1970, we multiply the reported value by .859 to obtain the residual standard error for 1980.

Next, a model constraining the totpercap coefficient to be the same across waves, dump the interaction:

coef(mtot <- gls(
  mrate ~ 0 + year.f + totpercap + year.f:spiritpercap, dat,
  corSymm(form = ~ 1 | state), varIdent(form = ~ 1 | year.f),
  method = "ML"))

             year.f1970              year.f1980               totpercap year.f1970:spiritpercap 
               56.38870                49.51440                10.64116               -14.55980 
year.f1980:spiritpercap 
              -16.62852 

Similarly for spiritpercap:

mspi <- gls(
  mrate ~ 0 + year.f + year.f:totpercap + spiritpercap, dat,
  corSymm(form = ~ 1 | state), varIdent(form = ~ 1 | year.f),
  method = "ML")

Finally, a model where all coefficients are constrained to be the same across years:

coef(mtotspi <- gls(
  mrate ~ 0 + year.f + totpercap + spiritpercap, dat,
  corSymm(form = ~ 1 | state), varIdent(form = ~ 1 | year.f),
  method = "ML"))

year.f1970   year.f1980    totpercap spiritpercap
59.240021    51.043122     7.910692   -11.070475 

We can then ask, should we constrain totpercap using a model comparison approach:

anova(m0, mtot)

     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
m0       1  9 889.7887 913.4135 -435.8944                        
mtot     2  8 890.4969 911.4967 -437.2484 1 vs 2 2.708178  0.0998

This permits comparison using IC and likelihood ratio tests. For what it's worth, analysis like this is preferably done with larger sample sizes.

We can do the same for spiritpercap using: anova(m0, mspi). We can ask if it is okay to constrain both: anova(m0, mtotspi).

This is how a regression-based approach would proceed. It is exactly equivalent to the SEM-based approach. I will run the equivalent SEM models and simply report the log likelihoods for comparison. First, transform the data from long to wide format for SEM:

dat.for.lavaan <- reshape(
  dat, idvar = "state", timevar = "year", direction = "wide")

And run the freest SEM model, meanstructure is required to estimate year intercepts as in regression:

logLik(
  sem("
mrate.1970 ~ totpercap.1970 + spiritpercap.1970
mrate.1980 ~ totpercap.1980 + spiritpercap.1980
mrate.1970 ~~ mrate.1980", dat.for.lavaan, meanstructure = TRUE)
)
'log Lik.' -435.8944 (df=9)

logLik(m0) # equivalent multilevel model
'log Lik.' -435.8944 (df=9)

And all constrained:

logLik(
  sem("
mrate.1970 ~ a * totpercap.1970 + b * spiritpercap.1970
mrate.1980 ~ a * totpercap.1980 + b * spiritpercap.1980
mrate.1970 ~~ mrate.1980", dat.for.lavaan, meanstructure = TRUE)
)
'log Lik.' -437.3729 (df=7)

logLik(mtotspi) # equivalent multilevel model
'log Lik.' -437.3729 (df=7)

So the SEM approach and regression approach are equivalent. You can compare parameter estimates. When they differ, it is usually a question of needing to standardize one. For example, the residual covariance between the 1970 outcome and 1980 in SEM if standardized will be the same as the GLS estimate of the correlation between the responses.

Answer 3

You may run a single regression with indicator variable (say $Z$ being 1 for wave 1 and 0 for wave 2) and its interactions with $X$'s.

Your model now becomes

$$Y=\beta_1X_1+\beta_2X_2+\beta_3Z+\gamma_1X_1Z+\gamma_2X_2Z$$

Significance of $\gamma_i$ now means that effect of $X_i$ on $Y$ is significantly different for wave 1 and wave 2.

Answer 4

The usual interpretation would be - one unit increase in variable X1 leads to Beta1 increase in Y. You could extend this to understand what the difference in impact of X1 is on Y across the two waves, assuming the explanatory and target variables remain the same. However I do not understand how you could possibly calculate the significance of this difference between them.