Poisson Point Process

A gente ja falou aqui sobre a distribuição dos pontos nos espaço, mas existe mais coisa a se pensar.

Eu comecei a ler o novo livro do Marc Kery, Applied hierarchical modeling in Ecology, eu ja gostava muito das coisas que ele escrevia desde os dois primeiros livros dele, Introduction to WinBUGS for Ecologists e Bayesian population analysis using WinBUGS, bem se o livro do Ben Bolker Ecological Models and Data in R abria os olhos para um monte de coisa, o Marc Kery detalhava bem modelos bayesianos desde de um começo simples, e me parece que com esse novo livro, vamos um pouquinho mais para a frente, ao mesmo estilo.

Bem uma coisa que comum no livro dele, que fala principalmente de modelos populacionais relacionados a abundância e ocupação, é o processo de pontos de poisson, ou Poisson Point Process. A gente ja falou da distribuição de poisson, que aliás é comum a muitas disciplinas, não so a ecologia. Mas vamos pensar um pouco mais sobre ela.

O Poisson point process (me perdoem usar o termo em inglês, mas as vezes é mais fácil ja ver o termo em inglês, tanto porque é mais fácil para procurar no google, achar na literatura, como muitas vezes eu traduzo as coisas erradas com meu inglês nórdico e depois vira uma confusão) tem a propriedade que cada ponto é estocasticamente independente de todos os outros pontos no processo. Apesar disso, ele é amplamente usado para modelar fenomemos representados por pontos, como fazemos direto ao usar o glm com distribuição de poisson, mas isso implica que ele nao descreve adequadamente fenomenos que tem interação suficientemente forte entre os pontos, o que é bem comum em biologia, e que explica porque a gente ve coisas como modelos quasi-poisson no glm e porque nos preocupamos com overdispersion.

O processo tem seu nome do matemático Siméon Denis Poisson e vem do fato que se uma coleção de pontos aleatórios no espaço formam um processo de poisson, então o número de pontos em uma região finita é uma variável aleatória com uma distribuição de poisson. O processo depende de um único parâmetro, que, dependendo do contexto, pode ser constante ou uma função localmente integrável (na maioria dos nossos casos a+bx ou uma função com várias médias, como numa anova). Esse parâmetro representa o valor esperado da densidade de pontos em um processo de Poisson, também conhecido como intensidade e normalmente representado pela letra grega lambda (\lambda).

A função massa de probabilidade é dada por:

p(n)=\frac{\lambda^n e^{-\lambda}}{n!}

Bem, a integral de p(n) da 1, como toda função massa de probabilidade (isso parece obvio hoje, mas eu demorei so me liguei disso depois de uma disciplina no coursera).

Para um \lambda grande, teremos o equivalente a uma distribuição normal.

Lembrando, a distribuição normal é:

p(x)=\frac{1}{\sqrt{2\sigma^2\pi}}e^{\frac{(x-\mu)^2}{2\sigma^2}}

Tudo isso basta dar ?rnorm e ?rpois, que você vai ver que ambas as distribuições são implementadas assim no R. Mas agora olha so que mágico.

Se você olhar uma parada chamada aproximação de stirling, que é uma formula fechada que mais ou menos da o resultado de um número fatorial para valores grandes, nos temos que

x! \approx \sqrt{2\pi x}e^{-x}x^x \ quando \ x \rightarrow \infty

Então, poisson é uma distribuição discreta e a normal continua, mas seja x=n=\lambda(1+\delta) onde \lambda \geq 10 (essa relação so vale para números grandes) e \delta \leq 1 (poisson é discreto, normal é contínuo). Dai substituindo temos…

p(x)=\frac{\lambda^{\lambda(1+\delta)} e^{-\lambda}}{\sqrt{2\pi \lambda(1+\delta)}e^{-\lambda(1+\delta)}(\lambda(1+\delta))^{\lambda(1+\delta)}}

Primeiro trazemos a potência do e do divisor para cima

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda(1+\delta)}(\lambda(1+\delta))^{\lambda(1+\delta) - \lambda(1+\delta)}}

Depois descemos a potência do lambda la de cima

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda(1+\delta)}(\lambda(1+\delta))}

simplificamos

p(x)=\frac{ e^{\lambda \delta}}{\sqrt{2\pi \lambda}(\lambda(1+\delta)) (1+\delta)^{1/2}}

ai tiramos o delta de dentro da raiz

p(x)=\frac{ e^{\lambda \delta} \lambda(1+\delta)^{-1/2}}{\sqrt{2\pi \lambda}(1+\delta)}

subimos ele

p(x)=\frac{ e^{\lambda \delta} \lambda^{-1/2}}{\sqrt{2\pi \lambda}}

cortamos o um mais delta

p(x)=\frac{ \lambda e^{-(\lambda \delta)/2}}{\sqrt{2\pi \lambda}}

e arrumamos as potências

E agora tem uma parte que eu não entendo perfeitamente, mas lembre-se que a média e a variância da distribuição de poisson são iguais, então \lambda=\mu=\sigma^2.
Assim, a parte debaixo já está pronta, e a parte de cima da equação, já é praticamente a pdf da distribuição normal, só não sei exatamente o que fazer agora, mas como vemos a relação dela com a distribuição normal nesse diagrama do blog do John Cook, é algo por ae, se alguém quiser comentar o que fazer daqui para frente, eu ficaria feliz 🙂

distribution_chart

Tudo bem que eu não sou muito bom em matemática, mas se a gente fizer um gráfico das pdfs, assim:

1
2
3
plot(1:30,dpois(1:30,10),type="l",lty=3,lwd=3,col="red",frame=F,xlab="x",ylab="Probabilidade")
points(1:30,dnorm(1:30,10,sqrt(10)),type="l",lty=3,lwd=3,col="blue")
legend("topright",legend=c("Poisson","Normal"),lty=3,lwd=3,col=c("red","blue"),bty="n")

A gente vai ver o seguinte:

01

Bem é isso ai, sem script hoje, mas olhe o repositório recologia de qualquer forma, e se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail, que para esse post em especial, eu tenho certeza que devo ter falado algo errado!!! Então correções são bem vindas.

Referência:

John D. Cook – Diagram of distribution relationships
Marc Kéry & Andy Royle 2015 Applied hierarchical modeling in Ecology Elsevier 183pp

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

Data Augmentation para dados de Captura e Recaptura com efeito do tempo como uma variável aleatória

A algum tempo atrás, a gente falou de data augmentation aqui e depois usando data augmentation num modelo com efeito do tempo aqui.

Muito legal, mas particularmente no modelo com efeito do tempo, a gente cogitou (no livro do Marc Kery é cogitado na verdade) a possibilidade de tratar esse efeito como uma variável aleatória, então vamos tentar fazer isso e uma pequena modificação depois para ver o que acontece.

