Regressão filogenética (PGLS)

Aula 5 · Métodos Filogenéticos Comparativos · PPGBA/UFMS

Diogo B. Provete

Hoje

1. Do OLS ao GLS
2. PGLS na prática
3. Diagnóstico
4. ANOVA e ANCOVA filogenéticas
5. Além do PGLS: GLM filogenético
6. Quando a correlação muda na árvore
7. Incerteza filogenética

Onde estamos

1. Do OLS ao GLS

Estão associados?

Antes, relembrando o modelo linear (OLS)

Os mínimos quadrados ordinários (ordinary least squares, OLS) escolhem a reta que minimiza a soma dos resíduos ao quadrado.

Pressupostos do OLS (Gauss-Markov)

Matemática	Tradução para língua humana
\(y = X\beta + \varepsilon\)	Relação linear entre \(y\) e \(X\)
\(X\) é uma matriz \(n \times k\) de posto cheio	As colunas de \(X\) são independentes
\(E(\varepsilon\varepsilon' \mid X) = \sigma^2 I\)	Homocedasticidade e ausência de autocorrelação
\(\varepsilon \mid X \sim N(0, \sigma^2 I)\)	Resíduos normais, com média zero e variância homogênea

A linha que ninguém desenha

Independência e homocedasticidade são duas afirmações sobre a mesma matriz:

\[ \operatorname{var}[\varepsilon_i \mid X] = \sigma^2 \qquad \operatorname{cov}[\varepsilon_i, \varepsilon_j \mid X] = 0 \;\; (i \neq j) \]

\[ E(\varepsilon\varepsilon' \mid X) = \begin{bmatrix} \sigma^2 & 0 & \cdots & 0 \\ 0 & \sigma^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sigma^2 \end{bmatrix} = \sigma^2 I \]

Reaproveitado do curso de Análise de Dados Univariados. Ver Zuur et al. (2009).

Toda violação é um buraco nessa matriz

\[ \underbrace{\sigma^2\begin{bmatrix}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{bmatrix}}_{\text{OLS}} \qquad \underbrace{\begin{bmatrix}\sigma_1^2&0&0&0\\0&\sigma_2^2&0&0\\0&0&\sigma_3^2&0\\0&0&0&\sigma_4^2\end{bmatrix}}_{\text{variância desigual}} \qquad \underbrace{\sigma^2\begin{bmatrix}1{,}44&1{,}26&0{,}20&0{,}20\\1{,}26&1{,}44&0{,}20&0{,}20\\0{,}20&0{,}20&1{,}44&1{,}42\\0{,}20&0{,}20&1{,}42&1{,}44\end{bmatrix}}_{\text{filogenia: } \sigma^2\mathbf{C}} \]

A filogenia preenche o que está fora da diagonal, com a matriz C da Aula 2 (aqui, as espécies 1 a 4 da árvore-guia).

O estimador GLS

\[\Large \hat\beta_{\text{GLS}} = \big(X' \color{#1c6e8c}{\mathbf{C}^{-1}} X\big)^{-1} X' \color{#1c6e8c}{\mathbf{C}^{-1}} y\]

Com \(\color{#1c6e8c}{\mathbf{C}} = I\) (uma estrela), volta a ser o OLS: \(\hat\beta = (X'X)^{-1}X'y\).

PGLS é o GLS com a matriz C da filogenia (ou uma transformação dela por um modelo evolutivo).

Filogenia, espaço e tempo são o mesmo problema estatístico, dados correlacionados: Ives & Zhu (2006), Ecol. Appl. 16: 20–32.

GLS na árvore-guia

OLS (C = I): inclinação 0,135; GLS com C: 1,398, a mesma dos contrastes independentes da Aula 2. Sob BM, a equivalência é exata.

Blomberg, Lefevre, Wells & Waterhouse (2012), Syst. Biol. 61: 382–391.

Quando usar o PGLS

• Normalmente, para testar a associação entre atributos
• Muito flexível: preditores contínuos, discretos (binários ou com vários estados) e combinações
• Resposta binária: regressão logística filogenética (phylolm::phyloglm)
• Pode incorporar diferentes modelos evolutivos, manipulando a estrutura de C

Grafen (1989), Phil. Trans. R. Soc. B 326: 119–157.

2. PGLS na prática

Estão associados?

A estimativa depende do modelo em C

Assume linearidade; a estimativa depende do modelo evolutivo (λ: Aula 3; OU: Aula 4). Nos Anolis, a inclinação muda pouco, mas em outros dados pode mudar muito.

Symonds & Blomberg (2014); ideia da figura: Mazel et al. (2016), Ecology 97: 286–293.

Erro comum num PGLS

Errado: testar o sinal filogenético na variável resposta ou na preditora para decidir se usa o PGLS.

Certo: o pressuposto é sobre os resíduos.

O melhor é ajustar um único modelo que estima, por máxima verossimilhança, o λ dos resíduos junto com os coeficientes da regressão.

Revell (2010), Methods Ecol. Evol. 1: 319–329.

Variáveis com sinal, resíduos sem sinal

X e Y simulados com sinal forte, mas Y = 1 + 1,2X + um erro sem estrutura filogenética: o λ dos resíduos é 0,00, e o OLS seria adequado.

Membros posteriores dos Anolis

O modelo em R

nlme::gls(hll ~ svl, data = dados, method = "ML",
          correlation = ape::corPagel(1, arvore, form = ~sp))

	Estimativa	EP	p
Intercepto	-0,499	0,194	= 0,012
log SVL	1,046	0,047	< 0,001

λ estimado: 1,013.

O corPagel pode dar λ um pouco acima de 1.

Quanto o modelo explica?

O R² do OLS não serve para modelos com resíduos correlacionados. O pacote rr2 calcula três R² para PGLS, phylolm e modelos mistos:

rr2::R2(pgls)

R²	Ideia	Anolis
R²lik	verossimilhança contra o modelo só com intercepto	0,957
R²resid	variância dos resíduos	0,961
R²pred	variância das predições, com a parte filogenética	0,961

Ives (2019), Syst. Biol. 68: 234–251; Ives & Li (2018), J. Open Source Softw. 3: 1028. R² do OLS nos mesmos dados: 0,867.

3. Diagnóstico

Estão associados?

Diagnosticando um PGLS

Use os resíduos normalizados (resid(modelo, type = "normalized")): os resíduos brutos de um GLS continuam correlacionados pela filogenia.

O grupo de resíduos negativos são os twig, de membros curtos: um preditor que falta no modelo (bloco 4).

Armadilhas comuns

• Nomes: as pontas da árvore e as linhas da tabela precisam bater (geiger::name.check())
• Ordem das linhas: no nlme::gls, diga qual coluna identifica a espécie (form = ~sp); sem isso, a correlação é aplicada na ordem das linhas
• Árvore: ultramétrica e com comprimentos de ramo (Aula 1)

4. ANOVA e ANCOVA filogenéticas

Estão associados?

Um preditor categórico: o ecomorfo

Membro posterior relativo: resíduo da regressão no SVL.

Ecomorfos dos Anolis: Losos (1990), Ecol. Monogr. 60: 369–388.

ANOVA filogenética

O valor crítico de F vem de simulações sob BM na árvore, e não da tabela F: F = 55,6, p = 0,001.

Garland, Dickerman, Janis & Jones (1993), Syst. Biol. 42: 265. Em R: phytools::phylANOVA.

Cuidados com a ANOVA filogenética

• Permutar os resíduos (RRPP) tem mais poder que simular dados sob BM na árvore, e serve para ANOVA e regressão, uni e multivariadas
• Grupos agregados na árvore (cada grupo é um clado) se confundem com a filogenia: o poder cai, e poucas transições entre grupos lembram a pseudorreplicação da Aula 6
• Em R: RRPP::lm.rrpp ou geomorph::procD.pgls

Adams & Collyer (2018), Evolution 72: 1204–1215.

ANCOVA filogenética

nlme::gls(hll ~ svl + eco, data = dados,
          correlation = ape::corBrownian(1, arvore, form = ~sp))

Termo	gl	F	p
log SVL	1	1275,9	< 0,001
ecomorfo	5	25,4	< 0,001

Os grupos podem diferir também na taxa de evolução, e não só no intercepto e na inclinação: a ANCOVA filogenética de Fuentes-G. et al. (2016) estima as duas coisas.

Revell & Harmon (2022), seção 3.5; Fuentes-G., Housworth, Weber & Martins (2016), Am. Nat. 188: 615–627.

5. Além do PGLS: GLM filogenético

Estão associados?

Respostas binárias e contagens

phylolm::phyloglm(tg ~ membro, data = dados, phy = arvore, method = "logistic_MPLE")
phylolm::phyloglm(lamelas ~ svl, data = dados, phy = arvore, method = "poisson_GEE")

Ives & Garland (2010); Paradis & Claude (2002); Ho & Ané (2014).

Lendo os dois modelos

Modelo	Coeficiente	Estimativa	p
Logística	membro posterior (DP)	2,42	= 0,003
Poisson	log SVL	0,49	< 0,001

DP: desvio padrão do membro posterior relativo. Trunk-ground têm membros posteriores relativamente longos; espécies maiores têm mais lamelas.

Lamelas: contagem reconstruída do log dos dados do phytools (round(exp(LAM)), de 13 a 42).

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

Estão associados?

Matrizes de covariância evolutiva por regime

E se a correlação entre dois atributos for diferente em partes da árvore (não estacionariedade)?

Revell & Collar (2009), Evolution 63: 1090–1100.

Versão bayesiana, com a incerteza das correlações: o pacote ratematrix.

Caetano & Harmon (2017), Methods Ecol. Evol. 8: 1920–1927.

Os peixes Centrarchidae

Os regimes (piscívoro ou não) são definidos a priori e mapeados nos ramos, como no BM com várias taxas (Aula 4).

Uma ou duas matrizes?

Modelo	var(abertura)	cov	var(bucal)	log L
Uma matriz	0,114	0,033	0,056	72,19
Não piscívoros	0,139	-0,002	0,012	81,24
Piscívoros	0,076	0,080	0,129

Razão de verossimilhança: p < 0,001. Nos piscívoros, os dois atributos evoluem juntos; nos outros, quase não covariam.

Calculado com phytools::evol.vcv nos dados sunfish do phytools. O log L de uma matriz bate com o do artigo (72,19); o de duas não (lá, 76,85), e as taxas estão em outra escala.

7. Incerteza filogenética

Estão associados?

A árvore também é uma estimativa

Refaça a análise numa amostra de árvores (uma posterior bayesiana, por exemplo) e relate a variação. Aqui, 100 árvores com ramos perturbados: inclinação de 1,023 a 1,083.

Ilustração simulada. Rangel et al. (2015), Evolution 69: 1301–1312.

Para fechar

Onde estamos

Leituras

• Revell & Harmon (2022), Phylogenetic Comparative Methods in R: cap. 3 (PGLS, ANOVA e ANCOVA filogenéticas)
• Symonds & Blomberg (2014): um guia do PGLS
• Revell (2010): sinal nos resíduos

Antes de usar qualquer método, leia o artigo que o propôs. Não fazer isso pode causar a rejeição do seu manuscrito.

Agora é com vocês

X e Y, um casal? Não: parentes

Belmiro de Almeida, Bom tempo ou Idílio campestre (1893), domínio público, via Wikimedia Commons. Prática: o walkthrough da Aula 5 e o desafio.