How to compare two regression slopes for one predictor on two different outcomes? 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 

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


*

*$APD=\beta_{0,1} + \beta_{1,1}\cdot NOL$

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


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


*

*Create a new dependent variable ($Y_{new}$) by just appending APD to HL

*Create a new independent variable by appending NOL to itself ($X_{new}$) (i.e. duplicating NOL)

*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).

*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.