Bem os dados serão os mesmo que para o post de data augumentation com efeito do tempo. Mas vamos mudar o modelo. Na verdade vamos mudar um pouco a geração dos dados, vamos pegar o efeito de tempo de uma distribuição normal, para condizer com nossa analise.

Mas no mais é tudo igual, agora o modelo na linguagem Bugs

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
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()

Basicamente ele continua o mesmo, só mudamos a parte do prior la em cima.

1
2
3
4
5
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo

Veja que agora temos ainda a geração dos efeitos do tempo, de cada tempo de amostragem, mas agora pegamos esses valores de uma distribuição normal, e não de uma distribuição uniforme como da outra vez.

Mas tentamos determinar agora, qual a média e o desvio padrão dessa distribuição.
Para a média partimos de uma distribuição normal com média zero e desvio 0.001, ou seja, um pico no zero, e para o desvio padrão saímos de uma distribuição uniforme de 0 a 5, lembre-se que o desvio padrão tem que ser um número maior que zero, e não pode ser negativo, senão o bugs vai reclamar.

Uma outra coisa da linguagem bugs é que a distribuição normal usa precisão, e não o desvio diretamente na função. Isso não tem nenhum impacto na estatística, tem haver com precisão numérica apenas eu acho. Até hoje não entendi esse conceito perfeitamente, mas isso não muda nada.

1
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)

No fim, tau vai ser 1 dividido pela variância, que é o desvio padrão ao quadrado, veja que o desvio só vai aumentar, e um sobre um valor entre zero e alguma coisa, vai ser cada vez mais um valor pequenininho, e bem é isso.

Fazemos umas cadeias maiores, adicionamos na função de geradores de valores iniciais nossos novos priors

1
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),sigma.tempo = runif(1,0, 5))

E é isso, agora é so rodar o modelo.

> print(out_tempo, dig = 3) Inference for Bugs model at "modelo_tempo.txt", Current: 3 chains, each with 5000 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 13500 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 99.825 0.995 99.000 99.000 100.000 100.000 102.000 1.001 7100 p[1] 0.656 0.047 0.562 0.625 0.658 0.689 0.743 1.001 14000 p[2] 0.140 0.034 0.080 0.115 0.138 0.162 0.212 1.002 3000 p[3] 0.070 0.025 0.029 0.051 0.067 0.085 0.124 1.002 2400 p[4] 0.559 0.049 0.464 0.525 0.559 0.593 0.653 1.001 14000 p[5] 0.931 0.027 0.869 0.916 0.935 0.951 0.974 1.002 2000 med.tempo 0.470 0.306 -0.123 0.325 0.469 0.612 1.073 1.001 14000 sigma.tempo 0.579 0.381 0.236 0.359 0.471 0.663 1.608 1.001 14000 omega 0.402 0.031 0.341 0.381 0.402 0.423 0.464 1.001 5100 deviance 436.569 10.169 424.600 428.000 435.400 441.800 461.100 1.001 4200 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = var(deviance)/2) pD = 51.7 and DIC = 488.3 DIC is an estimate of expected predictive error (lower deviance is better). >

E apesar da estimativa do tamanho populacional estar ok, as estimativas do tempo ficam bem ruins, confira os valores que geraram os dados com os que foram estimados.

> variacao_Tempo [1] 1.517343 -1.058690 -1.568697 1.357698 3.792000

Veja que a estimativa do desvio padrão ficou bem abaixo do real usado, que era 2. Ela ficou algo entre zero e 1.6, enquanto o desvio real era 2. A média ficou correta quanto ao usado para gerar os dados, englobando o zero, mas o desvio deixou a desejar.

Se a gente olhar a figura de das estimativas do tamanho da população.

01

Veja que a gente perdeu bastante em variabilidade, já que o tempo variou menos, porque ficou “preso” dentro do nosso desvio padrão subestimado.

Mas ainda existe uma possibilidade. Lembre-se que aqui temos poucas medidas de de tempo, apenas 5. Lembra de alguma distribuição que tentamos tratar do número de amostras para controlar o desvio da média? Se pensou na distribuição t, pensou certo.

Podemos modelar a variável aleatória que gera o tempo como uma distribuição t, isso vai adicionar um parâmetro a mais no modelo, mas no entanto, teremos como ter a possibilidade de ter os valores mais longe da média como valores mais comuns para o valor do tempo. Ja que se você ler o post sobre a distribuição t, vera que isso que o novo parâmetro no t faz, o grau de liberdade, quanto mais graus de liberdade, mais parecido com uma distribuição normal, mas quanto menos, deixando a distribuição mais “uniforme”, com uma maior chance de valor mais longe da média, e como temos apenas 5 valores, isso pode ser uma propriedade interessante. É meio maluco pensar que vamos estimar esse parâmetro, já que na distribuição t ele é simplesmente o valor de amostras, mas ok, podemos estimar um valor para ele, para tentar ajustar um modelo mais realista.

De todo aquele modelo, vamos mudar apenar a parte dos prior, mais especificamente onde está a variável aleatória que é o tempo.

1
2
3
4
5
6
7
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
        k.tempo ~ dunif(2, 30) # Hyperprior da variacao do tempo
 
        for (i in 1:tempo){
            p[i] ~ dt(med.tempo, tau.tempo, k.tempo)
        }

Veja que usamos o dt ao invés de dnorm dentro do loop, e temos um parâmetro chamado de k.tempo agora, que vamos tirar de uma distribuição uniforme de 2 a 30, quando ele for 30, temos uma distribuição normal, quanto menor, deixamos tudo mais uniforme. Começamos no 2 porque o valor de k tem que ser no mínimo 2 para a linguagem bugs, veja a definição aqui.
E é isso, agora é só colocar ele no inicializador.

1
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),sigma.tempo = runif(1,0, 5),k.tempo= runif(1,2, 30))

E pronto, podemos rodar. Não esquecendo que temos que colocar ele entre os parâmetros que temos que acompanhar.

> print(out_tempo, dig = 3) Inference for Bugs model at "modelo_tempo.txt", Current: 3 chains, each with 5000 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 13500 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 99.866 1.010 99.000 99.000 100.000 100.000 102.000 1.001 11000 med.tempo 0.470 0.294 -0.100 0.322 0.469 0.613 1.065 1.001 14000 sigma.tempo 0.558 0.361 0.217 0.345 0.462 0.648 1.519 1.001 5900 k.tempo 16.565 7.884 2.979 9.950 16.660 23.360 29.310 1.001 5100 omega 0.402 0.031 0.343 0.380 0.402 0.423 0.465 1.001 14000 deviance 437.026 10.301 424.700 428.200 435.800 442.900 461.900 1.001 14000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = var(deviance)/2) pD = 53.1 and DIC = 490.1 DIC is an estimate of expected predictive error (lower deviance is better).

E olhando para a figura temos a mesma coisa

02

