Phylogenetic Generalized Least Squares (PGLS)
Renata Brandt e Fabio Barros
July 20, 2015
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:
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): 
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-16plot(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: nlmepgls.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 residualplot(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 residualCom 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.4884E 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.9015pgls.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.601That’s all folks!