Does the variable order matter in linear regression I'm investigating interplay between two variables ($x_1$ and $x_2$). There is a great deal of linear correlation between these variables with $r>0.9$. From the nature of the problem I cannot say anything about the causation (whether $x_1$ causes $x_2$ or the other way around). I would like to study deviations from the regression line, in order to detect outliers. In order to do this I can either build a linear regression of $x_1$ as a function of $x_2$, or the other way around. Can my choice of variable order influence my results?
 A: It surely can (actually, it even matters with regard to the assumptions on your data - you only make assumptions about the distribution of the outcome given the covariate). In this light, you might look up a term like "inverse prediction variance". Either way, linear regression says nothing about causation! At best, you can say something about causation through careful design.
A: To make the case symmetrical, one may regress the difference between the two variables ($\Delta x$) vs their average value. 
A: Standard regression minimizes the vertical distance between the points and the line, so switching the 2 variables will now minimize the horizontal distance (given the same scatterplot).  Another option (which goes by several names) is to minimize the perpendicular distance, this can be done using principle components.
Here is some R code that shows the differences:
library(MASS)

tmp <- mvrnorm(100, c(0,0), rbind( c(1,.9),c(.9,1)) )

plot(tmp, asp=1)

fit1 <- lm(tmp[,1] ~ tmp[,2])  # horizontal residuals
segments( tmp[,1], tmp[,2], fitted(fit1),tmp[,2], col='blue' )
o <- order(tmp[,2])
lines( fitted(fit1)[o], tmp[o,2], col='blue' )

fit2 <- lm(tmp[,2] ~ tmp[,1])  # vertical residuals
segments( tmp[,1], tmp[,2], tmp[,1], fitted(fit2), col='green' )
o <- order(tmp[,1])
lines( tmp[o,1], fitted(fit2)[o], col='green' )

fit3 <- prcomp(tmp)
b <- -fit3$rotation[1,2]/fit3$rotation[2,2]
a <- fit3$center[2] - b*fit3$center[1]
abline(a,b, col='red')
segments(tmp[,1], tmp[,2], tmp[,1]-fit3$x[,2]*fit3$rotation[1,2], tmp[,2]-fit3$x[,2]*fit3$rotation[2,2], col='red')

legend('bottomright', legend=c('Horizontal','Vertical','Perpendicular'), lty=1, col=c('blue','green','red'))

To look for outliers you can just plot the results of the principle components analysis.
You may also want to look at:

Bland and Altman (1986), Statistical
  Methods for Assessing Agreement
  Between Two Methods of CLinical
  Measurement. Lancet, pp 307-310

A: Your x1 and x2 variables are collinear. In the presence of multicollinearity, your parameter estimates are still unbiased, but their variance is large, i.e., your inference on the significance of the parameter estimates is not valid, and your prediction will have large confidence intervals. 
Interpretation of the parameter estimates is also difficult. In the linear regression framework, the parameter estimate on x1 is the change in Y for a unit change in x1 given every other exogeneous variable in the model is held constant. In your case, x1 and x2 are highly correlated, and you cannot hold x2 constant when x1 is changing.