E no final das contas, o valor de k ficou entre 2 e 30, ou seja, não conseguimos estimar nada para ele, ja que de 2 a 30 era as possibilidades que abrimos. As estimativas de média para a distribuição de parâmetros continua ok, mas o desvio ainda está subestimado. Então apesar dessa ideia de usar a distribuição t seja legal, ela não contribuiu em nada praticamente, aumentamos um parâmetro e não melhoramos muito o modelo. Veja que o deviance continua a mesma coisa. Bom mas valeu a tentativa.

Eu imagino que dado o pequeno número de tempos, essa ideia de usar ele como uma variável aleatória não seja uma boa ideia, ja que estamos trocando 5 parâmetros por 3, ou 2 no primeiro exemplo. Além disso, precisamos lembrar que essa é uma afirmativa forte, pensar que que o tempo é uma variável aleatória, talvez num span pequeno de tempo, pode estar tudo ok, mas num span de tempo maior, dependendo da espécie que estamos trabalhando, podemos pensar o tanto que falamos em aquecimento global, em degradação do ambiente, em tantas mudanças que falar que nada acontece no tempo é algo bem “forte”. Mas apesar de tudo, talvez para mais medidas e um span pequeno de tempo, essa seja uma possibilidade bem legal.

Bem é isso ai, a gente tem que tentar fazer as coisas para aprender, 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:
Kéry e Shaub 2011 – Bayesian population analysis using WinBUGS. A hierarchical perspective – Academic Press 554pp

John K. Kruschke 2014 Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan.2nd Edition. Academic Press / Elsevier.

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
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
library(R2OpenBUGS)
 
##################################
## Gerando dados
##################################
set.seed(51)
n = 100
probabilidade_detecao = 0.25
tempo = 5
variacao_Tempo = rnorm(tempo,0,2)
 
yreal<- matrix(NA, nrow=n,ncol=tempo)
 
for (j in 1:tempo){
    p <- plogis(log(probabilidade_detecao / (1-probabilidade_detecao)) + variacao_Tempo[j])
    yreal[,j] <- rbinom(n = n, size = 1, prob = p)
}
 
detectado_pelo_menos_uma_vez <- apply(yreal, 1, max)
C <- sum(detectado_pelo_menos_uma_vez)
 
yobservado <- matrix(NA, nrow=n,ncol=tempo)
yobservado <- yreal[detectado_pelo_menos_uma_vez == 1,]
 
##################################
## Modelo com tempo
##################################
 
# Data Augumentation
nz <- 150
yaug <- rbind(yobservado, array(0, dim = c(nz, tempo)))
 
# Modelo na linguagem bugs
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dnorm(med.tempo, tau.tempo)
        }
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()
 
# Juntando os dados
bugs.data <- list(yaug = yaug, amostras = nrow(yaug), tempo = ncol(yaug))
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),
                  sigma.tempo = runif(1,0, 5))
 
# Parametros para saida
params <- c("N", "p","med.tempo","sigma.tempo","omega")
 
# Configurações do MCMC
ni <- 5000
nt <- 2
nb <- 500
nc <- 3
 
#Rodando o modelo
out_tempo <- bugs(bugs.data, inits, params, "modelo_tempo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out_tempo, dig = 3)
 
#Figura1
jpeg("01.jpg")
hist(out_tempo$sims.list$N,breaks=seq(90,110,by=1), col = "gray", main = "", xlab = "Tamanho Populacional", las = 1,xlim=c(95,105))
abline(v = C, col = "black", lwd = 3)
abline(v = n, col = "red", lwd = 3)
legend("topright",c("Detectado","Real"),lty=1,lwd=3,col=c("black","red"),bty="n")
dev.off()
#####################################
## Efeito aleatorio com distribuiçao t
#####################################
 
# Modelo na linguagem bugs
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
 
        med.tempo ~ dnorm(0,0.001) # Hyperprior da media tempo
        sigma.tempo ~ dunif(0, 5) # Hyperprior da variacao do tempo
        k.tempo ~ dunif(2, 30) # Hyperprior da variacao do tempo
 
        for (i in 1:tempo){
            p[i] ~ dt(med.tempo, tau.tempo, k.tempo)
        }
 
 
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        tau.tempo <- 1 / (sigma.tempo * sigma.tempo)
    }
    ",fill = TRUE)
sink()
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1),med.tempo = rnorm(1,0,0.001),
                         sigma.tempo = runif(1,0, 5),k.tempo= runif(1,2, 30))
inits()
# Parametros para saida
params <- c("N","med.tempo","sigma.tempo","k.tempo","omega")
 
#Rodando o modelo
out_tempo <- bugs(bugs.data, inits, params, "modelo_tempo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out_tempo, dig = 3)
 
#Figura2
jpeg("02.jpg")
hist(out_tempo$sims.list$N,breaks=seq(90,110,by=1), col = "gray", main = "", xlab = "Tamanho Populacional", las = 1,xlim=c(95,105))
abline(v = C, col = "black", lwd = 3)
abline(v = n, col = "red", lwd = 3)
legend("topright",c("Detectado","Real"),lty=1,lwd=3,col=c("black","red"),bty="n")
dev.off()

Data Augmentation para dados de Captura e Recaptura com efeito do tempo

Ano passado, a gente olhou aqui a ideia de data augumentation. Que é muito interessante. Ai eu queria fazer posts sobre todos esses modelos, que são muito legais, mas fiquei enrolando, mas agora sai mais um e sabe-se la quando o resto.

Mas vamos la, esse modelo assume que a probabilidade de detecção p varia por ocasião de coleta, se estivermos coletando passarinhos, talvez a chuva faça com que eles não saiam para voar e cair na rede, ou para pequenos mamíferos, as condições do clima estejam ruins, o dia para ficar na cama abaixo das cobertas, ou ainda diferentes armadilhas ou mecanismos de detecção foram usados.

Então vamos gerar dados, pensando em passarinhos caindo numa rede, que colocamos no campo, e conseguimos ir pegando os passarinhos e anilhando eles, conforme caem na rede, então identificamos eles perfeitamente.

Então imagine:

1
2
3
4
n = 100
probabilidade_detecao = 0.25
tempo = 5
variacao_Tempo = runif(tempo,-2,2)

Cem passarinhos (n), coletados em 5 campanha diferentes (tempo), mas além desses parâmetros, sabemos que… , bem nos não sabemos, mas a mãe natureza estabelece uma probabilidade desse passarinho cair na rede (probabilidade_detecao), mas essa probabilidade não é constante ao longo das campanhas, as vezes ela é mais um pouquinho, quando o clima ta bom, sem vento, e as vezes menos um pouquinho (esses valores estão no vetor variacao_Tempo).

Ok, agora vamos fazer a matriz de dados que a mãe natureza vê, ela vê todos dos 100 passarinhos que la residem, e ela observa as 5 campanhas também, ou seja uma matriz de n pelo tempo.

1
yreal<- matrix(NA, nrow=n,ncol=tempo)

Para simular os dados fazemos o seguinte

1
2
3
4
5
6
7
#Num loop para todas as campanha
for (j in 1:tempo){
    #Qual a probabilidade de detecção nessa campanha j?
    p <- plogis(log(probabilidade_detecao / (1-probabilidade_detecao)) + variacao_Tempo[j])
    #Ai simulamos todos os 100 individuos de uma vez
    yreal[,j] <- rbinom(n = n, size = 1, prob = p)
}

E pronto, a simulação é isso, gerar números a partir de uma distribuição binomial, na verdade.
Só que quem a gente nunca viu em nenhuma campanha, quem nunca foi anilhado não existe pra gente, só foi anilhado quem foi pego pelo menos uma vez

1
2
detectado_pelo_menos_uma_vez <- apply(yreal, 1, max)
C <- sum(detectado_pelo_menos_uma_vez)

E assim esses são os cara que observamos

1
2
yobservado <- matrix(NA, nrow=n,ncol=tempo)
yobservado <- yreal[detectado_pelo_menos_uma_vez == 1,]

Agora basta a gente estabelecer nosso modelo para a linguagem bugs

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dunif(0, 1)
        }
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        p.medio <- mean(p[])
    }

