Banco de Dados – Usando MySQL do R

Manga-Guide-to-Databases

Um banco de dados é uma coleção organizada de dados. Os dados são organizados em modelos que tentam representar a realidade de uma forma que o processamento(adição, remoção e consulta de dados) se torna mais eficiente, além de uma facilidade de guardar.

Para manter esses dados, a gente pode usar sistemas específicos para exercer todas essas funções da melhor forma possível, muitos deles são de uso geral, e existe a possibilidade de migrar entre eles sem começar tudo do zero. Acho que um bem comum para sistemas não tão grandes é o MySQL.

Mas primeiro, porque diabos estamos preocupados com isso? Porque dados em biologia tendem a desaparecer muito rápido, além de serem normalmente desorganizados. Nesse artigo aqui da revista CELL os autores falam como a cada ano da publicação de um artigo a chance de os dados sumirem aumenta em 17%. Na maioria dos editais de pesquisa, ou projetos, não existe uma preocupação com isso dentro da biologia. Alguma áreas, como em biologia molecular, isso ja está mais organizado, ao forçar todos os autores a publicarem seus dados numa base publica, o genbank. Mas em ecologia, as coisas não são bem assim.
Temos algumas revistas, como a Ecology da ESA que publica base de dados, mas eu acho que são casos a parte.
Mas além disso, eu sempre achei que uma maior preocupação com como os dados são guardados, pode ser muito útil mesmo localmente, ja que o ecossistema da biologia são normalmente laboratórios.
Muitos laboratórios, das mais diversas áreas na biologia, é comum você ver muitas pessoas coletando dados muitos similares(no sentido de todo mundo trabalhar no mesmo ambiente por exemplo), mas guardando em planilhas pessoais, cada um guardando de um jeito. Tudo bem que não existe nada de errado com isso, não é meu objetivo reclamar de como outras pessoas fazem seu trabalho. Mas existe algumas possibilidades que talvez sejam perdidas por isso.

Por exemplo, se cada pessoa guarda os dados do seu jeito, temos muitas dificuldades as vezes de conseguir unir esses dados depois, para fazer coisas maiores.
E continuando nesse pensamento, muitas pessoas vem e vão dos laboratórios, fazendo monografias, dissertações e teses, se todos os dados são guardados num banco de dados bem pensado e organizadinho, posteriormente o laboratório ganha uma base de dados, que pode por si só responder muitas perguntas, e além disso pode ajudar a começar a responder outras, já que dados ao longo do tempo são bem complexos de se conseguir.

Mas antes de se preocupar em montar um banco de dados, vamos ver como seria olhar para esses dados do R, so para ter um gostinho da idéia. Para tal vamos ver um pequeno exemplo usando o pacote RMySQL, que serve para acessar banco de dados administrasdos com o MySQL.

Primeiro, SQL é uma linguagem estruturada para “falar” com o banco de dados, que serve para criar eles, acessar os dados que você quer, alterar, enfim todas as atividades, e ela é diferente do programa que administra o banco de dados(mas todos tentam implementar SQL do mesmo jeito), mais ou menos isso.

O MySQL é um programa gerenciador de banco de dados que implementa essa linguagem, e podemos montar um servidor com vários banco de dados, cada banco de dados com seu próprio schema.

Mas vamos ver como funciona para acessar o banco de dados da UCSC de bioinformática.

Basicamente vão ser três atividades, a gente tem que se conectar ao banco de dados, pra isso a gente tem que ter o endereço dele, e um usuário, que pode ou não ter senha, e o gerenciador pode regular o que um usuário pode fazer. Nesse caso tudo esta explicado nessa pagina aqui para esse banco de dados em questão.

> ucscDb <- dbConnect(MySQL(),user="genome",host="genome-mysql.cse.ucsc.edu") > resultado <- dbGetQuery(ucscDb,"show databases;") > head(resultado) Database 1 information_schema 2 ailMel1 3 allMis1 4 anoCar1 5 anoCar2 6 anoGam1 > dbDisconnect(ucscDb) [1] TRUE

Aqui o primeiro comando serve para conseguir acesso ao banco de dados (usando a internet), veja que fornecemos um endereço, o usuário publico do banco de dados e o primeiro argumento é mais complicado, mas é o driver para fazer a conexão. Isso tudo nos retorna um objeto do R que é a conexão. Com ela agora a gente pode “fazer perguntas” para o banco de dados e recuperar os dados. Isso é o que fizemos com o dbGetQuery, a gente perguntou quais bancos de dados tem la, com o “show databases;”, e depois olhamos a saída, que é uma lista de bases que podemos acessar.

Certo, mas onde estão os dados, estão dentro dessas bases de dados, podemos escolher uma e acessar, por exemplo, se pegarmos a base de dados “hg19”, a gente acessa ele assim.

> hg19 <- dbConnect(MySQL(),user="genome", db="hg19",host="genome-mysql.cse.ucsc.edu") allTables <- dbListTables(hg19) > length(allTables) > [1] 11006 > allTables[1:5] [1] "HInv" "HInvGeneMrna" "acembly" "acemblyClass" "acemblyPep"

Como a gente desconectou la em cima, a gente conecta novamente, e podemos listar as tabelas que contém os dados. É um bocado de tabelas, com informações específicas. Podemos tentar olhar quais os atributos, os campos, colunas de uma tabela dessas assim:

> dbListFields(hg19,"affyU133Plus2") [1] "bin" "matches" "misMatches" "repMatches" "nCount" "qNumInsert" [7] "qBaseInsert" "tNumInsert" "tBaseInsert" "strand" "qName" "qSize" [13] "qStart" "qEnd" "tName" "tSize" "tStart" "tEnd" [19] "blockCount" "blockSizes" "qStarts" "tStarts"

Agora o mais legal é que o SQL é uma linguagem de busca, então podemos mais que olhar uma tabela, podemos pegar os dados de forma organizada, e até fazer buscas usando funções do SQL, por exemplo, podemos mandar um query para falar qual o tamanho de uma tabela, o número de linhas da seguinte forma.

> dbGetQuery(hg19, "select count(*) from affyU133Plus2") count(*) 1 58463

Bem aqui a gente ve que essa tabela tem realmente muitas linhas, isso implica que leva um tempo para baixar ela, mas temos um comando para baixar ela inteirinha como um dataframe no R, e tendo ela na memória no R, é so fazer o que queremos.

