Roteiro prático — Regressão filogenética (PGLS)

Aula 5 · Métodos Filogenéticos Comparativos

Este roteiro roda de verdade: os resultados abaixo saíram do código mostrado. Copie os blocos para um script e rode em ordem numa sessão limpa do R. No fim, o desafio pede o mesmo raciocínio com outro par de atributos.

0. Pacotes e dados

library(ape)
library(phytools)
library(nlme)
library(phylolm)
library(rr2)

Os Anolis das Grandes Antilhas (Mahler et al. 2010, no phytools). A árvore vem com os ecomorfos mapeados (classe simmap); para o gls(), basta a árvore comum.

data(anoletree, anole.data)
ecomorfo <- getStates(anoletree, "tips")
arv <- as.phylo(anoletree)
arv

Phylogenetic tree with 82 tips and 81 internal nodes.

Tip labels:
  ahli, allogus, rubribarbus, imias, sagrei, bremeri, ...

Rooted; includes branch length(s).
head(anole.data)
                SVL      HL     HLL     FLL     LAM      TL
ahli        4.03913 2.88266 3.96202 3.34498 2.86620 4.50400
allogus     4.04014 2.86103 3.94018 3.33829 2.80827 4.52189
rubribarbus 4.07847 2.89425 3.96135 3.35641 2.86751 4.56108
imias       4.09969 2.85293 3.98565 3.41402 2.94375 4.65242
sagrei      4.06716 2.83515 3.85786 3.24267 2.91872 4.77603
bremeri     4.11337 2.86044 3.90039 3.30585 2.97009 4.72996

Os atributos do anole.data já estão em logaritmo natural (SVL de 3,5 a 5,1 corresponde a 32–166 mm): não aplique log() de novo.

dados <- data.frame(sp = rownames(anole.data),
                    svl = anole.data$SVL, hll = anole.data$HLL,
                    eco = factor(ecomorfo[rownames(anole.data)]))
rownames(dados) <- dados$sp
head(dados)
                     sp     svl     hll eco
ahli               ahli 4.03913 3.96202  TG
allogus         allogus 4.04014 3.94018  TG
rubribarbus rubribarbus 4.07847 3.96135  TG
imias             imias 4.09969 3.98565  TG
sagrei           sagrei 4.06716 3.85786  TG
bremeri         bremeri 4.11337 3.90039  TG

1. O GLS à mão, na árvore-guia

O estimador GLS usa a matriz C (Aula 2) para pesar as observações: β = (X′C⁻¹X)⁻¹ X′C⁻¹y. Com os atributos A e B da árvore-guia:

guia <- read.tree(text = "(((sp1:0.18,sp2:0.18):1.06,(sp3:0.02,sp4:0.02):1.22):0.20,(sp5:1.37,(sp6:0.23,sp7:0.23):1.14):0.07);")
A <- c(2.51, 2.80, 2.23, 2.32, 1.64, 1.18, 0.58)
B <- c(0.38, 0.89, 0.88, 1.60, 0.33, 1.42, 0.31)
C <- vcv(guia)
X <- cbind(1, A)
Ci <- solve(C)
beta_gls <- solve(t(X) %*% Ci %*% X, t(X) %*% Ci %*% B)
beta_gls
       [,1]
  -1.796734
A  1.398177

A inclinação, 1,398, é a mesma da regressão dos contrastes pela origem (Aula 2). Com C = identidade (a árvore em estrela), o GLS vira o OLS:

solve(t(X) %*% X, t(X) %*% B)
       [,1]
  0.5749712
A 0.1346306
coef(lm(B ~ A))
(Intercept)           A 
  0.5749712   0.1346306

2. PGLS no R

Agora os Anolis: comprimento do membro posterior contra SVL. Três estruturas para os resíduos: nenhuma (OLS), BM e λ de Pagel. O form = ~sp amarra cada linha da tabela à ponta certa da árvore.

m_ols <- gls(hll ~ svl, data = dados, method = "ML")
m_bm  <- gls(hll ~ svl, data = dados, method = "ML",
             correlation = corBrownian(1, arv, form = ~sp))
