Manually transforming two dependent variables according to their correlation matrix in R? I have two continuous variables, X and Y, that are correlated - they are not independent. To correct for non-independence, I have a known correlation structure, a matrix S.
If one calls gls(Y ~ X, correlation = S), what I think happens is that, internally, gls() transforms X and Y in some way so that the regression ends up being S^(-1)*Y = S^(-1) * X.
How is this transformation actually performed? From the literature I've consulted, I've seen everything from:
X.transformed <- solve(chol(S)) %*% X 
#The inverse of the Choleski decomposition of S times the vertical vector X, 
#which in my case does nothing to the data

to
X.transformed <- chol(solve(S)) %*% X 
# which has negative values and gives meaningless values of X

Another method I've seen is transforming the dependent variable by 
chol(solve(S)) %*% Y 

and the independent variable by 
chol(solve(S)) %*% cbind(1,X) 

and doing the linear model using the transformed intercept terms in the first column of the X matrix: 
lm(Y ~ X - 1)

On a related note, is there any point to manually transforming the data in order to plot it? Do the transformed values have any meaning, or are they simply there to estimate regression coefficients? (In other words, if X is a variable of body mass figures, X values are not necessarily errant if they're negative since they're still linear?) I suppose it would follow from this that an $R^2$ statistic on transformed variables is also meaningless?
 A: When errors are dependent, for example with multi-level models, the matrix solution is:
$\beta=(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}Y$
Where $\Sigma$ is the covariance matrix.  For more information, see here.  The point of this procedure is not that $\hat\beta$ would otherwise be erroneous (i.e., biased), but rather that otherwise it would be inefficient.  That is, the sampling distribution of the OLS estimate of $\beta$ has greater variance than the GLS estimate.  In addition, the confidence intervals (e.g., p-values) would be incorrect.  (So, this is as much a problem of inference as estimation.)  
As for your question about how one could transform the data for plotting, I've never done that, but I will copy the information from the wikipaedia page:

$\Sigma = BB'$
$Y^* = B^{−1}Y$
$X^* = B^{−1}X$
$\epsilon^* = B^{−1}\epsilon$
$Y^*=X^*\beta+\epsilon^*$

As I say, I've never done that, but I don't see why it would be meaningless.  You could also just plot your individual data-lines, with a bolded regression line (using the mean slope and mean intercept values as parameters, for example) superimposed over them.  