E a principio o likelihood pode parecer idêntico ao do post anterior, mas veja que a diferença está em usar um vetor para as probabilidade de detecção da campanha

1
p.eff[i,j] <- z[i] * p[j]

Veja que no prior, inicializamos um vetor agora, e não um único valor, para contabilizar essa variabilidade entre campanha

1
2
3
for (i in 1:tempo){
     p[i] ~ dunif(0, 1)
}

Ai agora é só estimar os parâmetros, temos que colocar os dados no formado do openbugs que estamos usando aqui, lembrando das outras opções aqui, estabelecer valores iniciais para as cadeia do MCMC, seu parâmetros como tamanho, valores iniciais a serem queimados, número de cadeia e rodar.

Aqui o resultado fica o seguinte:

> print(out_tempo, dig = 3) Inference for Bugs model at "modelo_tempo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 99.604 4.646 92.000 96.000 99.000 102.000 110.000 1.001 6000 p[1] 0.542 0.055 0.436 0.505 0.542 0.579 0.649 1.001 6000 p[2] 0.108 0.031 0.055 0.086 0.105 0.127 0.179 1.001 6000 p[3] 0.138 0.034 0.078 0.113 0.135 0.160 0.211 1.001 6000 p[4] 0.672 0.056 0.560 0.635 0.673 0.711 0.776 1.001 6000 p[5] 0.099 0.030 0.049 0.078 0.096 0.117 0.164 1.001 6000 p.medio 0.312 0.022 0.268 0.297 0.312 0.327 0.354 1.001 6000 omega 0.418 0.037 0.349 0.392 0.417 0.443 0.495 1.002 2700 deviance 470.030 20.741 434.300 455.500 468.300 483.100 514.302 1.001 6000 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = var(deviance)/2) pD = 215.1 and DIC = 685.2 DIC is an estimate of expected predictive error (lower deviance is better).

Bem existem mais ou menos entre 92 a 110 indivíduos por ali, o intervalo de confiança, o que cobre a verdade (já que sabemos os valores usados pela mãe natureza foi 100), e temos um valor de detecção por campanha e duas quantidades derivadas, o p médio e o omega, que falamos depois.

Podemos tentar olhar esse resultado numa figura das nossas estimativas.

01

Beleza, mas será que valeu a pena o esforço, o que aconteceria se a gente ajusta-se o modelo considerando um único valor de detecção. Bem, podemos ajustar e ver o que acontece, por curiosidade.

E a saída é o seguinte

> print(out, dig = 3) Inference for Bugs model at "modelo.txt", Current: 3 chains, each with 2500 iterations (first 500 discarded), n.thin = 2 Cumulative: n.sims = 6000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 112.268 7.386 100.000 107.000 112.000 117.000 129.000 1.002 1700 p 0.274 0.026 0.224 0.256 0.274 0.292 0.326 1.002 1700 omega 0.471 0.044 0.390 0.440 0.469 0.499 0.561 1.002 1400 deviance 657.828 23.208 616.100 641.300 656.900 672.425 707.800 1.002 1700 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = var(deviance)/2) pD = 269.1 and DIC = 926.9 DIC is an estimate of expected predictive error (lower deviance is better).

Aqui o intervalo de confiança ainda pegou o valor real, que era de 100, o tamanho populacional teve uma estimativa de 100 a 129 indivíduos, o que está ok, com um p de 274, a probabilidade média de captura. Agora o que isso quer dizer?
Vamos olhar os valores da mãe natureza.

Veja que os valores de detecção foram

> plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo) [1] 0.50132242 0.09139408 0.12948952 0.69882099 0.09683113

Da média original, na primeira campanha foi bem mais, na segunda e terceira bem menos, na quarta a mais, e na quinta a menos. se fizermos a média disso da

> mean(plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo)) [1] 0.3035716

O valor que aquele p realmente queria chegar era 0.304, mas ele ficou em 0.274, isso quer dizer, que na segunda, terceira e quinta campanha, a gente deveria ter pensado em achar menos passarinhos, mas o modelo fez simulação com valores altos, sempre achando mais, por isso a gente super-estima muito, veja que comparando as duas distribuições, nesse exemplo o modelo sem levar o tempo foi la para frente.

figura2

Veja que isso é a quantidade derivada la do modelo com o tempo, que ficou próximo desse valor real.

p.medio 0.312 0.022 0.268 0.297 0.312 0.327 0.354 1.001 6000

Um pouquinho maior, mas próximo.

Veja que isso vai depender claro do número de campanhas, pensando que a não ser que haja uma tendência nesses valores de detecção, por exemplo, a não ser que eles sempre melhorem, ou piorem, se eles simplesmente forem um processo estocástico com média zero, com muitas campanhas, os modelo sem o tempo conseguiria um valor bom de estimativa populacional, mas basta pensar um pouco para saber que isso não vai acontecer, porque?

Porque campanhas de coletas são caras, dificilmente o dinheiro vai dar para 30 ou 90 campanhas, e com poucas campanhas, a chance de você considerar um valor fixo de detecção e dar errado é maior, chance no sentido de que o modelo sempre vai ser pior, então talvez começar com um modelo assim seja mais seguro, e é fácil você simplificar se for possível.

Veja o deviance do modelo sem o tempo é muito pior:

deviance 657.828 23.208 616.100 641.300 656.900 672.425 707.800 1.002 1700

Compare como o deviance do modelo com o tempo é muito melhor

deviance 470.030 20.741 434.300 455.500 468.300 483.100 514.302 1.001 6000

Veja o quanto, quando comparamos as estimativas, o desvio padrão das estimativas para o modelo com tempo é menor.

figura3