> affyData <- dbReadTable(hg19, "affyU133Plus2") > head(affyData) bin matches misMatches repMatches nCount qNumInsert qBaseInsert tNumInsert tBaseInsert strand 1 585 530 4 0 23 3 41 3 898 - 2 585 3355 17 0 109 9 67 9 11621 - 3 585 4156 14 0 83 16 18 2 93 - 4 585 4667 9 0 68 21 42 3 5743 - 5 585 5180 14 0 167 10 38 1 29 - 6 585 468 5 0 14 0 0 0 0 - qName qSize qStart qEnd tName tSize tStart tEnd blockCount 1 225995_x_at 637 5 603 chr1 249250621 14361 15816 5 2 225035_x_at 3635 0 3548 chr1 249250621 14381 29483 17 3 226340_x_at 4318 3 4274 chr1 249250621 14399 18745 18 4 1557034_s_at 4834 48 4834 chr1 249250621 14406 24893 23 5 231811_at 5399 0 5399 chr1 249250621 19688 25078 11 6 236841_at 487 0 487 chr1 249250621 27542 28029 1 . . . tStarts 1 14361,14454,14599,14968,15795, 2 14381,14454,14969,15075,15240,15543,15903,16104,16853,17054,17232,17492,17914,17988,18267,24736,29320,

Eu suprimi um pouco dos resultados, mas aqui a gente tem a tabela inteirinha num dataframe no R. Pronto para usar como quisermos. Mas o SQL vai muito além disso, ele permite comandos para selecionarmos exatamente o que queremos dos dados, além do que essas tabelas são montadas de forma estruturada, de forma que conseguimos recuperar diretamente os dados estruturadinhos.

Mas mais um exemplo colocando uma cláusula para pegar todas as colunas da tabela affyU133Plus2 onde a coluna misMatches estão entre 1 e 3.

> query <- dbSendQuery(hg19, "select * from affyU133Plus2 where misMatches between 1 and 3") > affyMis <- fetch(query) > affyMis[1:10,1:3] bin matches misMatches 1 585 723 3 2 586 740 2 3 73 986 3 4 587 741 1 5 587 985 1 6 587 938 1 7 588 784 3 8 588 1397 3 9 589 483 1 10 589 938 1

E daqui podemos analisar os dados

> quantile(affyMis$misMatches) 0% 25% 50% 75% 100% 1 1 2 2 3 Lembrando que como temos uma conexão com o banco de dados, sempre temos que fechar ela depois. > dbClearResult(query) [1] TRUE > dbDisconnect(hg19) [1] TRUE

Mas imaginem que legal, primeiro que teríamos uma base de dados disponível no laboratório que trabalhos direto, pessoas vem e vão, mas os dados sempre vão estar la. Provavelmente mais seguros, porque é mais fácil fazer uma cópia de segurança de um banco de dados que de mil planilha espalhadas nas mãos de muitas pessoas. Outra coisa é que para montar um banco de dados você precisa de um schema, uma estrutura de dados, isso pode nos forçar a pensar melhor sobre o que fazemos, sobre como coletamos dados, quais dados coletamos e porque coletamos, e pessoas que estiverem começando algum trabalho, podem se beneficiar sobre as caneladas passadas que foram corrigidas. Além do que assim, cedo ou tarde é fácil ter uma base de dados grande, se ela for alimentada por todas as pessoas que passam pelo laboratório. Além de que, quem passa pelo laboratório tem respaldo de dados para escrever projetos, o que facilitaria justificar quantidades, número de amostras ou pensar a quantidade de dados que vamos achar, ou se algum projeto é viável.

Outra coisa muito legal, é tornar mais fácil tudo ser reproduzível, só a gente guardar um script do R com os querys do banco de dados e analises realizadas (só controlar datas da busca, já que a base pode crescer sempre) e qualquer analise fica na mão.

Bem eu acho uma idéia bem legal, já que ninguém perde, porque todo mundo tem que coletar seus dados para concluir seja la o que for que estiver fazendo, mas a ciência ganha se tudo ficar guardadinho, e puder ser re-analisado conjugando dados dos vários trabalhos feitos anteriormente, e quem sabendo abrindo a possibilidade de responder coisas mais legais ainda.

Eu vi esse exemplo no curso de R do coursera chamado Getting and Cleaning Data, que mostra muitas outras formas de dados, e é legal para ver como a preocupação com esse tipo de coisa vem aumentando, com a tecnologia fazendo a gente cada vez conseguir mais dados, é importante conseguir guardar tudo direito, de forma estruturada para não se afogar em informação, mas banco de dados é uma área bem complexa, mas muito legal, e acho que não distante para começar empreitadas para guardar dados nos diversos laboratórios que cada um fica, mesmo disciplinas de pós-graduação no campo tendem a gerar muitos dados que ficariam muito legais se guardados num banco de dados bonitinho. 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 e até mais.

install.packages("RMySQL")
 
library(RMySQL)
 
ucscDb <- dbConnect(MySQL(),user="genome",host="genome-mysql.cse.ucsc.edu")
resultado <- dbGetQuery(ucscDb,"show databases;")
head(resultado)
dbDisconnect(ucscDb)
 
 
#Connecting to hg19 and listing tables
hg19 <- dbConnect(MySQL(),user="genome", db="hg19",host="genome-mysql.cse.ucsc.edu")
allTables <- dbListTables(hg19)
length(allTables)
allTables[1:5]
 
#Get dimensions of a specific table
dbListFields(hg19,"affyU133Plus2")
dbGetQuery(hg19, "select count(*) from affyU133Plus2")
 
#Read from the table
affyData <- dbReadTable(hg19, "affyU133Plus2")
head(affyData)
 
#Select a specific subset
query <- dbSendQuery(hg19, "select * from affyU133Plus2 where misMatches between 1 and 3")
affyMis <- fetch(query)
affyMis[1:10,1:3]
quantile(affyMis$misMatches)
 
dbClearResult(query)
dbDisconnect(hg19)

Data Augmentation para dados de Captura e Recaptura

Acho que um tipo de dados muito comum em ecologia são os de Captura e Recaptura. Desde aqueles feitos com armadilhas, onde elas são vistoriadas periodicamente num grid, até as armadilhas de câmeras, que tiram fotos por sensor de movimento, redes capturando passarinhos, peixes, borboletas, talvez mais algumas formas de conseguir dados assim.

A ideia é que capturando indivíduos de alguma espécie, e recapturando depois (basta dar um jeito de identificar o um indivíduo) a gente pode calcular uma serie parâmetros populacionais que temos muito interesse, taxas de crescimento, sobrevivência, migração e assim vai.

