Desenhando árvores com mr bayes

Mr Bayes é um software para inferência bayesiana e escolha de modelos para uma grande variedade de modelos filogenéticos e evolutivos. Ele usa Markov Chain Monte Carlo (MCMC) para estimar a distribuição posterior dos parâmetros do modelo.
Basicamente ele tem um interface de linha de comando, que é similar a sintaxe usada por outro programa muito comum, o Paup, mas o Paup é um software proprietário e o Mr é um Software Livre.

Certo, bem ele pode ser conseguido aqui, e depois de baixar, no linux é só rodar um ./configure, sendo que temos que ter a biblioteca beagle instalada, senão temos que avisar o ./configure para ele não pedir o beagle, feito isso é só compilar com o make. Pro windows ele vem como um instalador, o que torna as coisas mais fácil, mas uma vez com ele instalado, podemos começar a dar uma olhada na sintaxe básica, ja que para como compilar, tem um capítulo inteiro no manual do Mr Bayes so de como compilar do fonte.

Bem, como entrada, usamos um alinhamento no formato nexus, permite guardar mais informações além das bases em si.

A cara dele é assim:

Captura de tela de 2016-05-31 15-38-30

Beleza, vamos usar o arquivo primatemitDNA.nexus, eu vou deixar ele disponível no meu repositório do github, mas basicamente é um alinhamento de 5 sequências de dna mitocondrial de primatas. É o arquivo do exercício do curso de Evolução Molecular que fiz no coursera.

Para abrir o arquivo no Mr Bayes, a gente usa o comando execute.

MrBayes > execute primatemitDNA.nexus Executing file "primatemitDNA.nexus" UNIX line termination Longest line length = 911 Parsing file Expecting NEXUS formatted file Reading data block Allocated taxon set Allocated matrix Defining new matrix with 5 taxa and 896 characters Data is Dna Gaps coded as - Taxon 1 -> Human Taxon 2 -> Chimpanzee Taxon 3 -> Gorilla Taxon 4 -> Orangutan Taxon 5 -> Gibbon Successfully read matrix Setting default partition (does not divide up characters) Setting model defaults Seed (for generating default start values) = 1464723572 Setting output file names to "primatemitDNA.nexus.run." Exiting data block Reached end of file MrBayes >

Veja que a gente abre as sequencias, de 911 bases, da para ver o nome das sequencias e algumas outras informações. A gente pode conferir ainda se os dados estão ok usando o comando ShowMatrix.

MrBayes > Showmatrix 1 Human AAGCTTCACCGGCGCAGTCATTCTCATAATCGCCCACGGACTTACATCCT Chimpanzee AAGCTTCACCGGCGCAATTATCCTCATAATCGCCCACGGACTTACATCCT Gorilla AAGCTTCACCGGCGCAGTTGTTCTTATAATTGCCCACGGACTTACATCAT Orangutan AAGCTTCACCGGCGCAACCACCCTCATGATTGCCCATGGACTCACATCCT Gibbon AAGCTTTACAGGTGCAACCGTCCTCATAATCGCCCACGGACTAACCTCTT . . . 851 Human ACTGACACTGAGCCACAACCCAAACAACCCAGCTCTCCCTAAGCTT Chimpanzee ACTGGCACTGAGCAACAACCCAAACAACCCAGCTCTCCCTAAGCTT Gorilla GCTGACACTGAGCAACAACCCAAACAATTCAACTCTCCCTAAGCTT Orangutan ACTGATGCTGAACAACCACCCAGACACTACAACTCTCACTAAGCTT Gibbon ACTGACACTGAACTGCAACCCAAACGCTAGAACTCTCCCTAAGCTT

Eu suprimi parte da saída, mas assim a gente olha os dados que foram carregados.
Agora vamos definir um “outgroup” para a árvore, para termos uma árvore com raiz

MrBayes > outgroup Gibbon Setting outgroup to taxon "Gibbon" e vamos definor o modelo que vamos usar MrBayes > lset nst=2 rates=gamma Setting Nst to 2 Setting Rates to Gamma Successfully set likelihood model parameters

Aqui estamos especificando que modelo evolutivo vamos usar. Nst é o número de parâmetros do modelo evolutivo, dizer que queremos 2 parâmetros aqui, dois tipos de substituições, usando o modelo de kimura, eu acho pelo menos, que distinguimos entre transições e transversões. Além disso, rates=gamma significa que nos queremos usar a distribuição gamma para representar diferentes taxas em diferentes partes da sequência.

Podemos mudar outras coisas, mas por enquanto vamos tentar fazer o básico, que é desenhar a árvore, podemos ver como é o modelo que está definido até agora com Showmodel, mas antes de iniciar a análise essa informação é mostrada, então não precisamos nos preocupar.

Para iniciar o MCMC, usamos o comando mcmc, com os parâmetros

ngen=100000 : tamanho da cadeia
samplefreq=100 : de quanto em quanto vamos salvar uma árvore e todos os parâmetros (amostras do distribuição posterior)

nchains=3 : número de cadeias paralelas, mas isso significa que uma é “fria”, onde as amostras são feitas e duas são “quentes” (“heated”) que se movem no espaço do parâmetros mais rapidamente para char picos na distribuição de probabilidade.

diagnfreq=5000 tem haver com ver se o MrBayes está rodando corretamente. Resumidamente, MrBayes vai começar duas sessões independentes, começando de diferentes árvores aleatórios. No começo, as duas sessões vão ter amostras bem diferentes, de árvores diferentes, mas quando elas convergirem (a distribuição posterior começar a ficar igual), as duas árvores devem ser bem parecidas, toda diagonfreq gerações, o MrBayes mede o quão similar as amostras das árvores estão e como uma regra geral (rule of thumb), nos vamos esperar esse valor ficar menor que 0.01.

Mas vamos la, vamos rodar o modelo.