A importante contribuição da incerteza sobre o tamanho da população vem da incerteza sobre a probabilidade p. Só que a gente paga em parâmetros a mais que precisamos para estimar p. Mas ignorar esses parâmetros levam a erros, que podem ser muito ruins. Talvez se a variação fosse pequena, esse modelo nem vale-se a pena, mas na dúvida, compensa começar do modelo mais complexo, e ir simplificando, vendo se é possível simplificar, afinal, se todas as probabilidades de detecção fossem iguais, daria para ver nos parâmetros estimados e não precisaríamos de um vetor de probabilidade de detecção, mas pero que si, pero que no, compensa tentar todas as possibilidades.

Nesse modelo, nos assumimos que os efeitos do tempo são fixos, mas poderíamos ainda ajustar um modelo estocástico também, pensando que ao invés de valores fixos de p, eles viriam de uma distribuição como a normal ou uniforme, que é o caso dos dados que geramos aqui.

Bem é isso ai, o script vai estar la no repositório recologia, falta uns dois modelos a mais sobre esse tema no livro muito legal do Kery, mas um dia chegamos la, isso tem no capítulo 6. Se eu escrevi alguma bobeira, algo errado, deixe um comentário corrigindo ou mande um e-mail e até mais.

Referência:
Kéry e Shaub 2011 – Bayesian population analysis using WinBUGS. A hierarchical perspective – Academic Press 554pp

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
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
#
library(R2OpenBUGS)
 
##################################
## Gerando dados
##################################
set.seed(51)
n = 100
probabilidade_detecao = 0.25
tempo = 5
variacao_Tempo = runif(tempo,-2,2)
 
yreal<- matrix(NA, nrow=n,ncol=tempo)
 
for (j in 1:tempo){
    p <- plogis(log(probabilidade_detecao / (1-probabilidade_detecao)) + variacao_Tempo[j])
    yreal[,j] <- rbinom(n = n, size = 1, prob = p)
}
 
detectado_pelo_menos_uma_vez <- apply(yreal, 1, max)
C <- sum(detectado_pelo_menos_uma_vez)
 
yobservado <- matrix(NA, nrow=n,ncol=tempo)
yobservado <- yreal[detectado_pelo_menos_uma_vez == 1,]
 
##################################
## Modelo com tempo
##################################
 
# Data Augumentation
nz <- 150
yaug <- rbind(yobservado, array(0, dim = c(nz, tempo)))
 
# Modelo na linguagem bugs
sink("modelo_tempo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        for (i in 1:tempo){
            p[i] ~ dunif(0, 1)
        }
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p[j]
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
        p.medio <- mean(p[])
    }
    ",fill = TRUE)
sink()
 
# Juntando os dados
bugs.data <- list(yaug = yaug, amostras = nrow(yaug), tempo = ncol(yaug))
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(tempo, 0, 1))
 
# Parametros para saida
params <- c("N", "p","p.medio", "omega")
 
# Configurações do MCMC
ni <- 2500
nt <- 2
nb <- 500
nc <- 3
 
#Rodando o modelo
out_tempo <- bugs(bugs.data, inits, params, "modelo_tempo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out_tempo, dig = 3)
 
#Figura1
hist(out_tempo$sims.list$N,nclass=30, col = "gray", main = "", xlab = "Tamanho Populacional", las = 1, xlim = c(80, 120))
abline(v = C, col = "black", lwd = 3)
abline(v = n, col = "red", lwd = 3)
legend("topright",c("Detectado","Real"),lty=1,lwd=3,col=c("black","red"),bty="n")
 
##################################
## Modelo sem variação de tempo
##################################
 
# Modelo na linguagem bugs
sink("modelo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        p ~ dunif(0, 1)
        # Likelihood
        for (i in 1:amostras){
            z[i] ~ dbern(omega)
            for (j in 1:tempo){
                yaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p
                } #j
            } #i
        # Quantidade derivadas
        N <- sum(z[])
    }
    ",fill = TRUE)
sink()
 
# Valores iniciais
inits <- function() list(z = rep(1, nrow(yaug)), p = runif(1, 0, 1))
 
#
params <- c("N", "p", "omega")
 
#Rodando o modelo
out <- bugs(bugs.data, inits, params, "modelo.txt", n.chains = nc, n.thin = nt, n.iter = ni, n.burnin = nb)
 
# Resultados
print(out, dig = 3)
print(out_tempo, dig = 3)
 
#Figura 2
hist(out_tempo$sims.list$N, nclass = 40, col = "darkgray", main = "", xlab = "Tamanho Populacional", las = 1, xlim = c(70, 150))
hist(out$sims.list$N, nclass = 40, col = "lightgray",add=T)
abline(v = mean(out_tempo$sims.list$N), col = "blue", lwd = 3)
abline(v = mean(out$sims.list$N), col = "red", lwd = 3)
legend("topright",c("Modelo com tempo","Modelo sem tempo"),lty=1,lwd=3,col=c("blue","red"),bty="n")
 
#
plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo)
mean(plogis(log(probabilidade_detecao/(1-probabilidade_detecao))+variacao_Tempo))
 
#Figura 3
par(mfrow=c(2,1))
hist(out_tempo$sims.list$N, nclass = 40, col = "darkgray", main = "Modelo com variação no tempo",
     xlab = "Tamanho Populacional", las = 1, xlim = c(70, 150),ylim=c(0,600))
text(130,400,paste("Desvio padrão =",round(sd(out_tempo$sims.list$N),3)))
hist(out$sims.list$N, nclass = 40, col = "darkgray", main = "Modelo sem variação no tempo",
     xlab = "Tamanho Populacional", las = 1, xlim = c(70, 150),ylim=c(0,600))
text(130,400,paste("Desvio padrão =",round(sd(out$sims.list$N),3)))

Regressão linear ajustada com lm, bbmle, OpenBugs, Jags e Stan no R

No último post a gente falou do pacote sads do R. Muito legal, e ele faz vasto uso do pacote bbmle.

O pacote bbmle é bem legal, ele permite você escrever uma função de verosimilhança qualquer e então ele ajusta ela pra você. A gente tem um monte de post, por exemplo aqui, aqui e aqui onde o que a gente faz no fim das contas é escrever a função de verosimilhança e ajustar. Mas a gente sempre tinha ajustado usando estatística bayesiana, usando o OpenBugs pra ser mais preciso, e para conferir usando estatística frequentista usando o lme4, ou algum outro pacote que fazia o modelo pra gente.

Mas o bblme vai funcionar mais ou menos na mesma estrategia da linguagem bugs, escreve seu modelo e o bbmle ajusta pra você.

Vamos pegar um modelo simples aqui, a boa e velha regressão linear, e criar alguns dados.

01

Beleza, o mais simples de fazer é usar a função lm que ajusta modelos lineares no R a ajustar o modelo e olhar o resultado e pronto, simples assim.

> fit.lm <- lm(resposta~preditor) > summary(fit.lm) Call: lm(formula = resposta ~ preditor) Residuals: Min 1Q Median 3Q Max -1.0381 -0.4135 0.1270 0.3737 0.7836 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.41354 0.32834 7.351 5.28e-08 *** preditor 0.47595 0.01535 31.004 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.491 on 28 degrees of freedom Multiple R-squared: 0.9717, Adjusted R-squared: 0.9707 F-statistic: 961.3 on 1 and 28 DF, p-value: < 2.2e-16