Mas esses modelo sempre foram meio misticos pra mim, um monte de formula fechada que vem em um programa, a gente usa sem saber o que é, e usa da forma fechadinha. Mas se a gente conseguir ver esse tipo de modelo como um glm, as portas se abrem para um monte de possibilidades. Nesse post aqui a gente viu como estimar a colonização e extinção, usando o pacote unmarked e um modelo Bugs equivalente.

No modelo final, a gente vê que coisas como colonização, extinção são na verdade medidas derivadas, e que a resposta propriamente dita é encontrar ou não um individuo em uma amostra, observar ele ou não, a captura. Sendo que indivíduos podem ser recapturados. Resumindo, tentamos fazer um modelo assim.

logit(p_{i,j})=\alpha_i+\beta_j+\delta \cdot F_{i,j}+ \gamma \cdot X_{i,j}

E o \alpha_i \sim Normal(\mu_{\alpha},\sigma^2_{\alpha})

Vamos entender isso tudo, primeiro aqui, a gente ta lidando com dois índices, o i e o j, aqui o i diz respeito aos indivíduos e j ao evento de captura. Por exemplo pegando um exemplo em que a gente arma uma rede para pegar passarinhos, o i é referente a identidade do passarinho, o número da anilha dele, e o j é o dia que a rede for armada, no final isso resulta numa matriz de dados onde os indivíduos são linhas e os dias de rede armada as colunas, e dentro dessa matriz a gente marca as captura com um e zero para quando não houve captura, sendo que uma linha pode ter mais de uma captura, a gente pode recapturar um mesmo passarinho.

Agora seguindo esse exemplo, indivíduos, passarinhos podem diferir em capturabilidade, a gente pode tratar isso, tratar o fato de alguns passarinhos serem mais bobinhos e outros não assumindo que eles vem de uma distribuição normal para essa variabilidade, teremos então alguns bobinhos, alguns espertalhões mas a maioria vai ter uma determinada capturabilidade, isso é o \alpha_i da formula, veja que ele usa o índice i de individuo, e dentro desse \alpha podemos modelar outras coisas dos indivíduos, sei la, sexo por exemplo, se machos são mais curiosos, derrepente sejam mais capturados.

Agora a chance de uma captura não depende só de um individuo, por exemplo um dia ensolarado pode ter mais passarinho voando que um dia de chuva, ou dia depois da chuva, ou seja, o dia da amostra pode influenciar na capturabilidade. Isso ai a gente vei colocar dentro do \beta_j, lembrando que j é o índice dos eventos de captura. Então a gente pode modelar as coisas desse tipo aqui dentro. Lembrando que tanto no \alpha quanto o \beta, a gente ta colocando uma letrinha, mas cabe uma formulinha aqui dentro, com intercepto e inclinação de por exemplo, dados que você coleta dos passarinhos, tipo a morfometria, se você acha que ela tem a ver com a chance do passarinho ser capturado, podemos tentar observar isso.

Agora as outras duas coisas, o \delta \cdot F_{i,j} tem haver com o famoso ditado, gato escaldado tem medo de água, ou o contrario. Apos primeira comida, a chance de captura pode ser alterada, tanto para menos quanto pra mais, por exemplo, se ser capturado é muito estressante, o individuo pode começar a ficar mais esperto e difícil de capturar, ou se a gente ta lidando com uma gaiola com alguma isca, é fácil pensar que passa pela cabeça do individuo a frase “boca livre”, e ele pode começar a ser capturado com mais facilidade por procurar uma boca livre, ou sei la, veja que esse parâmetro depende tanto do individuo i como da captura j, já que esse efeito a partir da segunda captura, se alguém é capturado apenas uma vez, não tem como apresentar esse efeito, por isso leva os dois índices.

O último parâmetro que tem ali é o \gamma \cdot X_{i,j}, que ai são as questões ambientais, como é o ambiente, ou a disponibilidade de algo especial. Isso a gente enquadra ai nesse \gamma, veja que ele também depende de i e j.

Certo então com o modelo ali em cima a gente conseguiria estimar uma probabilidade de captura, mas qual o interesse nisso?
A gente se interessa em quantidades mais palpáveis, o tamanho da população, sobrevivência, imigração, mas essas quantidade pode ser derivada dessa chance de captura, sendo que com o modelo acima a gente tem a chance de considerar praticamente tudo que é do nosso interesse.

Agora vamos simular uns dados aqui para ver melhor como isso acontece. Vamos supor uma população de 100 indivíduos, ai vamos supor também uma chance de captura de 50% e que coletamos por três vezes, em três eventos de coleta.

Beleza, se tudo que colocamos no modelo la em cima for constante, teríamos o seguinte.

> N = 100 > p = 0.5 > T = 3 > Real<-array(NA, dim = c(N, T)) > for (j in 1:T){ + Real[,j] <- rbinom(n = N, size = 1, prob = p) + } > Real [,1] [,2] [,3] [1,] 0 1 1 [2,] 1 0 1 [3,] 0 1 0 [4,] 0 0 1 . . . [95,] 0 1 1 [96,] 1 0 1 [97,] 0 0 0 [98,] 1 0 1 [99,] 1 1 1 [100,] 0 1 1

Certo, colocamos nossos parâmetros e meramente sorteamos de uma distribuição binomial a com a chance de captura que definimos se capturamos ou não. Temos uma matriz de 100 linhas, que são 100 indivíduos e 3 colunas, que são os três eventos de coleta, agora da uma olhada o que aconteceu na linha 97, a gente nunca capturou esse individuo nos três eventos.

Se a gente somar as linhas dessa matriz, a gente vai ver que

> detectado <- apply(Real, 1, max) > sum(detectado) [1] 86

vimos somente 86 indivíduos dos 100 que existe, existem 14 sortudos que nunca caíram na rede, simples assim.
Agora a pergunta é a seguinte, existem 100 indivíduos, mas a gente não ve uma matriz como a acima, esse valor de 100 esta oculto pela natureza, aqui a gente sabe porque inventamos ele, mas com dados reais a gente não sabe.
Os dados reais que a gente ve vão ser assim:

> dados<-array(NA, dim = c(N, T)) > dados <- Real[detectado == 1,] > dados [,1] [,2] [,3] [1,] 0 1 1 [2,] 1 0 1 [3,] 0 1 0 . . . [83,] 1 0 1 [84,] 1 0 1 [85,] 1 1 1 [86,] 0 1 1

So vão até a linha 86, porque as 14 linhas que so tem zero a gente não viu, não ve e nunca verá. Nossas observações são sempre imperfeitas. Agora como faz para chegar aquele valor de 100 indivíduos, se a gente so tem nossas observações acima?