MrBayes > mcmc ngen=100000 samplefreq=100 nchains=3 diagnfreq=5000 Setting number of generations to 100000 Setting sample frequency to 100 Setting number of chains to 3 Setting diagnosing frequency to 5000 Running Markov chain MCMC stamp = 4955053902 Seed = 557265027 Swapseed = 1464723572 Model settings: Data not partitioned -- Datatype = DNA Nucmodel = 4by4 Nst = 2 Transition and transversion rates, expressed as proportions of the rate sum, have a Beta(1.00,1.00) prior Covarion = No # States = 4 State frequencies have a Dirichlet prior (1.00,1.00,1.00,1.00) Rates = Gamma The distribution is approximated using 4 categories. Likelihood summarized over all rate categories in each generation. Shape parameter is exponentially distributed with parameter (1.00). Active parameters: Parameters --------------------- Tratio 1 Statefreq 2 Shape 3 Ratemultiplier 4 Topology 5 Brlens 6 --------------------- 1 -- Parameter = Tratio Type = Transition and transversion rates Prior = Beta(1.00,1.00) 2 -- Parameter = Pi Type = Stationary state frequencies Prior = Dirichlet 3 -- Parameter = Alpha Type = Shape of scaled gamma distribution of site rates Prior = Exponential(1.00) 4 -- Parameter = Ratemultiplier Type = Partition-specific rate multiplier Prior = Fixed(1.0) 5 -- Parameter = Tau Type = Topology Prior = All topologies equally probable a priori Subparam. = V 6 -- Parameter = V Type = Branch lengths Prior = Unconstrained:GammaDir(1.0,0.1000,1.0,1.0) The MCMC sampler will use the following moves: With prob. Chain will use move 1.89 % Dirichlet(Tratio) 0.94 % Dirichlet(Pi) 0.94 % Slider(Pi) 1.89 % Multiplier(Alpha) 9.43 % ExtSPR(Tau,V) 9.43 % ExtTBR(Tau,V) 9.43 % NNI(Tau,V) 9.43 % ParsSPR(Tau,V) 37.74 % Multiplier(V) 13.21 % Nodeslider(V) 5.66 % TLMultiplier(V) Division 1 has 87 unique site patterns Initializing conditional likelihoods Using standard SSE likelihood calculator for division 1 (single-precision) Initial log likelihoods and log prior probs for run 1: Chain 1 -- -3213.775230 -- 14.143053 Chain 2 -- -3181.879527 -- 14.143053 Chain 3 -- -3086.072235 -- 14.143053 Initial log likelihoods and log prior probs for run 2: Chain 1 -- -3177.836224 -- 14.143053 Chain 2 -- -3222.902892 -- 14.143053 Chain 3 -- -3230.891466 -- 14.143053 There are results from a previous run saved using the same filename(s). Do you want to overwrite these results? (yes/no): yes Overwriting file "primatemitDNA.nexus.mcmc" Overwriting file "primatemitDNA.nexus.run1.p" Overwriting file "primatemitDNA.nexus.run1.t" Overwriting file "primatemitDNA.nexus.run2.p" Overwriting file "primatemitDNA.nexus.run2.t" Using a relative burnin of 25.0 % for diagnostics Chain results (100000 generations requested): 0 -- [-3213.775] (-3181.880) (-3086.072) * [-3177.836] (-3222.903) (-3230.891) 1000 -- (-2648.611) (-2710.304) [-2622.754] * (-2681.408) (-2664.530) [-2620.152] -- 0:00:00 2000 -- [-2617.298] (-2620.320) (-2615.762) * [-2617.340] (-2627.230) (-2617.470) -- 0:00:00 3000 -- (-2618.457) (-2620.569) [-2620.056] * (-2618.272) (-2618.647) [-2616.704] -- 0:00:00 4000 -- [-2623.450] (-2617.093) (-2620.776) * (-2624.657) (-2617.662) [-2618.030] -- 0:00:24 5000 -- (-2615.517) [-2616.710] (-2626.709) * [-2617.718] (-2615.441) (-2618.128) -- 0:00:19 . . . Average standard deviation of split frequencies: 0.005447 96000 -- [-2616.134] (-2616.792) (-2615.722) * (-2621.091) (-2624.519) [-2616.766] -- 0:00:00 97000 -- (-2617.755) [-2614.612] (-2617.002) * (-2621.728) [-2614.617] (-2616.776) -- 0:00:00 98000 -- [-2613.363] (-2624.533) (-2620.102) * (-2616.026) (-2619.449) [-2615.646] -- 0:00:00 99000 -- (-2617.593) [-2614.617] (-2617.432) * (-2624.591) [-2617.346] (-2620.014) -- 0:00:00 100000 -- [-2620.524] (-2620.783) (-2616.498) * (-2618.557) [-2616.113] (-2623.220) -- 0:00:00 Average standard deviation of split frequencies: 0.006120 Continue with analysis? (yes/no): no Analysis completed in 13 seconds Analysis used 13.26 seconds of CPU time Likelihood of best state for "cold" chain of run 1 was -2611.84 Likelihood of best state for "cold" chain of run 2 was -2611.87 Acceptance rates for the moves in the "cold" chain of run 1: With prob. (last 100) chain accepted proposals by move 22.7 % ( 18 %) Dirichlet(Tratio) 22.4 % ( 34 %) Dirichlet(Pi) 25.9 % ( 27 %) Slider(Pi) 35.1 % ( 34 %) Multiplier(Alpha) 1.2 % ( 0 %) ExtSPR(Tau,V) 1.8 % ( 3 %) ExtTBR(Tau,V) 1.8 % ( 3 %) NNI(Tau,V) 1.6 % ( 3 %) ParsSPR(Tau,V) 25.2 % ( 20 %) Multiplier(V) 26.4 % ( 22 %) Nodeslider(V) 25.7 % ( 26 %) TLMultiplier(V) Acceptance rates for the moves in the "cold" chain of run 2: With prob. (last 100) chain accepted proposals by move 22.7 % ( 23 %) Dirichlet(Tratio) 20.9 % ( 21 %) Dirichlet(Pi) 27.1 % ( 26 %) Slider(Pi) 34.2 % ( 28 %) Multiplier(Alpha) 1.4 % ( 0 %) ExtSPR(Tau,V) 1.3 % ( 4 %) ExtTBR(Tau,V) 1.7 % ( 0 %) NNI(Tau,V) 1.1 % ( 1 %) ParsSPR(Tau,V) 25.7 % ( 17 %) Multiplier(V) 25.8 % ( 25 %) Nodeslider(V) 25.1 % ( 28 %) TLMultiplier(V) Chain swap information for run 1: 1 2 3 ----------------------- 1 | 0.86 0.73 2 | 33001 0.87 3 | 33411 33588 Chain swap information for run 2: 1 2 3 ----------------------- 1 | 0.86 0.73 2 | 33298 0.87 3 | 33395 33307 Upper diagonal: Proportion of successful state exchanges between chains Lower diagonal: Number of attempted state exchanges between chains Chain information: ID -- Heat ----------- 1 -- 1.00 (cold chain) 2 -- 0.91 3 -- 0.83 Heat = 1 / (1 + T * (ID - 1)) (where T = 0.10 is the temperature and ID is the chain number) MrBayes >