Certo, agora vamos lembrar como é a formula da reta da regressão.
O que a gente fez é:

resposta = \alpha + \beta \cdot preditor + erro

E esse erro é a parte onde entra a distribuição normal, gaussiana ou sei la, mas o erro então é o seguinte.

erro = Normal(0,\sigma^2)

A distribuição com média zero e desvio sigma.

Daria na mesma o que a gente escreveu acima estar tudo junto:

resposta = \alpha + \beta \cdot preditor + Normal(0,\sigma^2)

Ou ainda:

resposta = Normal(\alpha + \beta \cdot preditor,\sigma^2)

Ou as vezes ao invés daquele escrito normal acho que a gente vê escrito p ou P, mas as pessoas usam p pra tanta coisa em estatística que não é a toa que tudo fica confuso, e não sou só eu que fica confuso, veja so aqui.

Mas a gente consegue perceber, pelo gráfico, que a resposta e o preditor são nossos dados, que temos valores, e o que sobra são os chamados parâmetros nessa formula toda. Na regressão, a gente está lidando com três deles, dois que determinam a reta, o alpha (\alpha) que é onde ela cruza o eixo de resposta e o beta (\beta) que é a inclinação. Além desses dois caras temos o desvio padrão, que é o sigma (\sigma), esse cara diz o quanto afastado daquela reta as nossas amostras podem estar, não exatamente isso mas mais ou menos isso, a dispersão dos pontos em relação a reta.

Ok, por mais confuso que isso possa parecer, se esse é o modelo, e resposta e preditor nos temos, sempre esses três caras tem que sair do modelo, que no caso a gente otimiza para conseguir o melhor ajuste.

Então o bbmle, assim como o OpenBugs e companhia vai precisar de um modelo construído, que diga pra ele como esta a qualidade do ajuste, verosimilhança do modelo, e ele te diz quais os melhores valores para esses três parâmetros.

Na saida da regressão ali em cima, o alpha é o intercept (2.41354), o beta é o que ele chamou de preditor(0.47595) e por último, o sigma é o Residual standard error (0.491). Os parâmetros alpha e beta não se tem como ter certeza absoluta, so quem tem certeza é quem gerou os dados, como na natureza é ela que gera os dados, então a gente nunca tem certeza, mas a gente acha que esses são os melhores valores (lembrando que cada valor tem seu erro aceitável também, ali do lado).

Mas então, o lm faz todo esse processo pra gente, e não tem porque não usar ele para modelos assim, mas só para treinar, ou curiosidade, vamos ver como ficaria esse modelo no bbmle.

No bbmle a gente tem que fazer nossa própria função de verosimilhança.

funcao.likelihood <- function(alpha=1,beta=1,desvio=1) {
    -sum(dnorm(resposta,mean=alpha+beta*preditor,sd=desvio,log=TRUE))
}

Antes de pensar que essa função faz muita coisa, vamos só usar ela, pra ver como é mais simples que parece, ela simplesmente calcula um número, que é a verosimilhança do modelo, na verdade o negative log-likelihood, veja que tem um menos antes do sum ali.

Então você fala um valor pra cada parâmetro e a função fala se o quão bom é o ajuste, em termos de verossimilhança.

> funcao.likelihood(alpha=1,beta=1,desvio=1) [1] 1488.255

Se a gente muda esses valores, o ajuste muda

> funcao.likelihood(alpha=1.5,beta=1,desvio=1) [1] 1632.549

Quanto menor melhor, quanto mais a gente aproxima todos os valores dos valores que usamos para gerar os dados, menor ele vai ficar.

> funcao.likelihood(alpha=1.5,beta=0.9,desvio=1) [1] 1038.321

E se a gente usar parâmetros muito nada a ver, esse número fica gigante.

> funcao.likelihood(alpha=8,beta=0.9,desvio=1) [1] 3195.411

E isso é tudo que o bbmle precisa. Agora para ajustar os parâmetros é so usar essa função de verosimilhança na função mle do pacote bbmle.

> fit.mle <- mle2(funcao.likelihood) Houve 11 avisos (use warnings() para vê-los) > summary(fit.mle) Maximum likelihood estimation Call: mle2(minuslogl = funcao.likelihood) Coefficients: Estimate Std. Error z value Pr(z) alpha 2.413540 0.317207 7.6087 2.768e-14 *** beta 0.475945 0.014830 32.0925 < 2.2e-16 *** desvio 0.474314 0.061234 7.7459 9.487e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -2 log L: 40.38304

E pimba, estimamos tudo do modelo. O desvio, ou sigma aqui fica menor que na função lm, porque o método usado é diferente, mas basicamente temos o mesmo resultado. Recebemos um aviso também, que acho que porque pode ser tentando valores negativos de desvio na função, mas não existe desvio negativo, dai a função produz um NaN, tipo assim.

> funcao.likelihood(alpha=8,beta=0.9,desvio=-1) [1] NaN Mensagens de aviso perdidas: In dnorm(resposta, mean = alpha + beta * preditor, sd = desvio, : NaNs produzidos

Mas o mle acha o caminho certo e ajusta. Ele usa o método do Nelder and Mead como default, a gente falou do método do Newton aqui, e vimos o gradient descent aqui, esse é outro método que serve para achar os melhores valores para os parâmetros.

Mas o bbmle é legal que da para usar outros métodos de otimização, só citar qual método você quer usar.

> fit.mle2<-mle2(funcao.likelihood,method="L-BFGS-B",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001)) > summary(fit.mle2) Maximum likelihood estimation Call: mle2(minuslogl = funcao.likelihood, method = "L-BFGS-B", lower = c(alpha = 1e-04, beta = 1e-04, desvio = 1e-04)) Coefficients: Estimate Std. Error z value Pr(z) alpha 2.413540 0.317207 7.6087 2.768e-14 *** beta 0.475945 0.014830 32.0924 < 2.2e-16 *** desvio 0.474314 0.061234 7.7459 9.488e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -2 log L: 40.38304

Ou ainda usar outra função para otimizar, que não o optmin que é o otimizador default do R, eu acho.

> fit.mle3<-mle2(funcao.likelihood,optimizer="nlminb",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001)) > summary(fit.mle3) Maximum likelihood estimation Call: mle2(minuslogl = funcao.likelihood, optimizer = "nlminb", lower = c(alpha = 1e-04, beta = 1e-04, desvio = 1e-04)) Coefficients: Estimate Std. Error z value Pr(z) alpha 2.413539 0.317206 7.6087 2.768e-14 *** beta 0.475945 0.014830 32.0925 < 2.2e-16 *** desvio 0.474313 0.061234 7.7460 9.486e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -2 log L: 40.38304

