# R: difference between Generalized Least Square and the Standard Least Squares with Cholesky

According to Wikipedia (source of all truth and knowledge...), http://en.wikipedia.org/wiki/Generalized_least_squares#Properties

a weighted least square regression is equivalent to a standard least square regression, if the variables have been previously "decorrelated" using a Cholesky decomposition.

I made up then a very simple example with the function pgls from the package CAPER to test it, where the correlation arises from a phylogeny tree:

tree.mod:
((A:0.2,(B:0.1,C:0.1):0.1):0.1,((E:0.1,F:0.1):0.1,D:0.2):0.1);


The two approaches are compared here:

library(caper)

## Data
species = c("A","B","C","D","E","F")
gene = c(0.1,0.2,0.3,0.5,0.6,0.7)
pheno = c( 0,0,0,1,1,1)
data=data.frame(species,gene,pheno)

## Phylogeny

## GLS regression
cat("\n     ===> GLS\n")
cdata   = comparative.data( phy = tree, data = data, names.col = "species" )
res = pgls( pheno~gene, cdata )
print(summary(res))

## Cholesky
cat("\n     ===> Cholesky\n")
corr = vcv( tree )
cholesky = chol( corr )
invCho = solve( cholesky )
data.gene =  invCho %*% as.vector( data$gene ) data.pheno = invCho %*% as.vector( data$pheno )
res=lm( data.pheno ~ data.gene )
print(summary(res))


and yield the outputs:

====> GLS
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.13000    0.27261 -0.4769  0.65834
gene         1.63333    0.59489  2.7456  0.05161 .

=====>Cholesky
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.02214    0.28551   0.078    0.942
data.gene    1.29188    0.35006   3.690    0.021 *


as you can see the results are different...

Does anyone have a clue why?

I think one of the issues is that you're only correcting one of the explanatory variables - there is an intercept in the model too and the transformation has to be applied to the entire model design matrix.

Explaining it to myself a while ago, I made the example below showing how the two models are calculated by the functions and through solving the linear equation Ax=b using the original and transformed data to generate the lm and pgls coefficients.

I can't work out how to tweak your example to replicate this. The vcv matrix is a covariance, not a correlation matrix and solve(cholesky) and chol2inv(cholesky) give different answers so it may be some combination of tweaking those bits of the code.

library(caper)

# a phylogeny

# a data frame of variables and a comparative dataset
data <- data.frame(species = tree$tip.label, A = 1:6, b = c(3.068, 2.213, 2.529, 5.681, 6.23, 6.05)) cdata <- comparative.data(tree, data=data, names='species') # create a design and response matrix for the model Amat <- model.matrix(b ~ A, data=data) b <- as.matrix(dat$b)

# OLS
# solving Ax=b for x:
# x = (A^TA)^{-1} A^Tb
x <- solve(t(Amat) %*% Amat) %*% t(Amat)%*%b

# PGLS

# same basic equation Ax=b, but correcting for phylogenetic signal
# get the inverse of the VCV matrix
V <- VCV.array(tree)
iV <- solve(V)

# definition of inverse: VV'=I
zapsmall(V %*% iV)

# apply phylogenetic correction to A and b
Amat.prime <- iV %*% Amat
b.prime <- iV %*% b

# solve Ax=b for x
# x = (A^TA')^{-1} A^Tb'
x.phy <- solve(t(Amat) %*% Amat.prime) %*% t(Amat)%*%b.prime

# CHECK USING LM AND PGLS
mod <- lm(b ~ A, data=dat)
pmod <- pgls(b ~ A, data=cdata)

# Visualise
par(mfrow=c(1,2), mar=c(3,3,1,1), mgp=c(2,1,0))
palette(c('khaki','salmon' ))
plot(tree)