Existiam outros estimadores, presentas em programas como o Capture e o Mark, eles estimavam tudo e depois faziam uma conta para estimar o tamanho da população, aquele valor de 100. Agora quando a gente ta usando estatística bayesiana, com MCMC e tal, existe um problema que a gente vai de iteração a iteração, e o vetor de N pode mudar de iteração a iteração, ou seja não cabe nos dados que a gente manda para o BUGS. Mas como problemas existem para serem solucionados, dois caras chamados Tanner e Wong inventaram uma técnica chamada de data augmentation, que o Royle e companhia trouxeram pros modelos de captura e recaptura usando MCMC.

Basicamente duas coisas são feitas, adicionamos um numero arbitrário de zeros a matriz de dados analisamos uma versão reparametrizada dos dados originais. Essencialmente convertendo um modelo de população fechada para um modelo de ocupação.

Então a gente está adicionando um grande numero de observações em potencial que não foram observadas. Os dados aumentados terão um número de linhas M e o mesmo numero de colunas T das outras observações, se o tamanho da população deveria ser o número de linhas, que deveria ser 100, e só temos 86 linhas, ao aumentarmos os dados teremos uma quantidade linhas agora muito muito maior que o tamanho da população, que é o que nos interessa. Nesses dados a gente ajusta um tipo de modelo de zero inflado se o tamanho da população fosse conhecido. Para isso a gente adiciona um indicador binario que indica se uma linha da matriz de dados aumentada representa um individuo real ou um que não existe na pratica. Esse indicador é dado por um prior da distribuição de Bernoulli, ou Binomial, e esse parâmetro, que podemos chamar de \omega (omega) é estimável dos dados. Uma probabilidade de inclusão. Com essa probabilidade de inclusão a gente tem o tamanho da população saindo dos dados.

Então criamos nosso modelo na linguagem bugs

# Modelo na linguagem Bugs
sink("modelo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        p ~ dunif(0, 1)
        # Likelihood
        for (i in 1:M){
            z[i] ~ dbern(omega)
            # Indicador de inclusão
            for (j in 1:T){
                dadosaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p # So pode ser detectado se z=1
                } #j
            } #i
        # Quantidades derivadas
        N <- sum(z[])
        }
    ",fill = TRUE)
sink()

Vamos observar alguma coisas, primeiro que temos dois priores

omega ~ dunif(0, 1)
p ~ dunif(0, 1)

Que são as chances daquela linha ser um indivíduo e a chance do indivíduo ter sido capturado.

Veja então que no modelo de verossimilhança, a gente tem uma chance do individuo ser existir