Veja que esses outros métodos nos permitiram fazer mais coisas, uma delas foi delimitar limites para os parâmetros. No caso colocamos que tudo tem que ser maior que 0.0001, isso impossibilita um valor negativo pro desvio e veja que depois disso a gente não recebe mais warnings.
Eu acho que a parte boa do bbmle, é que se amanha a gente lê um artigo, com um modelo muito louco que alguém acabou de inventar, se a gente conseguir passar ele para uma função de verossimilhança, a gente consegue ajustar valores pra ele.

Alias o cara que escreveu esse pacote, o Ben Bolker, que é um cara que parece ser, além de inteligente, bem legal, já que ele responde todo mundo nos fóruns de R, e em tudo que é fórum da internet, bem ele escreveu um livro chamado Ecological Models and Data in R, que é muito legal, pois ao invés de ensinar sobre testes e mais testes com seus nomes complicados, ele ensina que existe funções determinísticas, distribuições estocásticas,e você juntas essas coisas e descreve o mundo, ai depois de muitos capítulos ele mostra que você continua fazendo a mesma coisa que fazia nas anovas da vida, mas usando o modelo que você pensa e não pensando o que pode testar com os dados que tem.

Mas já que estamos aqui, veja que na estatística bayesiana, apesar de partir de pressupostos diferentes, paradigmas diferente e tal, no final a gente quer descrever a mesma coisa, e vai ser a mesma coisa na pratica, escrever um modelo, ajustar os dados, essa mesma regressão é feita da seguinte forma no OpenBugs

Escrevemos o modelo

sink("linreg.txt")
cat("
model {
# Prior
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 100)
 tau <- 1/ (sigma * sigma)
 
# Verossimilhança
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha + beta*x[i]
 }
}
",fill=TRUE)
sink()

Depois temos que preparar o terreno, preparar os dados para exportar, criar uma função para inicializar a cadeia, definir os parâmetros da cadeia e então otimizar.

> bugs.data <- list(x=preditor,y=resposta,n=length(resposta)) > inits <- function(){ list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))} > params <- c("alpha","beta","sigma") > nc = 3 #Numero de cadeias > ni=1200 #Tamanho da cadeira > nb=200 #Numero de simulação que serão descartadas > nt=1 #Thinning rate > modelo1.bugs <-bugs(data = bugs.data, inits = inits, parameters = params,model = "linreg.txt",n.thin = nt,n.chains = nc,n.burnin = nb, n.iter = ni) > print(modelo1.bugs, dig = 3) Inference for Bugs model at "linreg.txt", Current: 3 chains, each with 1200 iterations (first 200 discarded) Cumulative: n.sims = 3000 iterations saved mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff alpha 2.460 0.354 1.812 2.223 2.438 2.680 3.230 1.013 170 beta 0.474 0.017 0.436 0.464 0.475 0.485 0.504 1.017 140 sigma 0.514 0.071 0.398 0.462 0.506 0.559 0.674 1.002 1600 deviance 43.640 2.655 40.630 41.700 42.960 44.780 50.750 1.003 1600 For each parameter, n.eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor (at convergence, Rhat=1). DIC info (using the rule, pD = Dbar-Dhat) pD = 2.9 and DIC = 46.5 DIC is an estimate of expected predictive error (lower deviance is better).

Aqui o sigma fica um pouquinho maior, mas ainda temos basicamente o mesmo resultado. O bbmle é um pacote do R, aqui é meio diferente, porque o OpenBugs é um programa, de vida própria, que do R a gente um pacote chamado R R2OpenBUGS que tudo que faz é mandar todas essas informações que separamos no R para o OpenBugs, o OpenBugs processa elas e depois devolve bonitinho pro R, tudo organizadinho para então a gente interpretar, e fazer o que mais quiser como gráficos.

O OpenBugs, que na verdade vem da linguagem Bugs, de fazer inferência bayesiana, no caso Bugs significa Bayesian inference Using Gibbs Sampling, que usa um tipo de amostrador para tentar descobrir os valores reais dos parâmetros, e chegar ao melhor chute pelo menos. No final a gente queria saber o valor real do parâmetro, mas nunca saberemos, tem um post sobre Gibbs Sampler aqui, se a curiosidade bater.

Mas o OpenBugs não é o único a usar a linguagem Bugs para inferência, o jeitão de escrever o modelo etc, temos o WinBugs, que não esta mais em desenvolvimento, e o Jags, Just another Gibbs Sampler, que está a todo vapor eu acho. Na verdade ajustar modelos com o Jags, que é um programa externo também, é bem parecido, a gente pode ajustar o mesmo modelo denovo, da seguinte forma.

> jags <- jags.model(file='linreg.txt',data=bugs.data,inits=inits,n.chains = 3,n.adapt=100) Compiling model graph Resolving undeclared variables Allocating nodes Graph Size: 130 Initializing model |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% > update(jags, 1000) |**************************************************| 100% > jags.samples(jags,c('alpha', 'beta','sigma'),1000) |**************************************************| 100% $alpha mcarray: [1] 2.380524 Marginalizing over: iteration(1000),chain(3) $beta mcarray: [1] 0.4774447 Marginalizing over: iteration(1000),chain(3) $sigma mcarray: [1] 0.5097927 Marginalizing over: iteration(1000),chain(3)

E temos nossos três parâmetros denovo, veja que usamos tudo igual ao que usamos no OpenBugs. O modelo, o jeito de mandar os dados, numa lista e tal. A gente ouve muito falar do Jags, eu acho que uso mais o OpenBugs porque os livros que olho sempre o pessoal ta usando ele, mas migrar para os modelos pro Jags é relativamente fácil, como acabamos de ver.

Agora uma outra nova opção que surgiu foi o Stan. Stan é um programa para ajustar modelos bayesianos, e mais novo, acho que tem mais suporte a computação paralela, mas ao contrario do OpenBugs e do Jags que usam o Gibbs Sampler, ele usa um negocio chamado Hamiltonian Monte Carlo, esse eu não tenho nem ideia de como funciona, mas a gente pode ajustar a mesma regressão com ele, o que não é tão difícil, e ver que funciona igual.

O jeito de construir os modelos nele é bem diferente, essa nossa regressão fica algo assim.

linreg_stan <- '
data{
  int n;
  real y[n];
  real x[n];
}
 
parameters{
  real alpha;
  real beta;
  real<lower=0, upper=100> sigma;
}
 
transformed parameters{
    real mu[n];
    for(i in 1:n){
    mu[i] <- alpha + beta*x[i];
  }
}
 
model{
  //Priors
  alpha ~ normal(0, 100);
  beta ~ normal(0, 100);
  sigma ~ uniform(0, 100);
 
  //Likelihood
  for(i in 1:n){
    y[i] ~ normal(mu[i], sigma);
}
}
'

O que no final das contas nem é tão diferente assim do Bugs, é definir a função de verosimilhança, como é o framework bayesiano, a gente tem os prior, só que além disso la no inicio, a gente tem que inicializar todos os valores que vamos usar, como no c++, linguagem que ele é feito. Dai é so preparar os dados para mandar, numa lista igual ao bugs e mandar ajustar.