Primeiro temos um report de todo o modelo que estamos usando, inicializamos as cadeias, e particularmente aqui, como é a segunda vez que estava rodando, ja tem um resultado pronto, mas eu peço para rodar denovo, e temos os passos do mcmc, e depois do número de passos definidos, ele pergunta se queremos continuar, como estamos com uma pequena diferença entre as duas sessões falamos que nao precisa, e pronto, ja temos a cadeia do mcmc. Arquivos foram gerados na pasta que estamos trabalhando, mas ainda não temos que sumarizar os resultados. E podemos fazer isso com o comando sumt

MrBayes > sumt contype=halfcompat conformat=simple relburnin=yes burninfrac=0.25 showtreeprobs=yes Setting sumt contype to Halfcompat Setting sumt conformat to Simple Using relative burnin (a fraction of samples discarded). Setting burnin fraction to 0.25 Setting showtreeprobs to yes Summarizing trees in files "primatemitDNA.nexus.run1.t" and "primatemitDNA.nexus.run2.t" Using relative burnin ('relburnin=yes'), discarding the first 25 % of sampled trees Writing statistics to files primatemitDNA.nexus. Examining first file ... Found one tree block in file "primatemitDNA.nexus.run1.t" with 1001 trees in last block Expecting the same number of trees in the last tree block of all files Tree reading status: 0 10 20 30 40 50 60 70 80 90 100 v-------v-------v-------v-------v-------v-------v-------v-------v-------v-------v ********************************************************************************* Read a total of 2002 trees in 2 files (sampling 1502 of them) (Each file contained 1001 trees of which 751 were sampled) Overwriting file "primatemitDNA.nexus.parts" Overwriting file "primatemitDNA.nexus.tstat" Overwriting file "primatemitDNA.nexus.vstat" Overwriting file "primatemitDNA.nexus.con.tre" Overwriting file "primatemitDNA.nexus.trprobs" General explanation: In an unrooted tree, a taxon bipartition (split) is specified by removing a branch, thereby dividing the species into those to the left and those to the right of the branch. Here, taxa to one side of the removed branch are denoted '.' and those to the other side are denoted '*'. Specifically, the '.' symbol is used for the taxa on the same side as the outgroup. In a rooted or clock tree, the tree is rooted using the model and not by reference to an outgroup. Each bipartition therefore corresponds to a clade, that is, a group that includes all the descendants of a particular branch in the tree. Taxa that are included in each clade are denoted using '*', and taxa that are not included are denoted using the '.' symbol. The output first includes a key to all the bipartitions with frequency larger or equual to (Minpartfreq) in at least one run. Minpartfreq is a parameter to sumt command and currently it is set to 0.10. This is followed by a table with statistics for the informative bipartitions (those including at least two taxa), sorted from highest to lowest probability. For each bipartition, the table gives the number of times the partition or split was observed in all runs (#obs) and the posterior probability of the bipartition (Probab.), which is the same as the split frequency. If several runs are summarized, this is followed by the minimum split frequency (Min(s)), the maximum frequency (Max(s)), and the standard deviation of frequencies (Stddev(s)) across runs. The latter value should approach 0 for all bipartitions as MCMC runs converge. This is followed by a table summarizing branch lengths, node heights (if a clock model was used) and relaxed clock parameters (if a relaxed clock model was used). The mean, variance, and 95 % credible interval are given for each of these parameters. If several runs are summarized, the potential scale reduction factor (PSRF) is also given; it should approach 1 as runs converge. Node heights will take calibration points into account, if such points were used in the analysis. Note that Stddev may be unreliable if the partition is not present in all runs (the last column indicates the number of runs that sampled the partition if more than one run is summarized). The PSRF is not calculated at all if the partition is not present in all runs.The PSRF is also sensitive to small sample sizes and it should only be considered a rough guide to convergence since some of the assumptions allowing one to interpret it as a true potential scale reduction factor are violated in MrBayes. List of taxa in bipartitions: 1 -- Human 2 -- Chimpanzee 3 -- Gorilla 4 -- Orangutan 5 -- Gibbon Key to taxon bipartitions (saved to file "primatemitDNA.nexus.parts"): ID -- Partition ----------- 1 -- *.... 2 -- .*... 3 -- ..*.. 4 -- ...*. 5 -- ****. 6 -- ***.. 7 -- **... ----------- Summary statistics for informative taxon bipartitions (saved to file "primatemitDNA.nexus.tstat"): ID #obs Probab. Sd(s)+ Min(s) Max(s) Nruns ---------------------------------------------------------------- 6 1501 0.999334 0.000942 0.998668 1.000000 2 7 1464 0.974700 0.011299 0.966711 0.982690 2 ---------------------------------------------------------------- + Convergence diagnostic (standard deviation of split frequencies) should approach 0.0 as runs converge. Summary statistics for branch and node parameters (saved to file "primatemitDNA.nexus.vstat"): 95% HPD Interval -------------------- Parameter Mean Variance Lower Upper Median PSRF+ Nruns ------------------------------------------------------------------------------------- length[1] 0.059671 0.000347 0.026276 0.095918 0.058513 1.001 2 length[2] 0.080990 0.000482 0.041008 0.126153 0.078997 1.000 2 length[3] 0.075218 0.000727 0.027528 0.127487 0.071678 1.000 2 length[4] 0.344477 0.015647 0.154275 0.596846 0.322574 1.002 2 length[5] 0.524957 0.031039 0.246341 0.882645 0.490579 1.001 2 length[6] 0.175448 0.005446 0.054393 0.326580 0.163428 0.999 2 length[7] 0.043979 0.000544 0.005098 0.088742 0.040356 1.001 2 ------------------------------------------------------------------------------------- + Convergence diagnostic (PSRF = Potential Scale Reduction Factor; Gelman and Rubin, 1992) should approach 1.0 as runs converge. NA is reported when deviation of parameter values within all runs is 0 or when a parameter value (a branch length, for instance) is not sampled in all runs. Summary statistics for partitions with frequency >= 0.10 in at least one run: Average standard deviation of split frequencies = 0.006120 Maximum standard deviation of split frequencies = 0.011299 Average PSRF for parameter values (excluding NA and >10.0) = 1.001 Maximum PSRF for parameter values = 1.002 Clade credibility values: /---------------------------------------------------------------- Gibbon (5) | |---------------------------------------------------------------- Orangutan (4) | + /--------------------- Human (1) | /----------97---------+ | | \--------------------- Chimpanzee (2) \---------100--------+ \------------------------------------------- Gorilla (3) Phylogram (based on average branch lengths): /-------------------------------------------------------------------- Gibbon (5) | |--------------------------------------------- Orangutan (4) | + /-------- Human (1) | /----+ | | \----------- Chimpanzee (2) \----------------------+ \---------- Gorilla (3) |------------| 0.100 expected changes per site Calculating tree probabilities... Tree 1 (p = 0.974, P = 0.974): /---------------------------------------------------------- | | /------------------ | /-------------------+ +-------------------+ \------------------ | | | \-------------------------------------- | \---------------------------------------------------------- Tree 2 (p = 0.015, P = 0.989): /------------------ /-------------------+ /-------------------+ \------------------ | | | \-------------------------------------- | +---------------------------------------------------------- | \---------------------------------------------------------- Tree 3 (p = 0.010, P = 0.999): /---------------------------- /-----------------------------+ | \---------------------------- | | /---------------------------- +-----------------------------+ | \---------------------------- | \---------------------------------------------------------- Tree 4 (p = 0.001, P = 1.000): /---------------------------------------------------------- | | /-------------------------------------- +-------------------+ | | /------------------ | \-------------------+ | \------------------ | \---------------------------------------------------------- Credible sets of trees (4 trees sampled): 99 % credible set contains 3 trees MrBayes >

Esse comando da um sumario das árvores estimadas. A opção contype=halfcompat faz com que a a arvore consensus seja calculada a partir das árvores que não foram descartadas no burnin. A arvore consensus é primeiro plotada na tela, abaixo está o cladograma, um phylograma consensual. As distâncias das arestas é a media sobre as árvores que tiveram essas arestas presentes. O cladograma também tem um valor de credibilidade.

O que nos interessa é a lista de árvores que está depois do phylogram. Essas são nomeadas como Tree 1, Tree 2, etc, e estão ordenadas apartir da sua probabilidade posterior, que é indicada pelo p minisculo depois de cada árvore. O P maisculo é a probabilidade acumulada, diferente de desenhar árvores com ML, aqui temos uma lista de árvores e suas probabilidades associadas.

Bem é isso ai, arquivo do alinhamento vai estar la no repositório recologia na pasta dados, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail, e hoje não tem script, são só meia duzia de comandos.

Referência:

Anders Gorm Pedersen – Computational Molecular Evolution

Relação espécie área

A relação entre o número de espécies encontradas em amostras de diferentes áreas é um padrão a muito tempo conhecido, relativamente simples de pensar no processo que gera esse tipo de padrão, já que áreas maiores cabem mais coisas, as vezes meio confuso, principalmente quando a gente realiza que áreas maiores podem ter tanto maior abundancia de espécies assim como mais espécies. Mas areas maiores possem mais de muitos tipo de recursos, o que diminui a chance de extinção, aumentam colonização e muitas possibilidade de manter mais espécies, mas não estamos falando de biogeografia de ilhas, ainda.

Certo, mas como isso funciona? Geralmente é um padrão empírico do número de espécies encontrados em manchas de diferentes tamanhos, plotado contra como uma função dos seus respectivos tamanhos, podendo essas manchas serem ilhas oceânicas ou uma mancha de vegetação numa matriz de campo aberto.

Quantitativamente, essa relação é representada como uma potencia (power law).

R=cA^z

onde o R é o número de espécies em uma mancha de tamanho A, e os c e z são constantes ajustaveis. Mas normalmente a gente ve essa lei plotada como uma relação de log-log, que a torna linear.

log(R)=b+zA

So aplica log dos dois lados.

Podemos gerar um plot aqui, para ter uma ideia melhor.

1
2
3
4
5
6
7
8
A<-10^(1:10)
c<-1.5
z<-0.25
 
fn<-function(x,c=1.5,z=0.25) {
    return(c*x^z)
}
ggplot(data.frame(x=c(10, 10^10)), aes(x)) + stat_function(fun=fn) + xlab(label="Área(ha)") + ylab(label="Número de espécies")

E isso gera a seguinte figura.

01

Essa é a cara da relação, como uma lei de potência, e podemos linearizar, transformando para log os dois eixos.

02

E a parte boa, é que linearizando, tudo fica mais fácil, em quesito de modelagem, porque agora é so passar uma reta. Para tanto, vamos dar uma olhada no artigo do Johnson 1968, na table I do paper.

> exemplo Location Area Species 1 Tiburon Peninsula 5.9 370 2 San Francisco 45.0 640 3 Santa Barbara area 110.0 680 4 Santa Monica Mountains 320.0 640 5 Marin County 529.0 1060 6 Santa Cruz Mountains 1386.0 1200 7 Monterey County 3324.0 1400 8 San Diego County 4260.0 1450 9 California Coast 24520.0 2525

Temos a área e o número de espécies encontradas, em manchas de vegetação. Vamos fazer um plot dos dados e dar uma olhada.

03

Parece que existe ai um padrão de espécie área bem claro, então a gente pode ajustar um power law nesses dados e extrair os valores constantes.

1
modelo<-lm(log(Species,10)~log(Area,10),data=exemplo)

E olhamos o resultado.

> summary(modelo) Call: lm(formula = log(Species, 10) ~ log(Area, 10), data = exemplo) Residuals: Min 1Q Median 3Q Max -0.129947 -0.013386 0.003153 0.041204 0.057274 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.38560 0.05568 42.85 9.84e-10 *** log(Area, 10) 0.21976 0.01917 11.46 8.64e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.06025 on 7 degrees of freedom Multiple R-squared: 0.9494, Adjusted R-squared: 0.9422 F-statistic: 131.4 on 1 and 7 DF, p-value: 8.642e-06

Agora aqui é importante para interpretar, a gente pensar em algumas coisa. A gente usou uma base 10 nos logs, isso porque é mais fácil de pensar, só pensar que log de 10 é 1, de 100 é 2, de 1000 é 3, melhor que aqueles números quebrados do log natural, só facilita visualizar as coisas a base do log aqui. Além disso o que interessa é a inclinação, que representa o nosso valor de potência, o quão rápido essa curva cresce.

Podemos ajustar como um modelo não linear também, para ver como fica o resultado.

1
modelo_nls<-nls(Species~a*Area^z,start=list(a=1,z=0.25),data=exemplo)
> summary(modelo_nls) Formula: Species ~ a * Area^z Parameters: Estimate Std. Error t value Pr(>|t|) a 195.7664 31.4382 6.227 0.000434 *** z 0.2494 0.0186 13.407 3.01e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 122 on 7 degrees of freedom Number of iterations to convergence: 4 Achieved convergence tolerance: 3.028e-06

E antes de ficar olhando os números, vamos ver como essas curvas ficam em relação aos dados.

04

O ajuste fica bem bom, parece predizer bem a quantidade de espécies encontradas, apesar de que os valores não são exatamente iguais, mas olhando os intervalos de confiança.

> confint(modelo) 2.5 % 97.5 % (Intercept) 2.2539475 2.5172501 log(Area, 10) 0.1744247 0.2650924 > confint(modelo_nls) Waiting for profiling to be done... 2.5% 97.5% a 128.9454047 283.572614 z 0.2055227 0.296902

A gente ve que ambos pegam muitos valores semelhantes, apesar do modelo linear ter fica com a média menor.
Bem um nome recorrente quando falamos de espécie área é o Preston, que propos que dado dado a distribuição de espécies log-normal, então áreas maiores devem acumular mais espécies em taxas específicas.

Bem é isso ai, o script vai estar la no repositório recologia, os gráficos com ggplot ainda são meia boca, mas é treinando que se aprende e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Johnson, Mason Raven 1968 – Ecological Parameters and Plant Species DiversityThe American Naturalist Vol. 102(926) pp. 297-306

M.H.H. Stevens 2011, A Primer of Ecology with R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
library(ggplot2)
 
###Desenhando as curvas da relação espécie área
A<-10^(1:10)
c<-1.5
z<-0.25
 
fn<-function(x,c=1.5,z=0.25) {
    return(c*x^z)
}
ggplot(data.frame(x=c(10, 10^10)), aes(x)) + stat_function(fun=fn) + xlab(label="Área(ha)") + ylab(label="Número de espécies")
ggsave("01.jpg")
 
 
fn_log<-function(x,c=1.5,z=0.25) {
    return(log(c,10)+z*x)
}
ggplot(data.frame(x=c(1, 10)), aes(x)) + stat_function(fun=fn_log) + xlab(label="Área(ha)") + ylab(label="Número de espécies")
ggsave("02.jpg")
 
###Johnson, Mason, and Raven 1968
exemplo<-data.frame(Location = c("Tiburon Peninsula", "San Francisco","Santa Barbara area", "Santa Monica Mountains", "Marin County", 
                                 "Santa Cruz Mountains", "Monterey County", "San Diego County","California Coast"),
                    Area = c(5.9, 45, 110, 320, 529, 1386, 3324,4260, 24520),
                    Species = c(370L, 640L, 680L, 640L, 1060L, 1200L,1400L, 1450L, 2525L))
exemplo
 
 
ggplot(data=exemplo,aes(x=Area,y=Species))+ geom_point() + xlab(label="Área(ha)") + ylab(label="Número de espécies")
ggsave("03.jpg")
 
log(1000,10)
 
modelo<-lm(log(Species,10)~log(Area,10),data=exemplo)
summary(modelo)
 
modelo_nls<-nls(Species~a*Area^z,start=list(a=1,z=0.25),data=exemplo)
summary(modelo_nls)
 
linha<-data.frame(x=log10(seq(1,25000,100)),y=log10(predict(modelo_nls,newdata=data.frame(Area=seq(1,25000,100)))))
 
 
 
ggplot(data=exemplo,aes(x=log10(Area),y=log10(Species),linetype="solid") )  + geom_point() +
    geom_abline(intercept = modelo$coefficients[1], slope = modelo$coefficients[2]) +
    geom_line(data=linha,aes(x=x,y=y,linetype="dashed"))+        
    xlab(label="Área(ha)") + ylab(label="Número de espécies")+
    scale_linetype_manual(name='Modelo',values =c('dashed','solid'), labels = c('Não Linear','Linear'))
ggsave("04.jpg")
 
confint(modelo)
confint(modelo_nls)

Gráficos circulares usando ggplot2

A gente começou a olhar o sistema de gráficos ggplot2 aqui.

Não tem como aprender nada novo sem começar a usar, então vamos tentar ver como fazer alguns gráficos simples, e um jeito legal de transformá-los em gráficos circulares. Porque se a gente fizer gráficos o suficiente, uma hora a gente fica bom de ggplot2, e como eu disse no outro post, ele tem lá suas vantagens, apesar de que qualquer coisa que da para fazer com ggplot, da para fazer com o sistema de gráficos default do R.

Mas primeiro vamos gerar um conjunto de dados bem simples:

1
dados<-data.frame(preditor=1:20,resposta=runif(20,5,10),grupo=rep(1:4,each=5),linha=runif(20,8,12))

Então temos uma coluna com uma sequencia de um a vinte, uma coluna com números aleatórios, um coluna agrupando eles e outra coluna de números aleatórios.

> head(dados) preditor resposta grupo linha 1 1 8.501650 1 8.049639 2 2 7.122299 1 8.521307 3 3 8.113232 1 9.414598 4 4 5.568166 1 8.886906 5 5 8.067636 1 11.594324 6 6 6.692768 2 8.299302

Primeiro, vamos fazer uma figura bem simples, um gráfico de barras.

1
ggplot(data=dados, aes(x=preditor, y=resposta)) + geom_bar(stat="identity")

O comando acima diz o seguinte, não função ggplot, a gente diz onde estão os dados, e o aesthetic, ou quem é quem. Mas até ai não há representação nenhuma, ou seja, precisamos de uma forma para representar os dados, e para isso vamos usar o geom_bar, que são barras, então barras representam os dados, e vamos usar a função identidade, ou seja, não vamos alterar em nada a apresentação, isso vai nos render a seguinte figura.

01

Bem simples até de usar não, parecia mais complicado no post anterior. Então, mas o aesthetic não se limita a coordenadas, podemos definir o preenchimento das barras. Dependendo do tipo de forma que presentara os dados, teremos os aesthetics de interesse, formato de linha, formato de ponto, essas coisas, mas vamos continuar nas barras aqui.

1
ggplot(data=dados, aes(x=preditor, y=resposta,fill=grupo)) + geom_bar(stat="identity")

Veja que agora somente adicionamos o fill ao aesthetic e definimos ele como o grupo.

02

O que está dentro do aesthetic do ggplot, é o que será usado de forma geral, veja que dentro desse geom que estamos usando, podemos definir característica da sua representação também, por exemplo um contorno preto.

1
ggplot(data=dados, aes(x=preditor, y=resposta,fill=grupo)) + geom_bar(colour="black",stat="identity") + guides(fill=FALSE)

Além disso, alguma coisas, como a legenda, vem por padrão, mas podemos retirar ela, se desejarmos, que é o que fazemos ao dar o guides como false.

03

Legal, agora vamos ver como fazer um gráfico de linhas.

1
ggplot(data=dados, aes(x=preditor, y=linha)) + geom_line() + guides(fill=FALSE)

Linhas são apenas outra forma de representação dos dados, so trocamos o geom, a sintaxe do aesthetic anterior continua a mesma coisa.

04

Mas é legal, que se queremos usar linhas e pontos, por exemplo, o que a gente quer são duas representações geométricas dos mesmos valores, então basta adicionar os dois geom.

1
ggplot(data=dados, aes(x=preditor, y=linha)) + geom_line() + geom_point() + guides(fill=FALSE)

E temos a seguinte figura

05

Agora e se quisermos adicionar dados diferentes de acordo com o geom, por exemplo um gráfico com barras e linhas?
Como falamos até agora, isso implica em mais de uma aesthetic, mas podemos definir mais de um aesthetic, dentro de cada geom, podemos definir seu próprio aesthetic, se não definimos, ele usa o padrão definido dentro da função ggplot.

1
2
ggplot(data=dados, aes(x=preditor, y=resposta,fill=grupo)) + geom_bar(colour="black",stat="identity") + guides(fill=FALSE) +
    geom_line(data=dados, aes(x=preditor, y=linha)) + geom_point(data=dados, aes(x=preditor, y=linha))

E temos a seguinte figura.

06

Uma coisa legal, é que como tudo é interpretado de uma so vez, o ggplot vai tentar colocar tudo que você quer na mesma figura, por exemplo, ele tenta acertar a escala, para caber todo mundo, diferente do sistema gráfico padrão do R, que aqui exigiria dois passos, e nos teríamos que nos preocupar com a escala da figura também, que teria seu limites definido baseado em um dos aesthetics.

Mas como tornar essa figura circular? Aqui está a parte legal, basta trocar o sistema de coordenadas, nesse caso, por padrão, estamos usando o plano cartesiano, basta mudar para o plano circular e ta pronto, uma figura circular.

1
2
ggplot(data=dados, aes(x=preditor, y=resposta,fill=grupo)) + geom_bar(colour="black",stat="identity") + guides(fill=FALSE) +
    geom_line(data=dados, aes(x=preditor, y=linha)) + geom_point(data=dados, aes(x=preditor, y=linha))+coord_polar()

Veja que o que mudou, foi a adição do coord_polar, e ficamos com o seguinte.

07

Agora aqui é preciso atenção, veja que as barras, e os pontos, bem os pontos são o centro das barras, então quando fazemos uma figura de barras, os pontos não vão fechar, veja que so mudamos o plano, o que era a extensão do eixo x virou ângulos, de 0 a 360, e o eixo y virou raio.

Se tiramos as barras

1
ggplot(data=dados, aes(x=preditor, y=linha)) + geom_line() + geom_point() + guides(fill=FALSE) + coord_polar()

Temos que os pontos começam e terminam no “mesmo lugar”.

08.

Mas temos que ter atenção ainda, porque não da para ver direito o eixo y, no sentido de onde começa e termina em valores, mas ele nao esta começando do zero, temos que delimitar o inicio e fim dele.

1
ggplot(data=dados, aes(x=preditor, y=linha)) + geom_line() + geom_point() + guides(fill=FALSE) + coord_polar() + ylim(0,12)

Agora o ponto central é zero.

09

Bem é isso ai, esses são alguns teste para aprender alguma coisa de ggplot2, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
library("ggplot2")
 
dados<-data.frame(preditor=1:20,resposta=runif(20,5,10),grupo=rep(1:4,each=5),linha=runif(20,8,12))
head(dados)
 
## Gráfico de barras
ggplot(data=dados, aes(x=preditor, y=resposta)) + geom_bar(stat="identity")
ggsave("01.jpg")
 
## Separando por grupos
ggplot(data=dados, aes(x=preditor, y=resposta,fill=grupo)) + geom_bar(stat="identity")
##ggplot(data=dados, aes(x=preditor, y=resposta)) + geom_bar(aes(fill=grupo),stat="identity")
ggsave("02.jpg")
 
## Adicionando um contorno preto, e retirando a legenda
ggplot(data=dados, aes(x=preditor, y=resposta,fill=grupo)) + geom_bar(colour="black",stat="identity") + guides(fill=FALSE)
ggsave("03.jpg")
 
##Gráfico de linhas
ggplot(data=dados, aes(x=preditor, y=linha)) + geom_line() + guides(fill=FALSE)
ggsave("04.jpg")
ggplot(data=dados, aes(x=preditor, y=linha)) + geom_line() + geom_point() + guides(fill=FALSE)
ggsave("05.jpg")
 
##Combinando os dois
ggplot(data=dados, aes(x=preditor, y=resposta,fill=grupo)) + geom_bar(colour="black",stat="identity") + guides(fill=FALSE) +
    geom_line(data=dados, aes(x=preditor, y=linha)) + geom_point(data=dados, aes(x=preditor, y=linha))
ggsave("06.jpg")
 
 
##Transformando em uma figura circular
ggplot(data=dados, aes(x=preditor, y=resposta,fill=grupo)) + geom_bar(colour="black",stat="identity") + guides(fill=FALSE) +
    geom_line(data=dados, aes(x=preditor, y=linha)) + geom_point(data=dados, aes(x=preditor, y=linha))+coord_polar()
ggsave("07.jpg")
 
ggplot(data=dados, aes(x=preditor, y=linha)) + geom_line() + geom_point() + guides(fill=FALSE) + coord_polar()
ggsave("08.jpg")
ggplot(data=dados, aes(x=preditor, y=linha)) + geom_line() + geom_point() + guides(fill=FALSE) + coord_polar() + ylim(0,12)
ggsave("09.jpg")

Método de Monte Carlo

De forma geral, o método de Monte Carlo é um método estatístico, experimental de calcular integrais usando posições aleatórias, chamadas amostras, cuja distribuição tem que ser escolhidas com muito cuidado. Aleatória, ou “Random” vem de uma palavra francesa antiga randon (to run around, andar em volta) e “sample”, amostras é derivada do latin exemplum, exemplo, um pouco de cultura inútil para alegrar o dia.

Eu as vezes estava dando uma olhada no livro Statistical Mechanics. Algorithms and Computation, mais para olhar o primeiro capitulo, que tem uma introdução muito legal sobre MCMC.

Primeiro ele fala do exemplo que já vimos a alguma tempo atrás aqui.
Mas ele explica isso como um jogo de crianças.

pebbles

Então o joguinho tem um circulo desenhado dentro de um quadrado, e os guris vão jogando pedras, e vão anotando quantas pedras caem dentro do circulo. Ai basicamente ele traz o seguinte algorítimo.

mcalg.

Podemos implementar ele no R, para testar.

1
2
3
4
5
6
7
8
9
10
11
direct_pi<- function(n=4000) {
    n_acertos<-0
    for(i in 1:n){
        x<-runif(n=1,min=-1,max=1)
        y<-runif(n=1,min=-1,max=1)
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}

E fazer alguns experimentos.

> direct_pi() [1] 3152

Se usarmos nossa função, temos o número de acertos, relativos a 4 mil tentativas, 4 mil jogadas de pedras no chão. E se repetirmos vezes o suficiente, que podemos fazer com a função replicate do R, que replica coisas.

> repeticoes<-replicate(100,direct_pi()) > mean(4*(repeticoes/4000)) [1] 3.14076 > pi pi>min(4*(repeticoes/4000)) [1] TRUE

Vemos que fazendo isso 100 vezes, a média é muito próxima do valor de pi real, também que pi está contido dentro do menor e maior valores encontrados.

Massa, mas até aqui estamos falando de Monte Carlo, agora entra o Markov Chain, que é a tal da cadeia. Aqui, ele propõe outro jogo, pensando que os adultos foram jogar, só que eles foram jogar num heliporto. A primeira grande diferença aqui é que agora não da para jogar uma pedra, e cobrir toda a área, porque agora o espaço é muito grande. Mas o que da para fazer é ir andando pela área. Então a gente joga uma pedra, vai onde ela caiu, e joga outra pedra, de olhos fechados e assim vai, devagarinho, cobrindo toda a área do heliporto.

mcmcheli

Agora, o acerto é estar dentro do circulo, então todo mundo vai com um saquinho cheio de pedras, e toda lugar onde joga sua pedra, deixa uma pedrinha onde ela caiu para contar depois.
Mas existe um problema, se estamos jogando a pedra, indo onde ela está e jogando denovo, pode acontecer de eventualmente, chegarmos na borda do heliporto e jogar a pedra pra fora. E agora, a gente vai buscar a pedra? Continua jogando ela la fora até voltar pro heliporto? Bom, uma solução, é contar denovo, deixar 2 pedrinhas no mesmo lugar para contar depois, e jogar denovo do mesmo lugar, a sua pedra de arremessar.

Nisso a gente vai ter o seguinte

mcmcheli2

Agora a gente tem uma pilhas de pedras nas beiradas, e só pedras espalhadas pelo meio. Meio louco, mas esse joguinho funciona (alias, isso é o metropolis-hasting).
E temos um algorítimo para isso também.

markov_pi

E podemos implementar ele no R

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
markov_pi<-function(n=4000,passo=0.1){
    n_acertos<-0
    x<-1
    y<-1
    for(i in 1:n){
        delta_x<-runif(n=1,min=-passo,max=passo)
        delta_y<-runif(n=1,min=-passo,max=passo)
        if(abs(x+delta_x)<1 & abs(y+delta_y)<1){
            x<-x+delta_x
            y<-y+delta_y
        }
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}

Veja que agora a gente começa de um lugar (x=1,y=1), e vai andando no heliporto, andando na área, então, o antes, a gente jogava as pedrinhas, simplesmente sorteava números, agora a gente sorteia uma distância, para adicionar ao passo, que são os deltas.
Outra coisa, veja que se a gente olhar agora, antes de adicionar um acerto, se a pedra está dentro do heliporto, se estiver fora, a gente não soma os deltas, ou seja o local onde a pedra está continua igual, o que é adicionado aos acertos, ou não.

Podemos testar essa nova função, e vemos o seguinte.

> markov_pi() [1] 3397

Temos a contagem de passos.

> repeticoes<-replicate(100,markov_pi()) > mean(4*(repeticoes/4000)) [1] 3.11109

E o na média, se realizamos esse procedimento várias vezes, o valor é próximo de pi, sendo que o valor real está contido entre o maior e menor caso.

> pi pi>min(4*(repeticoes/4000)) [1] TRUE

Mas diferente do Monte Carlo, veja que aqui, nos não pegamos as amostras diretamente, como fizemos antes, aqui uma amostra depende da anterior, depende de onde você está no espaço amostral, e isso pode gerar uma certa dependência a curto prazo (que é que a gente tenta resolver quanto usa o thin no openbugs), mas ainda sim funciona.
Outra coisa é que veja que dependendo do tamanho do passo, a gente pode não explorar totalmente o espaço.

Se aumentarmos, ainda sim conseguimos andar bem em tudo

> repeticoes<-replicate(100,markov_pi(passo=0.5)) > mean(4*(repeticoes/4000)) [1] 3.13863

Mas se os passos forem muito pequenos, simplesmente podemos nunca chegar na borda, simplesmente demorar muito mais para percorrer todo o espaço.

> repeticoes<-replicate(100,markov_pi(passo=0.01)) > mean(4*(repeticoes/4000)) [1] 0.83486

Isso implica que para o MCMC, quando não retiramos nossas amostras diretamente da distribuição de interesse, temos que nos preocupar com esse ajuste do algorítimo, uma regra consensual, é que a rejeição, ou aceitação de novas amostras deve ficar em torno de 50%, o que é possível ser avaliado, é só contar o quanto estamos rejeitando amostras, para dar um parecer.

Bem é isso ai, so queria compartilhar essa explicação sobre MCMC, que achei muito legal, eu vi o cara do livro num curso do Coursera na verdade, aqui e depois eu consegui uma cópia do livro para dar uma olhada, o que é bem legal também. Alias o livro tem uma excelente página wiki, de acompanhamento, com os algorítimos em tudo que é linguagem, menos no R hehe, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:

Werner Krauth 2006 Statistical Mechanics – Algorithms and Computations

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
#Monte Carlo
direct_pi<- function(n=4000) {
    n_acertos<-0
    for(i in 1:n){
        x<-runif(n=1,min=-1,max=1)
        y<-runif(n=1,min=-1,max=1)
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}
 
#Testes
direct_pi()
 
repeticoes<-replicate(100,direct_pi())
mean(4*(repeticoes/4000))
pi<max(4*(repeticoes/4000))
pi>min(4*(repeticoes/4000))
 
#MCMC
markov_pi<-function(n=4000,passo=0.1){
    n_acertos<-0
    x<-1
    y<-1
    for(i in 1:n){
        delta_x<-runif(n=1,min=-passo,max=passo)
        delta_y<-runif(n=1,min=-passo,max=passo)
        if(abs(x+delta_x)<1 & abs(y+delta_y)<1){
            x<-x+delta_x
            y<-y+delta_y
        }
        if(x^2 + y^2 < 1){
            n_acertos<-n_acertos+1
        }
    }
    return(n_acertos)
}
 
 
#Teste
markov_pi()
 
repeticoes<-replicate(100,markov_pi())
mean(4*(repeticoes/4000))
pi<max(4*(repeticoes/4000))
pi>min(4*(repeticoes/4000))
 
 
repeticoes<-replicate(100,markov_pi(passo=0.5))
mean(4*(repeticoes/4000))
 
repeticoes<-replicate(100,markov_pi(passo=0.01))
mean(4*(repeticoes/4000))
pi<max(4*(repeticoes/4000))
pi>min(4*(repeticoes/4000))

Rosalind – Partial Permutations

Um exercício do Rosalind, para matar a saudade de fazer probleminhas.

Chama-se partial permutations.

Quando eu li esse nome, eu fui no google e fiz uma busca por partial permutations. E tem uma página do Wikipédia com esse título aqui. Veja que essa não é exatamente a solução. O que o problema realmente quer é o que está no fim dessa página, que é fazer os Restricted partial permutations.

O que a gente quer, por exemplo, é de um conjunto de 1 a 8, combinando de 3 em 3, quantas possibilidades temos, sem repetir os números.

Então alguns exemplos de conjuntos seriam:

(1,2,3) (1,2,4) (1,2,5) ... (7,8,9)

Então P(n,k) é quantos casos desses podemos escrever, onde n é as coisas que usamos para formar os conjuntos e k é quanto “espaço” temos. Isso serve para pensar por exemplo, quantos possibilidades temos de organizar os genes no genoma.

Bem, isso tem uma formula fechada, so consultar aqui

E assim perdeu a graça, em python, a função fatorial tem na biblioteca math, então:

1
2
3
4
import math
n=21
k=7
print (math.factorial(n)/math.factorial(n-k))%1000000

E temos

Python 2.7.10 (default, Oct 14 2015, 16:09:02) [GCC 5.2.1 20151010] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> >>> >>> >>> 51200

Bem é isso ai, o script vai estar la no repositório recologia, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail.

Referência:
Rosalind – Partial Permutations