m_lam <- gls(hll ~ svl, data = dados, method = "ML",
             correlation = corPagel(1, arv, form = ~sp))
summary(m_lam)$tTable
                 Value  Std.Error   t-value      p-value
(Intercept) -0.4991833 0.19408501 -2.571983 1.196221e-02
svl          1.0461523 0.04665884 22.421308 3.059259e-36
tab <- sapply(list(OLS = m_ols, BM = m_bm, lambda = m_lam), \(m)
  c(inclinacao = coef(m)[["svl"]], EP = sqrt(diag(vcov(m)))[["svl"]], AIC = AIC(m)))
round(tab, 3)
               OLS       BM   lambda
inclinacao   1.031    1.049    1.046
EP           0.045    0.047    0.047
AIC        -55.050 -147.709 -146.811
m_lam$modelStruct$corStruct   # o lambda estimado
Correlation structure of class corPagel representing
  lambda 
1.013439

Nos Anolis, as três inclinações são parecidas, mas o AIC prefere com folga os modelos com a árvore: o OLS subestima a dependência entre as espécies. O corPagel às vezes estima λ um pouco acima de 1; isso acontece e não quer dizer erro.

3. Quanto o modelo explica?

O R² do OLS não serve para resíduos correlacionados. O rr2 calcula três versões (Ives 2019):

R2(m_lam)
   R2_lik  R2_resid   R2_pred 
0.9574886 0.9612440 0.9610952

4. Diagnóstico

O pressuposto do PGLS é sobre os resíduos. Os resíduos “normalizados” do gls() já descontam a estrutura filogenética e devem parecer ruído.

par(mfrow = c(1, 2))
plot(fitted(m_lam), resid(m_lam, type = "normalized"),
     xlab = "ajustado", ylab = "resíduo normalizado"); abline(h = 0, lty = 2)
qqnorm(resid(m_lam, type = "normalized")); qqline(resid(m_lam, type = "normalized"))

Quanto sinal há nos resíduos do OLS? Se houver, o OLS está errado e o PGLS é necessário, independentemente do sinal em cada atributo.

res_ols <- setNames(resid(m_ols), dados$sp)
phylosig(arv, res_ols, method = "lambda", test = TRUE)

Phylogenetic signal lambda : 0.999927 
logL(lambda) : 76.7724 
LR(lambda=0) : 92.4945 
P-value (based on LR test) : 6.75106e-22

λ ≈ 1 nos resíduos do OLS: o modelo sem árvore viola o pressuposto de independência.

5. ANOVA e ANCOVA filogenéticas

O membro posterior relativo (o resíduo da regressão no SVL) difere entre ecomorfos? A ANOVA filogenética compara o F observado com o de dados simulados por BM na árvore.

dados$res <- resid(m_bm)
set.seed(43)
phylANOVA(arv, setNames(dados$eco, dados$sp), setNames(dados$res, dados$sp), nsim = 999)$Pf
[1] 0.001001001

A ANCOVA filogenética põe o tamanho e o ecomorfo no mesmo modelo:

m_ancova <- gls(hll ~ svl + eco, data = dados, method = "ML",
                correlation = corBrownian(1, arv, form = ~sp))
anova(m_ancova)
Denom. DF: 75 
            numDF   F-value p-value
(Intercept)     1 19396.198  <.0001
svl             1  1275.876  <.0001
eco             5    25.425  <.0001

Os ecomorfos diferem no membro posterior relativo mesmo considerando a filogenia; o anova() do gls() testa os termos em sequência (primeiro o SVL, depois o ecomorfo).

Adams & Collyer (2018) mostram que permutar os resíduos (RRPP) tem mais poder que simular dados sob BM, e alertam para grupos agregados na árvore (aqui, cada ecomorfo surgiu várias vezes, nas várias ilhas). O procD.pgls() do geomorph faz a ANCOVA com RRPP:

library(geomorph)
gdf <- geomorph.data.frame(hll = setNames(dados$hll, dados$sp),
                           svl = setNames(dados$svl, dados$sp),
                           eco = setNames(dados$eco, dados$sp), phy = arv)
set.seed(45)
summary(procD.pgls(hll ~ svl + eco, phy = phy, data = gdf, iter = 999, print.progress = FALSE))$table
          Df      SS      MS     Rsq        F       Z Pr(>F)    
