5
$\begingroup$

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?

$\endgroup$

1 Answer 1

4
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Thanks, that's very helpful. In the literature, covariance matrix and correlation matrix seems to be used interchangeably - but the R functions computing these give different matrices depending on which matrix you request. Am I only to use the covariance matrix? Also, is there a built-in function to save things like Total Sum of Squares from a gls() call, by any chance? gls() doesn't explicitly calculate $R^2$, but in principle it should be possible given the stuff gls() does internally. I seem to get different values doing manual transformation, so getting the internal info seems best. $\endgroup$ Jan 7, 2012 at 16:11
  • 1
    $\begingroup$ I'm not sure exactly what happens behind the scenes. Statistical software always requests correlation matrices as inputs, probably because they're intelligible. E.g., if you fit a model that nests students within classrooms, you can use an exchangeable matrix, as all students would be equally correlated with each other on average; or if you looked at factory output over a business cycle, you could use an autoregressive matrix, since outputs closer in time will be more strongly correlated. But these may be converted to covariance matrices for the algorithm. $\endgroup$ Jan 7, 2012 at 19:07
  • 1
    $\begingroup$ But the underlying theory is in terms of covariances. I don't do this manually, so it's harder to advise on that. As for getting sums of squares, you should be able to call anova(model) to get that. $\endgroup$ Jan 7, 2012 at 19:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.