Existe uma literatura extensa nos convencendo de que é necessário considerar a história evolutiva quando estamos lidando com dados específicos. Para simplificar não vamos tratar disso no nosso curso, já que ele é sobre implementação de análises. Mas quem estiver interessado pode checar o apêndice A de Lavin et al. 2008 para uma breve revisão do desenvolvimento dos métodos comparativos filogenéticos. Praticamente tudo foi desenvolvido depois da proposta dos contrastes independentes do Felsenstein (1985)

Para ilustrar o principal motivo de usar estatística filogenética quando nossos dados são de espécies, Ted Garland produziu o seguinte slide:

alt text

Bem, vamos ao que interessa. Quando você for começar sua análise, lembre-se de especificar o diretório de trabalho. Você vai abrir os dados usando read.table() ou read.csv(), abrir a arvore usando read.tree() ou read.nexus(). No caso agora para facilitar, vou simular dados contínuos que são dependentes e respeitam a estrutura hierárquica da árvore.

require(phytools)
## Loading required package: phytools
## Loading required package: ape
## Loading required package: maps
## simulando uma arvore de coalescencia com 100 tips
set.seed(1)
tree<-rcoal(n=100)
plotTree(tree,ftype="off")

## simulando um dado de desempenho (sprint speed) dependente do tamanho dos membros
hindlimb<-fastBM(tree, a=100)
forelimb<-fastBM(tree, a=70)
sprint<-0.99*hindlimb+ 0.99*forelimb+ 0.01*fastBM(tree, a=10)
habitat <- rep(c('florestas', 'desertos', 'pântanos', 'campos'), 25)

names(hindlimb) <- tree$tip.label
names(forelimb) <- tree$tip.label
names(sprint) <- tree$tip.label

data <- data.frame(sprint, hindlimb, forelimb, habitat)

# Vamos checar a relacao
plot(sprint~forelimb)

plot(sprint~hindlimb)

Uma das primeiras soluções para estudar a relação evolutiva entre dados contínuos foi a dos contrastes independentes do Felsenstein (1985): alt text

Vamos aplicá-lo aos nossos dados:

require(ape)
# Calculate PICs
FLPic <- pic(forelimb, tree)
sprintPic <- pic(sprint, tree)

# A regressao dos contrastes é feita sem intercepto,
# pois a direcao do calculo é arbitrária.
picFit <- lm(sprintPic ~ FLPic - 1)
summary(picFit)
## 
## Call:
## lm(formula = sprintPic ~ FLPic - 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.55808 -0.61373 -0.03062  0.65485  2.44227 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## FLPic  0.99225    0.09524   10.42   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9736 on 98 degrees of freedom
## Multiple R-squared:  0.5255, Adjusted R-squared:  0.5207 
## F-statistic: 108.5 on 1 and 98 DF,  p-value: < 2.2e-16
plot(sprintPic ~ FLPic)
abline(picFit)

Podemos ver que a velocidade de corrida das espécies que estamos estudando é dependente do tamanho do membro anterior. Essa análise também pode ser feita usando PGLS (phylogenetic generalized least squares), com menos passos:

require(nlme)
## Loading required package: nlme
pgls.fit1 <- gls(sprint~forelimb, correlation = corBrownian(1, tree), 
                method = 'ML', data = data)
summary(pgls.fit1)
## Generalized least squares fit by maximum likelihood
##   Model: sprint ~ forelimb 
##   Data: data 
##         AIC       BIC   logLik
##   -53.98179 -46.16628 29.99089
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 99.35596  6.630838 14.98392       0
## forelimb     0.99225  0.095240 10.41834       0
## 
##  Correlation: 
##          (Intr)
## forelimb -0.998
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.2040302 -0.4470672 -0.1799686  0.1310079  0.7452023 
## 
## Residual standard error: 0.7845377 
## Degrees of freedom: 100 total; 98 residual
plot(sprint~forelimb)
abline(pgls.fit1)

O PGLS tem ainda a vantagem de ser mais flexível, permitindo que se modele relações mais complexas. Vejam que minha sprint speed simulada é uma função também do membro posterior dos bichos.

pgls.fit2 <- gls(sprint~hindlimb + forelimb, correlation = corBrownian(1, tree), 
                method = 'ML', data = data)
summary(pgls.fit2)
## Generalized least squares fit by maximum likelihood
##   Model: sprint ~ hindlimb + forelimb 
##   Data: data 
##         AIC       BIC   logLik
##   -954.1334 -943.7127 481.0667
## 
## Correlation Structure: corBrownian
##  Formula: ~1 
##  Parameter estimate(s):
## numeric(0)
## 
## Coefficients:
##                  Value  Std.Error  t-value p-value
## (Intercept) -0.1765387 0.13305080  -1.3269  0.1877
## hindlimb     0.9917361 0.00110672 896.1009  0.0000
## forelimb     0.9913744 0.00105209 942.2938  0.0000
## 
##  Correlation: 
##          (Intr) hndlmb
## hindlimb -0.835       
## forelimb -0.549 -0.001
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.6693799 -1.0943025 -0.7312644 -0.1100457  1.4065997 
## 
## Residual standard error: 0.008622169 
## Degrees of freedom: 100 total; 97 residual

Com o PGLS é possível ainda incluir uma variável discreta.

pgls.fit3 <- gls(sprint~habitat, correlation = corBrownian(1, tree), 
                method = 'ML', data = data)
anova(pgls.fit3)
## Denom. DF: 96 
##             numDF  F-value p-value
## (Intercept)     1 82692.90  <.0001
## habitat         3     0.82  0.4884

E ajustar modelos muito mais complexos

pgls.fit4 <- gls(sprint~(hindlimb + forelimb)*habitat, 
                 correlation = corBrownian(1, tree), method = 'ML', data = data)
anova(pgls.fit4)
## Denom. DF: 88 
##                  numDF    F-value p-value
## (Intercept)          1 1242122425  <.0001
## hindlimb             1     675503  <.0001
## forelimb             1     745489  <.0001
## habitat              3          1  0.6405
## hindlimb:habitat     3          1  0.6311
## forelimb:habitat     3          0  0.9015
pgls.fit5 <- gls(sprint~hindlimb + forelimb + habitat, 
                 correlation = corBrownian(1, tree), method = 'ML', data = data)
anova(pgls.fit5)
## Denom. DF: 94 
##             numDF    F-value p-value
## (Intercept)     1 1376226949  <.0001
## hindlimb        1     748433  <.0001
## forelimb        1     825974  <.0001
## habitat         3          1   0.601

That’s all folks!