for (i in 1:M){
    z[i] ~ dbern(omega)

E a chance dele ter sido capturado num dado evento pra esse individuo

for (j in 1:T){
    dadosaug[i,j] ~ dbern(p.eff[i,j])
    p.eff[i,j] <- z[i] * p # So pode ser detectado se z=1

Mas veza que ele é multiplicado pelo valor de z, que vai ser zero ou um, com uma chance omega, que é a chance dele ser um individuo real.

N <- sum(z[])

Ou seja a soma do vetor de z vai ser exatamente a estimativa do tamanho populacional.

Agora é só juntar os dados, gerar valores iniciais, configurar as cadeias e rodar o esse modelo. E temos o seguinte resultado.

> 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 95.451 4.114 89.000 92.000 95.000 98.000 105.000 1.001 2700 p 0.542 0.037 0.466 0.517 0.542 0.567 0.613 1.001 4100 omega 0.405 0.036 0.337 0.380 0.404 0.429 0.478 1.001 4000 deviance 395.385 18.993 363.200 383.300 393.500 407.200 436.600 1.001 2800 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 = 180.3 and DIC = 575.7 DIC is an estimate of expected predictive error (lower deviance is better).

Sempre temos que fazer aquelas validações, que tivemos uma convergência dos parâmetros, mas podemos ver que funciona bem.

Primeiro vamos olhar o N que é o tamanho populacional

mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff N 95.451 4.114 89.000 92.000 95.000 98.000 105.000 1.001 2700

Veja que temos uma média de 95 para o parâmetro, que é perto de 100, o valor real do tamanho populacional, com um desvio padrão de 4.
Olhando o intervalo de confiança vemos que ele vai de 89 a 105, o que inclui o valor real de 100, ou seja, deu certo, esse negocio funciona mesmo para estimar o tamanho populacional.

E estimamos bem também a chance de captura, que na realidade era de 0.5 (50%), que esta dentro do intervalo de confiança.

mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff p 0.542 0.037 0.466 0.517 0.542 0.567 0.613 1.001 4100

Então vamos tentar fazer um gráfico para ver melhor essa saída.

01

Certo, nossa estimativa cobre o valor real, dentro do intervalo de confiança de 95%, não sei nem se posso falar intervalo de confiança, mas é o 95% ali das saídas, as vezes esses termos tem pequenos pressupostos que a gente não sabe, até escrever errado.
Mas agora podemos ver como o \omega (omega) faz, vamos olhar o valor dele.

mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff omega 0.405 0.036 0.337 0.380 0.404 0.429 0.478 1.001 4000

Média de 0.405 com desvio de 0.036, esse 40.5% é o número de linhas da matriz de dados aumentada que mandamos que são provavelmente realmente indivíduos. Como a gente tinha 86 observações, e aumentamos a matriz em 150 linhas, temos um total de 236 linhas. E dessas 40.5% consideramos como indivíduos.

> 0.405*(86+150) [1] 95.58

As outras não devem ser observações, apenas linhas. Mas o aumento deu a possibilidade de considerar esses individuos que não observamos, mas deviam estar la.
Mas uma pergunta razoavel é o quanto devemos aumentar a matriz. Porque a gente aumentou 150 aqui e deu certo, mas será que qualquer número funciona?
Bem sabemos que o tamanho da população vai ser igual ao número das observações, ou maior. Se aumentamos pouco a matriz, por exemplo, somente em 5 linhas, podemos ver o que acontece.

> print(out5, 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 90.004 1.093 87.000 89.000 90.000 91.000 91.000 1.001 6000 p 0.574 0.031 0.512 0.553 0.574 0.595 0.636 1.001 5300 omega 0.978 0.019 0.930 0.970 0.984 0.992 0.999 1.002 6000 deviance 369.337 5.821 354.500 365.200 369.700 373.700 376.900 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 = 16.9 and DIC = 386.3 DIC is an estimate of expected predictive error (lower deviance is better).

Veja que o valor de omega ficou espremido, grudado no 100%, olhando o grafico isso fica mais facil de visualizar.

02

Veja que a distribuição é bem truncada, como se tivesse uma barreira ali segurando a distribuição de ir mais para direita, ou seja a gente precisava de mais “espaço”. Então se a gente ver uma distribuição truncada assim, sabemos que temos que aumentar o tamanho da matriz mais um pouco. Com 150 deu legal, podemos olhar denovo isso.

> print(out150, 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 95.551 4.222 89.000 93.000 95.000 98.000 105.000 1.001 6000 p 0.542 0.038 0.467 0.517 0.543 0.568 0.615 1.002 1700 omega 0.406 0.037 0.339 0.381 0.405 0.430 0.481 1.001 5800 deviance 395.845 19.407 363.300 383.300 393.500 407.400 438.710 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 = 188.3 and DIC = 584.2 DIC is an estimate of expected predictive error (lower deviance is better).

03

Aqui a distribuição ficou legal, veja que o omega não ficou espremido, truncado nem para a esquerda nem para a direita. Então está tudo ok. Mas e se a gente aumentar demais, tem problema? Será que o precavido aqui vai se ferrar?

> print(out1500, 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 95.501 4.175 89.000 92.000 95.000 98.000 105.000 1.001 6000 p 0.542 0.037 0.468 0.517 0.542 0.568 0.614 1.001 2800 omega 0.061 0.006 0.049 0.056 0.061 0.065 0.074 1.001 3400 deviance 395.610 19.188 363.300 383.300 393.600 407.100 437.600 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 = 184.1 and DIC = 579.7 DIC is an estimate of expected predictive error (lower deviance is better).

04

E a resposta é não, não tem problema nenhum aumentar muito. O que acontece é que o omega vai ficar baixinho, mas não ficando truncado ta tudo certo. Veja que aumentando em 1500 linhas a matriz, temos exatamente a mesma estimativa de tamanho populacional, unica coisa é que aqui tivemos uma matriz gigante para processar, e aqui o 6.1% de 1586 linhas que é o total de indivíduos.

Eu achei essa solução muito mirabolante. Primeiro porque a gente cai em um glm, o que torna um modelo de captura e recaptura algo bem próximo dos glm que comumente olhamos e essa sacada de aumentar a matriz é muito foda. Melhor sobrar que faltar. Agora aqui conseguimos ter uma noção de que data augumentation funciona, mas em um modelo simples.

No inicio do post colocamos o seguinte modelo.
logit(p_{i,j})=\alpha_i+\beta_j+\delta \cdot F_{i,j}+ \gamma \cdot X_{i,j}
E o \alpha_i \sim Normal(\mu_{\alpha},\sigma^2_{\alpha})

Mas o que ajustamos foi algo como

logit(p_{i,j})=\alpha

Somente um parâmetro, so lembrar do modelo, que so tem uma linha aproximando nossas observações.

dadosaug[i,j] ~ dbern(p.eff[i,j])

Agora temos que começar a expandir esse parametro para o modelo completo. Sendo que não somos obrigados a tentar ajustar tudo, podemos usar somente aquilo que nos interessa, ou que temos dados coletados.

Bem é isso ai, hora de ir dormir. 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 boa noite.

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

#Gerando dados
N = 100
p = 0.5
T = 3
Real<-array(NA, dim = c(N, T))
for (j in 1:T){
    Real[,j] <- rbinom(n = N, size = 1, prob = p)
    }
#A visão da natureza
Real
 
detectado <- apply(Real, 1, max)
sum(detectado)
 
dados<-array(NA, dim = c(N, T))
dados <- Real[detectado == 1,]
#A nossa visão
dados
 
# Aumentando os dados com 150 individuos em potencial
nz <- 150
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
 
# Modelo na linguagem Bugs
sink("modelo.txt")
cat("
    model {
        # Priors
        omega ~ dunif(0, 1)
        p ~ dunif(0, 1)
        # Likelihood
        for (i in 1:M){
            z[i] ~ dbern(omega)
            # Indicador de inclusão
            for (j in 1:T){
                dadosaug[i,j] ~ dbern(p.eff[i,j])
                p.eff[i,j] <- z[i] * p # So pode ser detectado se z=1
                } #j
            } #i
        # Quantidades derivadas
        N <- sum(z[])
        }
    ",fill = TRUE)
sink()
 
# Juntando os dados para mandar para o bugs
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
# Função geradora de valores iniciais
inits <- function() list(z = rep(1, nrow(dadosaug)), p = runif(1, 0, 1))
# Parametros a serem monitorados
params <- c("N", "p", "omega")
 
# Configurando o MCMC
ni <- 2500
nt <- 2
nb <- 500
nc <- 3
# Chamando o openbugs do R
library(R2OpenBUGS)
out <- bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
 
 
# Distribuições posteriores
print(out, dig = 3)
 
str(out)
out$summary
 
#Figura 1
hist(out$sims.list$N,nclass=50,col="gray",main="",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))
 
 
nz <- 5
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
 
out5 <- bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
print(out5, dig = 3)
 
#Figura 2
hist(out5$sims.list$N,nclass=50,col="gray",main="Dados aumentados em 5",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out5$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))
 
 
nz <- 150
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
 
 
out150 <-bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
print(out150, dig = 3)
 
#Figura 3
hist(out150$sims.list$N,nclass=50,col="gray",main="Dados aumentados em 150",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out150$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))
 
 
nz <- 1500
dadosaug <- rbind(dados, array(0, dim = c(nz,T)))
bugs.data <- list(dadosaug = dadosaug, M = nrow(dadosaug), T = ncol(dadosaug))
 
out1500 <- bugs(data=bugs.data,inits=inits,parameters=params,model="modelo.txt",
            n.chains = nc,n.thin = nt, n.iter = ni, n.burnin = nb)
print(out1500, dig = 3)
 
#Figura 4
hist(out1500$sims.list$N,nclass=50,col="gray",main="Dados aumentados em 1500",xlab ="Tamanho Populacional",las = 1, xlim = c(80, 150),
     ylab="Frequência")
abline(v = sum(detectado),lwd = 3)
abline(v = 100,lwd=3,col="red")
abline(v = out1500$summary[1,c(3,7)],lwd = 1,col="blue",lty=2)
legend("topright",legend=c("Detectados","Tamanho real","Intervalo de 95% estimado"),lwd=c(3,3,1),lty=c(1,1,2),
       col=c("black","red","blue"))

Phase Plane

Quando falamos de equilíbrio de duas populações aqui a gente citou a palavra chave State Space, que são as possibilidades que o sistema pode assumir, sejam elas possíveis ou não.

Sistema é um sistema oras, um conjunto de elementos interconectados de modo a formar um todo organizado.

Quando a gente falava de competição entre espécies, o todo é o ecossistema de duas espécies, os elementos são as espécies, a elas são interconectadas pelo efeito que uma tem sobre a outra, a competição.

Lembrando nosso sistema

\begin{cases}  {dN_{1}\over{dt}} = r_{1}\cdot N_{1}\cdot (1-{N_{1}\cdot \alpha_{11}} - {N_{2}\cdot \alpha_{12}}) \\  {dN_{2}\over{dt}} = r_{2}\cdot N_{2}\cdot (1-{N_{2}\cdot \alpha_{22}} - {N_{1}\cdot \alpha_{21}})  \end{cases}

Que era mais fácil de visualizar com essa figura.

03

Beleza, comumente a gente resolvia esse sistema ao longo do tempo, para ver o que acontecia com as espécies, pra ver o que é possível e o que não é possível.

01

Nesse caso, elas iam crescendo, ai quando a população aumentava, elas começavam competir, as abundâncias iam mudando mas chegava uma hora que tudo estacionava num ponto, um ponto de equilíbrio, e la ficava a não ser que houve-se uma pertubação.

Mas ao invés de olhar dessa forma, podemos fazer um gráfico da abundância da espécie 1 em relação a abundância da espécie 2.

02

E esse que é o tal do Phase Plane, aqui em duas dimensões, porque temos duas espécies, mas poderia ser em três, n dimensões, claro que para mais de três dimensões não da para fazer exatamente um gráfico, mas da para tentar visualizar ele também.

Mais que fluflu para fazer charme, esse tipo de figura nos ajuda a determinar como funciona um sistema, se ele tem algum tipo de equilíbrio ou não, ciclo, essas coisas. A gente viu aqui os casos de equilíbrio para duas populações. Mas sempre tudo ia para algum ponto, meio diretão assim, Mas existem outros equilíbrios bem legais também.

Se a gente olhar o help da função ode a gente conseguir ver que a autora mandou um exemplo pronto do famoso sistema de predador e presa com uma presa com crescimento logístico.

\begin{cases}  {Ingestion} = rIng * Prey * Predator \\  {GrowthPrey} = rGrow * Prey * (1-\frac{Prey}{K}) \\  {MortPredator} = rMort * Predator \\  {dPrey} = GrowthPrey - Ingestion \\  {dPredator} = Ingestion * assEff - MortPredator \\  \end{cases}

O sistema é esse conjunto de equações, que fica mais simples talvez escrita na função do R para usar com o ode

LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    Ingestion    <- rIng  * Prey * Predator
    GrowthPrey   <- rGrow * Prey * (1 - Prey/K)
    MortPredator <- rMort * Predator
 
    dPrey        <- GrowthPrey - Ingestion
    dPredator    <- Ingestion * assEff - MortPredator
 
    return(list(c(dPrey, dPredator)))
  })
}

