Modelo animal bivariado
Diogo Melo
Alastair Wilson
2026-05-06
Base genética da relação agressividade-tamanho em esquilos
Agora sabemos que o tamanho dos esquilos é herdável, além de estar sob seleção positiva, de modo que a evolução para um tamanho maior é prevista sob modelos simples (equação do criador univariada). No entanto, a seleção no tamanho corporal pode ter consequências evolutivas para outros aspectos do fenótipo, dependendo de estarem ou não geneticamente correlacionados. Aqui, testaremos se existe uma base genética para a (hipotetizada) covariância entre tamanho corporal e agressividade. A priori, poderíamos esperar que essas características covariem positivamente em nível fenotípico se (i) indivíduos agressivos forem capazes de conquistar mais recursos alimentares em competição (e, assim, crescerem mais rápido) e/ou (ii) a expressão da agressividade em relação a conespecíficos for mediada pelo tamanho relativo (ou seja, indivíduos maiores podem se dar ao luxo de agir agressivamente contra rivais menores, pois correm menor risco de ferimentos caso um conflito se agrave).
Retornaremos ao conjunto de dados que já vimos, estendendo nossa análise anterior de modelo animal para o caso bivariado. Como os programas gratuitos de REML que utilizamos anteriormente não executam modelos multivariados, mudaremos para uma abordagem MCMC usando o pacote MCMCglmm. Note que estaremos estimando uma matriz \(\mathbf{G}\) aqui, mas apenas para duas características. A extensão para qualquer número de características é possível em princípio. Na prática, modelos com mais características tornam-se mais difíceis de ajustar e levam mais tempo para rodar.
Preparação.
Antes de qualquer coisa, precisamos carregar os arquivos com os dados que vamos usar no resto do tutorial.
#1. read in the main data file
squirels<-read.table("https://raw.githubusercontent.com/diogro/evofencom/refs/heads/main/Tutoriais/squirels2.txt", header=T)
squirels$ID<-as.factor(squirels$ID)
squirels$sex<-as.factor(squirels$sex)
squirels$SIRE<-as.factor(squirels$SIRE)
squirels$DAM<-as.factor(squirels$DAM)
#read in the pedigree file
squirelsPED<-read.table("https://raw.githubusercontent.com/diogro/evofencom/refs/heads/main/Tutoriais/squirelsPED.txt",header=T)Vamos também carregar os pacotes que vamos usar.
#2. Automatically installs packages in this list that user does not already have
list_pkgs <- c("lme4","lmerTest","MCMCglmm","pedigreemm")
new_pkgs <- list_pkgs[!(list_pkgs %in% installed.packages()[,"Package"])]
if(length(new_pkgs) > 0){ install.packages(new_pkgs) }
library(MCMCglmm)
library(lme4)
library(lmerTest)
library(pedigreemm)Existe uma relação entre tamanho e agressividade?
Vale a pena fazer um gráfico com os valores de agressividade e tamanhoe para verificar se existe uma relação simples entre os dois fenótipos.
#3. Test for size-aggression relationship
mod_regression<-lm(aggression~size, data=squirels)
plot(squirels$size,squirels$aggression, pch = 19)
abline(mod_regression)
Parece haver uma relação entre os dois. Realmente, olhando pro summary do modelo linear, vemos que existe uma relação significativa entre os dois fenótipos.
##
## Call:
## lm(formula = aggression ~ size, data = squirels)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3234 -1.2258 0.0623 1.2960 5.5054
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.26370 0.78236 0.337 0.736
## size 0.30567 0.03647 8.381 5.39e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.881 on 498 degrees of freedom
## Multiple R-squared: 0.1236, Adjusted R-squared: 0.1219
## F-statistic: 70.25 on 1 and 498 DF, p-value: 5.393e-16Poderíamos também olhar pra uma regressão múltipla, incluindo os efeitos de tamanho e forrageamento, que talvez expliquem a relação entre essas duas variáveis.
##
## Call:
## lm(formula = aggression ~ size + sex + forage, data = squirels)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.2266 -1.2194 0.0691 1.3164 5.5223
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.08987 0.78875 0.114 0.9093
## size 0.31095 0.03726 8.345 7.11e-16 ***
## sex2 0.13925 0.17218 0.809 0.4191
## forage -0.19103 0.08583 -2.226 0.0265 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.874 on 496 degrees of freedom
## Multiple R-squared: 0.1341, Adjusted R-squared: 0.1289
## F-statistic: 25.61 on 3 and 496 DF, p-value: 2.043e-15A partir desse modelo, você poderia concluir que forrageadores melhores são, na verdade, menos agressivos quando condicionamos por tamanho e sexo; e que existe uma relação positiva entre tamanho e agressividade, mesmo quando condicionamos por sexo. Mas será que essa relação fenotípica que observamos é devida a efeitos genéticos? Vamos testar essa hipótese usando um modelo animal bivariado no pacote MCMCglmm.
Modelo animal no MCMCglmm
O MCMCglmm é um método bayesiano — precisamos especificar as distribuições a priori (priors) ou permitir que o pacote utilize seus priors padrão. A especificação de priors é uma questão complexa e não entraremos em detalhes aqui. As notas do curso de MCMCglmm de Jarrod Hadfield fornecem muito mais informações e a explicação completa dos parametros e opções disponíveis (uma leitura OBRIGATÓRIA se você deseja trabalhar com este pacote). Se você quiser saber mais MESMO, paga uma cerveja pro Diogo.
Vamos executar um modelo bivariado, mas, apenas para ver como o MCMCglmm funciona, vamos rodar primeiro um modelo univariado para a agressividade para verificar se ela é herdável (se não for, faz pouco sentido estimar uma correlação genética com o tamanho).
#4. Univariate animal models in mcmcglmm
#first we'll set a prior for the variance components
prior_uni <- list(G = list(G1 = list(V = 1, nu = 0.002)),
R = list(V = 1, nu = 0.002))
#I think the package likes to reserve ID for an identity matrix so we’ll duplicate
#our identity column and call it “animal” in the data file
squirels$animal<-as.factor(squirels$ID)
#then run the model
model_uni <- MCMCglmm(aggression ~ 1+sex+forage, random = ~animal, pedigree = squirelsPED,
data = squirels, prior = prior_uni)##
## MCMC iteration = 0
##
## MCMC iteration = 1000
##
## MCMC iteration = 2000
##
## MCMC iteration = 3000
##
## MCMC iteration = 4000
##
## MCMC iteration = 5000
##
## MCMC iteration = 6000
##
## MCMC iteration = 7000
##
## MCMC iteration = 8000
##
## MCMC iteration = 9000
##
## MCMC iteration = 10000
##
## MCMC iteration = 11000
##
## MCMC iteration = 12000
##
## MCMC iteration = 13000Novamente, usamos a função summary para obter os resultados do modelo.
##
## Iterations = 3001:12991
## Thinning interval = 10
## Sample size = 1000
##
## DIC: 1998.099
##
## G-structure: ~animal
##
## post.mean l-95% CI u-95% CI eff.samp
## animal 1.548 0.8976 2.248 343.7
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 2.381 1.832 2.858 623
##
## Location effects: aggression ~ 1 + sex + forage
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 6.522638 6.195084 6.835678 913 <0.001 ***
## sex2 0.319847 -0.008055 0.664824 1000 0.060 .
## forage -0.183663 -0.343690 -0.031379 1182 0.032 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Acima, vemos as estimativas da variança associada aos efeitos genéticos no termo “animal”, a variança residual no termo “units”, e estimativas dos efeitos fixos associados ao efeito de sexo e forragamento. Todas essas estimativas acompanham intervalos de confiança que podem ser usados para teste de hipótese.
Estimando a herdabilidades para agressividade
Uma vantagem grande de modelos bayesianos é que temos uma distribuição posterior e completa para os parâmetros. Portanto, é trivial calcular tanto a herdabilidade quanto o intervalo de confiança da herdabilidade, usando o resultado desse modelo.
#5.Get a posterior for h2
posterior.heritability<-model_uni$VCV[,1]/(model_uni$VCV[,1]+ model_uni$VCV[,2])Você pode olhar pra essa distribuição a posteriori dos valores de herdabilidade usando um histograma, ou usar a média ou a mediana como uma estimativa da herdabilidade, e usar essa distribuição para obter intervalo de confiança. (Um cuidado necessário é que a herdabilidade neste modelo é restrita a valores positivos. Então, o intervalo de confiança nunca vai interceptar o zero e, portanto, não pode ser usado para um teste de hipótese clássico.)
## [1] 0.3918977## lower upper
## var1 0.2601502 0.5402656
## attr(,"Probability")
## [1] 0.95Complicações
Ao usar MCMC, realmente precisamos realizar um trabalho de pós-ajuste do modelo para verificar se devemos confiar em seus resultados. Para que nossa a posteriori seja robusta, precisamos garantir que a cadeia de MCMC tenha convergido e que as distribuições a posteriori de nossas estimativas estejam estacionárias. Caso contrário, provavelmente precisaremos executar a cadeia por mais tempo do que o padrão e/ou alterar outras opções que controlam o processo de ajuste do modelo.
Provavelmente rodaríamos o modelo várias vezes para garantir que nossos resultados permaneçam consistentes (note que eles não serão idênticos como seriam sob REML, pois há variação de amostragem do MCMC aqui, e não uma solução única de máxima verossimilhança). Como bons bayesianos, provavelmente também deveríamos testar diferentes especificações de a priori, apenas para garantir que nossa a priori não esteja influenciando excessivamente a a posteriori (ou seja, que o resultado venha realmente dos dados, e não de nossa crença prévia sobre o que esperávamos obter).
Embora os pontos acima sejam extremamente importantes para análises reais que se pretenda publicar, aqui vamos apenas plotar as distribuições a posteriori para fazer uma verificação visual rápida da execução do MCMC!
#6. Some diagnostic model plots
plot(model_uni$Sol) #plots traces and posterior distributions of fixed effects
Nos trace plots, o ideal é ver uma variação totalmente aleatória ao redor de uma média, e que não tenha nenhum ponto em que os valores fiquem estagnados em um único valor ou saltos entre valores médios diferentes. Da mesma forma, nos gráficos de densidade, o ideal é que tenhamos uma distribuição que se pareça com uma gaussiana, com apenas um valor máximo.
Se você obtiver médias a posteriori instáveis e/ou distribuições com aspecto estranho, o primeiro passo habitual seria alterar as opções de execução do MCMCglmm em relação aos valores padrão (comprimento total de 13.000 iterações, sendo as primeiras 3.000 descartadas como burn-in, e as 10.000 restantes amostradas a cada 10 iterações — o intervalo de thinning —, resultando em uma amostra a posteriori de tamanho 1.000).
Por exemplo, você poderia utilizar um período de burn-in mais longo e mais iterações, além de amostrar a parte não descartada da cadeia com menos frequência para reduzir a autocorrelação. Naturalmente, isso vai demorar mais pra rodar.
Ajustando o modelo bivariado no MCMCglmm
Agora que tivemos um gostinho dos modelos MCMCglmm, vamos estender o modelo e calcular a matriz G entre o tamanho e a agressividade. É possível que esse modelo demore alguns minutos pra rodar no computador de vocês.
#8. Fitting a bivariate model, first specify a prior
prior_bi <- list(G = list(G1 = list(V = diag(2), nu = 1.002)),
R = list(V = diag(2), nu = 1.002))
# Aqui estamos usando a função cbind para definir uma resposta bivariada.
# A variável "trait", que aparece do lado direito da fórmula, é importante para que os efeitos de sexo, forrageamento, e os interceptos sejam diferentes entre os dois fenótipos analisados.
model_bi <- MCMCglmm(cbind(size, aggression) ~ trait + trait:sex +trait:forage - 1,
random = ~us(trait):animal,
rcov = ~us(trait):units, family = c("gaussian", "gaussian"),
pedigree = squirelsPED, data = squirels, prior = prior_bi,
nitt = 65000, thin = 50, burnin = 15000, verbose = TRUE)##
## MCMC iteration = 0
##
## MCMC iteration = 1000
##
## MCMC iteration = 2000
##
## MCMC iteration = 3000
##
## MCMC iteration = 4000
##
## MCMC iteration = 5000
##
## MCMC iteration = 6000
##
## MCMC iteration = 7000
##
## MCMC iteration = 8000
##
## MCMC iteration = 9000
##
## MCMC iteration = 10000
##
## MCMC iteration = 11000
##
## MCMC iteration = 12000
##
## MCMC iteration = 13000
##
## MCMC iteration = 14000
##
## MCMC iteration = 15000
##
## MCMC iteration = 16000
##
## MCMC iteration = 17000
##
## MCMC iteration = 18000
##
## MCMC iteration = 19000
##
## MCMC iteration = 20000
##
## MCMC iteration = 21000
##
## MCMC iteration = 22000
##
## MCMC iteration = 23000
##
## MCMC iteration = 24000
##
## MCMC iteration = 25000
##
## MCMC iteration = 26000
##
## MCMC iteration = 27000
##
## MCMC iteration = 28000
##
## MCMC iteration = 29000
##
## MCMC iteration = 30000
##
## MCMC iteration = 31000
##
## MCMC iteration = 32000
##
## MCMC iteration = 33000
##
## MCMC iteration = 34000
##
## MCMC iteration = 35000
##
## MCMC iteration = 36000
##
## MCMC iteration = 37000
##
## MCMC iteration = 38000
##
## MCMC iteration = 39000
##
## MCMC iteration = 40000
##
## MCMC iteration = 41000
##
## MCMC iteration = 42000
##
## MCMC iteration = 43000
##
## MCMC iteration = 44000
##
## MCMC iteration = 45000
##
## MCMC iteration = 46000
##
## MCMC iteration = 47000
##
## MCMC iteration = 48000
##
## MCMC iteration = 49000
##
## MCMC iteration = 50000
##
## MCMC iteration = 51000
##
## MCMC iteration = 52000
##
## MCMC iteration = 53000
##
## MCMC iteration = 54000
##
## MCMC iteration = 55000
##
## MCMC iteration = 56000
##
## MCMC iteration = 57000
##
## MCMC iteration = 58000
##
## MCMC iteration = 59000
##
## MCMC iteration = 60000
##
## MCMC iteration = 61000
##
## MCMC iteration = 62000
##
## MCMC iteration = 63000
##
## MCMC iteration = 64000
##
## MCMC iteration = 65000UMA PEQUENA DIGRESSÃO — Esse modelo provavelmente levou um minuto (ou vários?) para rodar. Um programa rápido de REML teria encontrado a solução de máxima verossimilhança para o modelo em segundos. Isso ilustra que o tempo computacional é uma desvantagem relativa do MCMC. Para conjuntos de dados pequenos e modelos simples como este, raramente isso é um problema; já para grandes conjuntos de dados e modelos complexos (por exemplo, muitas características e/ou efeitos aleatórios), as análises de MCMC podem levar dias ou semanas. No entanto, você acaba obtendo distribuições a posteriori que são muito úteis para propagar a incerteza para análises subsequentes. O MCMCglmm também possui capacidades muito mais avançadas (e robustas) para ajustar modelos generalizados do que os métodos baseados em verossimilhança. Portanto, há vantagens e desvantagens a considerar!
Vamos observar os gráficos das estimativas dos componentes de variância para ver se consideramos que a execução gerou resultados confiáveis. Reiterando: isso não é uma verificação de modelo suficiente para uma análise formal que gere resultados publicáveis. Leia a vignette e as notas do curso de MCMCglmm de Hadfield para entender o que deve ser feito adequadamente. Estamos pegando atalhos consideráveis aqui para economizar tempo.
#9. Plot variance covariance estimates for diagnostics
plot(model_bi$VCV) #plots traces and posterior distributions of random effect (co)variancesSe os trace plots estiverem razoáveis, vamos continuar olhando o modelo.
##
## Iterations = 15001:64951
## Thinning interval = 50
## Sample size = 1000
##
## DIC: 3997.363
##
## G-structure: ~us(trait):animal
##
## post.mean l-95% CI u-95% CI eff.samp
## traitsize:traitsize.animal 2.341 1.461 3.147 993.6
## traitaggression:traitsize.animal 1.661 1.063 2.253 1000.0
## traitsize:traitaggression.animal 1.661 1.063 2.253 1000.0
## traitaggression:traitaggression.animal 1.608 1.021 2.331 876.8
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitsize:traitsize.units 2.5788 2.0171 3.2088 906.9
## traitaggression:traitsize.units -0.2476 -0.6083 0.1649 1000.0
## traitsize:traitaggression.units -0.2476 -0.6083 0.1649 1000.0
## traitaggression:traitaggression.units 2.3373 1.9002 2.8877 1000.0
##
## Location effects: cbind(size, aggression) ~ trait + trait:sex + trait:forage - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitsize 20.87616 20.50133 21.26506 1000 <0.001 ***
## traitaggression 6.53207 6.22367 6.80319 1000 <0.001 ***
## traitsize:sex2 0.77043 0.42431 1.12593 1000 <0.001 ***
## traitaggression:sex2 0.32054 -0.01434 0.63831 1000 0.052 .
## traitsize:forage 0.25317 0.07247 0.41559 1000 <0.001 ***
## traitaggression:forage -0.17796 -0.34075 -0.03300 1000 0.030 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This output contains estimates of the elements in G, the corresponding residual (environment) variance-covariance matrix (which we could call R), and the fixed effects. In all cases the posterior mean is provided as the point estimate.
Esse sumário contém estimativa de todos os elementos da matriz G, os elementos da matriz residual, das variações ambientais e todos os efeitos fixos. Em todos os casos, a média é a média da posteriori, que te dá uma estimativa do valor do parâmetro.
Fazendo uma matriz G e calculando a correlação genética.
Na parte de “G-Structure” da tabela de sumário, os elementos estão na ordem: variância aditiva de tamanho, covariância entre tamanho e agressividade, covariância entre agressividade e tamanho, e variância aditiva de agressividade. Escrever o termo de covariança duas vezes pode parecer desnecessário, mas essa duplicidade facilita o trabalho com os objetos nas análises posteriores.
#11. Make a G matrix and calcular rG
G_elements<-c(mean(model_bi$VCV[,1]), mean(model_bi$VCV[,2]),
mean(model_bi$VCV[,3]), mean(model_bi$VCV[,4]))
#model_bi$VCV contains the posteriors (each 1000 values) of the
#eight var-covar parameters (VA1,COVA12, COVA12, VA2, VR1, COVR12, COVR12, VR2).
#Here we get the means of columns 1:4 to get point estimates of genetic
#parameters then arrange them in a matrix
G<-matrix(G_elements, nrow=2)
#Turn G into the correspoding correlation matrix
(Gcor<-cov2cor(G))## [,1] [,2]
## [1,] 1.0000000 0.8561676
## [2,] 0.8561676 1.0000000Comentários finais
No final, chegamos a uma estimativa de mais ou menos 0,85 para a correlação genética entre agressividade e tamanho. Essa é uma correlação genética bastante alta, compatível com dados usados num tutorial simples para uma disciplina introdutória. Alguns comentários aqui. Primeiro, note como, apesar do nome “Covariância Genética”, nós não usamos nenhuma informação genética sobre os indivíduos, apenas as suas relações de parentesco. Nesse sentido, esse é um modelo puramente correlativo: estamos atribuindo a similaridade entre indivíduos aparentados a uma similaridade genética ao chamar essa correlação de “correlação genética”. É sempre possível, e precisamos considerar essa possibilidade, que a similaridade entre parentes venha de outras causas, por exemplo, ambientes correlacionados. Isso é fundamental no estudo de populações humanas, por exemplo, em que uma fração grande da similaridade entre parentes vem de ambientes compartilhados. Você herda os genes e a biblioteca dos seus pais.