Using pairwise differences as variable in regression I have a dataset consisting of 72 data points, with 1 dependent variable $y$ and 37 predictors $X_j$.
The typical way to perform a linear regression would be to model $y$ as
$$y = XB + \epsilon$$
where an intercept may or may not be included in the matrix of predictors. 
I thought about an alternative model, in which the difference in $y$ for pairs of cases is regressed on the difference in the predictor's values for those same cases, i.e., defining
$$ \{y_i - y_j\} \equiv y^*\\ \{X_i - X_j\} \equiv X^*$$
for all possible $(i,j)$ pairs (or a subset of all possible pairs), and then performing a linear regression on these variables
$$y^* = X^*B^* + \epsilon^*$$
I just wonder whether this is would be valid approach or if there is something inherently wrong in it. I could not find find any useful information on the web. I would be very grateful for any helpful comments or advice. 
Additional information
Simple tests with synthetic data suggest that the two models give identical results, as expected. 
However, there is a situation in which it seems to me that the second model could be preferable: suppose that the data was collected over several days, and that the measurements of $y$ are subject to a systematic bias due to, e.g., the temperature of the room which varies randomly from day to day. In this situation, taking the pairwise differences of all data points collected on the same day would remove the temperature bias, and thus yield a more accurate regression.
Again, a simple test with synthetic data showed that as expected, in the situation described above the pairwise-difference model produces smaller errors in the estimation of the regression weights than does the regular model. 
 A: It's an interesting idea, but it's a different model and there are inherent difficulties in fitting it.
You did well by explicitly including the error term $\epsilon$ in the original model.  This determines the error term in the differenced model; to wit,
$$\eqalign{
y^{*}_{ji} = y_j - y_i &= (x_iB + \epsilon_j) - (x_jB + \epsilon_i) \\
&= (x_i-x_j)B + (\epsilon_j-\epsilon_i) \\
&= x^{*}_{ji}B^{*} + \epsilon^{*}_{ji}.
}$$
Comparing the last two lines shows $B^{*}=B$ (as you had hoped) and $\epsilon^{*}_{ji} = \epsilon_j - \epsilon_i.$
The usual assumptions are that the $\epsilon_i$ are independent, have zero means, and have (at least approximately) the same common variance.  Although the latter two will be the case for the differenced model when they are true originally, the $\epsilon^{*}_{ji}$ are decidedly not independent!  Indeed, assuming the $\epsilon$ have a common (finite) variance $\sigma^2,$ their independence implies
$$\operatorname{Var}(\epsilon^{*}_{ji}) = \operatorname{Var}(\epsilon_j) + \operatorname{Var}(\epsilon_i) = \sigma^2 + \sigma^2 = 2\sigma^2,$$
showing that the errors in the differenced model have twice the variance of the original errors and for $k\ne i,$
$$\operatorname{Cov}(\epsilon^{*}_{ji}, \epsilon^{*}_{jk}) = \operatorname{Cov}(\epsilon_j-\epsilon_i, \epsilon_j-\epsilon_k) = \sigma^2 - 0 - 0 + 0 = \sigma^2;$$
$$\operatorname{Cov}(\epsilon^{*}_{ij}, \epsilon^{*}_{jk}) = \operatorname{Cov}(\epsilon_i-\epsilon_j, \epsilon_j-\epsilon_k) = 0 - \sigma^2 - 0 + 0 = -\sigma^2,$$
showing that when one $\epsilon^{*}$ shares a subscript with another $\epsilon^{*},$ their correlation is $\pm1/2.$  Although you could handle this circumstance with generalized least squares, the matrices involved will grow large rather quickly: for instance, with your tiny dataset of just 72 observation, the weights matrix is square with dimension $72(71)/2 = 2556,$ containing almost seven million entries (albeit somewhat sparsely).
What are the consequences of ignoring this lack of independence?  I'm afraid your synthetic data may have misled you in several ways into drawing false conclusions.


*

*When there is a constant term, all its differences are zero and therefore drop out of the model.  Thus, the differenced model should include no constant term (intercept) at all, and if it does, its estimate should be insignificantly different from zero--but usually won't be.  Regardless, by taking differences of the responses, you can no longer estimate the constant term: because it always multiplies $1-1=0,$ it is not identifiable.

*The strong correlations among the errors imply the effective amount of data--"degrees of freedom"--to use for computing standard errors and p-values is much smaller than the number of pairs.  Indeed, because subtracting the paired data introduces no new information at all, a correct calculation of the degrees of freedom is the same as before: $72-37=35$ (without an intercept), even though it appears you have $2556$ observations!

*As a result, if you do not use generalized least squares, all p-values will be (far) too small and many estimates may appear to be "significant" when they are the result of noise alone.

*If you include an intercept in the difference model, the parameter estimates will differ.  If you leave it out, you should obtain correct estimates: after all, you are estimating exactly the same parameters $B$ using exactly the same data.
Finally, although it is instructive to compare these approaches, one wonders what could possibly be the advantage of turning a problem involving computing with a $72\times 37$ matrix into one involving a $2556\times 37$ matrix: it takes (much) more computation to obtain the same estimates (ignoring the loss of ability to estimate the intercept) and far more work to obtain correct p-values, confidence intervals, hypothesis tests, and so on.

This R code produces sample data, fits both models, and compares their results.  It can be instructive to review the model summaries that it outputs.
n <- 72                   # Number of observations
k <- 37                   # Number of variables, apart from the constant
beta <- (-1)^seq_len(k+1) # The model coefficients, beginning with the intercept
sigma <- 0.5              # Error s.d.
#
# Create data.
#
X <- as.data.frame(matrix(rnorm(n*k), nrow=n, dimnames=list(NULL, paste0("x.", 1:k))))
X$y <- cbind(1, as.matrix(X)) %*% beta + rnorm(n, 0, sigma) # The response
#
# Base model.
#
fit <- lm(y ~ ., X)
summary(fit)
#
# Paired model.
#
xx <- as.matrix(X)
XX <- as.data.frame(t(apply(combn(1:n, 2), 2, function(ij) xx[ij[1], ] - xx[ij[2], ])))
names(XX) <- paste0("d", names(XX))
fit.d <- lm(dy ~ . - 1, XX)         # No intercept!
summary(fit.d)
#
# Check whether coefficients agree.
#
all.equal(coefficients(fit)[-1], coefficients(fit.d), check.attributes=FALSE)
#
# Plot the data.
#
# pairs(X)
# pairs(XX)

A: This second approach is often used in experimental situation where some data is collected before and after treatment for same subjects. If you would have all observation from two time periods then observations would not be independent and you should specify covariance structure for generalized regression model. It is easiear to use difference approach in this case. 
So I would say your approach is OK.
A: A quick simulation shows that doing this just gives you anti-conservative (rejects the null too often) normal regression
Example:
set.seed(123)
x <- seq(-3,3, by=0.1)
y <- rnorm(length(x), x, 1)

dx <- outer(x,x,`-`)
dy <- outer(y,y,`-`)

summary(lm(c(dy) ~ c(dx)))
summary(lm(y ~ x))

gives a test stat of 15 in the "normal" case and 122 in the pairwise difference case.
You might be interested to know that the linear regression coefficient is the weighted average of the pairwise slope. That is:
$$\hat{\beta}_{OLS} = \sum_{i=1}^n \sum_{j \ne i} w_{ij} \frac{Y_i - Y_j}{X_i - X_j} / \sum_{i=1}^n \sum_{j \ne i} w_{ij} $$
where $w_{ij} = (X_i - X_j)^2$.