Ao resolver esse sistemas a gente tem o seguinte comportamento de predadores e presas ao longo do tempo.

03

Agora se a gente fizer um Phase Plane, para olhar esse sistema (lembrando que definimos um certo set de parâmetros, como vimos para populações, dependendo dos valores dos parâmetros, os alphas la, víamos coisas diferentes).

04

Uma espiral, o início é onde tudo está oscilando mais, mas com o tempo as oscilações entre abundâncias de presas e predadores vai diminuindo até estabilizar em um ponto fixo.

Agora vamos olhar ainda mais um pouquinho o sistema de predadores e presas, mas sem a presa com crescimento logístico, vamos pegar o modelo que aparece no wikipedia nessa página aqui.

O modelo é esse aqui.

\begin{cases}  {\frac{dx}{dt}} = x(\alpha - \beta y) \\  {\frac{dy}{dt}} = - y (\gamma - \delta x) \\  \end{cases}

Que podemos escrever no R assim:

predadorpresa <- function(Time, State, Pars) {
    with(as.list(c(State, Pars)), {
        dpresa <- presa*(alpha-beta*predador)
        dpredador <- -predador*(gamma-delta*presa)
        list(c(dpresa, dpredador))
  } )
}

Claro que cada parâmetro tem um significado especial, mas não vamos ficar discutindo isso agora, a gente só quer ver os phase planes legais.
Mas então resolvendo esse sistema para um set de parâmetros e valores iniciais de predador e presa vemos o seguinte phase plane.

05

Uma bolinha, nesse caso aqui, a gente não tem um ponto de estabilidade, mas a gente tem um tipo de ciclo, esse trilho ai e o sistema nunca vai parar, ele vai ficar oscilando entre vários estados, mas os estados vão estar nessa figura ae, nesse circulo. Mas ainda diferente dos sistemas anteriores onde víamos ele “andando certinho”, a gente pode ligar linhas de um tempo para o próximo e ver algo legal.

06

A gente não esta indo suavemente de um ponto para o próximo, a gente ta indo de um canto do phase plane para o outro la longe. Melhor que isso, resolvemos o sistema para 20 tempos, então podemos numerar as bolinhas na ordem, para ver melhor onde tudo esta começando e acabando.

07

Bem mais bagunçado que a competição entre duas espécies não? Mas ainda sim existe um padrão ai, e o equilíbrio é mais algo como um ciclo que um ponto.

Veja só como a população esta oscilando ao longo do tempo, para gerar esse Phase Plane

08

Muito legal, agora a gente sabe que esse tipo de gráfico tem um nome, e pode ajudar, visualmente a entender melhor como as coisas funcionam, e na biologia a gente consegue descrever uma boa quantidade de coisas com sistemas de equações assim. Claro que é mais complicado que isso, entender, a partir formulações matemática de sistemas, como eles funcionam, o que podemos esperar e o que não devemos esperar, mas esse tipo de gráfico é bem legal, para visualizar algumas coisas que são talvez bem complicadas de entender somente olhando as equações.

Bem, acho que é isso ai por enquanto. O script está no repositório do github, além de aqui em baixo, que da pra olhar os parâmetros usados, que não citei la em cima, ainda fica faltando estudar as tal das matrizes jacobianas com mais detalhe, que é uma forma de achar pontos de equilíbrio, eu acho que serve pra isso pelo menos, na verdade falta estudar tanta coisa que nem sei por onde começar as vezes.

Referência:
M Henry H Stevens 2011 A primer of ecology with R. Springer

