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)