svl        1 1.64520 1.64520 0.86324 1275.876 10.1534  0.001 ***
eco        5 0.16393 0.03279 0.08601   25.425  6.5091  0.001 ***
Residuals 75 0.09671 0.00129 0.05074                            
Total     81 1.90583                                            
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se os grupos também diferirem na taxa de evolução, a ANCOVA filogenética de Fuentes-G. et al. (2016) estima isso junto com as diferenças de intercepto e inclinação.

6. Respostas binárias e contagens: o GLM filogenético

Ser trunk-ground (0 ou 1) depende do membro posterior relativo? E o número de lamelas (uma contagem) cresce com o tamanho?

dados$tg  <- as.integer(dados$eco == "TG")
dados$res_z <- as.numeric(scale(dados$res))       # padronizar ajuda o ajuste
dados$lam <- round(exp(anole.data$LAM))          # de volta para contagem
m_log <- phyloglm(tg ~ res_z, data = dados, phy = arv, method = "logistic_MPLE")
summary(m_log)$coefficients
             Estimate   StdErr   z.value     p.value
(Intercept) -1.925896 1.160165 -1.660019 0.096910725
res_z        2.419256 0.816681  2.962302 0.003053483
m_poi <- phyloglm(lam ~ svl, data = dados, phy = arv, method = "poisson_GEE")
summary(m_poi)$coefficients
             Estimate     StdErr   z.value      p.value
(Intercept) 1.0116451 0.21017338  4.813384 1.483959e-06
svl         0.4914623 0.04751698 10.342878 4.507929e-25

7. Quando a correlação muda na árvore

A covariância evolutiva entre dois atributos pode ser diferente entre regimes. Nos Centrarchidae (phytools), o modo de alimentação está mapeado na árvore: uma matriz de taxas para todos ou uma por regime?

data(sunfish.tree, sunfish.data)
vv <- evol.vcv(sunfish.tree, as.matrix(sunfish.data[, c("gape.width", "buccal.length")]))
vv
ML single-matrix model:
    R[1,1]  R[1,2]  R[2,2]  k   log(L)
fitted  0.114   0.033   0.0556  5   72.1893 

ML multi-matrix model:
    R[1,1]  R[1,2]  R[2,2]  k   log(L)
non 0.1387  -0.0022 0.012   8   81.2393 
pisc    0.0758  0.0796  0.1291  

P-value (based on X^2): 4e-04 

R thinks it has found the ML solution.

Duas matrizes ajustam melhor que uma: nos piscívoros, abertura da boca e comprimento bucal evoluem juntos (covariância positiva); nos não piscívoros, quase independentes.

8. Incerteza filogenética

A árvore também é uma estimativa. Uma forma simples de ver o efeito: refazer o PGLS em várias árvores (aqui, os mesmos ramos perturbados ao acaso, só para ilustrar; numa análise real, use uma amostra da posterior).

set.seed(44)
b <- replicate(100, {
  a2 <- arv
  a2$edge.length <- arv$edge.length * exp(rnorm(nrow(arv$edge), 0, .4))
  coef(gls(hll ~ svl, data = dados, correlation = corBrownian(1, a2, form = ~sp)))[["svl"]]
})
summary(b)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.022   1.044   1.049   1.050   1.056   1.087 
hist(b, main = "", xlab = "inclinação do PGLS em 100 árvores")
abline(v = coef(m_bm)[["svl"]], col = "firebrick", lwd = 2)

Resumo

  • O PGLS é um GLS com C; sob BM, dá a mesma inclinação que os contrastes.
  • Estime λ junto com a regressão e cheque os resíduos, não os atributos.
  • rr2::R2() para o R²; phylANOVA() e gls(... + fator) para ANOVA e ANCOVA; phyloglm() para respostas binárias e contagens.
  • evol.vcv() compara uma matriz de taxas com uma por regime; refaça a análise em várias árvores para ver a incerteza.

O PGLS e a ANOVA filogenética da edição de 2020 estão no roteiro da Aula 2, e o evol.vcv, no roteiro da Aula 4.