library(deSolve)
 
#Modelo para crescimento continuo
lvcomp2 <- function(t, n, parms) {
  with(as.list(parms), {
    dn1dt <- r1*n[1]*(1-a11*n[1] - a12*n[2])
    dn2dt <- r2*n[2]*(1-a22*n[2] - a21*n[1])
    list(c(dn1dt, dn2dt))
  } )
}
 
parms <- c(r1=0.8,r2=0.5,a11=0.010,a21=0.005,a22=0.010,a12=0.005);
initialN<-c(1,1)
out<-ode(y=initialN, times=1:50, func=lvcomp2, parms=parms)
 
#Espécies ao longo do tempo
#figura 1
matplot(out[,1],out[,-1],type='l',xlab="Tempo",ylab="Tamanho populacional",
        frame.plot=F,main="Crescimento logístico contínuo para duas espécies")
legend("right",c(expression("Espécie 1 "*(alpha[21]==0.005)),
                expression("Espécie 2 "*(alpha[12]==0.005))),
       lty=1:2,bty='n',col=c("black","red"))
 
#Representação das abundancias das duas espécies num plano
#figura 2
plot(1, 1, type = "n", ylim = c(0, 100), xlim = c(0,100),frame=F,
     ylab = expression("N"[2]),xlab=expression("N"[1]))
points(out[,2],out[,3],type="b",cex=0.5,pch=19)
 
 
##########################################
#  #
##########################################
LVmod <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    Ingestion    <- rIng  * Prey * Predator
    GrowthPrey   <- rGrow * Prey * (1 - Prey/K)
    MortPredator <- rMort * Predator
 
    dPrey        <- GrowthPrey - Ingestion
    dPredator    <- Ingestion * assEff - MortPredator
 
    return(list(c(dPrey, dPredator)))
  })
}
 
pars  <- c(rIng   = 0.2,    # /day, rate of ingestion
           rGrow  = 1.0,    # /day, growth rate of prey
           rMort  = 0.2 ,   # /day, mortality rate of predator
           assEff = 0.5,    # -, assimilation efficiency
           K      = 10)     # mmol/m3, carrying capacity
 
yini  <- c(Prey = 1, Predator = 2)
times <- seq(0, 200, by = 1)
out   <- ode(yini, times, LVmod, pars)
 
## User specified plotting
#figura 3
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "Tempo", ylab = "Abundância",
        main = "Lotka-Volterra", lwd = 2,frame=F)
legend("topright", c("Presa", "Predador"), col = 1:2, lty = 1:2,bty="n")
 
#figura 4
plot(out[,2],out[,3],type="b",pch=19,frame=F,xlab="Presa",ylab="Predador")
 
 
#########################################
 
predadorpresa <- function(Time, State, Pars) {
    with(as.list(c(State, Pars)), {
        dpresa <- presa*(alpha-beta*predador)
        dpredador <- -predador*(gamma-delta*presa)
        list(c(dpresa, dpredador))
  } )
}
 
pars <- c(alpha=2,beta=1.1,gamma=1,delta=0.4)
yini<-c(presa=5,predador=1)
times <- seq(0, 20, by = 1)
out   <- ode(yini, times, predadorpresa, pars)
 
#figura 5
plot(out[,2],out[,3],type="p",pch=19,xlab="Presa",ylab="Predador")
 
#figura 6
plot(out[,2],out[,3],type="b",pch=19,xlab="Presa",ylab="Predador")
 
#figura 7
plot(out[,2],out[,3],type="l",pch=19,lty=3,lwd=0.5,xlab="Presa",ylab="Predador")
text(out[,2],out[,3],0:20)
 
#figura 8
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "Tempo", ylab = "Abundância",
        main = "Predador-Presa", lwd = 2,frame=F)
legend("topright", c("Presa", "Predador"), col = 1:2, lty = 1:2,bty="n")

Mergesort em C++ usando o pacote Rcpp

A muito tempo atrás eu comecei a ver como dava para usar código em linguagem C do R neste post aqui. Depois disso, no comentário alguém mais esperto que eu falou sobre um tal de pacote Rcpp, que permitia de maneira mais simples, usar código em linguagem C++ do R, depois disso eu fiz mais alguns testes, tentando implementar alguns algorítimos de ordenação, o bubblesort aqui e o insertion sort aqui.

Esse é o tipo de coisa inútil de se fazer, mas é divertido. E esses algorítimos de ordenação tem em qualquer curso introdutório de programação, então são implementações básicas para qualquer linguagem. Depois do mergesort eu ainda tenho o quicksort, mas no mergesort tudo começa a ficar mais complexo.

Por curiosidade, o método padrão na função sort() do R é o Shellsort, criado por Donald Shell mas o que está implementado no R é uma variante do Robert Sedgewick, esse cara da quatro cursos sobre programação gratuitos no Coursera, so olhar aqui.

Mas na função sort do R também temos o Quicksort, mas ele vai ficar pra próxima.

Agora vamos ver o mergesort. Ele é um algorítimo recursivo, ou seja uma função que chama ela mesma, e a ideia básica é, melhor que ordenar todo um vetor, é pegar dois vetores ordenados e então simplesmente intercalar eles. Mas a gente precisa de dois vetores ordenados pra poder usar esse esquema certo?
Pra isso é so a gente quebrar o vetor no meio, e de novo no meio para as duas metades, e de novo para os quatro quartos, e de novo até que vai chegar um momento que so vamos ter mini vetores de um elemento, um vetor de um elemento está necessariamente ordenado, logo agora é so voltar intercalando esses vetores que fica tudo bem.

MergeSort

É até estranho como isso funciona, mas funciona. Bem o legal que diferente do InsertionSort e o BubbleSort, aqui a gente vai precisar de uma função auxiliar, então a gente vai escrever o código em um arquivo, que vamos chamar de merge.cpp, e depois vamos chamar esse arquivo do R usando a função sourceCpp() do pacote Rcpp.

Então o código do mergesort é assim:

#include <Rcpp.h>
using namespace Rcpp;
 
void intercala(int p, int q, int r, NumericVector v,  NumericVector w)
{
  int i, j, k;
   i = p;
   j = q;
   k = 0;
   while (i < q && j < r) {
      if (v[i] < v[j]) {
         w[k] = v[i];
         i++;
      }
      else {
         w[k] = v[j];
         j++;
      }
      k++;
   }
   while (i < q) {
      w[k] = v[i];
      i++;
      k++;
   }
   while (j < r) {
      w[k] = v[j];
      j++;
      k++;
   }
   for (i = p; i < r; i++)
      v[i] = w[i-p];
}
 
