What is the relationship between the ratio $y_{i}/x_{i}$ and the residuals of a linear regression $y_{i}$~$x_{i}$? My background in statistical mathematics is poor. I don't know why $y_{i}/x_{i}$ and the residuals of the linear regression $y_{i}$~$x_{i}$ are always linearly related (whether X and Y are related or not).  Can someone explain what is going on there (mathematically and intuitively)?  
This question is motivated by wanting to use a ratio in a regression (e.g., $y_{i}$~ $x_{i}/z_{i}$), and not quite understanding how this is different from $y_{i}$~ (residuals of $x_{i}$~$z_{i}$), and knowing that one shouldn't use residuals in a regression, rather just regress $y_{i}$ ~ $x_{i}$ + $z_{i}$).
Example in R, with no relationship:
x= rnorm(20,100,10)
y= rnorm(20,100,10)
plot(resid(lm(y~x)), y/x)

or with relationship:
x= rnorm(20,100,10)
y= x+rnorm(20,0,10)
plot(resid(lm(y~x)), y/x)

 A: Keep in mind that in a linear model, you are assuming that $Y$ is simply a linear combination of the predictors and an error term, $Y = \beta X + \epsilon$. If the linear model produced an unbiased estimate of the $\beta$s, residuals from lm is just the $\epsilon$ vector. Also, $y/x = (\beta x + \epsilon)/x = \beta  + \epsilon/x$.
Thus, when you run
plot(resid(lm(y~x)), y/x)

you are in fact plotting $\epsilon$ against $\beta + \epsilon / x$. If your $\epsilon$ is not much smaller than $\sigma(X)$, you should be able to clearly see a linear relationship. Compare:
par(mfrow=c(1,2))
x = rnorm(20,100,10)
e_small = rnorm(20,0,10)
e_large = rnorm(20,0,10)
y_small = x+e_small
y_large = x+e_large
plot(resid(lm(y_small~x)), y_small/x)
plot(resid(lm(y_large~x)), y_large/x)


As for the part about regressing y ~ residuals(x~z): it is perfectly legitimate, and equivalent to y ~ x + z, if the predictors are orthogonal. Regressing each predictor successively, like y ~ residuals(x[n-1] ~ residuals(x[n-2] ~ ... residuals(x[2] ~ x[1]) ... )), is also known as a Gram-Schmidt process. The main problem with it is that once the algorithm reaches a predictor that is collinear with the already included ones, the residuals become small and subject to rounding errors, thus producing numerically unstable $\beta$ estimates.
A: OK, lets see:
x  <-  rnorm(20,100,10)
y  <-  rnorm(20,100,10)

mod  <-  lm(y ~ x)
plot(resid(mod), y/x)  # Plot shown below
cor(resid(mod), y/x)   # about 0.67


so yes, there is some correlation, albeit long from a perfect one. Why is this?  If $y/x$ is large, then $y$ is larger than $x$, yes? But in that case we should expect the residual from the regression to be positive! If you look at the plot above, the relationship you can see is far from linear, it is as if all the points are above the diagonal, and this is also as expected. To see this, since $x$ and $y$ are simulated from the same distribution, independently, we expect a slope of 0 and an intercept close to the mean. Print the model object mod and you will see this is close!  If we use this expected values and not the estimated ones, we get the residual is basically $y-100$ which you plot against $y/x$. Got to be a relationship!