> stan.data <- list(n = length(resposta), x = preditor, y = resposta) > fit.stan<- stan(model_code=linreg_stan,data=stan.data,iter = 1000,chains = 3,thin = 1) TRANSLATING MODEL 'linreg_stan' FROM Stan CODE TO C++ CODE NOW. COMPILING THE C++ CODE FOR MODEL 'linreg_stan' NOW. SAMPLING FOR MODEL 'linreg_stan' NOW (CHAIN 1). Iteration: 1 / 1000 [ 0%] (Warmup) Iteration: 100 / 1000 [ 10%] (Warmup) Iteration: 200 / 1000 [ 20%] (Warmup) . . . Um monte de mensagens . . . Iteration: 900 / 1000 [ 90%] (Sampling) Iteration: 1000 / 1000 [100%] (Sampling) # Elapsed Time: 0.166985 seconds (Warm-up) # 0.131985 seconds (Sampling) # 0.29897 seconds (Total)

E podemos olhar nossos resultado de forma parecida ao output do OpenBugs

> print(fit.stan,digits=3) Inference for Stan model: linreg_stan. 3 chains, each with iter=1000; warmup=500; thin=1; post-warmup draws per chain=500, total post-warmup draws=1500. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat alpha 2.420 0.020 0.369 1.668 2.184 2.419 2.656 3.127 324 1.005 beta 0.476 0.001 0.017 0.442 0.465 0.475 0.487 0.511 331 1.003 sigma 0.520 0.004 0.077 0.395 0.464 0.513 0.567 0.693 454 1.003 lp__ 4.950 0.078 1.449 1.411 4.339 5.278 5.954 6.517 347 1.002 Samples were drawn using NUTS(diag_e) at Thu Jul 17 23:01:30 2014. For each parameter, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat=1).

E denovo os mesmo valores de alpha, beta e sigma.
Bem a gente começou olhando o bbmle, usando o lm pra ver se estávamos chegando a valores razoáveis para os parâmetros da regressão e acabamos olhando também como ajusta usando o OpenBugs, Jags e Stan.

A primeira pergunta agora seria, qual caminho seguir?

A minha melhor resposta é:
Eu não faço a menor ideia, mas é legal pra caralho chegar quase no mesmo valor de alpha, beta e sigma de todos os jeitos.
Na pior da hipóteses, se o modelo que estivermos tentando ajustar for certo, vai dar certo. Claro que para uma regressão linear a gente vai usar o lm na verdade, acho que não faz sentido adicionar toda a complexidade desses outros métodos para algo simples. Mas para muitos outros modelos de nosso interesse, que não da para ajustar com lm, ou mesmo glm, o jeito é começar a embrenhar nesses métodos.

Claro que se da para encontrar as coisas prontas, como os modelos prontinhos no pacote sads, vamos usar la e não tentar reinventar a roda. Mas as vezes a gente quer mudar algo um pouquinho, e ai saber como tudo funciona é legal. Além disso, como no livro do Bolker, um comentário que acho muito interessante, não interessa o método que você usa, é sempre bom saber um pouquinho de cada método para conseguir ler melhor a literatura, já que sempre cada pessoa vai fazer as coisas de um jeito, e não vamos deixar de ler um artigo porque ele usa estatística bayesiana ou frequentista.

Bem e esse post é legal, porque eu sempre quis ver alguém ajustando algo que eu entenda, como uma regressão linear com todos esses programas e pacotes, como ninguém fez isso pra mim ler, eu mesmo fiz ^^. Agora eu tenho uma referencia quando quiser usar algo diferente do OpenBugs.

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

set.seed(171)
preditor<-runif(30,10,30)
resposta<-rnorm(30,2+0.5*preditor,0.5)
 
plot(resposta~preditor)
 
############################
#lm do stats
############################
fit.lm <- lm(resposta~preditor)
summary(fit.lm)
 
############################
#bbmle
############################
library(bbmle)
 
funcao.likelihood <- function(alpha=1,beta=1,desvio=1) {
    -sum(dnorm(resposta,mean=alpha+beta*preditor,sd=desvio,log=TRUE))
}
 
fit.mle <- mle2(funcao.likelihood)
summary(fit.mle)
warnings()
 
fit.mle2<-mle2(funcao.likelihood,method="L-BFGS-B",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001))
summary(fit.mle2)
 
fit.mle3<-mle2(funcao.likelihood,optimizer="nlminb",lower=c(alpha=0.0001,beta=0.0001,desvio=0.0001))
summary(fit.mle3)
 
############################
#OpenBugs
############################
library(R2OpenBUGS)
 
sink("linreg.txt")
cat("
model {
# Prior
 alpha ~ dnorm(0,0.001)
 beta ~ dnorm(0,0.001)
 sigma ~ dunif(0, 100)
 tau <- 1/ (sigma * sigma)
 
# Verossimilhança
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha + beta*x[i]
 }
}
",fill=TRUE)
sink()
 
bugs.data <- list(x=preditor,y=resposta,n=length(resposta))
 
# Função de inicialização
inits <- function(){ list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))}
 
# Escolha os parametros que quer estimar
params <- c("alpha","beta","sigma")
 
# Caracteristicas do MCMC
nc = 3 #Numero de cadeias
ni=1200 #Tamanho da cadeira
nb=200 #Numero de simulação que serão descartadas
nt=1 #Thinning rate
 
# Inicie o Amostrador
 
modelo1.bugs <-bugs(data = bugs.data, inits = inits, parameters = params,model = "linreg.txt",n.thin = nt,
                     n.chains = nc,n.burnin = nb, n.iter = ni)
 
print(modelo1.bugs, dig = 3)
 
############################
#Jags
############################
library(rjags)
 
jags <- jags.model(file='linreg.txt',data=bugs.data,inits=inits,n.chains = 3,n.adapt=100)
update(jags, 1000)
jags.samples(jags,c('alpha', 'beta','sigma'),1000)
 
 
############################
#Stan
############################
library(rstan)
 
linreg_stan <- '
data{
  int n;
  real y[n];
  real x[n];
}
 
parameters{
  real alpha;
  real beta;
  real<lower=0, upper=100> sigma;
}
 
transformed parameters{
    real mu[n];
    for(i in 1:n){
    mu[i] <- alpha + beta*x[i];
  }
}
 
model{
  //Priors
  alpha ~ normal(0, 100);
  beta ~ normal(0, 100);
  sigma ~ uniform(0, 100);
 
  //Likelihood
  for(i in 1:n){
    y[i] ~ normal(mu[i], sigma);
}
}
'
 
stan.data <- list(n = length(resposta), x = preditor, y = resposta)
 
fit.stan<- stan(model_code=linreg_stan,data=stan.data,iter = 1000,chains = 3,thin = 1)
print(fit.stan,digits=3)