1
$\begingroup$

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
tree = read.tree( "small/tree_small.mod" )

## 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?

$\endgroup$

migrated from stackoverflow.com Oct 5 '14 at 3:52

This question came from our site for professional and enthusiast programmers.

2
$\begingroup$

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
tree <- read.tree(text="((t2:1.46,(t6:0.73,t1:0.73):0.733):1.18,(t3:0.119,(t5:0.0554,t4:0.0554):0.0635):2.53);")
group <- gl(2,3, labels=c('Clade1','Clade2'))

# 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)
tiplabels(pch=21, col='black', bg=group, adj=0.8)
plot(b ~ A, data=dat, pch=21, col='black', bg=group)
abline(x, col='grey', lwd=3)
abline(mod, col='white', lty=2, lwd=1.5)
abline(x.phy, col='cornflowerblue', lwd=3)
abline(pmod, col='white', lty=2, lwd=1.5)
$\endgroup$
  • $\begingroup$ Thanks, concerning the covariance matrix,I replaced it by a correlation matrix using vcv(tree,cor=TRUE), but the result of the Cholesky result is the same. What do you mean by "you're only correcting one of the explanatory variables" ? Both input variables (data.gene and data.pheno) have been multiplied by the inverse of the Cholesky matrix. Besides that, the discrepancy between the GLS and the Cholesky approaches diminishes when using chol2inv rather than solve(). $\endgroup$ – Xavier Prudent Sep 24 '14 at 14:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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