## Phylogeny (roughly) and data from 
## Evolution of biometric and life-history traits in lizards (Gallotia)
## from the Canary Islands by M. Molina-Borja nd A. Rodriguez-Doninguez
## J. Zool. Syst Evol. Research 42 (2004): 44-53


##load package ape FIRST

## Then read in tree

lizard_tree <- read.tree()
(Gst:1,
(((((Ggg:0.08, 
Gge:0.08):0.09, 
Ggp:0.17):0.16,
(Gcc:0.15, 
Gcg:0.15):0.18):0.25,
(Gsi:0.33, 
Gin:0.33):0.25):0.09, 
(Gaa:0.17, 
Gam:0.17):0.50
):0.33);

## you can look at the tree

plot(lizard_tree)

## here is some data, in the same order as the tree

SVL <- c(184.93, 107.36, 115.69, 103.29, 73.33, 96.34, 161.44, 146.8, 85.92, 62.76)
PileusWidth<- c(20.55, 13.17, 13.64, 12.11, 8.84, 11.6, 16.02, 16.7, 10.02, 7.41)

## Here is the OU model on one trait

compar.ou(SVL, lizard_tree, alpha=1)


## Here is the species relationship matrix

G <- vcv(lizard_tree)
G

## Here is the bivariate model in package ape

compar.gee(SVL~PileusWidth, phy=lizard_tree)


## Here is how to do Felsenstein Independent Contrasts, essentially
cG <- chol(G)
invG <- chol2inv(G)

ones = array(1, dim=c(length(SVL)))

SVLmean = ones%*%invG%*%SVL/(ones%*%invG%*%ones)
SVLmean

SVLvar = (SVL-SVLmean)%*%invG%*%(SVL-SVLmean)/length(SVL)
SVLvar

Pmean = ones%*%invG%*%PileusWidth/(ones%*%invG%*%ones)
Pmean

Pvar = (PileusWidth-Pmean)%*%invG%*%(PileusWidth-Pmean)/length(PileusWidth)
Pvar

corr = (SVL-SVLmean)%*%invG%*%(PileusWidth-Pmean)/(length(SVL)*sqrt(SVLvar*Pvar))
corr

## Here is a connectivity matrix for compar.cheverud

W <- 1/(array(1, dim=c(length(SVL), length(SVL))) - G)
for(i in 1:length(SVL)){
W[i,i]=0
}
W