void mergesort(int p, int r, NumericVector v, NumericVector aux)
{
   int q;
   if (p < r - 1) {
      q = (p + r) / 2;
      mergesort(p, q, v,aux);
      mergesort(q, r, v,aux);
      intercala(p, q, r, v,aux);
   }
}
 
// [[Rcpp::export]]
NumericVector mergesortC(NumericVector vetor) {
  Rcpp::NumericVector saida = Rcpp::clone(vetor);
  Rcpp::NumericVector aux = Rcpp::clone(vetor);
  int n = saida.size();
  mergesort(0,n,saida,aux);
  return saida;
}

Então, tudo isso está dentro de um arquivo que eu chamei de merge.cpp
Basicamente o código você acha em qualquer esquina da rede mundial de computadores, mas especificamente aqui, temos que importar a biblioteca Rcpp.h, já que estamos usando um objeto específico dela, que é o NumericVector. Outra coisa é que o R manda um ponteiro pra cá, que aponta aonde estão os dados, desse objeto. Se a gente tentar mexer direto nele, não vai dar certo, eu tentei e fiz o R travar com isso heheh, ai perguntando no forum aqui, o cara que criou o pacote me mandou ir estudar, ai lendo o livro dele eu vi que mexer direto nos dados, muitas vezes não é uma boa ideia, dai a gente usa o método clone de NumericVector, a gente clona ele para outro vetor, e então manipula esse novo vetor, pra não mexer na memoria que não devemos.
Outra coisa é que, talvez um ponto fraco do mergeSort, é que para intercalar dois vetores, precisamos de outro lugar para ir guardando eles em ordem, ou seja, diferente do insertionsort, não da para ir arrumando dentro do próprio vetor, por isso precisamos de um vetor auxiliar, aqui eu ja criei um vetor auxiliar chamado aux, para ser usado para esse propósito. Eu até tentei criar vetores dentro da função intercala, mas isso deixou tudo bem lento, desse jeito a velocidade ficou um pouco melhor. Por último, mas bem legal, é que nesse caso, veja que temos três funções, uma que recebe e retorna tudo para o R, e outras duas que fazem o trabalho propriamente dito. Mas quando a gente olha no R, desse jeito a gente só vai ver a mergesortC, isso porque em cima dela tem escrito

// [[Rcpp::export]]

Isso faz com que a gente só veja ela, e não as funções auxiliares.

Concluído a construção da função em C++, agora é simples no R

> library(Rcpp) > sourceCpp("merge.cpp") > mergesortC(sample(1:10)) [1] 1 2 3 4 5 6 7 8 9 10

Basicamente, só damos um sourceCpp no arquivo com o código, podemos usar o argumento verbose=TRUE se quisermos acompanhar a compilação e ta tudo pronto para usar, e podemos ver que o mergesortC funciona certinho.

Agora a gente pode até medir o tempo que ele precisa para organizar, por exemplo, 1000 números inteiros, usando a função system.time.

> system.time(mergesortC(sample(1:1000))) usuário sistema decorrido 0.000 0.000 0.001

Alias, para ter uma noção de porque as pessoas continuaram inventando algorítimos de ordenação, a gente pode fazer uma comparação, entre o mergesort e o insertionsort por exemplo, e ver qual deles é o mais rápido, para organizar um vetor de um tamanho qualquer, é so a gente ficar bagunçando um vetor, com sample, e registar o tempo que ele leva para organizar, e repetir isso varias vezes.

dados<-data.frame(Tempo=rep(NA,1000),Algoritimo=rep(c("InsertionSort","MergeSort"),each=500)) for(i in 1:1000){ if(i<=500){ dados[i,1]<-system.time(insertionsortC(sample(1:10000)))[3] }else{ dados[i,1]<-system.time(mergesortC(sample(1:10000)))[3] } print(i) }

A gente cria então um dataframe para receber os dados, e faz o experimento, pegando o tempo total gasto no processo, o insertionsortC eu peguei do post passado.

O resultado disso é o seguinte.

01

É, o mergesort pode ser mais complicado de entender, e implementar, mas é muito mais rápido.

> aggregate(dados$Tempo, list(dados$Algoritimo), mean) Group.1 x 1 InsertionSort 0.043494 2 MergeSort 0.005362 > aggregate(dados$Tempo, list(dados$Algoritimo), sd) Group.1 x 1 InsertionSort 0.0044236978 2 MergeSort 0.0004933993

Ele é muitas vezes mais rápido, e varia muito menos no tempo de execução, não é a toa que o sort do R não usa o InsertionSort, mas ele também não usa o mergesort, ou seja, tem como melhorar ainda, mas essa vai ficar para a próxima porque ja está na hora de dormir.

Bem o script está aqui embaixo além do repositório recologia no github. Esse foi o primeiro código no Rcpp que envolvia mais de uma função, e no final muitas coisas envolvem varias funções, então acho que assim vamos abrindo portas para possibilidades legais. Outra coisa é que esse tipo de estratégia para resolver problemas, de partir o problema em partes menores, e resolver essas partes antes de resolver um problema maior, pode parecer algo longe para biologia, mas a gente vai pela mesma estrada em problemas como o de alinhamento de sequências. E vamos dormir.

library(Rcpp)
 
sourceCpp("merge.cpp",verbose=T)
 
mergesortC(sample(1:10))
 
system.time(mergesortC(sample(1:1000)))
 
cppFunction("
    NumericVector insertionsortC(NumericVector vetor) {
        int n = vetor.size();
 
        double aux;
        int i , j;
 
        for(i=1;i<n;i++) {
            aux=vetor[i];
            j=i-1;
            while(j>=0 && vetor[j]>aux) {
                vetor[j+1]=vetor[j];
                j=j-1;
                }
            vetor[j+1]=aux;
            }
        return vetor;
        }
")
 
dados<-data.frame(Tempo=rep(NA,1000),Algoritimo=rep(c("InsertionSort","MergeSort"),each=500))
 
for(i in 1:1000){
    if(i<=500){
        dados[i,1]<-system.time(insertionsortC(sample(1:10000)))[3]
    }else{
        dados[i,1]<-system.time(mergesortC(sample(1:10000)))[3]
    }
    print(i)
}
 
#Figura 1
boxplot(Tempo~Algoritimo,data=dados)
 
aggregate(dados$Tempo, list(dados$Algoritimo), mean)
aggregate(dados$Tempo, list(dados$Algoritimo), sd)