The question asks about "identifying underlying [linear] relationships" among variables.
The quick and easy way to detect relationships is to regress any other variable (use a constant, even) against those variables using your favorite software: any good regression procedure will detect and diagnose collinearity. (You will not even bother to look at the regression results: we're just relying on a useful side-effect of setting up and analyzing the regression matrix.)
Assuming collinearity is detected, though, what next? Principal Components Analysis (PCA) is exactly what is needed: its smallest components correspond to near-linear relations. These relations can be read directly off the "loadings," which are linear combinations of the original variables. Small loadings (that is, those associated with small eigenvalues) correspond to near-collinearities. An eigenvalue of $0$ would correspond to a perfect linear relation. Slightly larger eigenvalues that are still much smaller than the largest would correspond to approximate linear relations.
(There is an art and quite a lot of literature associated with identifying what a "small" loading is. For modeling a dependent variable, I would suggest including it within the independent variables in the PCA in order to identify the components--regardless of their sizes--in which the dependent variable plays an important role. From this point of view, "small" means much smaller than any such component.)
Let's look at some examples. (These use R for the calculations and plotting.) Begin with a function to perform PCA, look for small components, plot them, and return the linear relations among them.
pca <- function(x, threshold, ...) {
fit <- princomp(x)
#
# Compute the relations among "small" components.
#
if(missing(threshold)) threshold <- max(fit$sdev) / ncol(x)
i <- which(fit$sdev < threshold)
relations <- fit$loadings[, i, drop=FALSE]
relations <- round(t(t(relations) / apply(relations, 2, max)), digits=2)
#
# Plot the loadings, highlighting those for the small components.
#
matplot(x, pch=1, cex=.8, col="Gray", xlab="Observation", ylab="Value", ...)
suppressWarnings(matplot(x %*% relations, pch=19, col="#e0404080", add=TRUE))
return(t(relations))
}
Let's apply this to some random data. These are built on four variables (the $B,C,D,$ and $E$ of the question). Here is a little function to compute $A$ as a given linear combination of the others. It then adds i.i.d. Normally-distributed values to all five variables (to see how well the procedure performs when multicollinearity is only approximate and not exact).
process <- function(z, beta, sd, ...) {
x <- z %*% beta; colnames(x) <- "A"
pca(cbind(x, z + rnorm(length(x), sd=sd)), ...)
}
We're all set to go: it remains only to generate $B, \ldots, E$ and apply these procedures. I use the two scenarios described in the question: $A=B+C+D+E$ (plus some error in each) and $A=B+(C+D)/2+E$ (plus some error in each). First, however, note that PCA is almost always applied to centered data, so these simulated data are centered (but not otherwise rescaled) using sweep.
n.obs <- 80 # Number of cases
n.vars <- 4 # Number of independent variables
set.seed(17)
z <- matrix(rnorm(n.obs*(n.vars)), ncol=n.vars)
z.mean <- apply(z, 2, mean)
z <- sweep(z, 2, z.mean)
colnames(z) <- c("B","C","D","E") # Optional; modify to match `n.vars` in length
Here we go with two scenarios and three levels of error applied to each. The original variables $B, \ldots, E$ are retained throughout without change: only $A$ and the error terms vary.
The output associated with the upper left panel was
A B C D E
Comp.5 1 -1 -1 -1 -1
This says that the row of red dots--which is constantly at $0$, demonstrating a perfect multicollinearity--consists of the combination $0 \approx A -B-C-D-E$: exactly what was specified.
The output for the upper middle panel was
A B C D E
Comp.5 1 -0.95 -1.03 -0.98 -1.02
The coefficients are still close to what we expected, but they are not quite the same due to the error introduced. It thickened the four-dimensional hyperplane within the five-dimensional space implied by $(A,B,C,D,E)$ and that tilted the estimated direction just a little. With more error, the thickening becomes comparable to the original spread of the points, making the hyperplane almost impossible to estimate. Now (in the upper right panel) the coefficients are
A B C D E
Comp.5 1 -1.33 -0.77 -0.74 -1.07
They have changed quite a bit but still reflect the basic underlying relationship $A' = B' + C' + D' + E'$ where the primes denote the values with the (unknown) error removed.
The bottom row is interpreted the same way and its output similarly reflects the coefficients $1, 1/2, 1/2, 1$.
In practice, it is often not the case that one variable is singled out as an obvious combination of the others: all coefficients may be of comparable sizes and of varying signs. Moreover, when there is more than one dimension of relations, there is no unique way to specify them: further analysis (such as row reduction) is needed to identify a useful basis for those relations. That's how the world works: all you can say is that these particular combinations that are output by PCA correspond to almost no variation in the data. To cope with this, some people use the largest ("principal") components directly as the independent variables in the regression or the subsequent analysis, whatever form it might take. If you do this, do not forget first to remove the dependent variable from the set of variables and redo the PCA!
Here is the code to reproduce this figure:
par(mfrow=c(2,3))
beta <- c(1,1,1,1) # Also can be a matrix with `n.obs` rows: try it!
process(z, beta, sd=0, main="A=B+C+D+E; No error")
process(z, beta, sd=1/10, main="A=B+C+D+E; Small error")
process(z, beta, sd=1/3, threshold=2/3, main="A=B+C+D+E; Large error")
beta <- c(1,1/2,1/2,1)
process(z, beta, sd=0, main="A=B+(C+D)/2+E; No error")
process(z, beta, sd=1/10, main="A=B+(C+D)/2+E; Small error")
process(z, beta, sd=1/3, threshold=2/3, main="A=B+(C+D)/2+E; Large error")
(I had to fiddle with the threshold in the large-error cases in order to display just a single component: that's the reason for supplying this value as a parameter to process.)
User ttnphns has kindly directed our attention to a closely related thread. One of its answers (by J.M.) suggests the approach described here.
