10
$\begingroup$

I need to compare two regression slopes where:

$
y_1 ~ a + b_1x
y_2 ~ a + b_2x
$

How can I compare b1 and b2?

Or in the language of my specific example in rodents, I want to compare

antero-posterior diameter ~  a + b1 * humeral length   
de naso-occipital length  ~  a + b2 * humeral length 
$\endgroup$
  • 2
    $\begingroup$ Calculate a regression model with both variables PLUS the interaction of the two variables (humeral length $\times$ antero-posterior diameter). The interaction tests the assumption of parallelism of the slopes of the two variables. If the interaction term is significant, the slopes are different. $\endgroup$ – COOLSerdash May 17 '13 at 18:42
  • $\begingroup$ Thanks!! But humeral length and antero-posterior diameter of the humerus are DVs and naso-occipital length is the IV. Can i run the analysis as you suggest? $\endgroup$ – Dra. Alejandra Echeverria May 17 '13 at 18:55
  • 1
    $\begingroup$ @Dra.AlejandraEcheverria Do you mean that you have one linear regression model with two independent variables and that you want to test the equality of the two coefficients on the independent variables, or, you have two simple linear regression models and you want to compare the coefficients across the two models? $\endgroup$ – Graeme Walsh May 17 '13 at 19:21
  • 1
    $\begingroup$ Dear @Graeme Walsh, I've two simple linear regression models and I want to compare the coefficients across the two models. $\endgroup$ – Dra. Alejandra Echeverria May 20 '13 at 17:28
11
$\begingroup$

Okay, let's look at your situation. You have basically two regressions (APD = antero-posterior diameter, NOL = naso-occipital length, HL = humeral length):

  1. $APD=\beta_{0,1} + \beta_{1,1}\cdot NOL$
  2. $HL=\beta_{0,2} + \beta_{1,2}\cdot NOL$

To test the hypothesis $\beta_{1,1}=\beta_{1,2}$, you can do the following:

  1. Create a new dependent variable ($Y_{new}$) by just appending APD to HL
  2. Create a new independent variable by appending NOL to itself ($X_{new}$) (i.e. duplicating NOL)
  3. Create a dummy variable ($D$) that is 1 if the data came from the second dataset (with HL) and 0 if the data came from the first dataset (APD).
  4. Calculate the regression with $Y_{new}$ as dependent variable, and the main effects and the interaction between $X_{new}$ and the dummy variable $D$ as explanatory variables. EDIT @Jake Westfall pointed out that the residual standard error could be different for the two regressions for each DV. Jake provided the answer which is to fit a generalized least squares model (GLS) which allows the residual standard error to differ between the two regressions.

Let's look at an example with made-up data (in R):

# Create artificial data

library(nlme) # needed for the generalized least squares

set.seed(1500)

NOL <- rnorm(10000,100,12)
APD <- 10 + 15*NOL+ rnorm(10000,0,2)
HL <- - 2  - 5*NOL+ rnorm(10000,0,3) 

mod1 <- lm(APD~NOL)
mod1

Coefficients:
(Intercept)          NOL
      10.11        15.00

mod2 <- lm(HL~NOL)
mod2

Coefficients:
(Intercept)          NOL
      -1.96        -5.00

# Combine the dependent variables and duplicate the independent variable

y.new <- c(APD, HL)
x.new <- c(NOL, NOL)

# Create a dummy variable that is 0 if the data are from the first data set (APD) and 1 if they are from the second dataset (HL)

dummy.var <- c(rep(0, length(APD)), rep(1, length(HL)))

# Generalized least squares model allowing for differend residual SDs for each regression (strata of dummy.var)

gls.mod3 <- gls(y.new~x.new*dummy.var, weights=varIdent(form=~1|dummy.var))

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | dummy.var 
 Parameter estimates:
       0        1 
1.000000 1.481274 

Coefficients:
                    Value  Std.Error   t-value p-value
(Intercept)      10.10886 0.17049120    59.293       0
x.new            14.99877 0.00169164  8866.430       0
dummy.var       -12.06858 0.30470618   -39.607       0
x.new:dummy.var -19.99917 0.00302333 -6614.939       0

Note: The intercept and the slope for $X_{new}$ are exactly the same as in the first regression (mod1). The coefficient of dummy.var denotes the difference between the intercept of the two regressions. Further: the residual standard deviation of the second regression was estimated larger than the SD of the first (about 1.5 times larger). This is exactly what we have specified in the generation of the data (2 vs. 3). We're nearly there: The coefficient of the interaction term (x.new:dummy.var) tests the equality of the slopes. Here the slope of the second regression (mod2) is about $\beta_{x.new} - \beta_{x.new\times dummy.var}$ or about $15-20=-5$. The difference of $20$ is exactly what we've specified when we generated the data. If you work in Stata, there is a nice explanation here.

Warning: This does only work if the antero-posterior diameter and naso-occipital length (the two dependent variables) are independent. Otherwise it can get very complicated.

EDIT

These two posts on the site deal with the same question: First and second.

$\endgroup$
  • $\begingroup$ Just to avoid confusion, it looks like you've got NOL and HL mixed up. HL was the predictor, NOL was the second DV (and APD was the first DV, as you pointed out). Though I just noticed that the poster herself changed the status of her variables in a comment... $\endgroup$ – Patrick Coulombe May 18 '13 at 6:51
  • $\begingroup$ @Patrick Coulombe Thanks for pointing it out. It wasn't clear from her comment yesterday. $\endgroup$ – COOLSerdash May 18 '13 at 6:57
  • $\begingroup$ @PatrickCoulombe On a second thought: I think Jeromy Anglim misunderstood Alejandra's comment and exchanged the variables. $\endgroup$ – COOLSerdash May 18 '13 at 7:42
  • 1
    $\begingroup$ This solution seems reasonable, but I am slightly concerned about the fact that in your combined/interactive model, the residual variance is assumed to be equal at both levels of dummy.var, that is, for both of the DVs. Depending on what the DVs are in the original context, it's possible that the residual variances are radically different in the separate regressions of each DV. I wonder if it would be better to use the same basic approach you proposed, but with a gls model where we estimate different residual variances for each DV. Any thoughts about this? $\endgroup$ – Jake Westfall May 18 '13 at 8:02
  • 1
    $\begingroup$ @COOLSerdash Sure, it would look something like this: library(nlme); mod4 <- gls(y.new~x.new*dummy.var, weights=varIdent(form= ~1 | dummy.var)) $\endgroup$ – Jake Westfall May 18 '13 at 8:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.